Can one keep manifolds when merging triangulations?

[nfehn]

GridGenerator::merge_triangulations() forgets about the manifolds, see

https://www.dealii.org/current/doxygen/deal.II/namespaceGridGenerator.html#a7cd88e7eacd46697dee80ad2b8438d54

Unlike most GridGenerator functions, this function does not attach any manifolds to result, nor does it set any manifold ids.

Is there a way to keep manifolds or simplify usability as it can become very complicated to reset manifolds after merge_triangulation()?

What happens with manifolds when applying functions like GridTools::shift(), GridTools::rotate()?

[masterleinad]

That's not quite true. The documentation also says

The function copies the material ids of the cells of the two input triangulations into the output triangulation. If copy_manifold_ids is set to true, manifold ids will be copied. [...]

The manifolds assigned to objects stay the same if the position in space is altered.

[masterleinad]

As long as the Manifold objects are invariant under the transformation, you should not have any problems after moving vertices.

[nfehn]

Thanks for the explanations. I want to give further information on the second point mentioned above (moving the grid).

Let's consider the GridTools::rotate() function by the example of GridGenerator::torus<3,3>():
The code is the following

  double const r = 0.5, R = 1.5;
  GridGenerator::torus(*triangulation, R, r);
  GridTools::rotate(0.5*numbers::PI, 1 /* rotate around y */, *triangulation);

This works fine and rotate() does not pose problems as the torus geometry is rotationally symmetric w.r.t. the y-axis. The result looks like this

Now, consider rotation around the z-axis

  double const r = 0.5, R = 1.5;
  GridGenerator::torus(*triangulation, R, r);
  GridTools::rotate(0.5*numbers::PI, 2 /* rotate around z */, *triangulation);

The code throws an exception in TransfiniteInterpolationManifold, probably because the manifolds are no longer valid after the geometry has been rotated.

In my opinion, this is clearly unexpected behavior. The user does not see the manifolds, as torus() is applying them under-the-hood. So currently my point of view is the following: A general function like rotate() needs to explicitly forget about the manifolds (and this needs to be highlighted in the documentation), as the manifolds can not be preserved in a general setting. That manifolds "are invariant under the transformation" is a special case but a function like rotate() should rely on the most general case and specialized functions should be written in addition. For the torus example, I think there is no practical interest to rotate around the y-axis, as the geometry itself is rotationally symmetric w.r.t. y.

To obtain a code that at least does not throw an exception, I used the following modifcation:

  GridGenerator::torus(*triangulation, R, r);
  triangulation->reset_all_manifolds();
  GridTools::rotate(0.5*numbers::PI, 2 /* rotate around z */, *triangulation);

with the following result

I think the second line riangulation->reset_all_manifolds(); should be integrated into rotate() (and similarily for other functions where the problem is the same), as the current version is not really transparent (why should I reset manifolds after the first line of code?). The mesh above is obviously not the optimal result, but the code is stable. The user might then ask himself, where has the manifold gone? But he will then look into the documentation of rotate() and find an explanation that this is what one can guarantee in general. However, I think this needs to be discussed in a broader context. Please share your opinion!

[masterleinad]

The history for this is as follows:
We had the functions in namespace GridGenerator a long time before introducing the Manifold classes. When introducing these classes manifolds were not set by the GridGenerator functions. Instead, users had to deal with refining completely on their own. Sometime later, we noticed that it is annoying for the user to always set the same manifold objects for the same generated mesh. Hence, we decided to attach the "correct" manifold at the same time. At that point, we also talked about adapting attached manifold objects when modifying the underlying geometry, but that would not be possible with most of the manifold objects.
That being said, I agree that simply removing the attached manifolds in GridTools::transform() is a good idea.

[bangerth]

Same issue as #8920.

What if we replace the manifolds in the triangulation with ones where we "wrap" the original manifold with one level of transformation, using the same transformation as use in GridTools::transform?

[nfehn]

I like that idea. @bangerth What do you mean with "one level of transformation"? It should also work when applying an arbitrary sequence of shift/rotate operations. For a general GridTools::transform() operation this is probably not possible?

Nevertheless, I could make sense to now explicitly deactivate manifolds when applying grid manipulations and adapt the documentation in a first step (to have a basis where the code is stable and the documentation clear regarding the limitations), and to extend/improve functionality in a second step?

[nfehn]

@masterleinad Regarding the original question related to merge_triangulations():

I think I did not completely get your point. What are the current capabilities and limitations of merge_triangulations()? You talked about manifold_ids but I also mean the manifolds, for example the question, whether the grid will be refined according to the manifolds when calling a refine_global() after merge_triangulations()?

[masterleinad]

Ah, I see, I misread your question. Manifolds are indeed not copied since it was not obvious what to do if multiple triangulations to be merged use the same manifold id, but have different Manifold objects attached.

[nfehn]

If I understand you correctly it would be possible and it is a question of how to update the data structures. For example, can one increment the manifold_ids of one triangulation by the maximum manifold_id of the other one, and then merge the two triangulations?

[masterleinad]

From a user perspective, you can just loop over the manifold indicators of the input triangulations, obtain the correct object via Triangulation::get_manifold() and attach it to the new triangulation via Triangulation::set_manifold. This does not require any changes to merge_triangulations.

For me, it would be confusing if merge_triangulations would shift manifold ids on its own, though. Hence, I would rather propose something like the above in the documentation.

[nfehn]

