Refactored for single-valued objective function

Author: ericneiva

src/AgFEM/CellAggregation.jl

@@ -78,7 +78,7 @@ function aggregate(
 
   facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut.bgmodel,geo)
   lid_to_gid = 1:num_cells(get_background_model(cut))
-  cell_to_cellin,_,_,all_aggregated =
+  cell_to_cellin,_,all_aggregated =
     _aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut,lid_to_gid)
   @assert all_aggregated "Not all cells were aggregated"
   cell_to_cellin
@@ -92,14 +92,14 @@ function aggregate(
   in_or_out)
 
   facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut.bgmodel,geo)
-  cell_to_cellin, cell_to_path_length, cell_to_bbox_diam, _ =
+  cell_to_cellin, cell_to_value, _ =
     _aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut,lid_to_gid)
   for i in eachindex(cell_to_cellin)
     if !iszero(cell_to_cellin[i])
       cell_to_cellin[i] = lid_to_gid[cell_to_cellin[i]]
     end
   end
-  cell_to_cellin, cell_to_path_length, cell_to_bbox_diam
+  cell_to_cellin, cell_to_value
 end
 
 function aggregate(
@@ -111,7 +111,7 @@ function aggregate(
 
   facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut_facets,geo)
   lid_to_gid = 1:num_cells(get_background_model(cut))
-  cell_to_cellin,_,_,all_aggregated =
+  cell_to_cellin,_,all_aggregated =
     _aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut,lid_to_gid)
   @assert all_aggregated "Not all cells were aggregated"
   cell_to_cellin
@@ -125,7 +125,7 @@ function aggregate(
   facet_to_inoutcut::AbstractVector)
 
   lid_to_gid = 1:num_cells(get_background_model(cut))
-  cell_to_cellin,_,_,all_aggregated =
+  cell_to_cellin,_,all_aggregated =
     _aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut,lid_to_gid)
   @assert all_aggregated "Not all cells were aggregated"
   cell_to_cellin
