Behavior of GridGenerator::extract_boundary_mesh

[Sebastian Stark]

Hello all,

I'm dealing with problems which involve unknowns in the (three dimemsional) domain and on interfaces separating several portions of this domain. To generate the triangulation of the interfaces from the triangulation of the domain, I decided to use a modified version of  GridGenerator::extract_boundary_mesh. While going over this function, I observed the following, which might be a bug: In order to get consistent orientations of the cells on the boundary, vertices 1 and 2 are swapped for cells of the boundary triangulation corresponding to faces 1, 3, 5 of the underlying cells of the domain. However, this does also change the ordering of the children of these cells; and I think this is not taken into account when later generating the map between the cells of the boundary triangulation and the underlying faces of the cells of the domain. I have attached an example demonstrating the problem.

Another related behavior of GridGenerator::extract_boundary_mesh, which is not quite as I would expect it, is the following: If I'm not missing something, GridGenerator::extract_boundary_mesh produces a mesh with the normal field generally pointing into the underlying domain in the three dimensional case. In 2 dimensions, no vertex swapping is done in GridGenerator::extract_boundary_mesh and the normal field points either into or out of the domain, depending on the exact situation (what happens depends also on the algorithm setting the direction_flag in create_triangulation). I think in regard to dimension independent programming this is not really what one wants to have. In my modified version of  extract_boundary_mesh I currently swap vertices 1, 2 of faces 2, 4, 6 in 3 dimensions and vertices 0, 1 of edges 1, 2 in two dimensions. This way I always get normals pointing outward to the underlying domain. However, as I am not familiar with the internals of the library, I'm unsure whether this is a reasonable approach.

Best,

Sebastian

[Wolfgang Bangerth]

Sebastian,

> I'm dealing with problems which involve unknowns in the (three dimemsional)
> domain and on interfaces separating several portions of this domain. To
> generate the triangulation of the interfaces from the triangulation of the
> domain, I decided to use a modified version of
> GridGenerator::extract_boundary_mesh. While going over this function, I
> observed the following, which might be a bug: In order to get consistent
> orientations of the cells on the boundary, vertices 1 and 2 are swapped for
> cells of the boundary triangulation corresponding to faces 1, 3, 5 of the
> underlying cells of the domain. However, this does also change the ordering of
> the children of these cells; and I think this is not taken into account when
> later generating the map between the cells of the boundary triangulation and
> the underlying faces of the cells of the domain. I have attached an example
> demonstrating the problem.

Yes, that seems entirely plausible and is a bug.

Since you already read through the implementation of the function, would you
be interested in writing a patch to fix this? There is a description of how to
contribute here
https://github.com/dealii/dealii/wiki/Contributing
but we'll also be happy to help you with the process!

> Another related behavior of GridGenerator::extract_boundary_mesh, which is not
> quite as I would expect it, is the following: If I'm not missing something,
> GridGenerator::extract_boundary_mesh produces a mesh with the normal field
> generally pointing into the underlying domain in the three dimensional case.

Yes, this was recently pointed out to me in a private mail as well. In some
sense, it doesn't really matter very much -- all that's necessary is that one
ends up with a consistent normal vector to the manifold. But it's true that if
we were to implement this function again, we'd probably choose the normal
vector to point outward, as it usually does.

> In 2 dimensions, no vertex swapping is done in
> GridGenerator::extract_boundary_mesh and the normal field points either into
> or out of the domain, depending on the exact situation (what happens depends
> also on the algorithm setting the direction_flag in create_triangulation). I
> think in regard to dimension independent programming this is not really what
> one wants to have. In my modified version of  extract_boundary_mesh I
> currently swap vertices 1, 2 of faces 2, 4, 6 in 3 dimensions and vertices 0,
> 1 of edges 1, 2 in two dimensions. This way I always get normals pointing
> outward to the underlying domain. However, as I am not familiar with the
> internals of the library, I'm unsure whether this is a reasonable approach.

I've been trying to think about the trade-off between breaking backward
compatibility and providing a consistent result. I think it's possible to call
the current behavior a bug. So if you wanted to submit patches for the things
you mention, that would be much appreciated as well!

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[stark.s...@gmx.de]

Wolfgang,

Yes, that seems entirely plausible and is a bug.

 

Since you already read through the implementation of the function, would you 

be interested in writing a patch to fix this? There is a description of how to
contribute here
   https://github.com/dealii/dealii/wiki/Contributing
but we'll also be happy to help you with the process! 

I'll use this as a learning example of how to contribute. Just to be sure that I'm not doing it overly complicated: If I only have an iterator to a face of a cell, there is no simple way of figuring out which face it is (in terms of the face number within the underlying cell), right?

 

Yes, this was recently pointed out to me in a private mail as well. In some
sense, it doesn't really matter very much -- all that's necessary is that one
ends up with a consistent normal vector to the manifold. But it's true that if
we were to implement this function again, we'd probably choose the normal
vector to point outward, as it usually does.

I'm actually not concerned too much about whether the normal vector points inward or outward (although my preferred choice would be outward) since this is just a matter of definition. In my opinion, the main issue is that the behavior depends on the spatial dimension. So, I think the best trade-off between backward compatibility and consistency would be to retain the choice that the normal vector points inward and write a patch ensuring that this is also the case in 2D. Does this make sense?

Best,

Sebastian

[Wolfgang Bangerth]

Sebastian,

> I'll use this as a learning example of how to contribute.

Great! Please do let us know if you need help with anything!

> Just to be sure that
> I'm not doing it overly complicated: If I only have an iterator to a face of a
> cell, there is no simple way of figuring out which face it is (in terms of the
> face number within the underlying cell), right?

Correct. You can ask a cell for its faces, and a face for its bounding edges
and vertices. But you can't ask questions in the other direction of this chain.

> I'm actually not concerned too much about whether the normal vector points
> inward or outward (although my preferred choice would be outward) since this
> is just a matter of definition. In my opinion, the main issue is that the
> behavior depends on the spatial dimension. So, I think the best trade-off
> between backward compatibility and consistency would be to retain the choice
> that the normal vector points inward and write a patch ensuring that this is
> also the case in 2D. Does this make sense?

Yes!

Best
W.