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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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 From 5af2a2b9e7b762576dbb0dca05a64cb2bfcb165d Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 11 Nov 2024 15:32:35 +1100 Subject: [PATCH 15/16] Revised some of the changes --- .../DistributedDiscreteGeometries.jl | 57 +++++++------------ src/Distributed/DistributedDiscretizations.jl | 43 +++++++------- src/LevelSetCutters/DiscreteGeometries.jl | 38 +++---------- 3 files changed, 51 insertions(+), 87 deletions(-) diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 935e3b9..0e1b602 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -4,41 +4,28 @@ end local_views(a::DistributedDiscreteGeometry) = a.geometries +# TODO: Is this really necessary? function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) where {Dc,Dp} @assert DomainStyle(φh) == ReferenceDomain() gids = get_cell_gids(model) 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 = 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) - # 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) - T = return_type(testitem(CellData.get_data(φh)),testitem(testitem(cell_node_coords))) + own_cells = get_face_to_parent_face(own_model,Dc) + + trian = get_triangulation(φh) + cell_points = get_cell_points(trian) + cell_ids = get_cell_node_ids(own_model) + cell_values = φh(cell_points) + + T = eltype(testitem(cell_values)) 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 + cell_ids_cache = array_cache(cell_ids) + cell_values_cache = array_cache(cell_values) + for (ocell,cell) in enumerate(own_cells) + ids = getindex!(cell_ids_cache,cell_ids,ocell) + vals = getindex!(cell_values_cache,cell_values,cell) + values[ids] .= vals end return values end @@ -56,7 +43,7 @@ function DiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::S DistributedDiscreteGeometry(geometries) end -function get_geometry(a::AbstractArray{<:DiscreteGeometry}) +function distributed_geometry(a::AbstractArray{<:DiscreteGeometry}) DistributedDiscreteGeometry(a) end @@ -104,8 +91,8 @@ function distributed_embedded_triangulation( T, cutgeo::DistributedEmbeddedDiscretization, cutinorout, - geom::DistributedDiscreteGeometry) - + geom::DistributedDiscreteGeometry +) trians = map(local_views(cutgeo),local_views(geom)) do lcutgeo,lgeom T(lcutgeo,cutinorout,lgeom) end @@ -117,8 +104,8 @@ function distributed_aggregate( strategy::AggregateCutCellsByThreshold, cut::DistributedEmbeddedDiscretization, geo::DistributedDiscreteGeometry, - in_or_out=IN) - + 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) @@ -133,8 +120,8 @@ end function compute_bgfacet_to_inoutcut( cutter::Cutter, bgmodel::DistributedDiscreteModel, - geo::DistributedDiscreteGeometry) - + 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) diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl index 9ae968c..8bb2d6e 100644 --- a/src/Distributed/DistributedDiscretizations.jl +++ b/src/Distributed/DistributedDiscretizations.jl @@ -3,11 +3,12 @@ struct DistributedEmbeddedDiscretization{A,B} <: GridapType discretizations::A model::B function DistributedEmbeddedDiscretization( - d::AbstractArray{<:AbstractEmbeddedDiscretization}, - model::DistributedDiscreteModel) - A = typeof(d) + discretizations::AbstractArray{<:AbstractEmbeddedDiscretization}, + model::DistributedDiscreteModel + ) + A = typeof(discretizations) B = typeof(model) - new{A,B}(d,model) + new{A,B}(discretizations,model) end end @@ -16,12 +17,13 @@ local_views(a::DistributedEmbeddedDiscretization) = a.discretizations get_background_model(a::DistributedEmbeddedDiscretization) = a.model function get_geometry(a::DistributedEmbeddedDiscretization) - loc_geometries = map(get_geometry,local_views(a)) - get_geometry(loc_geometries) + geometries = map(get_geometry,local_views(a)) + distributed_geometry(geometries) end -function get_geometry(a::AbstractArray{<:CSG.Geometry}) - PartitionedArrays.getany(a) +# Neeed for dispatching between analytical geometries and distributed geometries +function distributed_geometry(geometries::AbstractArray{<:CSG.Geometry}) + PartitionedArrays.getany(geometries) end function cut(bgmodel::DistributedDiscreteModel,args...) @@ -127,8 +129,8 @@ end function compute_bgfacet_to_inoutcut( cutter::Cutter, bgmodel::DistributedDiscreteModel, - geo) - + geo +) gids = get_cell_gids(bgmodel) bgf_to_ioc = map(local_views(bgmodel),local_views(gids)) do model,gids ownmodel = remove_ghost_cells(model,gids) @@ -243,12 +245,11 @@ function change_bgmodel( DistributedEmbeddedDiscretization(cuts,model) end - function _change_bgmodels( cutgeo::DistributedEmbeddedDiscretization, model::DistributedDiscreteModel, - cell_to_newcell) - + cell_to_newcell +) map(local_views(cutgeo),local_views(model),cell_to_newcell) do c,m,c_to_nc change_bgmodel(c,m,c_to_nc) end @@ -256,8 +257,8 @@ end function _change_bgmodels( cutgeo::DistributedEmbeddedDiscretization, - model::DistributedDiscreteModel) - + model::DistributedDiscreteModel +) map(local_views(cutgeo),local_views(model)) do c,m change_bgmodel(c,m) end @@ -266,8 +267,8 @@ end function change_bgmodel( cut::EmbeddedDiscretization, newmodel::DiscreteModel, - cell_to_newcell=1:num_cells(get_background_model(cut))) - + cell_to_newcell=1:num_cells(get_background_model(cut)) +) ls_to_bgc_to_ioc = map(cut.ls_to_bgcell_to_inoutcut) do bgc_to_ioc new_bgc_to_ioc = Vector{Int8}(undef,num_cells(newmodel)) new_bgc_to_ioc[cell_to_newcell] = bgc_to_ioc @@ -289,10 +290,9 @@ end function change_bgmodel( cut::EmbeddedFacetDiscretization, newmodel::DiscreteModel, - facet_to_newfacet=1:num_facets(get_background_model(cut))) - + facet_to_newfacet=1:num_facets(get_background_model(cut)) +) nfacets = num_facets(newmodel) - ls_to_bgf_to_ioc = map(cut.ls_to_facet_to_inoutcut) do bgf_to_ioc new_bgf_to_ioc = Vector{Int8}(undef,nfacets) new_bgf_to_ioc[facet_to_newfacet] = bgf_to_ioc @@ -305,7 +305,8 @@ function change_bgmodel( subfacets, cut.ls_to_subfacet_to_inout, cut.oid_to_ls, - cut.geo) + cut.geo + ) end function change_bgmodel(cells::SubCellData,cell_to_newcell) diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index 23d8cbe..472fa27 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -34,7 +34,6 @@ function discretize(a::AnalyticalGeometry,point_to_coords::AbstractArray{<:Point end function discretize(a::AnalyticalGeometry,point_to_coords::Vector{<:Point}) - tree = get_tree(a) j_to_fun, oid_to_j = _find_unique_leaves(tree) j_to_ls = [ fun.(point_to_coords) for fun in j_to_fun ] @@ -48,13 +47,10 @@ function discretize(a::AnalyticalGeometry,point_to_coords::Vector{<:Point}) end newtree = replace_data(identity,conversion,tree) - DiscreteGeometry(newtree,point_to_coords) - end function _find_unique_leaves(tree) - i_to_fun = map(n->first(n.data),collect(Leaves(tree))) i_to_oid = map(objectid,i_to_fun) j_to_oid = unique(i_to_oid) @@ -65,31 +61,6 @@ 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) - 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 - 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 - return values -end - function DiscreteGeometry( point_to_value::AbstractVector,point_to_coords::AbstractVector;name::String="") data = (point_to_value,name,nothing) @@ -97,9 +68,14 @@ function DiscreteGeometry( DiscreteGeometry(tree,point_to_coords) end +# TODO: This assumes that the level set φh is 1st order, i.e that there is a 1-to-1 correspondence +# between nodes in the mesh and dofs in φh. +# Even if we allowed higher order, the cuts are always linear. Not only it would be a waste +# of time to use higher order, but cuts could actually be wrong. +# This might be developped in the future. function DiscreteGeometry( φh::CellField,model::DiscreteModel;name::String="") - point_to_value = _get_value_at_coords(φh,model) + point_to_value = get_free_dof_values(φh) point_to_coords = collect1d(get_node_coordinates(model)) DiscreteGeometry(point_to_value,point_to_coords;name) -end \ No newline at end of file +end From f04cbcbb90995be46d8b45aa5364917e45bc9e7e Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 11 Nov 2024 15:38:24 +1100 Subject: [PATCH 16/16] Minor --- src/Distributed/DistributedDiscretizations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl index 8bb2d6e..e7a8aa1 100644 --- a/src/Distributed/DistributedDiscretizations.jl +++ b/src/Distributed/DistributedDiscretizations.jl @@ -21,7 +21,7 @@ function get_geometry(a::DistributedEmbeddedDiscretization) distributed_geometry(geometries) end -# Neeed for dispatching between analytical geometries and distributed geometries +# Needed for dispatching between analytical geometries and discrete geometries function distributed_geometry(geometries::AbstractArray{<:CSG.Geometry}) PartitionedArrays.getany(geometries) end