Merge pull request #90 from zjwegert/distributed_discrete_geometry

Support for distributed discrete geometries

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,144 @@
+struct DistributedDiscreteGeometry{A} <: GridapType
+  geometries::A
+end
+
+local_views(a::DistributedDiscreteGeometry) = a.geometries
+
+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)))
+    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(loc_φ,point_to_coords;name)
+  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

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

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

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