Merge branch 'master' of https://github.com/gridap/GridapEmbedded.jl …

…into distributed_agfem_with_algoim_quads

NEWS.md

@@ -8,8 +8,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added support for distributed level-set geometries. Since PR[#99](https://github.com/gridap/GridapEmbedded.jl/pull/99).
+
+## [0.9.5] - 2024-10-18
+
+### Added
+
 - Adding `compute_redistribute_weights` and `compute_adaptive_flags` functions for load balancing and adaptive mesh refinement, respectively. Since PR [#95](https://github.com/gridap/GridapEmbedded.jl/pull/95).
 
+### Changed
+
+- Updated to Algoim v0.2.2, which runs on Julia 1.11. Since PR [#97](https://github.com/gridap/GridapEmbedded.jl/pull/97).
 
 ## [0.9.4] - 2024-07-09
 

Project.toml

@@ -1,7 +1,7 @@
 name = "GridapEmbedded"
 uuid = "8838a6a3-0006-4405-b874-385995508d5d"
 authors = ["Francesc Verdugo <[email protected]>", "Eric Neiva <[email protected]>", "Pere Antoni Martorell <[email protected]>", "Santiago Badia <[email protected]>"]
-version = "0.9.4"
+version = "0.9.5"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -22,9 +22,9 @@ algoimWrapper_jll = "3c43aa7b-5398-51f3-8d75-8f051e6faa4d"
 
 [compat]
 AbstractTrees = "0.3.3, 0.4"
-Algoim = "0.2"
+Algoim = "0.2.2"
 Combinatorics = "1"
-CxxWrap = "0.14"
+CxxWrap = "0.16"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 Gridap = "0.17, 0.18"
 GridapDistributed = "0.3, 0.4"

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)

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)
+# 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...)
@@ -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,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
@@ -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)

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,41 +61,21 @@ 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
 
+# 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
+end