diff --git a/NEWS.md b/NEWS.md index a2c41f59b..507214249 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,18 @@ 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). -## [Unreleased] +## [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 + + +### Fixed + +- BUG in `FineToCoarseFields.jl`. Since PR[#1074](https://github.com/gridap/Gridap.jl/pull/1074) ### Added diff --git a/Project.toml b/Project.toml index 7a5793270..b91eb271d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gridap" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" authors = ["Santiago Badia ", "Francesc Verdugo ", "Alberto F. Martin "] -version = "0.18.8" +version = "0.18.9" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/src/Adaptivity/FineToCoarseFields.jl b/src/Adaptivity/FineToCoarseFields.jl index 648b88600..9be4a063b 100644 --- a/src/Adaptivity/FineToCoarseFields.jl +++ b/src/Adaptivity/FineToCoarseFields.jl @@ -108,7 +108,7 @@ function Arrays.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point}) field_id = id_map[child_id] if field_id != 0 fi = getindex!(fi_cache,fields,field_id) - mi = getindex!(mi_cache,cmaps,field_id) + mi = getindex!(mi_cache,cmaps,child_id) zi = Fields.evaluate!(zi_cache,mi,xi) y[i] = Fields.evaluate!(yi_cache,fi,zi) end 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 274ae27c2..f00c65389 100644 --- 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 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 fc8566acc..a18a36902 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -237,6 +237,24 @@ function get_facet_normal(trian::BoundaryTriangulation, face_glue::FaceToCellGlu 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 index 19a4fd017..e85618a25 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -73,6 +73,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 @@ -110,6 +111,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/test/AdaptivityTests/FineToCoarseFieldsTests.jl b/test/AdaptivityTests/FineToCoarseFieldsTests.jl index 89beb196b..2e51b6948 100644 --- a/test/AdaptivityTests/FineToCoarseFieldsTests.jl +++ b/test/AdaptivityTests/FineToCoarseFieldsTests.jl @@ -100,4 +100,91 @@ xHh = change_domain(xH,get_triangulation(modelh),ReferenceDomain()) evaluate(Gridap.CellData.get_data(xHh)[1],[Point(0.0,0.0),Point(0.5,0.5)]) evaluate(Gridap.CellData.get_data(xHh)[1],Point(0.5,0.5)) +function test_PR_1074() + + function setup_coarse_discrete_model() + ptr = [ 1, 5 ] + data = [ 1,2,3,4 ] + cell_vertex_lids = Gridap.Arrays.Table(data,ptr) + node_coordinates = Vector{Point{2,Float64}}(undef,4) + node_coordinates[1]=Point{2,Float64}(0.0,0.0) + node_coordinates[2]=Point{2,Float64}(1.0,0.0) + node_coordinates[3]=Point{2,Float64}(0.0,1.0) + node_coordinates[4]=Point{2,Float64}(1.0,1.0) + + polytope=QUAD + scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1) + cell_types=collect(Fill(1,length(cell_vertex_lids))) + cell_reffes=[scalar_reffe] + grid = Gridap.Geometry.UnstructuredGrid(node_coordinates, + cell_vertex_lids, + cell_reffes, + cell_types, + Gridap.Geometry.NonOriented()) + Gridap.Geometry.UnstructuredDiscreteModel(grid) + end + + function setup_adapted_discrete_model(parent) + ptr = [ 1, 5, 9 ] + data = [ 1,2,3,4, 2,5,4,6 ] + cell_vertex_lids = Gridap.Arrays.Table(data,ptr) + node_coordinates = Vector{Point{2,Float64}}(undef,6) + node_coordinates[1]=Point{2,Float64}(0.0,0.5) + node_coordinates[2]=Point{2,Float64}(0.5,0.5) + node_coordinates[3]=Point{2,Float64}(0.0,1.0) + node_coordinates[4]=Point{2,Float64}(0.5,1.0) + node_coordinates[5]=Point{2,Float64}(1.0,0.5) + node_coordinates[6]=Point{2,Float64}(1.0,1.0) + + polytope=QUAD + scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1) + cell_types=collect(Fill(1,length(cell_vertex_lids))) + cell_reffes=[scalar_reffe] + grid = Gridap.Geometry.UnstructuredGrid(node_coordinates, + cell_vertex_lids, + cell_reffes, + cell_types, + Gridap.Geometry.NonOriented()) + model=Gridap.Geometry.UnstructuredDiscreteModel(grid) + + n2o_faces_map=Vector{Vector{Int64}}(undef,3) + n2o_faces_map[3]=[1,1] + n2o_cell_to_child_id=[2,3] + + ref_rules = [Gridap.Adaptivity.RefinementRule(Gridap.ReferenceFEs.LagrangianRefFE(Float64,QUAD,1),2)] + + glue=Gridap.Adaptivity.AdaptivityGlue(Gridap.Adaptivity.RefinementGlue(), + n2o_faces_map, + n2o_cell_to_child_id, + ref_rules) + Gridap.Adaptivity.AdaptedDiscreteModel(model,parent,glue) + end + + coarse_model=setup_coarse_discrete_model() + fine_model=setup_adapted_discrete_model(coarse_model) + + order=0 + + VH = FESpace(coarse_model, + ReferenceFE(raviart_thomas,Float64,order), + conformity=:Hdiv) + + UH = TrialFESpace(VH) + + Vh = FESpace(fine_model, + ReferenceFE(raviart_thomas,Float64,order), + conformity=:Hdiv) + + Uh = TrialFESpace(Vh) + + u(x) = VectorValue(x[1],x[2]) + uh = interpolate(u,Uh) + uH = interpolate(u,UH) + uhH = interpolate(uh,UH) + + @test get_free_dof_values(uhH)[2] ≈ 1.0 +end + +test_PR_1074() + end \ No newline at end of file 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 d4b6bd793..d2c971c73 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 @@ -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])