Setting all cells to artificial prior to refinement

[dougshidong]

I am currently using a Manifold that depends on a MappingFEField in the refinement process. However, in the current implementation of the parallel::distributed::Triangulation, all the cells are set to artificial. As a a result, the MappingFEField cannot reference to those cells since they are not locally_owned anymore.

dealii/source/distributed/tria.cc

Lines 3526 to 3532 in 7751fe9

	// set all cells to artificial. we will later set it to the correct
	// subdomain in match_tree_recursively
	for (typename Triangulation<dim, spacedim>::cell_iterator cell =
	this->begin(0);
	cell != this->end(0);
	++cell)
	cell->recursively_set_subdomain_id(numbers::artificial_subdomain_id);

Is setting all those cells to artificial necessary? Shouldn't the subdomain assignment after the refinement process go through all the cells? The way I see it, this line of code is like resetting an entire vector to 0, even though we will assign a new value to it later. We simply lose the possibly relevant information when resetting it. I have tried to comment this out and my simple case (global uniform refinement) and it works fine.

If not necessary, can we remove this part of the code since it would allow the Manifold to depend on the current Triangulation? Otherwise, I can still maintain this in my own fork.

[peterrum]

@dougshidong Could you add an example code that demonstrates your issue?

[dougshidong]

This is related with my current effort to basically have a high order mesh that I can deform and adapt as described by my post. It is rather hard to provide a succinct code without going into some details. I think this issue can be explained just by showing the class I'm trying to derive from Manifold.

Below is a snippet of the code I am using:

https://github.com/dougshidong/PHiLiP/blob/7e2dbb45b1f770ed067f6b5d49e914e2ef3f04da/tests/unit_tests/grid/mapping_manifold.cpp#L23-L109

The important parts are line 28-35, which shows that the Manifold has a Triangulation and a MappingFEField. During the refinement process, this Manifold will then use get_intermediate_point() to generate new point locations. This function will then use MappingFEField to evaluate the current cell geometry and determine the exact midpoints on the curved cell. However, line 74, 84, 100, 101, and 105 will hit an assertion because all cells have been set to artificial as described in the first post.

[tjhei]

The important parts are line 28-35, which shows that the Manifold has a Triangulation and a MappingFEField. During the refinement process, this Manifold will then use get_intermediate_point() to generate new point locations.

Ah, that makes sense. You need access to the DoFs on the old mesh, while generating the new mesh during refinement.

Ideally, we would keep the old mesh intact while creating the new mesh, but we can't easily do that. We also need to first set all ids to invalid before setting exactly the ones we need. Otherwise, we might end up with invalid ids set in parts that are no longer on this machine (when repartitioning or coarsening). So your change would possibly break refinement in some situations.

I think the problem is more complicated though: The MappingFEField assumes that the refinement is done on a processor that had knowledge and ownership of the parent before the refinement. In general, this assumption is false as the set of owned cells can change with an adaptive refinement step.

Sorry, but I don't think there is an easy solution to this.

[dougshidong]

Ideally, we would keep the old mesh intact while creating the new mesh, but we can't easily do that. We also need to first set all ids to invalid before setting exactly the ones we need. Otherwise, we might end up with invalid ids set in parts that are no longer on this machine (when repartitioning or coarsening). So your change would possibly break refinement in some situations.

Can we go through those cells after and set them to artificial? I can hold a list of visited cells that have been tagged as locally owned or ghost cells, then go through the other cells and tag them as artificial.

I think the problem is more complicated though: The MappingFEField assumes that the refinement is done on a processor that had knowledge and ownership of the parent before the refinement. In general, this assumption is false as the set of owned cells can change with an adaptive refinement step.

Yes, that is true. Currently, my solution to this is to perform a SolutionTransfer on the degrees of freedom euler_vector. Not 100% sure that would work, but I don't currently see a reason it wouldn't.

Sorry, but I don't think there is an easy solution to this.

Please let me know if either solution does not seem viable. I'll have to go back to the drawing board.

[tjhei]

Can we go through those cells after and set them to artificial? I can hold a list of visited cells that have been tagged as locally owned or ghost cells, then go through the other cells and tag them as artificial.

yes, we could do that, but it doesn't fix the problem I mention below.

Currently, my solution to this is to perform a SolutionTransfer on the degrees of freedom euler_vector. Not 100% sure that would work, but I don't currently see a reason it wouldn't.

No, this won't work in general. It is likely that you have to refine a cell on a processor that did not have the parent before. So your solution won't work. The SolutionTransfer will be done after the mesh has been refined, but you will need the information during the refinement to place new points.

Have you considered using MappingQEulerian or MappingQ1Eulerian instead? This way, the mesh will stay in the reference configuration and refinement does not require access to entries that don't exist locally.

Alternatively, you could try removing the manifold before doing refinement and then moving points of the newly created cells around. It requires writing some fixup code, though.

[dougshidong]

Oh I see where the issue lies, the processor doesn't necessarily own the cell it is refining. Had to re-read your first reply.

Have you considered using MappingQEulerian or MappingQ1Eulerian instead? This way, the mesh will stay in the reference configuration and refinement does not require access to entries that don't exist locally.

Wouldn't the same issue arise with MappingQEulerian? From my understanding, it is basically identical to MappingFEField, except that it takes relative positions instead of absolute. During the refinement process, I would still need to use transform_real_to_unit_cell and transform_unit_to_real_cell, which requires me to access those non-local entries.

Alternatively, you could try removing the manifold before doing refinement and then moving points of the newly created cells around. It requires writing some fixup code, though.

I guess that in this scenario, I would simply never need a Manifold. Simply use the FlatManifold at all times. However, even in that scenario, I need to find a way to "refine" my Mapping. It kind of becomes this circular issue where I don't have the correct Triangulation to transfer my Mapping solution, and I don't have the correct refined Mapping to move the Triangulation.

[tjhei]

which requires me to access those non-local entries.

wait, now I am confused. Why do you need anything but a flat manifold when using Eulerian or FEField mapping? They are made in such a way that the mesh stays in reference configuration. SolutionTransfer takes care of interpolation of your displacement vector.

[dougshidong]

Oh wow, ok yeah. That completely solves the issue. We don't need to move the code to change the artificial cells.

I guess that I did not fully understand how MappingFEField works. I thought that the cells vertices had to match between the Triangulation and the MappingFEField, or at least that there was a link from the Triangulation to the MappingFEField.

I now realize that when using the Triangulation with a MappingFEField, the Triangulation vertices' locations are now "useless" after the initial get_position_vector(). They only become used for connectivity purposes. For all we care, I could move all my Triangulation's vertices to a single point at the origin, and I would still have a valid domain to perform computations (I just tried it and it works).

I've been tackling and strategizing this non-issue for a week. Thank you so much for the help, it's greatly appreciated.