Implementation of get_tangent_vector with similar behavior as get_normal_vector

src/CellData/CellData.jl

@@ -61,6 +61,7 @@ export Integrand
 export ∫
 export CellDof
 export get_normal_vector
+export get_tangent_vector
 export get_cell_measure
 export Interpolable
 export KDTreeSearch

src/CellData/CellFields.jl

@@ -113,11 +113,20 @@
 
 function get_normal_vector(trian::Triangulation)
   cell_normal = get_facet_normal(trian)
-  get_normal_vector(trian,cell_normal)
+  get_normal_vector(trian, cell_normal)
 end
 
-function get_normal_vector(trian::Triangulation,cell_normal::AbstractArray)
-  GenericCellField(cell_normal,trian,ReferenceDomain())
+function get_tangent_vector(trian::Triangulation)
+  cell_tangent = get_edge_tangent(trian)
+  get_tangent_vector(trian, cell_tangent)
+end
+
+function get_normal_vector(trian::Triangulation,cell_vectors::AbstractArray)
+  GenericCellField(cell_vectors,trian,ReferenceDomain())
+end
+
+function get_tangent_vector(trian::Triangulation,cell_vectors::AbstractArray)
+  GenericCellField(cell_vectors,trian,ReferenceDomain())
 end
 
 evaluate!(cache,f::Function,x::CellPoint) = CellField(f,get_triangulation(x))(x)
@@ -712,6 +721,12 @@
   SkeletonPair(plus,minus)
 end
 
+function get_tangent_vector(trian::Triangulation,cell_tangent::SkeletonPair)
+  plus = get_normal_vector(trian,cell_tangent.plus)
+  minus = get_normal_vector(trian,cell_tangent.minus)
+  SkeletonPair(plus,minus)
+end
+
 for op in (:outer,:*,:dot)
   @eval begin
     ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)

src/Exports.jl

@@ -141,6 +141,7 @@ using Gridap.TensorValues: ⊗; export ⊗
 @publish CellData mean
 @publish CellData update_state!
 @publish CellData get_normal_vector
+@publish CellData get_tangent_vector
 using Gridap.CellData: ∫; export ∫
 @publish CellData get_cell_measure
 @publish CellData get_physical_coordinate

src/Geometry/BoundaryTriangulations.jl

@@ -238,10 +238,32 @@ function get_facet_normal(trian::BoundaryTriangulation, boundary_trian_glue::Fac
   Fields.MemoArray(face_s_n)
 end
 
+function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue)
+
+  face_s_n = get_facet_normal(trian, boundary_trian_glue)
+
+  function rotate_normal_2d(n)
+    t = VectorValue(n[2], -n[1])
+    m = sqrt(t[1]^2 + t[2]^2)
+    return m < eps() ? VectorValue(0.0, 0.0) : t
+  end
+
+  face_s_t = lazy_map(Operation(rotate_normal_2d), face_s_n)
+
+  return Fields.MemoArray(face_s_t)
+end
+
+function get_edge_tangent(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A}
+  Dp != 2 && error("get_edge_tangent only implemented for 2D")
+  glue = trian.glue
+  get_edge_tangent(trian, glue)
+end
+
 function get_facet_normal(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A}
   glue = trian.glue
   get_facet_normal(trian, glue)
 end
+
 function push_normal(invJt,n)
   v = invJt⋅n
   m = sqrt(inner(v,v))

src/Geometry/Geometry.jl

@@ -70,6 +70,7 @@ import Gridap.ReferenceFEs: num_cell_dims
 import Gridap.ReferenceFEs: num_point_dims
 import Gridap.ReferenceFEs: simplexify
 import Gridap.ReferenceFEs: get_facet_normal
+import Gridap.ReferenceFEs: get_edge_tangent
 import Gridap.ReferenceFEs: Quadrature
 
 export GridTopology
@@ -107,6 +108,7 @@ export get_cell_ref_coordinates
 export get_cell_reffe
 export get_cell_shapefuns
 export get_facet_normal
+export get_edge_tangent
 export test_triangulation
 export get_cell_map
 

src/Geometry/SkeletonTriangulations.jl

