Implementation of get_tangent_vector with similar behavior as get_normal_vector

[aerappa]

My motivation for the PR came from the desire to implement a DG scheme for vector value problems with jumping tangential components, e.g.

where the tangent and normal are define schematically via

I followed the outline I found for get_normal_vector that works for both SkeletonTrianguations and BoundaryTriangulations but instead of using get_facet_normal I used get_edge_tangent for the underlying Polytope.

Here are the resulting normal (resp. tangent) cell fields for a BoundaryTriangulation

And idem for the plus component of a SkeletonTriangulation

I still need to add some tests in the spirit of those in SkeletonTriangulationTests and BoundaryTriangulationTests but I wanted to make this PR first and get some feedback. I also think there could probably be a more elegant way to unify the logic of get_normal_vector and get_tangent_vector to reduce code duplication but this would also probably obfuscate the code a bit.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 89.06%. Comparing base (b568864) to head (e7726f4).
Report is 8 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1071   +/-   ##
=======================================
  Coverage   89.05%   89.06%           
=======================================
  Files         197      197           
  Lines       23905    23927   +22     
=======================================
+ Hits        21289    21310   +21     
- Misses       2616     2617    +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[JordiManyer]

Hi @aerappa , thanks for the contribution!

Some comments:

• What is the purpose of get_normal_or_tangen_vector(...)? I don't think that should stay, right?
• Don't worry about code duplication, I think it's manageable. Maybe the part pushing the vectors to the physical space could be reused, but I think it's fine.

The big question/doubt I have is the following:

• When computing the normal, it is guaranteed that the reference normal faces outward. This means that on a shared facet the left and right normals are the opposite of eachother. We use that fact while computing the jump.
• I am not sure that the same can be said for the tangent. In a 2D Cartesian mesh, for instance, I am quite sure that the tangents are consistent (i.e they are the same on the left and right cells) if you take them oriented like the nodes in the reference edge.
  So I would be careful and check that the drawing you posted is indeed correct, I am not convinced every single case will be like that. In general, I am not sure we can guarantee a property like the one the normal has, unless we compute the reference tangent such that it has a a fixed orientation wrt the reference normal.

Let me know how the testing goes.

[aerappa]

Hey @JordiManyer

What is the purpose of get_normal_or_tangen_vector(...)? I don't think that should stay, right?

I will delete this and just duplicate. I agree it's not useful

The big question/doubt I have is the following:

• When computing the normal, it is guaranteed that the reference normal faces outward. This means that on a shared facet the left and right normals are the opposite of eachother. We use that fact while computing the jump.

Yes, you are right. The quick and dirty solution I can think of is the following

function get_edge_tangent(trian::SkeletonTriangulation)
  plus = get_edge_tangent(trian.plus)
  minus = lazy_map(x -> -x, plus)
  SkeletonPair(plus,minus)
end

although I'm not sure what the implications are for "ignoring" the information from minus. Otherwise, there are more involved approaches I can think of:

1. In 2D as someone mentioned on gitter, we can just rotate all the normals 90 degrees
2. Forget the reference element and define a canonical orientation for each edge (e.g. low to high node GID) and set plus to be the canonical orientation and minus the opposite

Let me know what you think

[JordiManyer]

Here are my thoughts:

• In 2D, each edge is also a facet (shared by 2 cells). Then it makes sense to compute a jump over an edge, so we indeed have the image you show above.
• In higher dimension this does not really make sense, does it? Each edge might be shared by an arbitrary number of faces/cells depending on the mesh, and I do not believe we even support skeleton-like terms on the edges (do they even make sense ?). So do we even want to support this?

I believe you want this for AMR aposteriori estimators, right? I have recently dipped into the subject, but I don't have much experience with the state of the art. Are edge-based tangencial jumps ever needed in 3D?

I would start by implementing this in 2D by taking a rotation of the normals. I would explicitly disable the 3D case for now. We might want to implement it in the future.

How do you see it?

[aerappa]

How do you see it?

Sounds good to me. Another benefit of the rotation idea is that the tangent vector has a consistent direction on the physical boundary.

I have implemented this idea (adding an error when called on anything but 2D) and added corresponding tests. Let me know what you think.

[JordiManyer]

Some minor changes, but I think it is mostly done. Thank you!

[JordiManyer]

@aerappa could you update NEWS.md with the changes? I'll merge after that is done.