Merge pull request #1071 from aerappa/tangent_vector

Implementation of `get_tangent_vector` with similar behavior as `get_normal_vector`
gridap · Jan 13, 2025 · ea403aa · ea403aa
2 parents b568864 + e7726f4
commit ea403aa
Show file tree

Hide file tree

Showing 9 changed files with 83 additions and 3 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.18.10] - 2025-01-13
+
+### Added
+
+- Added corresponding function `get_tangent_vector` to `get_normal_vector`. This method calculates the (unique up to sign) tangential unit vector to edges in 2D meshes, by rotating the normal (nx, ny) -> (ny, -nx). Since PR[#1071](https://github.com/gridap/Gridap.jl/pull/1071).
+
 ## [0.18.9] - 2025-01-13
 
 

diff --git a/src/CellData/CellData.jl b/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

diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl
@@ -113,11 +113,20 @@ end
 
 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 @@ function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair)
   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)

diff --git a/src/Exports.jl 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

diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl
@@ -238,10 +238,28 @@ 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)
+
+  rotate_normal_2d(n) = VectorValue(n[2], -n[1])
+
+  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))

diff --git a/src/Geometry/Geometry.jl b/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
 

diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl
@@ -94,6 +94,12 @@ function get_facet_normal(trian::SkeletonTriangulation)
   SkeletonPair(plus,minus)
 end
 
+function get_edge_tangent(trian::SkeletonTriangulation)
+  plus = get_edge_tangent(trian.plus)
+  minus = get_edge_tangent(trian.minus)
+  SkeletonPair(plus,minus)
+end
+
 # Related with CompositeTriangulation
 function _compose_glues(rglue::FaceToFaceGlue,dglue::SkeletonPair)
   plus = _compose_glues(rglue,dglue.plus)

diff --git a/test/GeometryTests/BoundaryTriangulationsTests.jl b/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)

diff --git a/test/GeometryTests/SkeletonTriangulationsTests.jl b/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))
@@ -110,13 +112,27 @@ itrian = InterfaceTriangulation(Ω_in,Ω_out)
 #writevtk(trian,"trian",celldata=["inout"=>cell_to_inout])
 #writevtk(itrian,"itrian",nsubcells=10,cellfields=["ni"=>ni,"nl"=>nl,"nr"=>nr])
 
+ti = get_tangent_vector(itrian)
+tl = get_tangent_vector(ltrian)
+tr = get_tangent_vector(rtrian)
+
+@test ti isa SkeletonPair
+@test tl isa Gridap.CellData.GenericCellField
+@test tr isa Gridap.CellData.GenericCellField 
+
 reffe = LagrangianRefFE(Float64,QUAD,(2,2))
 conf = CDConformity((CONT,DISC))
 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])