diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index c68af9a..ee2ea62 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 PartitionedArrays: VectorFromDict @@ -27,6 +28,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, testitem, return_type using Gridap.FESpaces: FESpaceWithLinearConstraints using Gridap.FESpaces: _dof_to_DOF, _DOF_to_dof using GridapDistributed: DistributedDiscreteModel @@ -47,6 +49,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, DiscreteGeometry import Gridap.Geometry: Triangulation import Gridap.Geometry: SkeletonTriangulation import Gridap.Geometry: BoundaryTriangulation @@ -56,6 +59,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..0e1b602 --- /dev/null +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -0,0 +1,131 @@ +struct DistributedDiscreteGeometry{A} <: GridapType + geometries::A +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 + own_model = remove_ghost_cells(model,gids) + 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)) + + 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 + 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(loc_φ,point_to_coords;name) + end + DistributedDiscreteGeometry(geometries) +end + +function distributed_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 08816f3..e7a8aa1 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,8 +17,13 @@ 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) + geometries = map(get_geometry,local_views(a)) + distributed_geometry(geometries) +end + +# Needed for dispatching between analytical geometries and discrete geometries +function distributed_geometry(geometries::AbstractArray{<:CSG.Geometry}) + PartitionedArrays.getany(geometries) end function cut(bgmodel::DistributedDiscreteModel,args...) @@ -123,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) @@ -239,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 @@ -252,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 @@ -262,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 @@ -285,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 @@ -301,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 5289299..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) @@ -71,3 +67,15 @@ function DiscreteGeometry( tree = Leaf(data) 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_free_dof_values(φh) + point_to_coords = collect1d(get_node_coordinates(model)) + DiscreteGeometry(point_to_value,point_to_coords;name) +end diff --git a/src/LevelSetCutters/LevelSetCutters.jl b/src/LevelSetCutters/LevelSetCutters.jl index aa5b147..ac7caf5 100644 --- a/src/LevelSetCutters/LevelSetCutters.jl +++ b/src/LevelSetCutters/LevelSetCutters.jl @@ -22,12 +22,14 @@ 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 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..26b34ba --- /dev/null +++ b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl @@ -0,0 +1,67 @@ +module PeriodicAgFEMSpacesTests + +using Test +using Gridap +using GridapEmbedded +using Gridap.Geometry: get_active_model + +const R = 0.55 +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 = 1 + +# 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-7 +@test sum( ∫(abs2(v-vhagg))dΩ ) < tol +@test sum( ∫(abs2(v-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-7 +@test sum( ∫(abs2(v-vhagg))dΩ ) < tol +@test sum( ∫(abs2(v-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/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..10a9959 --- /dev/null +++ b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,143 @@ +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 + +end # module \ No newline at end of file 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/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 fe600b7..d0510c0 100644 --- a/test/DistributedTests/mpi/runtests_body.jl +++ b/test/DistributedTests/mpi/runtests_body.jl @@ -6,6 +6,9 @@ using MPI include("../PoissonTests.jl") include("../AggregatesTests.jl") +include("../DistributedDiscreteGeometryPoissonTest.jl") +include("../DistributedLSDiscreteGeometryPoissonTest.jl") +include("../PeriodicDistributedDiscreteGeometryPoissonTest.jl") if ! MPI.Initialized() MPI.Init() @@ -21,6 +24,21 @@ 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") + + PArrays.tic!(t) + 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/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl new file mode 100644 index 0000000..080d4b9 --- /dev/null +++ b/test/DistributedTests/sequential/DistributedDiscreteGeometryPoissonTest.jl @@ -0,0 +1,8 @@ +module DistributedDiscreteGeometryPoissonTestsSeq +using PartitionedArrays +include("../DistributedDiscreteGeometryPoissonTest.jl") +with_debug() do distribute + DistributedDiscreteGeometryPoissonTest.main(distribute,(2,2)) + DistributedDiscreteGeometryPoissonTest.main(distribute,(4,1),cells=(12,12),geometry=:remotes) +end +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/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 e68831c..57910a8 100644 --- a/test/DistributedTests/sequential/runtests.jl +++ b/test/DistributedTests/sequential/runtests.jl @@ -3,5 +3,10 @@ module SequentialTests 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 + +## 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 1576bc4..f6d52df 100644 --- a/test/GridapEmbeddedTests/runtests.jl +++ b/test/GridapEmbeddedTests/runtests.jl @@ -6,6 +6,8 @@ using Test @time @testset "PoissonAgFEM" begin include("PoissonAgFEMTests.jl") end +@time @testset "PeriodicPoissonAgFEM" begin include("PeriodicPoissonAgFEMTests.jl") end + @time @testset "PoissonModalC0AgFEM" begin include("PoissonModalC0AgFEMTests.jl") end @time @testset "BimaterialPoissonCutFEM" begin include("BimaterialPoissonCutFEMTests.jl") 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