@@ -161,8 +161,7 @@ function _aggregate_by_threshold_barrier(
 
   n_cells = length(cell_to_unit_cut_meas)
   cell_to_cellin = zeros(Int32,n_cells)
-  cell_to_path_length = zeros(Int8,n_cells)
-  cell_to_bbox_diam = zeros(Float64,n_cells)
+  cell_to_value = zeros(Int128,n_cells)
   cell_to_touched = fill(false,n_cells)
 
   for cell in 1:n_cells
@@ -184,7 +183,7 @@ function _aggregate_by_threshold_barrier(
     all_aggregated = true
     for cell in 1:n_cells
       if ! cell_to_touched[cell] && cell_to_inoutcut[cell] == CUT
-        neigh_cell, bbox_diam = _find_best_neighbor(
+        neigh_cell, value = _find_best_neighbor(
           c1,c2,c3,c4,cell,
           cell_to_faces,
           face_to_cells,
@@ -193,12 +192,12 @@ function _aggregate_by_threshold_barrier(
           cell_to_cellin,
           lid_to_gid,
           facet_to_inoutcut,
+          iter,
           loc)
         if neigh_cell > 0
           cellin = cell_to_cellin[neigh_cell]
           cell_to_cellin[cell] = cellin
-          cell_to_path_length[cell] = iter
-          cell_to_bbox_diam[cell] = bbox_diam
+          cell_to_value[cell] = value
         else
           all_aggregated = false
         end
@@ -210,7 +209,7 @@ function _aggregate_by_threshold_barrier(
     _touch_aggregated_cells!(cell_to_touched,cell_to_cellin)
   end
 
-  cell_to_cellin, cell_to_path_length, cell_to_bbox_diam, all_aggregated
+  cell_to_cellin, cell_to_value, all_aggregated
 end
 
 function _find_best_neighbor(
@@ -222,13 +221,13 @@ function _find_best_neighbor(
   cell_to_cellin,
   lid_to_gid,
   facet_to_inoutcut,
+  iter,
   loc)
 
   faces = getindex!(c1,cell_to_faces,cell)
-  min_bb_diam = Inf
   T = eltype(eltype(face_to_cells))
   best_neigh_cell = zero(T)
-  min_cellin = zero(T)
+  min_value = typemax(Int128)
   i_to_coords = getindex!(c3,cell_to_coords,cell)
   for face in faces
     inoutcut = facet_to_inoutcut[face]
@@ -247,19 +246,19 @@ function _find_best_neighbor(
             bb_diam = max(bb_diam,Float64(norm(p-q)))
           end
         end
-        if ( min_bb_diam == Inf ) |
-           ( ( bb_diam<min_bb_diam ) &
-             !isapprox(bb_diam,min_bb_diam,atol=1.0e-9) ) |
-           ( isapprox(bb_diam,min_bb_diam,atol=1.0e-9) &
-             (lid_to_gid[cellin] < min_cellin ) )
-          min_bb_diam = bb_diam
+        # Evaluate objective function and update best values
+        bit_iter    = bitstring(Int8(iter))
+        bit_bb_diam = bitstring(Float32(round(bb_diam,sigdigits=6)))
+        bit_cellin  = bitstring(lid_to_gid[cellin]) # Assuming ≤ 64-bit
+        value = parse(Int128,bit_iter*bit_bb_diam*bit_cellin;base=2)
+        if value < min_value
           best_neigh_cell = neigh_cell
-          min_cellin = lid_to_gid[cellin]
+          min_value = value
         end
       end
     end
   end
-  best_neigh_cell, min_bb_diam
+  best_neigh_cell, min_value
 end
 
 function _touch_aggregated_cells!(cell_to_touched,cell_to_cellin)
@@ -370,7 +369,7 @@ function aggregate(bgtrian,cell_to_is_active,cell_to_is_cut,in_or_out)
 
   threshold = 1.0
   lid_to_gid = 1:num_cells(model)
-  cell_to_cellin,_,_,all_aggregated = _aggregate_by_threshold_barrier(
+  cell_to_cellin,_,all_aggregated = _aggregate_by_threshold_barrier(
     threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
     in_or_out,cell_to_coords,cell_to_faces,face_to_cells,lid_to_gid)
   @assert all_aggregated "Not all cells were aggregated"
@@ -401,7 +400,7 @@ function aggregate_narrow_band(bgtrian,cell_to_is_in_narrow,cell_to_is_active,ce
 
   threshold = 1.0
   lid_to_gid = 1:num_cells(model)
-  cell_to_cellin,_,_,all_aggregated = _aggregate_by_threshold_barrier(
+  cell_to_cellin,_,all_aggregated = _aggregate_by_threshold_barrier(
     threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
     in_or_out,cell_to_coords,cell_to_faces,face_to_cells,lid_to_gid)
   @assert all_aggregated "Not all cells were aggregated"

test/dev/distributed_aggregation.jl

@@ -6,7 +6,29 @@ using PartitionedArrays
 
 using GridapEmbedded: aggregate
 
-function generate_dumbell(ranks, np, nc, ng)
+function generate_asymmetric_kettlebell(ranks, np, nc, ng)
+  L = 1
+  p0 = Point(0.0,0.0)
+  pmin = Point(-L/4,-L/8)
+  pmax = Point(L/4,3*L/8)
+
+  R = 0.2
+  x = VectorValue(1.0,0.0)
+  h = L/(2*nc)
+
+  geo1 = disk(R,x0=p0)
+  geo2 = quadrilateral(;x0=Point(-L/4+h/4,5*h/4),
+                        d1=VectorValue(6*h+h/2,0.0),
+                        d2=VectorValue(0.0,4*h+h/2))
+  geo3 = quadrilateral(;x0=Point(-L/4+3*h/4,7*h/4),
+                        d1=VectorValue(6*h-h/2,0.0),
+                        d2=VectorValue(0.0,4*h-h/2))
+  geo = union(geo1,intersect(geo2,!geo3))
+  bgmodel = CartesianDiscreteModel(ranks,np,pmin,pmax,(nc,nc);ghost=(ng,ng))
+  return bgmodel, geo
+end
+
+function generate_symmetric_kettlebell(ranks, np, nc, ng)
   L = 1
   p0 = Point(0.0,0.0)
   pmin = Point(-L/4,-L/8)
@@ -43,59 +65,42 @@ function exchange_impl!(vector_partition,cache)
   return t
 end
 
-function find_optimal_roots!(lcell_to_root,lcell_to_path_length,lcell_to_bbox_diam,cell_indices)
+function find_optimal_roots!(lcell_to_root,lcell_to_value,cell_indices)
   # Bring all root candidates to the owner of the cut cell
   roots_cache = PartitionedArrays.p_vector_cache(lcell_to_root,cell_indices)
   t1 = exchange_impl!(lcell_to_root,roots_cache)
-  lengths_cache = PartitionedArrays.p_vector_cache(lcell_to_path_length,cell_indices)
-  t2 = exchange_impl!(lcell_to_path_length,lengths_cache)
-  bb_diams_cache = PartitionedArrays.p_vector_cache(lcell_to_bbox_diam,cell_indices)
-  t3 = exchange_impl!(lcell_to_bbox_diam,bb_diams_cache)
+  values_cache = PartitionedArrays.p_vector_cache(lcell_to_value,cell_indices)
+  t2 = exchange_impl!(lcell_to_value,values_cache)
   lcell_to_owner = map(copy∘local_to_owner,cell_indices)
   owners_cache = PartitionedArrays.p_vector_cache(lcell_to_owner,cell_indices)
   wait(t1)
   wait(t2)
-  wait(t3)
 
   # Select the optimal root for each local cut cell
   map(
-    lcell_to_root,lcell_to_path_length,lcell_to_bbox_diam,
-      lcell_to_owner,roots_cache,lengths_cache,bb_diams_cache
-  ) do lcell_to_root,lcell_to_path_length,lcell_to_bbox_diam,
-      lcell_to_owner,roots_cache,lengths_cache,bb_diams_cache
-
+    lcell_to_root,lcell_to_value,lcell_to_owner,roots_cache,values_cache
+  ) do lcell_to_root, lcell_to_value, lcell_to_owner, roots_cache, values_cache
     lids_rcv = roots_cache.local_indices_rcv
-    lengths_rcv = lengths_cache.buffer_rcv
-    bb_diams_rcv = bb_diams_cache.buffer_rcv
+    values_rcv = values_cache.buffer_rcv
     roots_rcv = roots_cache.buffer_rcv
-
-    for (k,nbor) in enumerate(roots_cache.neighbors_rcv)
-      for (lcell,root,len,bb_diam) in zip(lids_rcv[k],roots_rcv[k],lengths_rcv[k],bb_diams_rcv[k])
+    for (k, nbor) in enumerate(roots_cache.neighbors_rcv)
+      for (lcell,root,value) in zip(lids_rcv[k], roots_rcv[k], values_rcv[k])
         root == 0 && continue
-        if lcell_to_path_length[lcell] > len
+        if lcell_to_value[lcell] > value # Take the minimum
           lcell_to_root[lcell] = root
+          lcell_to_value[lcell] = value
           lcell_to_owner[lcell] = nbor
-        elseif lcell_to_path_length[lcell] == len
-          if ( lcell_to_bbox_diam[lcell] > bb_diam ) &
-              !isapprox(lcell_to_bbox_diam[lcell],bb_diam,atol=1.0e-9)
-            lcell_to_root[lcell] = root
-            lcell_to_owner[lcell] = nbor
-          elseif isapprox(lcell_to_bbox_diam[lcell],bb_diam,atol=1.0e-9) &
-              ( lcell_to_root[lcell] > root )
-            lcell_to_root[lcell] = root
-            lcell_to_owner[lcell] = nbor
-          end
         end
       end
     end
-
   end
 
   # Scatter the optimal roots and values
   # Technically, we do not need to scatter the values, it can be removed after 
   # we are done debugging
   t1 = consistent!(PVector(lcell_to_owner,cell_indices,owners_cache))
   t2 = consistent!(PVector(lcell_to_root,cell_indices,roots_cache))
+  t3 = consistent!(PVector(lcell_to_value,cell_indices,values_cache))
 
   wait(t1)
 
@@ -104,29 +109,30 @@ function find_optimal_roots!(lcell_to_root,lcell_to_path_length,lcell_to_bbox_di
   end
 
   wait(t2)
+  wait(t3)
 
-  return lcell_to_root, lcell_to_owner, agg_cell_indices
+  return lcell_to_root, lcell_to_owner, lcell_to_value, agg_cell_indices
 end
 
 distribute = PartitionedArrays.DebugArray
 np = (3,1)
 ranks = distribute(LinearIndices((prod(np),)))
-bgmodel, geo = generate_dumbell(ranks, np, 9, 2)
+bgmodel, geo = generate_symmetric_kettlebell(ranks, np, 9, 2)
 cutgeo = cut(bgmodel, geo)
 
 cell_indices = partition(get_cell_gids(bgmodel))
 
 strategy = AggregateCutCellsByThreshold(1.0)
-lcell_to_root, lcell_to_path_length, lcell_to_bbox_diam =
+lcell_to_root, lcell_to_value =
   map(local_views(cutgeo),cell_indices) do cutgeo,cell_indices
     lid_to_gid = local_to_global(cell_indices)
     aggregate(strategy,cutgeo,geo,lid_to_gid,IN)
   end |> tuple_of_arrays
 
 lcell_to_inconsistent_root = map(copy,lcell_to_root)
 
-lcell_to_root, lcell_to_owner, agg_cell_indices =
-  find_optimal_roots!(lcell_to_root,lcell_to_path_length,lcell_to_bbox_diam,cell_indices);
+lcell_to_root, lcell_to_owner, lcell_to_value, agg_cell_indices =
+  find_optimal_roots!(lcell_to_root,lcell_to_value,cell_indices);
 
 # Creating aggregate-conforming cell partition
 
@@ -137,28 +143,22 @@ end
 
 writevtk(EmbeddedBoundary(cutgeo),"data/bnd");
 writevtk(
-  Triangulation(bgmodel), "data/dumbell_aggregates", 
+  Triangulation(bgmodel), "data/kettlebell_aggregates", 
   celldata = ["aggregate" => ocell_to_root],
 );
 
 map(ranks,
     local_views(bgmodel),
     lcell_to_root,
     lcell_to_inconsistent_root,
-    lcell_to_owner,
-    lcell_to_path_length,
-    lcell_to_bbox_diam) do r,bgmodel,
-                          lcell_to_root,
-                          lcell_to_iroot,
-                          lcell_to_owner,
-                          lcell_to_path_length,
-                          lcell_to_bbox_diam
+    lcell_to_owner) do r,bgmodel,
+                       lcell_to_root,
+                       lcell_to_iroot,
+                       lcell_to_owner
   writevtk(
-    Triangulation(bgmodel), "data/dumbell_aggregates_$(r)", 
+    Triangulation(bgmodel), "data/kettlebell_aggregates_$(r)", 
     celldata = [ "roots"              => lcell_to_root,
                  "inconsistent roots" => lcell_to_iroot,
-                 "owners"             => lcell_to_owner,
-                 "length"             => lcell_to_path_length,
-                 "diameter"           => lcell_to_bbox_diam ],
+                 "owners"             => lcell_to_owner ],
   );
 end