From 5bb7d4af9fc13766f3deebae48b080daa776e4ba Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 7 Jun 2024 17:14:20 +1000 Subject: [PATCH 01/14] DistributedDiscreteGeometry implementation --- src/Distributed/Distributed.jl | 3 + .../DistributedDiscreteGeometries.jl | 102 ++++++++++++ src/Distributed/DistributedDiscretizations.jl | 8 +- .../DistributedDiscreteGeometryPoissonTest.jl | 148 ++++++++++++++++++ .../DistributedDiscreteGeometryPoissonTest.jl | 8 + test/DistributedTests/sequential/runtests.jl | 1 + 6 files changed, 268 insertions(+), 2 deletions(-) create mode 100644 src/Distributed/DistributedDiscreteGeometries.jl create mode 100644 test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index aadb4d5..4fc3e95 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -43,6 +43,7 @@ import GridapEmbedded.Interfaces: EmbeddedBoundary import GridapEmbedded.Interfaces: compute_bgfacet_to_inoutcut import GridapEmbedded.Interfaces: compute_bgcell_to_inoutcut import GridapEmbedded.CSG: get_geometry +import GridapEmbedded.LevelSetCutters: discretize import Gridap.Geometry: Triangulation import Gridap.Geometry: SkeletonTriangulation import Gridap.Geometry: BoundaryTriangulation @@ -52,6 +53,8 @@ import GridapDistributed: remove_ghost_cells include("DistributedDiscretizations.jl") +include("DistributedDiscreteGeometries.jl") + include("DistributedAgFEM.jl") include("DistributedQuadratures.jl") diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl new file mode 100644 index 0000000..ce3c8be --- /dev/null +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -0,0 +1,102 @@ +struct DistributedDiscreteGeometry{A} <: GridapType + geometries::A +end + +local_views(a::DistributedDiscreteGeometry) = a.geometries + +function DistributedDiscreteGeometry(φh::CellField,model::DistributedDiscreteModel) + gids = get_cell_gids(model) + geometries = map(local_views(model),local_views(gids),local_views(φh)) do model,gids,φh + ownmodel = remove_ghost_cells(model,gids) + point_to_coords = collect1d(get_node_coordinates(ownmodel)) + DiscreteGeometry(φh(point_to_coords),point_to_coords) + end + DistributedDiscreteGeometry(geometries) +end + +function get_geometry(a::AbstractArray{<:DiscreteGeometry}) + DistributedDiscreteGeometry(a) +end + +function discretize(a::AnalyticalGeometry,model::DistributedDiscreteModel) + gids = get_cell_gids(model) + geometries = map(local_views(model),local_views(gids)) do model,gids + ownmodel = remove_ghost_cells(model,gids) + point_to_coords = collect1d(get_node_coordinates(ownmodel)) + discretize(a,point_to_coords) + end + DistributedDiscreteGeometry(geometries) +end + +function cut(cutter::Cutter,bgmodel::DistributedDiscreteModel,geom::DistributedDiscreteGeometry) + gids = get_cell_gids(bgmodel) + cuts = map(local_views(bgmodel),local_views(gids),local_views(geom)) do bgmodel,gids,geom + ownmodel = remove_ghost_cells(bgmodel,gids) + cutgeo = cut(cutter,ownmodel,geom) + change_bgmodel(cutgeo,bgmodel,own_to_local(gids)) + end + consistent_bgcell_to_inoutcut!(cuts,gids) + DistributedEmbeddedDiscretization(cuts,bgmodel) +end + +function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel,geom::DistributedDiscreteGeometry) + D = map(num_dims,local_views(bgmodel)) |> PartitionedArrays.getany + cell_gids = get_cell_gids(bgmodel) + facet_gids = get_face_gids(bgmodel,D-1) + cuts = map( + local_views(bgmodel), + local_views(cell_gids), + local_views(facet_gids), + local_views(geom)) do bgmodel,cell_gids,facet_gids,geom + ownmodel = remove_ghost_cells(bgmodel,cell_gids) + facet_to_pfacet = get_face_to_parent_face(ownmodel,D-1) + cutfacets = cut_facets(cutter,ownmodel,geom) + cutfacets = change_bgmodel(cutfacets,bgmodel,facet_to_pfacet) + remove_ghost_subfacets(cutfacets,facet_gids) + end + consistent_bgfacet_to_inoutcut!(cuts,facet_gids) + DistributedEmbeddedDiscretization(cuts,bgmodel) +end + +function distributed_embedded_triangulation( + T, + cutgeo::DistributedEmbeddedDiscretization, + cutinorout, + geom::DistributedDiscreteGeometry) + + trians = map(local_views(cutgeo),local_views(geom)) do lcutgeo,lgeom + T(lcutgeo,cutinorout,lgeom) + end + bgmodel = get_background_model(cutgeo) + DistributedTriangulation(trians,bgmodel) +end + +function distributed_aggregate( + strategy::AggregateCutCellsByThreshold, + cut::DistributedEmbeddedDiscretization, + geo::DistributedDiscreteGeometry, + in_or_out=IN) + + bgmodel = get_background_model(cut) + facet_to_inoutcut = compute_bgfacet_to_inoutcut(bgmodel,geo) + _distributed_aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut) +end + +function compute_bgcell_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,geo::DistributedDiscreteGeometry) + map(local_views(cutgeo),local_views(geo)) do cutgeo,geo + compute_bgcell_to_inoutcut(cutgeo,geo) + end +end + +function compute_bgfacet_to_inoutcut( + cutter::Cutter, + bgmodel::DistributedDiscreteModel, + geo::DistributedDiscreteGeometry) + + gids = get_cell_gids(bgmodel) + bgf_to_ioc = map(local_views(bgmodel),local_views(gids),local_views(geo)) do model,gids,geo + ownmodel = remove_ghost_cells(model,gids) + compute_bgfacet_to_inoutcut(cutter,ownmodel,geo) + end + compute_bgfacet_to_inoutcut(bgmodel,bgf_to_ioc) +end \ No newline at end of file diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl index 850eb18..742a233 100644 --- a/src/Distributed/DistributedDiscretizations.jl +++ b/src/Distributed/DistributedDiscretizations.jl @@ -16,8 +16,12 @@ local_views(a::DistributedEmbeddedDiscretization) = a.discretizations get_background_model(a::DistributedEmbeddedDiscretization) = a.model function get_geometry(a::DistributedEmbeddedDiscretization) - cut = local_views(a) |> PartitionedArrays.getany - get_geometry(cut) + loc_geometries = map(get_geometry,local_views(a)) + get_geometry(loc_geometries) +end + +function get_geometry(a::AbstractArray{<:CSG.Geometry}) + PartitionedArrays.getany(a) end function cut(bgmodel::DistributedDiscreteModel,args...) diff --git a/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..ffe8f67 --- /dev/null +++ b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,148 @@ +module DistributedDiscreteGeometryPoissonTest + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapEmbedded.CSG +using GridapEmbedded.LevelSetCutters + +function main(distribute,parts; + threshold=1, + n=8, + cells=(n,n), + geometry=:circle) + + ranks = distribute(LinearIndices((prod(parts),))) + + u(x) = x[1] - x[2] + f(x) = -Δ(u)(x) + ud(x) = u(x) + + geometries = Dict( + :circle => circle_geometry, + :remotes => remotes_geometry, + ) + + bgmodel,_geo = geometries[geometry](ranks,parts,cells) + geo = discretize(_geo,bgmodel) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_bg = Triangulation(bgmodel) + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 1 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + reffe = ReferenceFE(lagrangian,Float64,order) + + Vstd = FESpace(Ω_act,reffe) + + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V) + + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +function circle_geometry(ranks,parts,cells) + L = 1 + p0 = Point(0.0,0.0) + pmin = p0-L/2 + pmax = p0+L/2 + R = 0.35 + geo = disk(R,x0=p0) + bgmodel = CartesianDiscreteModel(ranks,parts,pmin,pmax,cells) + bgmodel,geo +end + +function remotes_geometry(ranks,parts,cells) + x0 = Point(0.05,0.05) + d1 = VectorValue(0.9,0.0) + d2 = VectorValue(0.0,0.1) + geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) + + x0 = Point(0.15,0.1) + d1 = VectorValue(0.25,0.0) + d2 = VectorValue(0.0,0.6) + geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) + geo = union(geo1,geo2) + + domain = (0, 1, 0, 1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells) + bgmodel,geo +end + +with_debug() do distribute + main(distribute,(2,2)) + main(distribute,(4,1),cells=(12,12),geometry=:remotes) +end + +end # module \ No newline at end of file diff --git a/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..9e29dee --- /dev/null +++ b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,8 @@ +module DistributedDiscreteGeometryPoissonTestsSeq +using PartitionedArrays +include("../DistributedDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + PoissonTests.main(distribute,(2,2)) + PoissonTests.main(distribute,(4,1),cells=(12,12),geometry=:remotes) +end +end diff --git a/test/DistributedTests/sequential/runtests.jl b/test/DistributedTests/sequential/runtests.jl index e68831c..5819b69 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -3,5 +3,6 @@ module SequentialTests using Test @time @testset "PoissonSeq" begin include("PoissonTests.jl") end +@time @testset "PoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end end From 1ac42f2b57551c94c9ae2132fa3a0bccfe07b8e2 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 7 Jun 2024 17:24:37 +1000 Subject: [PATCH 02/14] Fix test --- .../sequential/DistributedDiscreteGeometryPoissonTest.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl index 9e29dee..080d4b9 100644 --- a/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl +++ b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl @@ -2,7 +2,7 @@ module DistributedDiscreteGeometryPoissonTestsSeq using PartitionedArrays include("../DistributedDiscreteGeometryPoissonTest.jl") with_debug() do distribute - PoissonTests.main(distribute,(2,2)) - PoissonTests.main(distribute,(4,1),cells=(12,12),geometry=:remotes) + DistributedDiscreteGeometryPoissonTest.main(distribute,(2,2)) + DistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(12,12),geometry=:remotes) end end From 643efc205a13c8a7a9b6a3e0352d19cb2f6cac3d Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Sat, 8 Jun 2024 10:41:18 +1000 Subject: [PATCH 03/14] Missing MPI test call --- test/DistributedTests/mpi/runtests_body.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl index fe600b7..7c44d8b 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -6,6 +6,7 @@ using MPI include("../PoissonTests.jl") include("../AggregatesTests.jl") +include("../DistributedDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -21,6 +22,11 @@ function all_tests(distribute,parts) PoissonTests.main(distribute,(prod(parts),1),cells=(12,12),geometry=:remotes) PArrays.toc!(t,"Poisson") + PArrays.tic!(t) + DistributedDiscreteGeometryPoissonTest.main(distribute,parts) + DistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(12,12),geometry=:remotes) + PArrays.toc!(t,"DistributedDiscreteGeometryPoisson") + if prod(parts) == 4 DistributedAggregatesTests.main(distribute,parts) end From 2c1d835b91ac94ffe14f4b5cb814800869310c88 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Mon, 10 Jun 2024 13:21:12 +1000 Subject: [PATCH 04/14] Include name in DistributedDiscreteGeometry --- src/Distributed/DistributedDiscreteGeometries.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index ce3c8be..1f3d32b 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -4,12 +4,12 @@ end local_views(a::DistributedDiscreteGeometry) = a.geometries -function DistributedDiscreteGeometry(φh::CellField,model::DistributedDiscreteModel) +function DistributedDiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::String="") gids = get_cell_gids(model) geometries = map(local_views(model),local_views(gids),local_views(φh)) do model,gids,φh ownmodel = remove_ghost_cells(model,gids) point_to_coords = collect1d(get_node_coordinates(ownmodel)) - DiscreteGeometry(φh(point_to_coords),point_to_coords) + DiscreteGeometry(φh(point_to_coords),point_to_coords;name) end DistributedDiscreteGeometry(geometries) end From 1b5a49014db5568b088a1c0a2a32961c68b27a2a Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 19 Jun 2024 00:59:43 +1000 Subject: [PATCH 05/14] Added _get_value_at_coords to avoid interpolation --- src/Distributed/Distributed.jl | 2 +- .../DistributedDiscreteGeometries.jl | 49 +++++++++++++++++-- src/LevelSetCutters/DiscreteGeometries.jl | 32 ++++++++++++ 3 files changed, 79 insertions(+), 4 deletions(-) diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index 4fc3e95..a08527e 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -43,7 +43,7 @@ import GridapEmbedded.Interfaces: EmbeddedBoundary import GridapEmbedded.Interfaces: compute_bgfacet_to_inoutcut import GridapEmbedded.Interfaces: compute_bgcell_to_inoutcut import GridapEmbedded.CSG: get_geometry -import GridapEmbedded.LevelSetCutters: discretize +import GridapEmbedded.LevelSetCutters: discretize, DiscreteGeometry import Gridap.Geometry: Triangulation import Gridap.Geometry: SkeletonTriangulation import Gridap.Geometry: BoundaryTriangulation diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 1f3d32b..616c712 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -4,12 +4,55 @@ end local_views(a::DistributedDiscreteGeometry) = a.geometries -function DistributedDiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::String="") +function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) where {Dc,Dp} + @assert DomainStyle(φh) == ReferenceDomain() gids = get_cell_gids(model) - geometries = map(local_views(model),local_views(gids),local_views(φh)) do model,gids,φh + values = map(local_views(φh),local_views(model),local_views(gids)) do φh, model, gids + # Maps from the no-ghost model to the original model + own_model = remove_ghost_cells(model,gids) + own_to_local_node = Geometry.get_face_to_parent_face(own_model,0) + local_to_own_node = Arrays.find_inverse_index_map(own_to_local_node,num_nodes(model)) + own_to_local_cell = Geometry.get_face_to_parent_face(own_model,Dc) + + # Cell-to-node map for the original model + # topo = get_grid_topology(model) + # c2n_map = get_faces(topo,Dc,0) + c2n_map = collect1d(get_cell_node_ids(model)) + + # Cell-wise node coordinates (in ReferenceDomain coordinates) + cell_reffe = get_cell_reffe(model) + cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) + + φh_data = CellData.get_data(φh) + space = get_fe_space(φh) + T = get_dof_value_type(space) + values = Vector{T}(undef,num_nodes(own_model)) + touched = fill(false,num_nodes(model)) + + cell_node_coords_cache = array_cache(cell_node_coords) + for cell in own_to_local_cell # For each owned cell + field = φh_data[cell] + node_coords = getindex!(cell_node_coords_cache,cell_node_coords,cell) + for (iN,node) in enumerate(c2n_map[cell]) # Go over local nodes + own_node = local_to_own_node[node] + if (own_node != 0) && !touched[node] # Compute value if suitable + values[own_node] = field(node_coords[iN]) + touched[node] = true + end + end + end + return values + end + return values +end + +function DiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::String="") + φ_values = _get_values_at_owned_coords(φh,model) + gids = get_cell_gids(model) + geometries = map(local_views(model),local_views(gids),local_views(φ_values)) do model,gids,loc_φ ownmodel = remove_ghost_cells(model,gids) point_to_coords = collect1d(get_node_coordinates(ownmodel)) - DiscreteGeometry(φh(point_to_coords),point_to_coords;name) + DiscreteGeometry(loc_φ,point_to_coords;name) end DistributedDiscreteGeometry(geometries) end diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index 5289299..d2e15ff 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -65,9 +65,41 @@ function _find_unique_leaves(tree) j_to_fun, oid_to_j end +function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where {Dc,Dp} + @assert DomainStyle(φh) == ReferenceDomain() + # Cell-to-node map for the original model + c2n_map = collect1d(get_cell_node_ids(model)) + + # Cell-wise node coordinates (in ReferenceDomain coordinates) + cell_reffe = get_cell_reffe(model) + cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) + + # Get cell data + φh_data = CellData.get_data(φh) + space = get_fe_space(φh) + T = get_dof_value_type(space) + values = Vector{T}(undef,num_nodes(model)) + cell_node_coords_cache = array_cache(cell_node_coords) + # Loop over cells + for cell in eachindex(c2n_map) + field = φh_data[cell] + node_coords = getindex!(cell_node_coords_cache,cell_node_coords,cell) + for (iN,node) in enumerate(c2n_map[cell]) + values[node] = field(node_coords[iN]) + end + end +end + function DiscreteGeometry( point_to_value::AbstractVector,point_to_coords::AbstractVector;name::String="") data = (point_to_value,name,nothing) tree = Leaf(data) DiscreteGeometry(tree,point_to_coords) end + +function DiscreteGeometry( + φh::CellField,model::DiscreteModel;name::String="") + point_to_value = _get_value_at_coords(φh,model) + point_to_coords = collect1d(get_node_coordinates(model)) + DiscreteGeometry(point_to_value,point_to_coords;name) +end \ No newline at end of file From e80c763c717ee28cdd8b215c7a0e191da300f1db Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 19 Jun 2024 15:12:25 +1000 Subject: [PATCH 06/14] Remove unused test --- .../DistributedDiscreteGeometryPoissonTest.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl index ffe8f67..10a9959 100644 --- a/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl +++ b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl @@ -140,9 +140,4 @@ function remotes_geometry(ranks,parts,cells) bgmodel,geo end -with_debug() do distribute - main(distribute,(2,2)) - main(distribute,(4,1),cells=(12,12),geometry=:remotes) -end - end # module \ No newline at end of file From 5195d51dcf189b3e401051114cbf80c6c26a93da Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 19 Jun 2024 16:50:21 +1000 Subject: [PATCH 07/14] added periodic testing, suggests that weak method is likely required here. --- src/LevelSetCutters/DiscreteGeometries.jl | 3 +- src/LevelSetCutters/LevelSetCutters.jl | 1 + test/AgFEMTests/PeriodicAgFEMSpacesTests.jl | 75 +++++++++ test/AgFEMTests/runtests.jl | 2 + ...cDistributedDiscreteGeometryPoissonTest.jl | 143 ++++++++++++++++++ ...istributedLSDiscreteGeometryPoissonTest.jl | 110 ++++++++++++++ test/DistributedTests/PeriodicPoissonTests.jl | 141 +++++++++++++++++ test/DistributedTests/mpi/runtests_body.jl | 17 +++ ...cDistributedDiscreteGeometryPoissonTest.jl | 8 + ...istributedLSDiscreteGeometryPoissonTest.jl | 7 + .../sequential/PeriodicPoissonTests.jl | 8 + test/DistributedTests/sequential/runtests.jl | 4 +- .../PeriodicDiscreteGeoPoissonAgFEMTests.jl | 76 ++++++++++ .../PeriodicPoissonAgFEMTests.jl | 74 +++++++++ test/GridapEmbeddedTests/runtests.jl | 4 + 15 files changed, 671 insertions(+), 2 deletions(-) create mode 100644 test/AgFEMTests/PeriodicAgFEMSpacesTests.jl create mode 100644 test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/PeriodicPoissonTests.jl create mode 100644 test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/sequential/PeriodicPoissonTests.jl create mode 100644 test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl create mode 100644 test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index d2e15ff..3eb9705 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -76,7 +76,7 @@ function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where # Get cell data φh_data = CellData.get_data(φh) - space = get_fe_space(φh) + space = FESpaces.get_fe_space(φh) T = get_dof_value_type(space) values = Vector{T}(undef,num_nodes(model)) cell_node_coords_cache = array_cache(cell_node_coords) @@ -88,6 +88,7 @@ function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where values[node] = field(node_coords[iN]) end end + return values end function DiscreteGeometry( diff --git a/src/LevelSetCutters/LevelSetCutters.jl b/src/LevelSetCutters/LevelSetCutters.jl index aa5b147..75296ea 100644 --- a/src/LevelSetCutters/LevelSetCutters.jl +++ b/src/LevelSetCutters/LevelSetCutters.jl @@ -28,6 +28,7 @@ using Gridap.Geometry using Gridap.CellData using Gridap.Polynomials using Gridap.Visualization +using Gridap.FESpaces export LevelSetCutter export AnalyticalGeometry diff --git a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl new file mode 100644 index 0000000..cc1845f --- /dev/null +++ b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl @@ -0,0 +1,75 @@ +module PeriodicAgFEMSpacesTests + +using Test +using Gridap +using GridapEmbedded +using Gridap.Geometry: get_active_model + +const R = 0.5 +geom = disk(R,x0=Point(0.5,0.5)) +n = 21 +partition = (n,n) + +domain = (0,1,0,1) +bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) + +cutdisc = cut(bgmodel,geom) + +strategy = AggregateCutCellsByThreshold(1) +aggregates = aggregate(strategy,cutdisc) + +Ω_bg = Triangulation(bgmodel) +Ω_ac = Triangulation(cutdisc,ACTIVE) +Ω = Triangulation(cutdisc,PHYSICAL) +Ω_in = Triangulation(cutdisc,IN) + +dΩ_bg = Measure(Ω_bg,2) +dΩ = Measure(Ω,2) +dΩ_in = Measure(Ω_in,2) + +model = get_active_model(Ω_ac) +order = 2 + +# In the physical domain +cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order) +Vstd = FESpace(Ω_ac,cell_fe) + +Vagg = AgFEMSpace(Vstd,aggregates) +U = TrialFESpace(Vagg) + +v(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 +vhagg = interpolate(v,Vagg) + +writevtk(Ω_ac,"test",cellfields=["v"=>vhagg]) + +tol = 10e-9 +@test sum( ∫(abs2(v-vhagg))dΩ ) < tol +@test sum( ∫(abs2(v-vhagg))dΩ_in ) < tol + +vh = FEFunction(V,rand(num_free_dofs(V))) +vhagg = interpolate(vh,Vagg) +@test sum( ∫(abs2(vh-vhagg))dΩ_in ) < tol + +# In the reference space + +reffe = ReferenceFE(lagrangian,Float64,order) +V = FESpace(Ω_ac,reffe) +Vagg = AgFEMSpace(V,aggregates) + +v(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 +vhagg = interpolate(v,Vagg) + +tol = 10e-9 +@test sum( ∫(abs2(v-vhagg))dΩ ) < tol +@test sum( ∫(abs2(v-vhagg))dΩ_in ) < tol + +vh = FEFunction(V,rand(num_free_dofs(V))) +vhagg = interpolate(vh,Vagg) +@test sum( ∫(abs2(vh-vhagg))dΩ_in ) < tol + +#cellfields = ["vh"=>vh,"vhagg"=>vhagg,"e"=>vh-vhagg] +#writevtk(Ω_bg,"trian_bg",nsubcells=10,cellfields=cellfields) +#writevtk(Ω_in,"trian_in",nsubcells=10,cellfields=cellfields) +#writevtk(Ω,"trian_phys",cellfields=cellfields) + +end # module diff --git a/test/AgFEMTests/runtests.jl b/test/AgFEMTests/runtests.jl index 9d03917..66d7f07 100644 --- a/test/AgFEMTests/runtests.jl +++ b/test/AgFEMTests/runtests.jl @@ -6,4 +6,6 @@ using Test @testset "AgFEMSpaces" begin include("AgFEMSpacesTests.jl") end +@testset "PeriodicAgFEMSpaces" begin include("PeriodicAgFEMSpacesTests.jl") end + end # module diff --git a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..8bd3e67 --- /dev/null +++ b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,143 @@ +module PeriodicDistributedDiscreteGeometryPoissonTest + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapEmbedded.CSG +using GridapEmbedded.LevelSetCutters + +function main(distribute,parts; + threshold=1, + n=8, + cells=(n,n), + geometry=:circle) + + ranks = distribute(LinearIndices((prod(parts),))) + + u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 + f(x) = -Δ(u)(x) + ud(x) = u(x) + + geometries = Dict( + :circle => circle_geometry, + :remotes => remotes_geometry, + ) + + bgmodel,_geo = geometries[geometry](ranks,parts,cells) + geo = discretize(_geo,bgmodel) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_bg = Triangulation(bgmodel) + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 2 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + reffe = ReferenceFE(lagrangian,Float64,order) + + Vstd = FESpace(Ω_act,reffe) + + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V) + + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +function circle_geometry(ranks,parts,cells) + L = 1 + p0 = Point(0.5,0.5) + pmin = p0-L/2 + pmax = p0+L/2 + R = 0.55 + geo = disk(R,x0=p0) + bgmodel = CartesianDiscreteModel(ranks,parts,pmin,pmax,cells;isperiodic=(true,true)) + bgmodel,geo +end + +function remotes_geometry(ranks,parts,cells) + x0 = Point(0.0,0.4) + d1 = VectorValue(1.0,0.0) + d2 = VectorValue(0.0,0.2) + geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) + + x0 = Point(0.4,0.0) + d1 = VectorValue(0.2,0.0) + d2 = VectorValue(0.0,1.0) + geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) + geo = union(geo1,geo2) + + domain = (0, 1, 0, 1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) + bgmodel,geo +end + +end # module \ No newline at end of file diff --git a/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..f0b659e --- /dev/null +++ b/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl @@ -0,0 +1,110 @@ +module PeriodicDistributedLSDiscreteGeometryPoissonTest + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapEmbedded.CSG +using GridapEmbedded.LevelSetCutters + +function main(distribute,parts; + threshold=1, + n=21, + cells=(n,n)) + + order = 2 + ranks = distribute(LinearIndices((prod(parts),))) + domain = (0,1,0,1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) + + u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 + f(x) = -Δ(u)(x) + ud(x) = u(x) + + reffe = ReferenceFE(lagrangian,Float64,order) + Ω_bg = Triangulation(bgmodel) + V_bg = FESpace(Ω_bg,reffe) + φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.55,V_bg) + geo = DiscreteGeometry(φh,bgmodel) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 2 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + Vstd = FESpace(Ω_act,reffe) + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V) + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +end # module \ No newline at end of file diff --git a/test/DistributedTests/PeriodicPoissonTests.jl b/test/DistributedTests/PeriodicPoissonTests.jl new file mode 100644 index 0000000..130ad41 --- /dev/null +++ b/test/DistributedTests/PeriodicPoissonTests.jl @@ -0,0 +1,141 @@ +module PeriodicPoissonTests + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapEmbedded.CSG + +function main(distribute,parts; + threshold=1, + n=8, + cells=(n,n), + geometry=:circle) + + ranks = distribute(LinearIndices((prod(parts),))) + + u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 + f(x) = -Δ(u)(x) + ud(x) = u(x) + + geometries = Dict( + :circle => circle_geometry, + :remotes => remotes_geometry, + ) + + bgmodel,geo = geometries[geometry](ranks,parts,cells) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_bg = Triangulation(bgmodel) + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 2 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + reffe = ReferenceFE(lagrangian,Float64,order) + + Vstd = FESpace(Ω_act,reffe) + + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V) + + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +function circle_geometry(ranks,parts,cells) + L = 1 + p0 = Point(0.5,0.5) + pmin = p0-L/2 + pmax = p0+L/2 + R = 0.55 + geo = disk(R,x0=p0) + bgmodel = CartesianDiscreteModel(ranks,parts,pmin,pmax,cells;isperiodic=(true,true)) + bgmodel,geo +end + +function remotes_geometry(ranks,parts,cells) + x0 = Point(0.0,0.4) + d1 = VectorValue(1.0,0.0) + d2 = VectorValue(0.0,0.2) + geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) + + x0 = Point(0.4,0.0) + d1 = VectorValue(0.2,0.0) + d2 = VectorValue(0.0,1.0) + geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) + geo = union(geo1,geo2) + + domain = (0, 1, 0, 1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) + bgmodel,geo +end + +end # module diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl index 7c44d8b..cdbb6dc 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -7,6 +7,9 @@ using MPI include("../PoissonTests.jl") include("../AggregatesTests.jl") include("../DistributedDiscreteGeometryPoissonTest.jl") +include("../PeriodicPoissonTests.jl") +include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") +include("../PeriodicDistributedLSDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -27,6 +30,20 @@ function all_tests(distribute,parts) DistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(12,12),geometry=:remotes) PArrays.toc!(t,"DistributedDiscreteGeometryPoisson") + PArrays.tic!(t) + PeriodicPoissonTests.main(distribute,parts,cells=(21,21)) + PeriodicPoissonTests.main(distribute,(prod(parts),1),cells=(21,21),geometry=:remotes) + PArrays.toc!(t,"PeriodicPoissonTests") + + PArrays.tic!(t) + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(21,21),geometry=:remotes) + PArrays.toc!(t,"PeriodicDistributedDiscreteGeometryPoissonTest") + + PArrays.tic!(t) + PeriodicDistributedLSDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) + PArrays.toc!(t,"PeriodicDistributedLSDiscreteGeometryPoissonTest") + if prod(parts) == 4 DistributedAggregatesTests.main(distribute,parts) end diff --git a/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..87bf08e --- /dev/null +++ b/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,8 @@ +module PeriodicDistributedDiscreteGeometryPoissonTestsSeq +using PartitionedArrays +include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(2,2),cells=(21,21)) + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(21,21),geometry=:remotes) +end +end diff --git a/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..52a24da --- /dev/null +++ b/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl @@ -0,0 +1,7 @@ +module PeriodicDistributedLSDiscreteGeometryPoissonTestsSeq +using PartitionedArrays +include("../PeriodicDistributedLSDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + PeriodicDistributedLSDiscreteGeometryPoissonTest.main(distribute,(2,2),cells=(21,21)) +end +end diff --git a/test/DistributedTests/sequential/PeriodicPoissonTests.jl b/test/DistributedTests/sequential/PeriodicPoissonTests.jl new file mode 100644 index 0000000..e0d8ebb --- /dev/null +++ b/test/DistributedTests/sequential/PeriodicPoissonTests.jl @@ -0,0 +1,8 @@ +module PeriodicPoissonTestsSeq +using PartitionedArrays +include("../PeriodicPoissonTests.jl") +with_debug() do distribute + PeriodicPoissonTests.main(distribute,(2,2),cells=(21,21)) + PeriodicPoissonTests.main(distribute,(4,1),cells=(21,21),geometry=:remotes) +end +end diff --git a/test/DistributedTests/sequential/runtests.jl b/test/DistributedTests/sequential/runtests.jl index 5819b69..a9ac993 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -3,6 +3,8 @@ module SequentialTests using Test @time @testset "PoissonSeq" begin include("PoissonTests.jl") end -@time @testset "PoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end +@time @testset "DiscreteGeoPoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end +@time @testset "PeriodicPoissonSeq" begin include("PeriodicPoissonTests.jl") end +@time @testset "PeriodicDiscreteGeoPoissonSeq" begin include("PeriodicDistributedDiscreteGeometryPoissonTest.jl") end end diff --git a/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl b/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl new file mode 100644 index 0000000..4e79136 --- /dev/null +++ b/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl @@ -0,0 +1,76 @@ +# module PoissonAgFEMTests + +using Gridap +using GridapEmbedded +using Test + +u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 +f(x) = -Δ(u)(x) +ud(x) = u(x) + +order = 2 +n = 31 +partition = (n,n) +domain = (0,1,0,1) +bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) +dp = 1 +const h = dp/n + +reffe = ReferenceFE(lagrangian,Float64,order) +Ω_bg = Triangulation(bgmodel) +V_bg = FESpace(Ω_bg,reffe) +φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.55,V_bg) +geo = DiscreteGeometry(φh,bgmodel) +cutgeo = cut(bgmodel,geo) + +strategy = AggregateAllCutCells() +aggregates = aggregate(strategy,cutgeo) + +Ω_bg = Triangulation(bgmodel) +Ω_act = Triangulation(cutgeo,ACTIVE) +Ω = Triangulation(cutgeo,PHYSICAL) +Γ = EmbeddedBoundary(cutgeo) + +n_Γ = get_normal_vector(Γ) + +degree = 2*order +dΩ = Measure(Ω,degree) +dΓ = Measure(Γ,degree) + +model = get_active_model(Ω_act) +Vstd = FESpace(Ω_act,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)) + +V = AgFEMSpace(Vstd,aggregates) +U = TrialFESpace(V) + +const γd = 10.0 + +a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + +l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + +op = AffineFEOperator(a,l,U,V) +uh = solve(op) + +e = u - uh + +l2(u) = sqrt(sum( ∫( u*u )*dΩ )) +h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + +el2 = l2(e) +eh1 = h1(e) +ul2 = l2(uh) +uh1 = h1(uh) + +colors = color_aggregates(aggregates,bgmodel) +writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh]) +writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) +writevtk(Γ,"trian_G") +@test el2/ul2 < 1.e-8 +@test eh1/uh1 < 1.e-7 + +# end # module diff --git a/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl b/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl new file mode 100644 index 0000000..c45c38d --- /dev/null +++ b/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl @@ -0,0 +1,74 @@ +module PoissonAgFEMTests + +using Gridap +using GridapEmbedded +using Test + +u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 +f(x) = -Δ(u)(x) +ud(x) = u(x) + +const R = 0.55 +geom = disk(R,x0=Point(0.5,0.5)) + +n = 31 +partition = (n,n) +domain = (0,1,0,1) +bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) +dp = 1 +const h = dp/n + +cutgeo = cut(bgmodel,geom) + +strategy = AggregateAllCutCells() +aggregates = aggregate(strategy,cutgeo) + +Ω_bg = Triangulation(bgmodel) +Ω_act = Triangulation(cutgeo,ACTIVE) +Ω = Triangulation(cutgeo,PHYSICAL) +Γ = EmbeddedBoundary(cutgeo) + +n_Γ = get_normal_vector(Γ) + +order = 2 +degree = 2*order +dΩ = Measure(Ω,degree) +dΓ = Measure(Γ,degree) + +model = get_active_model(Ω_act) +Vstd = FESpace(Ω_act,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)) + +V = AgFEMSpace(Vstd,aggregates) +U = TrialFESpace(V) + +const γd = 10.0 + +a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + +l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + +op = AffineFEOperator(a,l,U,V) +uh = solve(op) + +e = u - uh + +l2(u) = sqrt(sum( ∫( u*u )*dΩ )) +h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + +el2 = l2(e) +eh1 = h1(e) +ul2 = l2(uh) +uh1 = h1(uh) + +#colors = color_aggregates(aggregates,bgmodel) +#writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh]) +#writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) +#writevtk(Γ,"trian_G") +@test el2/ul2 < 1.e-8 +@test eh1/uh1 < 1.e-7 + +end # module diff --git a/test/GridapEmbeddedTests/runtests.jl b/test/GridapEmbeddedTests/runtests.jl index 1576bc4..a7e3d78 100644 --- a/test/GridapEmbeddedTests/runtests.jl +++ b/test/GridapEmbeddedTests/runtests.jl @@ -6,6 +6,10 @@ using Test @time @testset "PoissonAgFEM" begin include("PoissonAgFEMTests.jl") end +@time @testset "PeriodicPoissonAgFEM" begin include("PeriodicPoissonAgFEMTests.jl") end + +@time @testset "PeriodicDiscreteGeoPoissonAgFEM" begin include("PeriodicDiscreteGeoPoissonAgFEMTests.jl") end + @time @testset "PoissonModalC0AgFEM" begin include("PoissonModalC0AgFEMTests.jl") end @time @testset "BimaterialPoissonCutFEM" begin include("BimaterialPoissonCutFEMTests.jl") end From e0fb29835305190ab4fd49b0b65d25c549eaf3d0 Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 19 Jun 2024 16:51:10 +1000 Subject: [PATCH 08/14] add dep --- src/Distributed/Distributed.jl | 1 + src/Distributed/DistributedDiscreteGeometries.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index a08527e..3c84570 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -10,6 +10,7 @@ using Gridap.CellData using Gridap.Geometry using Gridap.Helpers using Gridap.ReferenceFEs +using Gridap.FESpaces using GridapEmbedded.CSG using GridapEmbedded.LevelSetCutters diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 616c712..2ea3788 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -24,7 +24,7 @@ function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) φh_data = CellData.get_data(φh) - space = get_fe_space(φh) + space = FESpaces.get_fe_space(φh) T = get_dof_value_type(space) values = Vector{T}(undef,num_nodes(own_model)) touched = fill(false,num_nodes(model)) From 67faf8399d598b91e2bfbfb26d3f8a46c64e52dc Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 19 Jun 2024 19:21:15 +1000 Subject: [PATCH 09/14] fix test --- test/AgFEMTests/PeriodicAgFEMSpacesTests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl index cc1845f..ec84d4e 100644 --- a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl +++ b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl @@ -5,7 +5,7 @@ using Gridap using GridapEmbedded using Gridap.Geometry: get_active_model -const R = 0.5 +const R = 0.55 geom = disk(R,x0=Point(0.5,0.5)) n = 21 partition = (n,n) @@ -46,7 +46,7 @@ tol = 10e-9 @test sum( ∫(abs2(v-vhagg))dΩ ) < tol @test sum( ∫(abs2(v-vhagg))dΩ_in ) < tol -vh = FEFunction(V,rand(num_free_dofs(V))) +vh = FEFunction(Vstd,rand(num_free_dofs(Vstd))) vhagg = interpolate(vh,Vagg) @test sum( ∫(abs2(vh-vhagg))dΩ_in ) < tol From 82d459791f0ac827678af6d185da54dd7cfc476c Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Thu, 20 Jun 2024 10:03:02 +1000 Subject: [PATCH 10/14] Revert changes --- test/AgFEMTests/PeriodicAgFEMSpacesTests.jl | 14 +- ...cDistributedDiscreteGeometryPoissonTest.jl | 143 ------------------ ...istributedLSDiscreteGeometryPoissonTest.jl | 110 -------------- test/DistributedTests/PeriodicPoissonTests.jl | 141 ----------------- test/DistributedTests/mpi/runtests_body.jl | 17 --- ...cDistributedDiscreteGeometryPoissonTest.jl | 8 - ...istributedLSDiscreteGeometryPoissonTest.jl | 7 - .../sequential/PeriodicPoissonTests.jl | 8 - test/DistributedTests/sequential/runtests.jl | 2 - .../PeriodicDiscreteGeoPoissonAgFEMTests.jl | 76 ---------- .../PeriodicPoissonAgFEMTests.jl | 74 --------- 11 files changed, 3 insertions(+), 597 deletions(-) delete mode 100644 test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl delete mode 100644 test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl delete mode 100644 test/DistributedTests/PeriodicPoissonTests.jl delete mode 100644 test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl delete mode 100644 test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl delete mode 100644 test/DistributedTests/sequential/PeriodicPoissonTests.jl delete mode 100644 test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl delete mode 100644 test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl diff --git a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl index ec84d4e..26b34ba 100644 --- a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl +++ b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl @@ -28,7 +28,7 @@ dΩ = Measure(Ω,2) dΩ_in = Measure(Ω_in,2) model = get_active_model(Ω_ac) -order = 2 +order = 1 # In the physical domain cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order) @@ -42,14 +42,10 @@ vhagg = interpolate(v,Vagg) writevtk(Ω_ac,"test",cellfields=["v"=>vhagg]) -tol = 10e-9 +tol = 10e-7 @test sum( ∫(abs2(v-vhagg))dΩ ) < tol @test sum( ∫(abs2(v-vhagg))dΩ_in ) < tol -vh = FEFunction(Vstd,rand(num_free_dofs(Vstd))) -vhagg = interpolate(vh,Vagg) -@test sum( ∫(abs2(vh-vhagg))dΩ_in ) < tol - # In the reference space reffe = ReferenceFE(lagrangian,Float64,order) @@ -59,14 +55,10 @@ Vagg = AgFEMSpace(V,aggregates) v(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 vhagg = interpolate(v,Vagg) -tol = 10e-9 +tol = 10e-7 @test sum( ∫(abs2(v-vhagg))dΩ ) < tol @test sum( ∫(abs2(v-vhagg))dΩ_in ) < tol -vh = FEFunction(V,rand(num_free_dofs(V))) -vhagg = interpolate(vh,Vagg) -@test sum( ∫(abs2(vh-vhagg))dΩ_in ) < tol - #cellfields = ["vh"=>vh,"vhagg"=>vhagg,"e"=>vh-vhagg] #writevtk(Ω_bg,"trian_bg",nsubcells=10,cellfields=cellfields) #writevtk(Ω_in,"trian_in",nsubcells=10,cellfields=cellfields) diff --git a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl deleted file mode 100644 index 8bd3e67..0000000 --- a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl +++ /dev/null @@ -1,143 +0,0 @@ -module PeriodicDistributedDiscreteGeometryPoissonTest - -using Gridap -using GridapEmbedded -using GridapDistributed -using PartitionedArrays -using Test - -using GridapEmbedded.CSG -using GridapEmbedded.LevelSetCutters - -function main(distribute,parts; - threshold=1, - n=8, - cells=(n,n), - geometry=:circle) - - ranks = distribute(LinearIndices((prod(parts),))) - - u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 - f(x) = -Δ(u)(x) - ud(x) = u(x) - - geometries = Dict( - :circle => circle_geometry, - :remotes => remotes_geometry, - ) - - bgmodel,_geo = geometries[geometry](ranks,parts,cells) - geo = discretize(_geo,bgmodel) - - D = 2 - cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) - meas = map(first,cell_meas) |> PartitionedArrays.getany - h = meas^(1/D) - - cutgeo = cut(bgmodel,geo) - cutgeo_facets = cut_facets(bgmodel,geo) - - strategy = AggregateCutCellsByThreshold(threshold) - bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) - - Ω_bg = Triangulation(bgmodel) - Ω_act = Triangulation(cutgeo,ACTIVE) - Ω = Triangulation(cutgeo,PHYSICAL) - Γ = EmbeddedBoundary(cutgeo) - - n_Γ = get_normal_vector(Γ) - - order = 2 - degree = 2*order - dΩ = Measure(Ω,degree) - dΓ = Measure(Γ,degree) - - reffe = ReferenceFE(lagrangian,Float64,order) - - Vstd = FESpace(Ω_act,reffe) - - V = AgFEMSpace(bgmodel,Vstd,aggregates) - U = TrialFESpace(V) - - - γd = 10.0 - - a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - - l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - - op = AffineFEOperator(a,l,U,V) - uh = solve(op) - - e = u - uh - - l2(u) = sqrt(sum( ∫( u*u )*dΩ )) - h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - - el2 = l2(e) - eh1 = h1(e) - ul2 = l2(uh) - uh1 = h1(uh) - - # - colors = map(color_aggregates,aggregates,local_views(bgmodel)) - gids = get_cell_gids(bgmodel) - - - global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid - map(i-> i==0 ? 0 : gid[i],agg) - end - own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid - map(Reindex(agg),oid) - end - own_colors = map(colors,own_to_local(gids)) do col,oid - map(Reindex(col),oid) - end - - writevtk(Ω_bg,"trian", - celldata=[ - "aggregate"=>own_aggregates, - "color"=>own_colors, - "gid"=>own_to_global(gids)])#, - # cellfields=["uh"=>uh]) - - writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) - writevtk(Γ,"trian_G") - @test el2/ul2 < 1.e-8 - @test eh1/uh1 < 1.e-7 - -end - -function circle_geometry(ranks,parts,cells) - L = 1 - p0 = Point(0.5,0.5) - pmin = p0-L/2 - pmax = p0+L/2 - R = 0.55 - geo = disk(R,x0=p0) - bgmodel = CartesianDiscreteModel(ranks,parts,pmin,pmax,cells;isperiodic=(true,true)) - bgmodel,geo -end - -function remotes_geometry(ranks,parts,cells) - x0 = Point(0.0,0.4) - d1 = VectorValue(1.0,0.0) - d2 = VectorValue(0.0,0.2) - geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) - - x0 = Point(0.4,0.0) - d1 = VectorValue(0.2,0.0) - d2 = VectorValue(0.0,1.0) - geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) - geo = union(geo1,geo2) - - domain = (0, 1, 0, 1) - bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) - bgmodel,geo -end - -end # module \ No newline at end of file diff --git a/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl deleted file mode 100644 index f0b659e..0000000 --- a/test/DistributedTests/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl +++ /dev/null @@ -1,110 +0,0 @@ -module PeriodicDistributedLSDiscreteGeometryPoissonTest - -using Gridap -using GridapEmbedded -using GridapDistributed -using PartitionedArrays -using Test - -using GridapEmbedded.CSG -using GridapEmbedded.LevelSetCutters - -function main(distribute,parts; - threshold=1, - n=21, - cells=(n,n)) - - order = 2 - ranks = distribute(LinearIndices((prod(parts),))) - domain = (0,1,0,1) - bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) - - u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 - f(x) = -Δ(u)(x) - ud(x) = u(x) - - reffe = ReferenceFE(lagrangian,Float64,order) - Ω_bg = Triangulation(bgmodel) - V_bg = FESpace(Ω_bg,reffe) - φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.55,V_bg) - geo = DiscreteGeometry(φh,bgmodel) - - D = 2 - cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) - meas = map(first,cell_meas) |> PartitionedArrays.getany - h = meas^(1/D) - - cutgeo = cut(bgmodel,geo) - cutgeo_facets = cut_facets(bgmodel,geo) - - strategy = AggregateCutCellsByThreshold(threshold) - bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) - - Ω_act = Triangulation(cutgeo,ACTIVE) - Ω = Triangulation(cutgeo,PHYSICAL) - Γ = EmbeddedBoundary(cutgeo) - - n_Γ = get_normal_vector(Γ) - - order = 2 - degree = 2*order - dΩ = Measure(Ω,degree) - dΓ = Measure(Γ,degree) - - Vstd = FESpace(Ω_act,reffe) - V = AgFEMSpace(bgmodel,Vstd,aggregates) - U = TrialFESpace(V) - - γd = 10.0 - - a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - - l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - - op = AffineFEOperator(a,l,U,V) - uh = solve(op) - - e = u - uh - - l2(u) = sqrt(sum( ∫( u*u )*dΩ )) - h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - - el2 = l2(e) - eh1 = h1(e) - ul2 = l2(uh) - uh1 = h1(uh) - - # - colors = map(color_aggregates,aggregates,local_views(bgmodel)) - gids = get_cell_gids(bgmodel) - - - global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid - map(i-> i==0 ? 0 : gid[i],agg) - end - own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid - map(Reindex(agg),oid) - end - own_colors = map(colors,own_to_local(gids)) do col,oid - map(Reindex(col),oid) - end - - writevtk(Ω_bg,"trian", - celldata=[ - "aggregate"=>own_aggregates, - "color"=>own_colors, - "gid"=>own_to_global(gids)])#, - # cellfields=["uh"=>uh]) - - writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) - writevtk(Γ,"trian_G") - @test el2/ul2 < 1.e-8 - @test eh1/uh1 < 1.e-7 - -end - -end # module \ No newline at end of file diff --git a/test/DistributedTests/PeriodicPoissonTests.jl b/test/DistributedTests/PeriodicPoissonTests.jl deleted file mode 100644 index 130ad41..0000000 --- a/test/DistributedTests/PeriodicPoissonTests.jl +++ /dev/null @@ -1,141 +0,0 @@ -module PeriodicPoissonTests - -using Gridap -using GridapEmbedded -using GridapDistributed -using PartitionedArrays -using Test - -using GridapEmbedded.CSG - -function main(distribute,parts; - threshold=1, - n=8, - cells=(n,n), - geometry=:circle) - - ranks = distribute(LinearIndices((prod(parts),))) - - u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 - f(x) = -Δ(u)(x) - ud(x) = u(x) - - geometries = Dict( - :circle => circle_geometry, - :remotes => remotes_geometry, - ) - - bgmodel,geo = geometries[geometry](ranks,parts,cells) - - D = 2 - cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) - meas = map(first,cell_meas) |> PartitionedArrays.getany - h = meas^(1/D) - - cutgeo = cut(bgmodel,geo) - cutgeo_facets = cut_facets(bgmodel,geo) - - strategy = AggregateCutCellsByThreshold(threshold) - bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) - - Ω_bg = Triangulation(bgmodel) - Ω_act = Triangulation(cutgeo,ACTIVE) - Ω = Triangulation(cutgeo,PHYSICAL) - Γ = EmbeddedBoundary(cutgeo) - - n_Γ = get_normal_vector(Γ) - - order = 2 - degree = 2*order - dΩ = Measure(Ω,degree) - dΓ = Measure(Γ,degree) - - reffe = ReferenceFE(lagrangian,Float64,order) - - Vstd = FESpace(Ω_act,reffe) - - V = AgFEMSpace(bgmodel,Vstd,aggregates) - U = TrialFESpace(V) - - - γd = 10.0 - - a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - - l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - - op = AffineFEOperator(a,l,U,V) - uh = solve(op) - - e = u - uh - - l2(u) = sqrt(sum( ∫( u*u )*dΩ )) - h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - - el2 = l2(e) - eh1 = h1(e) - ul2 = l2(uh) - uh1 = h1(uh) - - # - colors = map(color_aggregates,aggregates,local_views(bgmodel)) - gids = get_cell_gids(bgmodel) - - - global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid - map(i-> i==0 ? 0 : gid[i],agg) - end - own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid - map(Reindex(agg),oid) - end - own_colors = map(colors,own_to_local(gids)) do col,oid - map(Reindex(col),oid) - end - - writevtk(Ω_bg,"trian", - celldata=[ - "aggregate"=>own_aggregates, - "color"=>own_colors, - "gid"=>own_to_global(gids)])#, - # cellfields=["uh"=>uh]) - - writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) - writevtk(Γ,"trian_G") - @test el2/ul2 < 1.e-8 - @test eh1/uh1 < 1.e-7 - -end - -function circle_geometry(ranks,parts,cells) - L = 1 - p0 = Point(0.5,0.5) - pmin = p0-L/2 - pmax = p0+L/2 - R = 0.55 - geo = disk(R,x0=p0) - bgmodel = CartesianDiscreteModel(ranks,parts,pmin,pmax,cells;isperiodic=(true,true)) - bgmodel,geo -end - -function remotes_geometry(ranks,parts,cells) - x0 = Point(0.0,0.4) - d1 = VectorValue(1.0,0.0) - d2 = VectorValue(0.0,0.2) - geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) - - x0 = Point(0.4,0.0) - d1 = VectorValue(0.2,0.0) - d2 = VectorValue(0.0,1.0) - geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) - geo = union(geo1,geo2) - - domain = (0, 1, 0, 1) - bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) - bgmodel,geo -end - -end # module diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl index cdbb6dc..7c44d8b 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -7,9 +7,6 @@ using MPI include("../PoissonTests.jl") include("../AggregatesTests.jl") include("../DistributedDiscreteGeometryPoissonTest.jl") -include("../PeriodicPoissonTests.jl") -include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") -include("../PeriodicDistributedLSDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -30,20 +27,6 @@ function all_tests(distribute,parts) DistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(12,12),geometry=:remotes) PArrays.toc!(t,"DistributedDiscreteGeometryPoisson") - PArrays.tic!(t) - PeriodicPoissonTests.main(distribute,parts,cells=(21,21)) - PeriodicPoissonTests.main(distribute,(prod(parts),1),cells=(21,21),geometry=:remotes) - PArrays.toc!(t,"PeriodicPoissonTests") - - PArrays.tic!(t) - PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) - PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(21,21),geometry=:remotes) - PArrays.toc!(t,"PeriodicDistributedDiscreteGeometryPoissonTest") - - PArrays.tic!(t) - PeriodicDistributedLSDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) - PArrays.toc!(t,"PeriodicDistributedLSDiscreteGeometryPoissonTest") - if prod(parts) == 4 DistributedAggregatesTests.main(distribute,parts) end diff --git a/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl deleted file mode 100644 index 87bf08e..0000000 --- a/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl +++ /dev/null @@ -1,8 +0,0 @@ -module PeriodicDistributedDiscreteGeometryPoissonTestsSeq -using PartitionedArrays -include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") -with_debug() do distribute - PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(2,2),cells=(21,21)) - PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(21,21),geometry=:remotes) -end -end diff --git a/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl deleted file mode 100644 index 52a24da..0000000 --- a/test/DistributedTests/sequential/PeriodicDistributedLSDiscreteGeometryPoissonTest.jl +++ /dev/null @@ -1,7 +0,0 @@ -module PeriodicDistributedLSDiscreteGeometryPoissonTestsSeq -using PartitionedArrays -include("../PeriodicDistributedLSDiscreteGeometryPoissonTest.jl") -with_debug() do distribute - PeriodicDistributedLSDiscreteGeometryPoissonTest.main(distribute,(2,2),cells=(21,21)) -end -end diff --git a/test/DistributedTests/sequential/PeriodicPoissonTests.jl b/test/DistributedTests/sequential/PeriodicPoissonTests.jl deleted file mode 100644 index e0d8ebb..0000000 --- a/test/DistributedTests/sequential/PeriodicPoissonTests.jl +++ /dev/null @@ -1,8 +0,0 @@ -module PeriodicPoissonTestsSeq -using PartitionedArrays -include("../PeriodicPoissonTests.jl") -with_debug() do distribute - PeriodicPoissonTests.main(distribute,(2,2),cells=(21,21)) - PeriodicPoissonTests.main(distribute,(4,1),cells=(21,21),geometry=:remotes) -end -end diff --git a/test/DistributedTests/sequential/runtests.jl b/test/DistributedTests/sequential/runtests.jl index a9ac993..bf7521b 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -4,7 +4,5 @@ using Test @time @testset "PoissonSeq" begin include("PoissonTests.jl") end @time @testset "DiscreteGeoPoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end -@time @testset "PeriodicPoissonSeq" begin include("PeriodicPoissonTests.jl") end -@time @testset "PeriodicDiscreteGeoPoissonSeq" begin include("PeriodicDistributedDiscreteGeometryPoissonTest.jl") end end diff --git a/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl b/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl deleted file mode 100644 index 4e79136..0000000 --- a/test/GridapEmbeddedTests/PeriodicDiscreteGeoPoissonAgFEMTests.jl +++ /dev/null @@ -1,76 +0,0 @@ -# module PoissonAgFEMTests - -using Gridap -using GridapEmbedded -using Test - -u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 -f(x) = -Δ(u)(x) -ud(x) = u(x) - -order = 2 -n = 31 -partition = (n,n) -domain = (0,1,0,1) -bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) -dp = 1 -const h = dp/n - -reffe = ReferenceFE(lagrangian,Float64,order) -Ω_bg = Triangulation(bgmodel) -V_bg = FESpace(Ω_bg,reffe) -φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.55,V_bg) -geo = DiscreteGeometry(φh,bgmodel) -cutgeo = cut(bgmodel,geo) - -strategy = AggregateAllCutCells() -aggregates = aggregate(strategy,cutgeo) - -Ω_bg = Triangulation(bgmodel) -Ω_act = Triangulation(cutgeo,ACTIVE) -Ω = Triangulation(cutgeo,PHYSICAL) -Γ = EmbeddedBoundary(cutgeo) - -n_Γ = get_normal_vector(Γ) - -degree = 2*order -dΩ = Measure(Ω,degree) -dΓ = Measure(Γ,degree) - -model = get_active_model(Ω_act) -Vstd = FESpace(Ω_act,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)) - -V = AgFEMSpace(Vstd,aggregates) -U = TrialFESpace(V) - -const γd = 10.0 - -a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - -l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - -op = AffineFEOperator(a,l,U,V) -uh = solve(op) - -e = u - uh - -l2(u) = sqrt(sum( ∫( u*u )*dΩ )) -h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - -el2 = l2(e) -eh1 = h1(e) -ul2 = l2(uh) -uh1 = h1(uh) - -colors = color_aggregates(aggregates,bgmodel) -writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh]) -writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) -writevtk(Γ,"trian_G") -@test el2/ul2 < 1.e-8 -@test eh1/uh1 < 1.e-7 - -# end # module diff --git a/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl b/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl deleted file mode 100644 index c45c38d..0000000 --- a/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl +++ /dev/null @@ -1,74 +0,0 @@ -module PoissonAgFEMTests - -using Gridap -using GridapEmbedded -using Test - -u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 -f(x) = -Δ(u)(x) -ud(x) = u(x) - -const R = 0.55 -geom = disk(R,x0=Point(0.5,0.5)) - -n = 31 -partition = (n,n) -domain = (0,1,0,1) -bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) -dp = 1 -const h = dp/n - -cutgeo = cut(bgmodel,geom) - -strategy = AggregateAllCutCells() -aggregates = aggregate(strategy,cutgeo) - -Ω_bg = Triangulation(bgmodel) -Ω_act = Triangulation(cutgeo,ACTIVE) -Ω = Triangulation(cutgeo,PHYSICAL) -Γ = EmbeddedBoundary(cutgeo) - -n_Γ = get_normal_vector(Γ) - -order = 2 -degree = 2*order -dΩ = Measure(Ω,degree) -dΓ = Measure(Γ,degree) - -model = get_active_model(Ω_act) -Vstd = FESpace(Ω_act,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)) - -V = AgFEMSpace(Vstd,aggregates) -U = TrialFESpace(V) - -const γd = 10.0 - -a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - -l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - -op = AffineFEOperator(a,l,U,V) -uh = solve(op) - -e = u - uh - -l2(u) = sqrt(sum( ∫( u*u )*dΩ )) -h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - -el2 = l2(e) -eh1 = h1(e) -ul2 = l2(uh) -uh1 = h1(uh) - -#colors = color_aggregates(aggregates,bgmodel) -#writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh]) -#writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) -#writevtk(Γ,"trian_G") -@test el2/ul2 < 1.e-8 -@test eh1/uh1 < 1.e-7 - -end # module From de2d2bbe87ce6399fc5fc306284449bf9a6a6026 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Thu, 20 Jun 2024 12:03:54 +1000 Subject: [PATCH 11/14] tests passing --- src/Distributed/Distributed.jl | 1 + .../DistributedDiscreteGeometries.jl | 6 +- ...istributedLSDiscreteGeometryPoissonTest.jl | 110 ++++++++++++++++++ test/DistributedTests/mpi/runtests_body.jl | 5 + ...istributedLSDiscreteGeometryPoissonTest.jl | 7 ++ test/DistributedTests/sequential/runtests.jl | 1 + .../DiscreteGeometriesTests.jl | 12 ++ 7 files changed, 139 insertions(+), 3 deletions(-) create mode 100644 test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/sequential/DistributedLSDiscreteGeometryPoissonTest.jl diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index 3c84570..dec113c 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -26,6 +26,7 @@ using GridapEmbedded.AgFEM: AggregateCutCellsByThreshold using GridapEmbedded.MomentFittedQuadratures: MomentFitted using Gridap.Geometry: AppendedTriangulation using Gridap.Geometry: get_face_to_parent_face +using Gridap.Arrays: find_inverse_index_map using GridapDistributed: DistributedDiscreteModel using GridapDistributed: DistributedTriangulation using GridapDistributed: DistributedFESpace diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 2ea3788..c8a949e 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -10,9 +10,9 @@ function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) values = map(local_views(φh),local_views(model),local_views(gids)) do φh, model, gids # Maps from the no-ghost model to the original model own_model = remove_ghost_cells(model,gids) - own_to_local_node = Geometry.get_face_to_parent_face(own_model,0) - local_to_own_node = Arrays.find_inverse_index_map(own_to_local_node,num_nodes(model)) - own_to_local_cell = Geometry.get_face_to_parent_face(own_model,Dc) + own_to_local_node = get_face_to_parent_face(own_model,0) + local_to_own_node = find_inverse_index_map(own_to_local_node,num_nodes(model)) + own_to_local_cell = get_face_to_parent_face(own_model,Dc) # Cell-to-node map for the original model # topo = get_grid_topology(model) diff --git a/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..68bbf39 --- /dev/null +++ b/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl @@ -0,0 +1,110 @@ +module DistributedLSDiscreteGeometryPoissonTest + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapEmbedded.CSG +using GridapEmbedded.LevelSetCutters + +function main(distribute,parts; + threshold=1, + n=21, + cells=(n,n)) + + order = 2 + ranks = distribute(LinearIndices((prod(parts),))) + domain = (0,1,0,1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells) + + u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 + f(x) = -Δ(u)(x) + ud(x) = u(x) + + reffe = ReferenceFE(lagrangian,Float64,order) + Ω_bg = Triangulation(bgmodel) + V_bg = FESpace(Ω_bg,reffe) + φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.35,V_bg) + geo = DiscreteGeometry(φh,bgmodel) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 2 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + Vstd = FESpace(Ω_act,reffe) + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V) + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +end # module \ No newline at end of file diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl index 7c44d8b..048144f 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -7,6 +7,7 @@ using MPI include("../PoissonTests.jl") include("../AggregatesTests.jl") include("../DistributedDiscreteGeometryPoissonTest.jl") +include("../DistributedLSDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -27,6 +28,10 @@ function all_tests(distribute,parts) DistributedDiscreteGeometryPoissonTest.main(distribute,(prod(parts),1),cells=(12,12),geometry=:remotes) PArrays.toc!(t,"DistributedDiscreteGeometryPoisson") + PArrays.tic!(t) + DistributedLSDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) + PArrays.toc!(t,"DistributedLSDiscreteGeometryPoissonTest") + if prod(parts) == 4 DistributedAggregatesTests.main(distribute,parts) end diff --git a/test/DistributedTests/sequential/DistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/DistributedLSDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..afdba27 --- /dev/null +++ b/test/DistributedTests/sequential/DistributedLSDiscreteGeometryPoissonTest.jl @@ -0,0 +1,7 @@ +module DistributedLSDiscreteGeometryPoissonTestsSeq +using PartitionedArrays +include("../DistributedLSDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + DistributedLSDiscreteGeometryPoissonTest.main(distribute,(2,2),cells=(21,21)) +end +end diff --git a/test/DistributedTests/sequential/runtests.jl b/test/DistributedTests/sequential/runtests.jl index bf7521b..0e8f21f 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -4,5 +4,6 @@ using Test @time @testset "PoissonSeq" begin include("PoissonTests.jl") end @time @testset "DiscreteGeoPoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end +@time @testset "LSDiscreteGeoPoissonSeq" begin include("DistributedLSDiscreteGeometryPoissonTest.jl") end end diff --git a/test/LevelSetCuttersTests/DiscreteGeometriesTests.jl b/test/LevelSetCuttersTests/DiscreteGeometriesTests.jl index 5a7479e..f246c2e 100644 --- a/test/LevelSetCuttersTests/DiscreteGeometriesTests.jl +++ b/test/LevelSetCuttersTests/DiscreteGeometriesTests.jl @@ -38,4 +38,16 @@ test_geometry(geo4_y) #print_tree(stdout,get_tree(geo4_x)) #print_tree(stdout,replace_data(d->objectid(first(d)),get_tree(geo4_x))) +#using a cellfield +domain = (0,1,0,1) +n = 10 +bgmodel = CartesianDiscreteModel(domain,(n,n)) + +reffe = ReferenceFE(lagrangian,Float64,1) +Ω_bg = Triangulation(bgmodel) +V_bg = FESpace(Ω_bg,reffe) +φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.55,V_bg) +geo = DiscreteGeometry(φh,bgmodel) +test_geometry(geo) + end # module From c9da43dedc9981cfdbe5fed169db517f9d8ec37b Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Thu, 20 Jun 2024 14:24:36 +1000 Subject: [PATCH 12/14] Added periodic tests --- ...cDistributedDiscreteGeometryPoissonTest.jl | 193 ++++++++++++++++++ test/DistributedTests/mpi/runtests_body.jl | 7 + ...cDistributedDiscreteGeometryPoissonTest.jl | 8 + test/DistributedTests/sequential/runtests.jl | 3 + .../PeriodicPoissonAgFEMTests.jl | 111 ++++++++++ test/GridapEmbeddedTests/runtests.jl | 2 - 6 files changed, 322 insertions(+), 2 deletions(-) create mode 100644 test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl create mode 100644 test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl create mode 100644 test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl diff --git a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..26cde8c --- /dev/null +++ b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,193 @@ +module PeriodicDistributedDiscreteGeometryPoissonTest + +using Gridap +using GridapEmbedded +using GridapDistributed +using PartitionedArrays +using Test + +using GridapDistributed: DistributedDiscreteModel + +using Gridap.Geometry +using Gridap.Geometry: get_vertex_coordinates,get_faces + +using GridapEmbedded.CSG +using GridapEmbedded.LevelSetCutters + +function update_labels!(e::Integer,model::DistributedDiscreteModel,f_Γ::Function,name::String) + mask = mark_nodes(f_Γ,model) + cell_to_entity = map(local_views(model),local_views(mask)) do model,mask + _update_labels_locally!(e,model,mask,name) + end + cell_gids=get_cell_gids(model) + cache=GridapDistributed.fetch_vector_ghost_values_cache(cell_to_entity,partition(cell_gids)) + GridapDistributed.fetch_vector_ghost_values!(cell_to_entity,cache) + nothing +end + +function _update_labels_locally!(e,model::CartesianDiscreteModel{2},mask,name) + topo = get_grid_topology(model) + labels = get_face_labeling(model) + cell_to_entity = labels.d_to_dface_to_entity[end] + entity = maximum(cell_to_entity) + e + # Vertices + vtxs_Γ = findall(mask) + vtx_edge_connectivity = Array(get_faces(topo,0,1)[vtxs_Γ]) + # Edges + edge_entries = [findall(x->any(x .∈ vtx_edge_connectivity[1:end.!=j]), + vtx_edge_connectivity[j]) for j = 1:length(vtx_edge_connectivity)] + edge_Γ = unique(reduce(vcat,getindex.(vtx_edge_connectivity,edge_entries),init=[])) + labels.d_to_dface_to_entity[1][vtxs_Γ] .= entity + labels.d_to_dface_to_entity[2][edge_Γ] .= entity + add_tag!(labels,name,[entity]) + return cell_to_entity +end + +function mark_nodes(f,model::DistributedDiscreteModel) + local_masks = map(local_views(model)) do model + mark_nodes(f,model) + end + gids = get_face_gids(model,0) + mask = PVector(local_masks,partition(gids)) + assemble!(|,mask) |> fetch # Ghosts -> Owned with `or` applied + consistent!(mask) |> fetch # Owned -> Ghost + return mask +end + +function mark_nodes(f,model::DiscreteModel) + topo = get_grid_topology(model) + coords = get_vertex_coordinates(topo) + mask = map(f,coords) + return mask +end + +function main(distribute,parts; + threshold=1, + n=8, + cells=(n,n), + geometry=:circle) + + ranks = distribute(LinearIndices((prod(parts),))) + + u(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2 + f(x) = -Δ(u)(x) + ud(x) = u(x) + + geometries = Dict( + :circle => circle_geometry, + :remotes => remotes_geometry, + ) + + bgmodel,_geo = geometries[geometry](ranks,parts,cells) + update_labels!(1,bgmodel,x->iszero(x[1]) || iszero(x[2]),"outer_boundary") + + geo = discretize(_geo,bgmodel) + + D = 2 + cell_meas = map(get_cell_measure∘Triangulation,local_views(bgmodel)) + meas = map(first,cell_meas) |> PartitionedArrays.getany + h = meas^(1/D) + + cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) + + strategy = AggregateCutCellsByThreshold(threshold) + bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo) + + Ω_bg = Triangulation(bgmodel) + Ω_act = Triangulation(cutgeo,ACTIVE) + Ω = Triangulation(cutgeo,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo) + + n_Γ = get_normal_vector(Γ) + + order = 2 + degree = 2*order + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + + reffe = ReferenceFE(lagrangian,Float64,order) + + Vstd = FESpace(Ω_act,reffe;dirichlet_tags="outer_boundary") + + V = AgFEMSpace(bgmodel,Vstd,aggregates) + U = TrialFESpace(V,ud) + + γd = 10.0 + + a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + + l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + + op = AffineFEOperator(a,l,U,V) + uh = solve(op) + + e = u - uh + + l2(u) = sqrt(sum( ∫( u*u )*dΩ )) + h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + + el2 = l2(e) + eh1 = h1(e) + ul2 = l2(uh) + uh1 = h1(uh) + + # + colors = map(color_aggregates,aggregates,local_views(bgmodel)) + gids = get_cell_gids(bgmodel) + + + global_aggregates = map(aggregates,local_to_global(gids)) do agg,gid + map(i-> i==0 ? 0 : gid[i],agg) + end + own_aggregates = map(global_aggregates,own_to_local(gids)) do agg,oid + map(Reindex(agg),oid) + end + own_colors = map(colors,own_to_local(gids)) do col,oid + map(Reindex(col),oid) + end + + writevtk(Ω_bg,"trian", + celldata=[ + "aggregate"=>own_aggregates, + "color"=>own_colors, + "gid"=>own_to_global(gids)])#, + # cellfields=["uh"=>uh]) + + writevtk(Ω,"trian_O",cellfields=["uh"=>uh]) + writevtk(Γ,"trian_G") + @test el2/ul2 < 1.e-8 + @test eh1/uh1 < 1.e-7 + +end + +function circle_geometry(ranks,parts,cells) + p0 = Point(0.5,0.5) + R = 0.55 + geo = disk(R,x0=p0) + bgmodel = CartesianDiscreteModel(ranks,parts,(0,1,0,1),cells;isperiodic=(true,true)) + bgmodel,geo +end + +function remotes_geometry(ranks,parts,cells) + x0 = Point(0.0,0.4) + d1 = VectorValue(1.0,0.0) + d2 = VectorValue(0.0,0.2) + geo1 = quadrilateral(;x0=x0,d1=d1,d2=d2) + + x0 = Point(0.4,0.0) + d1 = VectorValue(0.2,0.0) + d2 = VectorValue(0.0,1.0) + geo2 = quadrilateral(;x0=x0,d1=d1,d2=d2) + geo = union(geo1,geo2) + + domain = (0, 1, 0, 1) + bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells;isperiodic=(true,true)) + bgmodel,geo +end + +end # module \ No newline at end of file diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl index 048144f..d0510c0 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -8,6 +8,7 @@ include("../PoissonTests.jl") include("../AggregatesTests.jl") include("../DistributedDiscreteGeometryPoissonTest.jl") include("../DistributedLSDiscreteGeometryPoissonTest.jl") +include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -32,6 +33,12 @@ function all_tests(distribute,parts) DistributedLSDiscreteGeometryPoissonTest.main(distribute,parts,cells=(21,21)) PArrays.toc!(t,"DistributedLSDiscreteGeometryPoissonTest") + ## Disabled due to Issue #87 + # PArrays.tic!(t) + # PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,4),cells=(31,31),geometry=:circle) + # PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(31,31),geometry=:remotes) + # PArrays.toc!(t,"PeriodicDistributedDiscreteGeometryPoisson") + if prod(parts) == 4 DistributedAggregatesTests.main(distribute,parts) end diff --git a/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..7802ba3 --- /dev/null +++ b/test/DistributedTests/sequential/PeriodicDistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,8 @@ +module PeriodicDistributedDiscreteGeometryPoissonTest +using PartitionedArrays +include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,4),cells=(31,31),geometry=:circle) + PeriodicDistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(31,31),geometry=:remotes) +end +end diff --git a/test/DistributedTests/sequential/runtests.jl b/test/DistributedTests/sequential/runtests.jl index 0e8f21f..57910a8 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -6,4 +6,7 @@ using Test @time @testset "DiscreteGeoPoissonSeq" begin include("DistributedDiscreteGeometryPoissonTest.jl") end @time @testset "LSDiscreteGeoPoissonSeq" begin include("DistributedLSDiscreteGeometryPoissonTest.jl") end +## Disabled due to Issue #87 +# @time @testset "PeriodicDiscreteGeoPoissonSeq" begin include("PeriodicDistributedDiscreteGeometryPoissonTest.jl") end + end diff --git a/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl b/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl new file mode 100644 index 0000000..0ce321e --- /dev/null +++ b/test/GridapEmbeddedTests/PeriodicPoissonAgFEMTests.jl @@ -0,0 +1,111 @@ +module PeriodicPoissonAgFEMTests + +using Gridap +using GridapEmbedded +using Test + +function update_labels!(e::Integer,model::CartesianDiscreteModel,f_Γ::Function,name::String) + mask = mark_nodes(f_Γ,model) + _update_labels_locally!(e,model,mask,name) + nothing +end + +function _update_labels_locally!(e,model::CartesianDiscreteModel{2},mask,name) + topo = Gridap.Geometry.get_grid_topology(model) + labels = Gridap.Geometry.get_face_labeling(model) + cell_to_entity = labels.d_to_dface_to_entity[end] + entity = maximum(cell_to_entity) + e + # Vertices + vtxs_Γ = findall(mask) + vtx_edge_connectivity = Array(Gridap.Geometry.get_faces(topo,0,1)[vtxs_Γ]) + # Edges + edge_entries = [findall(x->any(x .∈ vtx_edge_connectivity[1:end.!=j]), + vtx_edge_connectivity[j]) for j = 1:length(vtx_edge_connectivity)] + edge_Γ = unique(reduce(vcat,getindex.(vtx_edge_connectivity,edge_entries),init=[])) + labels.d_to_dface_to_entity[1][vtxs_Γ] .= entity + labels.d_to_dface_to_entity[2][edge_Γ] .= entity + add_tag!(labels,name,[entity]) + return cell_to_entity +end + +function mark_nodes(f,model::DiscreteModel) + topo = Gridap.Geometry.get_grid_topology(model) + coords = Gridap.Geometry.get_vertex_coordinates(topo) + mask = map(f,coords) + return mask +end + +u(x) = (x[1]-0.5)^2 + 2(x[2]-0.5)^2 +f(x) = -Δ(u)(x) +ud(x) = u(x) + +# R = 0.3 +# geo1 = square(;L=2) +# geo2 = disk(R,x0=Point(0.5,0.5)) + +geom = disk(0.55,x0=Point(0.5,0.5)) +# geom = setdiff(geo1,geo2) + +n = 31 +partition = (n,n) +domain = (0,1,0,1) +bgmodel = CartesianDiscreteModel(domain,partition;isperiodic=(true,true)) +update_labels!(1,bgmodel,x->iszero(x[1]) || iszero(x[2]),"outer_boundary") + +dp = 1 +const h = dp/n + +cutgeo = cut(bgmodel,geom) + +strategy = AggregateAllCutCells() +aggregates = aggregate(strategy,cutgeo) + +Ω_bg = Triangulation(bgmodel) +Ω_act = Triangulation(cutgeo,ACTIVE) +Ω = Triangulation(cutgeo,PHYSICAL) +Γ = EmbeddedBoundary(cutgeo) + +n_Γ = get_normal_vector(Γ) + +order = 2 +degree = 2*order +dΩ = Measure(Ω,degree) +dΓ = Measure(Γ,degree) + +model = get_active_model(Ω_act) +Vstd = FESpace(Ω_act,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order);dirichlet_tags="outer_boundary") + +V = AgFEMSpace(Vstd,aggregates) +U = TrialFESpace(V,ud) + +const γd = 10.0 + +a(u,v) = + ∫( ∇(v)⋅∇(u) ) * dΩ + + ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ + +l(v) = + ∫( v*f ) * dΩ + + ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ + +op = AffineFEOperator(a,l,U,V) +uh = solve(op) + +e = u - uh + +l2(u) = sqrt(sum( ∫( u*u )*dΩ )) +h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) + +el2 = l2(e) +eh1 = h1(e) +ul2 = l2(uh) +uh1 = h1(uh) + +#colors = color_aggregates(aggregates,bgmodel) +#writevtk(Ω_bg,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh]) +# writevtk(Ω,"trian_O",cellfields=["uh"=>uh,"u"=>u]) +# writevtk(Γ,"trian_G") +@test el2/ul2 < 1.e-8 +@test eh1/uh1 < 1.e-7 + +end # module diff --git a/test/GridapEmbeddedTests/runtests.jl b/test/GridapEmbeddedTests/runtests.jl index a7e3d78..f6d52df 100644 --- a/test/GridapEmbeddedTests/runtests.jl +++ b/test/GridapEmbeddedTests/runtests.jl @@ -8,8 +8,6 @@ using Test @time @testset "PeriodicPoissonAgFEM" begin include("PeriodicPoissonAgFEMTests.jl") end -@time @testset "PeriodicDiscreteGeoPoissonAgFEM" begin include("PeriodicDiscreteGeoPoissonAgFEMTests.jl") end - @time @testset "PoissonModalC0AgFEM" begin include("PoissonModalC0AgFEMTests.jl") end @time @testset "BimaterialPoissonCutFEM" begin include("BimaterialPoissonCutFEMTests.jl") end From ef1fc9ab14d48328703872901ceb5dc060547de8 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 21 Jun 2024 10:51:28 +1000 Subject: [PATCH 13/14] Allow propagation of dual numbers in _get_value_... --- src/Distributed/DistributedDiscreteGeometries.jl | 3 +-- src/LevelSetCutters/DiscreteGeometries.jl | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index c8a949e..5346328 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -24,8 +24,7 @@ function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) φh_data = CellData.get_data(φh) - space = FESpaces.get_fe_space(φh) - T = get_dof_value_type(space) + T = typeof(first(φh_data)(first(first(cell_node_coords)))) values = Vector{T}(undef,num_nodes(own_model)) touched = fill(false,num_nodes(model)) diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index 3eb9705..15824e8 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -76,8 +76,7 @@ function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where # Get cell data φh_data = CellData.get_data(φh) - space = FESpaces.get_fe_space(φh) - T = get_dof_value_type(space) + T = typeof(first(φh_data)(first(first(cell_node_coords)))) # Allow propogation of dual numbers values = Vector{T}(undef,num_nodes(model)) cell_node_coords_cache = array_cache(cell_node_coords) # Loop over cells From 0fcedea8b1993caf00f381d9379502e56f4dbf9e Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 21 Jun 2024 11:44:10 +1000 Subject: [PATCH 14/14] Adjust to use Gridap.Arrays API --- src/Distributed/Distributed.jl | 2 +- src/Distributed/DistributedDiscreteGeometries.jl | 2 +- src/LevelSetCutters/DiscreteGeometries.jl | 2 +- src/LevelSetCutters/LevelSetCutters.jl | 1 + 4 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index dec113c..71d34cd 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -26,7 +26,7 @@ using GridapEmbedded.AgFEM: AggregateCutCellsByThreshold using GridapEmbedded.MomentFittedQuadratures: MomentFitted using Gridap.Geometry: AppendedTriangulation using Gridap.Geometry: get_face_to_parent_face -using Gridap.Arrays: find_inverse_index_map +using Gridap.Arrays: find_inverse_index_map, testitem, return_type using GridapDistributed: DistributedDiscreteModel using GridapDistributed: DistributedTriangulation using GridapDistributed: DistributedFESpace diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 5346328..935e3b9 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -24,7 +24,7 @@ function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) φh_data = CellData.get_data(φh) - T = typeof(first(φh_data)(first(first(cell_node_coords)))) + T = return_type(testitem(CellData.get_data(φh)),testitem(testitem(cell_node_coords))) values = Vector{T}(undef,num_nodes(own_model)) touched = fill(false,num_nodes(model)) diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index 15824e8..23d8cbe 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -76,7 +76,7 @@ function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where # Get cell data φh_data = CellData.get_data(φh) - T = typeof(first(φh_data)(first(first(cell_node_coords)))) # Allow propogation of dual numbers + T = return_type(testitem(CellData.get_data(φh)),testitem(testitem(cell_node_coords))) values = Vector{T}(undef,num_nodes(model)) cell_node_coords_cache = array_cache(cell_node_coords) # Loop over cells diff --git a/src/LevelSetCutters/LevelSetCutters.jl b/src/LevelSetCutters/LevelSetCutters.jl index 75296ea..ac7caf5 100644 --- a/src/LevelSetCutters/LevelSetCutters.jl +++ b/src/LevelSetCutters/LevelSetCutters.jl @@ -22,6 +22,7 @@ using MiniQhull using Gridap.TensorValues using Gridap.ReferenceFEs using Gridap.Arrays +using Gridap.Arrays: testitem, return_type using Gridap.Fields using Gridap.Helpers using Gridap.Geometry