From 9fd6eb4766c4766c3afb6f4a0c305dd1e12f8a90 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Wed, 1 Jan 2025 10:39:34 -0700 Subject: [PATCH 1/6] Tangent vector running but not tested --- src/CellData/CellFields.jl | 19 ++++-- src/Geometry/BoundaryTriangulations.jl | 61 +++++++++++++++++++ src/Geometry/Geometry.jl | 2 + src/Geometry/SkeletonTriangulations.jl | 6 ++ src/Geometry/Triangulations.jl | 1 + .../SkeletonTriangulationsTests.jl | 3 + 6 files changed, 86 insertions(+), 6 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 1158c4b9e..b60474a08 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -112,14 +112,21 @@ function CellField(f,trian::Triangulation) end function get_normal_vector(trian::Triangulation) + println("Calling get_normal_vector") cell_normal = get_facet_normal(trian) - get_normal_vector(trian,cell_normal) + get_vector_skeleton_pair(trian,cell_normal, get_normal_vector) 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_vector_skeleton_pair(trian,cell_tangent, get_edge_tangent) end +function get_normal_vector(trian::Triangulation,cell_vectors::AbstractArray) + GenericCellField(cell_vectors,trian,ReferenceDomain()) +end + + evaluate!(cache,f::Function,x::CellPoint) = CellField(f,get_triangulation(x))(x) function change_domain(a::CellField,::ReferenceDomain,::PhysicalDomain) @@ -706,9 +713,9 @@ function CellFieldAt{T}(parent::OperationCellField) where T OperationCellField(parent.op,args...) end -function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair) - plus = get_normal_vector(trian,cell_normal.plus) - minus = get_normal_vector(trian,cell_normal.minus) +function get_vector_skeleton_pair(trian::Triangulation,cell_normal::SkeletonPair, f::Function) + plus = f(trian,cell_normal.plus) + minus = f(trian,cell_normal.minus) SkeletonPair(plus,minus) end diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl index 53431876f..141593d5c 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -238,10 +238,71 @@ 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) + # 1) Retrieve the cell grid from the boundary triangulation + cell_grid = get_grid(get_background_model(trian.trian)) + + # 2) Define how we get reference tangents for each local face + function f(r) + p = get_polytope(r) + lface_to_t = get_edge_tangent(p) + + lface_to_pindex_to_perm = get_face_vertex_permutations(p, num_cell_dims(p)-1) + nlfaces = length(lface_to_t) + + # Repeat tangents according to the permutations, just like normals + lface_pindex_to_t = [ + fill(lface_to_t[lface], length(lface_to_pindex_to_perm[lface])) + for lface in 1:nlfaces + ] + return lface_pindex_to_t + end + + # 3) Map the function f over all reference polytope types + ctype_lface_pindex_to_tref = map(f, get_reffes(cell_grid)) + + # 4) Convert into a FaceCompressedVector + face_to_tref = FaceCompressedVector(ctype_lface_pindex_to_tref, boundary_trian_glue) + + # 5) Mark it as a constant field on each facet + face_s_tref = lazy_map(constant_field, face_to_tref) + + # 6) Geometric transformation from reference to real domain + cell_q_x = get_cell_map(cell_grid) + cell_q_Jt = lazy_map(∇, cell_q_x) + face_q_Jt = lazy_map(Reindex(cell_q_Jt), boundary_trian_glue.face_to_cell) + face_q_J = lazy_map(transpose, face_q_Jt) + + # 7) Get dimension, get glue + D = num_cell_dims(cell_grid) + boundary_trian_glue_D = get_glue(trian, Val(D)) + face_s_q = boundary_trian_glue_D.tface_to_mface_map + face_s_J = lazy_map(∘, face_q_J, face_s_q) + + # 8) Push the reference tangent to physical domain + face_s_t = lazy_map(Broadcasting(Operation(push_tangent)), face_s_J, face_s_tref) + + # 9) Memoize + return Fields.MemoArray(face_s_t) +end + +function get_edge_tangent(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A} + 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 + + +# Identitcal to push_normal, just for emphasis that we are pushing forward the +# Jacobian and not the inverse transpose +function push_tangent(J, t) + return push_normal(J, t) +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 index 471b420ff..d7d34ae53 100644 --- 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 index 5313e040b..936834956 100644 --- 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/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl index 4c6b956a3..5a0f30018 100644 --- a/src/Geometry/Triangulations.jl +++ b/src/Geometry/Triangulations.jl @@ -41,6 +41,7 @@ get_cell_node_ids(trian::Triangulation) = get_cell_node_ids(get_grid(trian)) 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 diff --git a/test/GeometryTests/SkeletonTriangulationsTests.jl b/test/GeometryTests/SkeletonTriangulationsTests.jl index d4b6bd793..29e0379b8 100644 --- 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 @@ -45,6 +46,7 @@ vglue = get_glue(vtrian,Val(3)) @test vglue.plus.tface_to_mface == sglue.plus.tface_to_mface[ids] @test vglue.minus.tface_to_mface == sglue.minus.tface_to_mface[ids] vn = get_facet_normal(vtrian) +vt = Gridap.Geometry.get_edge_tangent(strian) @test isa(vn,SkeletonPair) @test isa(vn.plus,AbstractArray) @test isa(vn.minus,AbstractArray) @@ -52,6 +54,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)) From fa50c430ba4b4fd52613e8b6910482a7c571a582 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Wed, 1 Jan 2025 17:22:39 -0700 Subject: [PATCH 2/6] Visually inspecting tangents looks correct. Need to add tests --- src/CellData/CellData.jl | 1 + src/CellData/CellFields.jl | 22 +++++++++++----- src/Exports.jl | 1 + src/Geometry/BoundaryTriangulations.jl | 25 +++++-------------- .../SkeletonTriangulationsTests.jl | 2 +- 5 files changed, 25 insertions(+), 26 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index d96ecca48..19efed72e 100644 --- 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 index b60474a08..5d3459d61 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -112,20 +112,22 @@ function CellField(f,trian::Triangulation) end function get_normal_vector(trian::Triangulation) - println("Calling get_normal_vector") cell_normal = get_facet_normal(trian) - get_vector_skeleton_pair(trian,cell_normal, get_normal_vector) + get_normal_vector(trian, cell_normal) end function get_tangent_vector(trian::Triangulation) cell_tangent = get_edge_tangent(trian) - get_vector_skeleton_pair(trian,cell_tangent, get_edge_tangent) + 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) @@ -713,9 +715,17 @@ function CellFieldAt{T}(parent::OperationCellField) where T OperationCellField(parent.op,args...) end -function get_vector_skeleton_pair(trian::Triangulation,cell_normal::SkeletonPair, f::Function) - plus = f(trian,cell_normal.plus) - minus = f(trian,cell_normal.minus) +function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair) + get_normal_or_tangent_vector(trian, cell_normal) +end + +function get_tangent_vector(trian::Triangulation,cell_tangent::SkeletonPair) + get_normal_or_tangent_vector(trian, cell_tangent) +end + +function get_normal_or_tangent_vector(trian::Triangulation, cell::SkeletonPair) + plus = get_normal_vector(trian,cell.plus) + minus = get_normal_vector(trian,cell.minus) SkeletonPair(plus,minus) end diff --git a/src/Exports.jl b/src/Exports.jl index 27dbe467d..345576e8a 100644 --- 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 index 141593d5c..36aceff34 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -239,17 +239,14 @@ function get_facet_normal(trian::BoundaryTriangulation, boundary_trian_glue::Fac end function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue) - # 1) Retrieve the cell grid from the boundary triangulation cell_grid = get_grid(get_background_model(trian.trian)) - # 2) Define how we get reference tangents for each local face + ## Reference tangent function f(r) p = get_polytope(r) lface_to_t = get_edge_tangent(p) - lface_to_pindex_to_perm = get_face_vertex_permutations(p, num_cell_dims(p)-1) nlfaces = length(lface_to_t) - # Repeat tangents according to the permutations, just like normals lface_pindex_to_t = [ fill(lface_to_t[lface], length(lface_to_pindex_to_perm[lface])) @@ -257,32 +254,22 @@ function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::Fac ] return lface_pindex_to_t end - - # 3) Map the function f over all reference polytope types ctype_lface_pindex_to_tref = map(f, get_reffes(cell_grid)) - - # 4) Convert into a FaceCompressedVector face_to_tref = FaceCompressedVector(ctype_lface_pindex_to_tref, boundary_trian_glue) - - # 5) Mark it as a constant field on each facet face_s_tref = lazy_map(constant_field, face_to_tref) - # 6) Geometric transformation from reference to real domain + # Tangent vectors transform like the Jacobian cell_q_x = get_cell_map(cell_grid) cell_q_Jt = lazy_map(∇, cell_q_x) - face_q_Jt = lazy_map(Reindex(cell_q_Jt), boundary_trian_glue.face_to_cell) - face_q_J = lazy_map(transpose, face_q_Jt) - - # 7) Get dimension, get glue + cell_q_J = lazy_map(Operation(transpose), cell_q_Jt) + face_q_J = lazy_map(Reindex(cell_q_J), boundary_trian_glue.face_to_cell) + + # Change of domain D = num_cell_dims(cell_grid) boundary_trian_glue_D = get_glue(trian, Val(D)) face_s_q = boundary_trian_glue_D.tface_to_mface_map face_s_J = lazy_map(∘, face_q_J, face_s_q) - - # 8) Push the reference tangent to physical domain face_s_t = lazy_map(Broadcasting(Operation(push_tangent)), face_s_J, face_s_tref) - - # 9) Memoize return Fields.MemoArray(face_s_t) end diff --git a/test/GeometryTests/SkeletonTriangulationsTests.jl b/test/GeometryTests/SkeletonTriangulationsTests.jl index 29e0379b8..ff1beba65 100644 --- a/test/GeometryTests/SkeletonTriangulationsTests.jl +++ b/test/GeometryTests/SkeletonTriangulationsTests.jl @@ -46,7 +46,7 @@ vglue = get_glue(vtrian,Val(3)) @test vglue.plus.tface_to_mface == sglue.plus.tface_to_mface[ids] @test vglue.minus.tface_to_mface == sglue.minus.tface_to_mface[ids] vn = get_facet_normal(vtrian) -vt = Gridap.Geometry.get_edge_tangent(strian) +# vt = Gridap.Geometry.get_edge_tangent(strian) @test isa(vn,SkeletonPair) @test isa(vn.plus,AbstractArray) @test isa(vn.minus,AbstractArray) From 3962b629f8eb0790e2a0d429df532db5b82c20e8 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Sun, 5 Jan 2025 10:25:37 +0100 Subject: [PATCH 3/6] Using rotated normal in 2D and add tests --- src/CellData/CellFields.jl | 12 +++-- src/Geometry/BoundaryTriangulations.jl | 44 ++++--------------- src/Geometry/SkeletonTriangulations.jl | 3 +- .../BoundaryTriangulationsTests.jl | 15 +++++++ .../SkeletonTriangulationsTests.jl | 7 ++- 5 files changed, 37 insertions(+), 44 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 5d3459d61..b2f0add7f 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -716,16 +716,14 @@ function CellFieldAt{T}(parent::OperationCellField) where T end function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair) - get_normal_or_tangent_vector(trian, cell_normal) + plus = get_normal_vector(trian,cell_normal.plus) + minus = get_normal_vector(trian,cell_normal.minus) + SkeletonPair(plus,minus) end function get_tangent_vector(trian::Triangulation,cell_tangent::SkeletonPair) - get_normal_or_tangent_vector(trian, cell_tangent) -end - -function get_normal_or_tangent_vector(trian::Triangulation, cell::SkeletonPair) - plus = get_normal_vector(trian,cell.plus) - minus = get_normal_vector(trian,cell.minus) + plus = get_normal_vector(trian,cell_tangent.plus) + minus = get_normal_vector(trian,cell_tangent.minus) SkeletonPair(plus,minus) end diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl index 36aceff34..a756e37de 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -239,41 +239,22 @@ function get_facet_normal(trian::BoundaryTriangulation, boundary_trian_glue::Fac end function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue) - cell_grid = get_grid(get_background_model(trian.trian)) - ## Reference tangent - function f(r) - p = get_polytope(r) - lface_to_t = get_edge_tangent(p) - lface_to_pindex_to_perm = get_face_vertex_permutations(p, num_cell_dims(p)-1) - nlfaces = length(lface_to_t) - # Repeat tangents according to the permutations, just like normals - lface_pindex_to_t = [ - fill(lface_to_t[lface], length(lface_to_pindex_to_perm[lface])) - for lface in 1:nlfaces - ] - return lface_pindex_to_t + 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 - ctype_lface_pindex_to_tref = map(f, get_reffes(cell_grid)) - face_to_tref = FaceCompressedVector(ctype_lface_pindex_to_tref, boundary_trian_glue) - face_s_tref = lazy_map(constant_field, face_to_tref) - # Tangent vectors transform like the Jacobian - cell_q_x = get_cell_map(cell_grid) - cell_q_Jt = lazy_map(∇, cell_q_x) - cell_q_J = lazy_map(Operation(transpose), cell_q_Jt) - face_q_J = lazy_map(Reindex(cell_q_J), boundary_trian_glue.face_to_cell) - - # Change of domain - D = num_cell_dims(cell_grid) - boundary_trian_glue_D = get_glue(trian, Val(D)) - face_s_q = boundary_trian_glue_D.tface_to_mface_map - face_s_J = lazy_map(∘, face_q_J, face_s_q) - face_s_t = lazy_map(Broadcasting(Operation(push_tangent)), face_s_J, face_s_tref) + 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 @@ -283,13 +264,6 @@ function get_facet_normal(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue} get_facet_normal(trian, glue) end - -# Identitcal to push_normal, just for emphasis that we are pushing forward the -# Jacobian and not the inverse transpose -function push_tangent(J, t) - return push_normal(J, t) -end - function push_normal(invJt,n) v = invJt⋅n m = sqrt(inner(v,v)) diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl index 936834956..dc0c05515 100644 --- a/src/Geometry/SkeletonTriangulations.jl +++ b/src/Geometry/SkeletonTriangulations.jl @@ -96,7 +96,8 @@ end function get_edge_tangent(trian::SkeletonTriangulation) plus = get_edge_tangent(trian.plus) - minus = get_edge_tangent(trian.minus) + # Flip the minus side + minus = lazy_map(x -> -x, plus) SkeletonPair(plus,minus) end diff --git a/test/GeometryTests/BoundaryTriangulationsTests.jl b/test/GeometryTests/BoundaryTriangulationsTests.jl index b775798ef..2b47c286d 100644 --- 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 index ff1beba65..3aa73e507 100644 --- a/test/GeometryTests/SkeletonTriangulationsTests.jl +++ b/test/GeometryTests/SkeletonTriangulationsTests.jl @@ -46,7 +46,6 @@ vglue = get_glue(vtrian,Val(3)) @test vglue.plus.tface_to_mface == sglue.plus.tface_to_mface[ids] @test vglue.minus.tface_to_mface == sglue.minus.tface_to_mface[ids] vn = get_facet_normal(vtrian) -# vt = Gridap.Geometry.get_edge_tangent(strian) @test isa(vn,SkeletonPair) @test isa(vn.plus,AbstractArray) @test isa(vn.minus,AbstractArray) @@ -119,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]) From 2ad1695a83c1009dfb4d2fb994e60b67790f6c07 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Tue, 7 Jan 2025 10:19:00 +0100 Subject: [PATCH 4/6] Adress requested changes --- src/Geometry/BoundaryTriangulations.jl | 6 +----- src/Geometry/SkeletonTriangulations.jl | 3 +-- test/GeometryTests/SkeletonTriangulationsTests.jl | 8 ++++++++ 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl index a756e37de..26d783882 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -242,11 +242,7 @@ function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::Fac 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 + rotate_normal_2d(n) = VectorValue(n[2], -n[1]) face_s_t = lazy_map(Operation(rotate_normal_2d), face_s_n) diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl index dc0c05515..936834956 100644 --- a/src/Geometry/SkeletonTriangulations.jl +++ b/src/Geometry/SkeletonTriangulations.jl @@ -96,8 +96,7 @@ end function get_edge_tangent(trian::SkeletonTriangulation) plus = get_edge_tangent(trian.plus) - # Flip the minus side - minus = lazy_map(x -> -x, plus) + minus = get_edge_tangent(trian.minus) SkeletonPair(plus,minus) end diff --git a/test/GeometryTests/SkeletonTriangulationsTests.jl b/test/GeometryTests/SkeletonTriangulationsTests.jl index 3aa73e507..d2c971c73 100644 --- a/test/GeometryTests/SkeletonTriangulationsTests.jl +++ b/test/GeometryTests/SkeletonTriangulationsTests.jl @@ -112,6 +112,14 @@ 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) From ba6e68562eae150e9c2a1632e7273d4b87348703 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Tue, 7 Jan 2025 11:05:48 +0100 Subject: [PATCH 5/6] Remove Grid interface for get_edge_tangent --- src/Geometry/Triangulations.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl index 5a0f30018..4c6b956a3 100644 --- a/src/Geometry/Triangulations.jl +++ b/src/Geometry/Triangulations.jl @@ -41,7 +41,6 @@ get_cell_node_ids(trian::Triangulation) = get_cell_node_ids(get_grid(trian)) 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 From e7726f42919e820c8cf4d536c89ecdae9e031bf7 Mon Sep 17 00:00:00 2001 From: RAPPAPORT Ari Date: Sun, 12 Jan 2025 21:37:54 +0100 Subject: [PATCH 6/6] Update NEWS.md --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index 8ebb57948..507214249 100644 --- 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