Provide a general concept of orientation that will be general for whatever dimension or polytope about gridap.jl

[fverdugo]

@santiagobadia, continuing the discussion about the normals.

my goal is to construct CellField objects that represent the normal vector to a given boundary. I want that the normal vector is the outwards one in order the divergence theorem to hold.

I want to achieve an API like this one:

using Gridap
import Gridap: ∇

ufun(x) = VectorValue(x[1]*x[2], x[2])
grad_ufun(x) = TensorValue(x[2],x[1],0.0,1.0)
∇(::typeof(ufun)) = grad_ufun

model = CartesianDiscreteModel(partition=(4,4,4))
trian = Triangulation(model)
trian_Γ = BoundaryTriangulation(model, "boundary") # An integration mesh of the entire boundary

quad = CellQuadrature(trian, order=2)
quad_Γ = CellQuadrature(trian_Γ, order=2)

u = CellField(trian, ufun)
u_Γ = restrict(u, trian_Γ)

n = NormalVector(trian_Γ) # This is the major part missing

# Check that divergence theorem holds
@assert sum( integrate( div(u) , trian, quad) )  ≈ sum( integrate( u_Γ * n , trian_Γ, quad_Γ) ) 

The function NormalVector is the major part missing (the divergence operator div is also missing, but this is minor).

My first guess is to construct the normal from the Jacobian of the CellGeomap of the boundary triangulation. I.e.:

phi = CellGeomap(trian_Γ)
jac = ∇(phi)

The jac object is a CellField whose value at a given point represents the vectors "tangent" to the boundary. From these tangent vectors I can construct a normal (e.g., with the cross product in 3D). However, this normal is not guaranteed to point outwards... Thus, I need a way to decide if I need to multiply the normal by -1 or not. That's why I need to extract this info from the polytope...

Do you understand what I mean? Do you have perhaps a better way to compute the normals?

[santiagobadia]

Hi @fverdugo ,

I think that the general way to define outward normals is: 1) Define the outward normal in the reference FE, using the nullspace of the squared matrix of tangent vectors, etc; 2) Since the tangent vectors transform with the local change of basis with the Jacobian matrix J, the normal vector in the physical space transforms with J^-T, in order to keep orthogonality with the tangent vectors in the physical space. We already have the reference normal N and J, so we can readily compute the normal vector at every cell as J^-T N. This expression is general for whatever dimension, polytope, and geometrical map (e.g., non-affine maps for curved FEs). Check, e.g., this article for more details.

What do you think?

[fverdugo]

hi @santiagobadia

Yes, we can follow this approach. It is simple, and we already have all ingredients.

I have just a question. This approach assumes that each face in the surface mesh touches a volume cell in some volumetric mesh (this is actually the case of our BoundaryTriangulation. Thus, everything OK). However, in some cases one can have surface meshes that do not have an associated volume mesh (e.g., an STL). In that case, how do you would compute the normals? In any case, solving this is not critical for the moment.

[santiagobadia]

@fverdugo, I think that if you talk about general manifolds in arbitrary dimension, the situation is more complicated. First, you need to define the concept of orientable (or oriented, not sure) manifold. Not all manifolds are orientable. Assuming that you have an orientable manifold, you could probably define the normal component of a facet in the reference domain (in dimension+1, normal = (0,...0,1) and use the tranpose of the Jacobian to determine the normal component of the manifold. Again, given a oriented manifold, the normal would always point to the same side if some restrictions apply. It seems that one way to define an orientable map is via triangulation, enforcing all triangles be clockwise (or counterclockwise). If we assume that the manifold is the boundary of domain, the manifold is orientable and we can make the normal to point outwards using what I proposed. By the way, STL meshes should provide a consistent normal, so in this case, I think that there is no problem, it is just data. They are (as far as I can say) a representation of an orientable manifold. Not sure about STL for Mobius strips.

Surely life is simpler in 3D using the vector product (and it can also be used in 2D by appending the Z-component). But there is no extension of the vector product in higher dimension. We can probably consider a particular case for 2D/3D manifolds in 3D/4D when needed. Or study a little bit more of differential geometry and try to provide an abstract definition. My differential geometry knowledge is quite basic.

[fverdugo]

@santiagobadia work done!

using Gridap

tags = [23,24,25]
model = CartesianDiscreteModel(partition=(3,4,2))

btrian = BoundaryTriangulation(model,tags)

n = NormalVector(btrian) # New function!

bquad = CellQuadrature(btrian,order=2)
q = coordinates(bquad)

n_q = evaluate(n,q)

bphi = CellGeomap(btrian)

x = evaluate(bphi,q)

writevtk(btrian,"btrian")
writevtk(x,"x",pointdata=["n"=>n_q])

[fverdugo]

Closed via PR #60