Distributed DiscreteGeometries

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")

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

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,21 +245,20 @@ 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
 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)

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

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

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

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