@@ -94,6 +94,13 @@ function get_facet_normal(trian::SkeletonTriangulation)
   SkeletonPair(plus,minus)
 end
 
+function get_edge_tangent(trian::SkeletonTriangulation)
+  plus = get_edge_tangent(trian.plus)
+  # Flip the minus side
+  minus = lazy_map(x -> -x, plus)
+  SkeletonPair(plus,minus)
+end
+
 # Related with CompositeTriangulation
 function _compose_glues(rglue::FaceToFaceGlue,dglue::SkeletonPair)
   plus = _compose_glues(rglue,dglue.plus)

src/Geometry/Triangulations.jl

@@ -41,6 +41,7 @@
 get_reffes(trian::Triangulation) = get_reffes(get_grid(trian))
 get_cell_type(trian::Triangulation) = get_cell_type(get_grid(trian))
 get_facet_normal(trian::Triangulation) = get_facet_normal(get_grid(trian))
+get_edge_tangent(trian::Triangulation) = get_edge_tangent(get_grid(trian))
 
 # The following are not strictly needed, sine there is a default implementation for them.
 # In any case, we delegate just in case the underlying grid defines more performant versions

test/GeometryTests/BoundaryTriangulationsTests.jl

@@ -110,6 +110,11 @@ face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
 test_array(face_to_nvec_s,collect(face_to_nvec_s))
 @test isa(face_to_nvec_s,Geometry.FaceCompressedVector)
 
+face_to_tvec = get_edge_tangent(btrian)
+face_to_tvec_s = lazy_map(evaluate,face_to_tvec,face_to_s)
+test_array(face_to_tvec_s,collect(face_to_tvec_s))
+@test isa(face_to_tvec_s,Geometry.FaceCompressedVector)
+
 #print_op_tree(face_shapefuns_q)
 #print_op_tree(face_grad_shapefuns_q)
 #print_op_tree(face_to_nvec_s)
@@ -147,6 +152,16 @@ r = Vector{Point{2,Float64}}[
   [(0.0,1.0),(0.0,1.0)],[(1.0,-0.0),(1.0,-0.0)]]
 test_array(nvec_s,r)
 
+tvec = get_edge_tangent(btrian)
+
+tvec_s = lazy_map(evaluate,tvec,s)
+
+rotate_90(x) = [Point(x[1][2], -x[1][1]), Point(x[2][2], -x[2][1])]
+
+rr = rotate_90.(r)
+
+test_array(rr, tvec_s)
+
 cellids = collect(1:num_cells(model))
 
 face_to_cellid = lazy_map(Reindex(cellids),glue.tface_to_mface)

test/GeometryTests/SkeletonTriangulationsTests.jl

@@ -1,6 +1,7 @@
 module SkeletonTriangulationsTests
 
 using Test
+using Gridap
 using Gridap.TensorValues
 using Gridap.Arrays
 using Gridap.Fields
@@ -52,6 +53,7 @@ vn = get_facet_normal(vtrian)
 Ω = Triangulation(model)
 Γ = BoundaryTriangulation(model)
 Λ = SkeletonTriangulation(Γ)
+normal = get_normal_vector(Λ)
 @test Λ.rtrian === Γ
 @test isa(Λ.dtrian,SkeletonTriangulation)
 glue = get_glue(Λ,Val(3))
@@ -116,7 +118,13 @@ face_own_dofs = get_face_own_dofs(reffe,conf)
 strian = SkeletonTriangulation(model,reffe,face_own_dofs)
 test_triangulation(strian)
 ns = get_facet_normal(strian)
+ts = get_edge_tangent(strian)
 @test length(ns.⁺) == num_cells(strian)
+@test length(ts.⁺) == num_cells(strian)
+# Test orthogonality
+should_be_zero = lazy_map(Broadcasting(Operation(dot)), ns.plus, ts.plus)
+ndot_t_evaluated = evaluate(should_be_zero, VectorValue.(0:0.05:1))
+@test maximum(ndot_t_evaluated) < 1e-15
 
 #using Gridap.Visualization
 #writevtk(strian,"strian",cellfields=["normal"=>ns])