Current approach to manifold mesh generation

[amartinhuertas]

Hi @fverdugo, @santiagobadia,

I have spent some time looking at the approach we are following to generate discretizations of a manifold, to represent them in Gridap.Geometry, and to compute the manifold mesh outward unit normal (I think @fverdugo has mentioned that the current approach is temporary, not sure if the following is among the causes.)

While this approach is clearly functional (e.g., it has let us solve PDE problems on the cubed sphere), I am not sure about the following:

• does it faithfully represent the underlying mathematical abstraction?
• is it general enough?
• it it efficient enough?

if the answer to any of these questions is "not enough", this may require some more profound thinking, that perhaps we can incorporate in the future to the very much needed Gridap.Geometry refactorization that we are anticipating (for different reasons, though).

(The following is inferred from my understanding of what I have read in the code, if I am not right, confused, etc., I am sorry ... let me please know if this is the case)

The current approach relies on BoundaryDiscreteModel. This data type is created from a volume discrete model. Taking profit of this, Triangulation(::BoundaryDiscreteModel) returns a BoundaryTriangulation instance. Using this latter instance you can compute the manifold's outward unit normal using the function get_normal_vector(::BoundaryTriangulation), which uses the inverse of the volume mesh cell Jacobian, the normal in reference space, etc. So far so good.

However, I see the following issues:

1. With BoundaryDiscreteModel/BoundaryTriangulation, by conception, there is a volume mesh in the background. However, my understanding is that manifolds can exist without a volume in the background, e.g., you might read a 2D manifold mesh embedded in 3D directly from data files, e.g., in STL or GeoMesh format.
2. Even if in the particular case of BoundaryDiscreteModel, there is a volume mesh in the background, we may compute the normal in a cheaper way from the computation POV (e.g., computation of inverse of volume mesh map no longer necessary).

I think that Triangulation(::BoundaryDiscreteModel) should return a Triangulation instance which is decoupled from the volume mesh. For example, UnstructuredGrid{2,3,Tp,Orientation} is a Triangulation instance which exists and has sense without a volume mesh in the background. (@fverdugo, is there any other data type that also fulfills this?) We can provide a new function, say, get_manifold_normal_vector, defined for Triangulation{Dc,Dp}, with Dc+1==Dp, which returns the unit outward normal to the manifold.

This normal vector can be computed without using the volume mesh cell mapping more efficiently. For example, this would look like as (for 2D and 3D, generalizable for any D?):

"""
    get_manifold_normal(trian::Triangulation{Dc,Dp}) where {Dc,Dp}
"""
function get_manifold_normal(trian::Triangulation{Dc,Dp}) where {Dc,Dp}
  @assert Dc+1==Dp
  function _unit_outward_normal(v::MultiValue{Tuple{2,3}})
    n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
    n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
    n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
    n = VectorValue(n1,n2,n3)
    n/norm(n)
  end
  function _unit_outward_normal(v::MultiValue{Tuple{1,2}})
    n1 = v[2,1]
    n2 = -v[1,1]
    n = VectorValue(n1,n2)
    n/norm(n)
  end
  map = get_cell_map(trian)
  Jt = lazy_map(∇,map)
  lazy_map(Operation(_unit_outward_normal),Jt)
end

The only caveat with this approach is that the cells in the input Triangulation (e.g., quadrilaterals/triangles in 3D, segments in 2D) have to be oriented accordingly to the global manifold, otherwise one gets the normal with the wrong sign.

For example, the following pseudocode:

model=CartesianDiscreteModel((0,1,0,1),(3,3))
trian = Triangulation(ReferenceFE{1},model)
manifold_trian = RestrictedTriangulation(bgface_trian,"boundary")

which returns a Triangulation of the boundary of a box, does not fullfilling this pre-requisite ... We should "orient" the edges in manifold_trian ... but I guess that using different semantics for the Orientation trait we use for volume meshes. (If I am not wrong, for volume meshes Oriented means that, for all faces, the local orientation of each face within each cell matches from the perspective of either cell around, and NonOriented, the opposite.) Besides the Orientation trait is not common for all Triangulations .....

Does this sound reasonable? What is your view? Sure I am missing a lot of things, there are smarter solutions, etc.

[fverdugo]

BoundaryDiscreteModel is likely to be removed in the future, hopefully near future.

About the computation of normals. I remember that we discuss this quite in detail with @santiagobadia and ended up with the current solution. As long as I remember the reasons include:

• In general, it is not possible to "orient" a manifold that is not the boundary of a volume, e.g., the möbius strip.
• When you read a manifold, eg from a STL, it also contains the normal vectors. The most efficient way would be to reuse these.
• The current approach is dimension independent, works also for D>3.

[amartinhuertas]

BoundaryDiscreteModel is likely to be removed in the future, hopefully near future.

Ok, but if you remove it, you have to develop something alternative to provide similar functionality, right?

I remember that we discuss this quite in detail with @santiagobadia and ended up with the current solution.

Ok. I see which are the benefits of the current solution. Using the inverse of the Jacobian of the cell combined with the normal in reference space you always get the outward unit normal, no matter whether the determinant of the Jacobian is positive or negative.

In general, it is not possible to "orient" a manifold that is not the boundary of a volume, e.g., the möbius strip.

Contravariant Piola elements are not available for non-orientable manifolds.

When you read a manifold, eg from a STL, it also contains the normal vectors. The most efficient way would be to reuse these.

Ok, that makes sense.

The current approach is dimension independent, works also for D>3.

Yest, at the cost of storing a volume mesh which is not actually needed.

[amartinhuertas]

Anyway, I can live with the current solution, I just wanted to contribute to the discussion ...

[amartinhuertas]

BTW, not 100% related to this discussion, but I found a geometry-related paper that has been recently accepted in ACM TOMS here: https://arxiv.org/abs/1910.09824

[fverdugo]

I think that the problem/missing part now is that we are not able to read eg from gmsh a discrete model on a manifold. In particular, we do not have a clear place to store the normal vector that gmsh would provide. This can/needs to be included in the refactoring of the geometry module I am planning.

[amartinhuertas]

I think that the problem/missing part now is that we are not able to read eg from gmsh a discrete model on a manifold. In particular, we do not have a clear place to store the normal vector that gmsh would provide. This can/needs to be included in the refactoring of the geometry module I am planning.

Yes, agreed. But I think we are also interested into generating manifold meshes in-core, i..e, without relying on data files, etc. E.g., a cubed mesh sphere, or an icosohedral mesh grid. (https://en.wikipedia.org/wiki/Geodesic_polyhedron)

With regards to the refactoring, it might help that we sketch a plan before starting to code (not now, just whenever we get into this). I might be able to help in this regard ...

[amartinhuertas]

BTW ... I forgot ... the following line is NOT correct for manifolds:

Gridap.jl/src/Geometry/SkeletonTriangulations.jl

Line 256 in f1ad186

minus = lazy_map(Broadcasting(Operation(-)),plus)

Two adjacent cells in a manifold might not be co-planar ...