The concerns that I have regarding the user perspective or the documentation is that some GridGenerator functions implicitly attach manifolds, so the user does not get into contact with the manifolds, i.e., we consider manifolds as part of a triangulation. This is why I would vote for a solution where merge_triangulation() cares about that. Otherwise, responsibilities are no longer clear. For me, the solution described above is rather a "developer perspective" from someone who exactly knows what one has to do.
Maybe the underlying question is whether we understand a manifold as part of a triangulation or as a standalone object. Considering it as a standalone object might also be valid but I think one can find several places in the code where this is not the case (see example mentioned above). Why not incrementing the manifolds, the new triangulation also does not have two cells with ID = 0?

[masterleinad]

Changing the manifold ids (in case they should be copied from the input triangulations) would not be backwards-compatible. Furthermore, the manifold ids would then depend on the order in which triangulations are merged. This makes it also pretty hard to argue about the manifold ids and which manifold object should be attached to which id afterwards.

In any way, I am happy to change the function signature (to allow copying manifolds) in a way that preserves the current function's behavior by default. For the other reasons above, I still like preserving the manifold ids better than incrementing, though. Let's see what other people think.

[bangerth]

I like the idea of re-attaching manifolds (by default, or by option). We've made an effort to move attaching the correct manifold to triangulations from user space into the library, and it's a bit of a shame that we don't do it here. It might be nice if we could do it here as well. That's particularly true because the function internally knows the mapping from old to new cells and can set manifold ids without too much trouble -- that's more complicated in user space.

As for shifting IDs: How about the function getting a manifold_id_shift_offset argument. If that argument is numbers::invalid_unsigned_int, then no setting of manifold ids happens; if it is a nonzero value, then that's by how much the manifold ids of the second triangulation are shifted (making things predictable to users; there then ought to be an assertion that we are shifting into a range that is not used by the first triangulation); and if it is zero, then the shift happens by one plus the largest manifold id used by the first triangulation (=user doesn't actually care about the numbers).

[masterleinad]

I like the idea of re-attaching manifolds (by default, or by option). We've made an effort to move attaching the correct manifold to triangulations from user space into the library, and it's a bit of a shame that we don't do it here. It might be nice if we could do it here as well. That's particularly true because the function internally knows the mapping from old to new cells and can set manifold ids without too much trouble -- that's more complicated in user space.

I am against changing default behavior when we can't deprecate. Hence, we likely need two additional parameters with your suggestion and the number of options seems to be getting a little bit out of hand. If you really want this functionality, I would rather propose to have a separate function that always shifts manifold ids and always re-attaches manifolds.

[nfehn]

From my perspective, the current implementation is a bug somehow (since properties that belong to the triangulation get lost), or at least an inconsistent behavior. In this sense, I think one can accept that the default behavior changes. I would be more happy with the additional parameter variant (with a default value to not change the interface), since a new function merge_triangulations_and_keep_manifolds() further complicates things and after some time users will ask themselves how it came to these names. For the case of a new function, I would prefer a new name like fuse_triangulations(), and get rid of merge_triangulations() in the long run.

@kronbichler Didn't we have changes in the past from surface manifolds to volume manifolds in order to obtain truly high-order geometry representations? In these cases the default behavior also changed?

[bangerth]

I tend to agree with @nfehn here. @dealii/dealii -- anyone else want to join in?

[Rombur]

I agree with @nfehn. I think we should change the default behavior.

[luca-heltai]

I think @nfehn has a point. The current implementation is not functional, and it is confusing to users = it is a bug. We should fix bugs, not maintain backward compatibility with them. I would be in favor of simply documenting what the behavior is, and fixing it according to what we expect it to do.

Merge triangulation should shift the manifold ids of the second triangulation by the maximum id of the first triangulation, just as it does with cell ids, and attach the manifolds accordingly. Any other behavior is in my opinion user dependent, and should be left to users.

As far as transformations are concerned, I'd be in favor of removing manifolds before any transformation. In a subsequent PR we could make sure that we can actually attach a transformation to a manifold, and manipulate it afterwards, but this is going to be very expensive, and case dependent, so I'm not sure we'll be able to it effectively (i.e.: transfinite manifold needs not be transformed: only the boundaries of the manifolds need to be transformed. How would we make this work?).

[drwells]

This is more-or-less (the discussion is different but the problem is the same) the same issue as that in #7052. I'll take a stab at finishing this for the release.

[drwells]

I haven't had time to work on this and its not critical so lets move this off 9.2.

[nfehn]

Anything new here?

[kronbichler]

It seems you have this higher on the priority list than most other developers, so I would invite you to start looking into it. I would be more than happy to assist on the way, but I do not have the resources to do it alone. Thus, we need to work towards this goal as a community. Keep in mind that deal.II is a voluntary project with no developer getting paid for these tasks (at least not more than you or me).

[nfehn]

I don't understand why I get so many accusations when asking after approximately one year? I just wanted to know whether there are already plans, or if someone else already did something in this direction or that this issue is just left for some reasons (because I do not follow everything in deal.II), etc., ...

[drwells]

I plan on getting this done in time for 9.3: I tagged myself and added it to the milestone for this reason. If you need it as soon as possible I can help you implement it!

[peterrum]

I plan on getting this done in time for 9.3: I tagged myself and added it to the milestone for this reason. If you need it as soon as possible I can help you implement it!

@drwells Thanks for offering your help here. However, I think any changes made in the context of this issue made in the next days will be too close to the release. I am adjusting the milestone.