Implemented compute FE distance functions for algoim distributed fields

Author: ericneiva

.gitignore

@@ -8,5 +8,6 @@
 Manifest.toml
 LocalPreferences.toml
 .vscode/
+data/
 *.swp
 data

src/AlgoimUtils/AlgoimUtils.jl

@@ -590,8 +590,7 @@ function compute_closest_point_projections(model::OctreeDistributedDiscreteModel
   fgrid = vec(map(i->xmin+Point((Tuple(i).-1).*fsizes),
               CartesianIndices((fcells...,))))
   fvals = φ.(fgrid)
-  @show length(fvals)
-
+  
   # Coordinates of free DoFs on each local partition
   coords = interpolate_everywhere(identity,V)
   coords_vals = get_free_dof_values(coords)
@@ -604,14 +603,15 @@ function compute_closest_point_projections(model::OctreeDistributedDiscreteModel
   FEFunction(V,PVector(cpps,partition(V.gids)))
 end
 
-function compute_closest_point_projections(model::OctreeDistributedDiscreteModel,
-                                           V::DistributedFESpace,
-                                           φ::DistributedAlgoimCallLevelSetFunction,
-                                           order::Int,
-                                           max_refinement_level::Int;
-                                           cppdegree::Int=2,
-                                           trim::Bool=false,
-                                           limitstol::Float64=1.0e-8)
+function compute_closest_point_projections(
+    model::OctreeDistributedDiscreteModel,
+    V::DistributedFESpace,
+    φ::DistributedAlgoimCallLevelSetFunction,
+    order::Int,
+    max_refinement_level::Int;
+    cppdegree::Int=2,
+    trim::Bool=false,
+    limitstol::Float64=1.0e-8)
 
   # Uniform partition at the maximum refinement level
   coarse_desc = get_cartesian_descriptor(model.coarse_model)
@@ -624,11 +624,10 @@ function compute_closest_point_projections(model::OctreeDistributedDiscreteModel
   # The level set values at the maximum refinement level
   fsizes = 
     coarse_desc.sizes ./ ( Int32(2^max_refinement_level) * Int32(order) )
-  fcells = fine_partition.+1
+  fcells = fine_partition .+ 1
   fgrid = vec(map(i->xmin+Point((Tuple(i).-1).*fsizes),
               CartesianIndices((fcells...,))))
-
-  # This is a hack to evaluate the LS function on 2D  
+  # A hack follows to evaluate the LS function on 2D  
   # OctreeDistributedDiscreteModels, as the current
   # implementation of Interpolable is not using the
   # NonConformingGridapTopology, so the list of cells
@@ -1127,6 +1126,34 @@ function compute_distance_fe_function(
   FEFunction(fespace_scalar_type,dists)
 end
 
+function compute_distance_fe_function(
+    bgmodel::OctreeDistributedDiscreteModel,
+    fespace_scalar_type::DistributedFESpace,
+    fespace_vector_type::DistributedFESpace,
+    φ::DistributedAlgoimCallLevelSetFunction,
+    order::Int,
+    max_refinement_level::Int;
+    cppdegree::Int=2)
+  
+  cps = compute_closest_point_projections(
+    bgmodel,fespace_vector_type,φ,order,
+    max_refinement_level,cppdegree=cppdegree)
+  cps_vals = get_free_dof_values(cps)
+  
+  coords = interpolate_everywhere(identity,fespace_vector_type)
+  coords_vals = get_free_dof_values(coords)
+
+  _dists = map(local_views(fespace_scalar_type),
+               local_views(cps_vals),
+               local_views(coords_vals),
+               local_views(φ)) do fs,cp,cos,φl
+    _compute_signed_distance(φl,cp,cos,num_free_dofs(fs),
+      num_dims(bgmodel.coarse_model))
+  end
+  dists = PVector(_dists,partition(fespace_scalar_type.gids))
+  FEFunction(fespace_scalar_type,dists)
+end
+
 function compute_distance_fe_function(
     fespace_scalar_type::FESpace,
     fespace_vector_type::FESpace,

test/dev/cpp_dev.jl

@@ -57,32 +57,37 @@ with_mpi() do distribute
 
   max_refinement_level = num_levels_initial_refinement
 
-  # # RMK: Algoim CPP algorithms work on uniform grids
-  # # In order to use them on non-uniform grids, we work
-  # # on an upper bound grid, corresponding to the uniform 
-  # # mesh obtained at maximum refinement level. On this 
-  # # maximal grid, the working arrays for the CPP are
-  # # computed using the coordinates of the grid points 
-  # # and its level set values.
-  # cpps = compute_closest_point_projections(
-  #   fmodel,Vₕ,φ_fun,order,max_refinement_level,cppdegree=3)
-  # dists = compute_distance_fe_function(
-  #   fmodel,Qₕ,Vₕ,φ_fun,order,max_refinement_level,cppdegree=3)
-
-  # writevtk(Ω,"Ω",cellfields=["cpp"=>cpps,"dist"=>dists],nsubcells=3)
+  # RMK: Algoim CPP algorithms work on uniform grids
+  # In order to use them on non-uniform grids, we work
+  # on an upper bound grid, corresponding to the uniform 
+  # mesh obtained at maximum refinement level. On this 
+  # maximal grid, the working arrays for the CPP are
+  # computed using the coordinates of the grid points 
+  # and its level set values.
+  cpps = compute_closest_point_projections(
+    fmodel,Vₕ,φ_fun,order,max_refinement_level)
+  dists = compute_distance_fe_function(
+    fmodel,Qₕ,Vₕ,φ_fun,order,max_refinement_level)
+
+  writevtk(Ω,"Ω_fun",cellfields=["cpp"=>cpps,"dist"=>dists],nsubcells=3)
 
   _φ = interpolate_everywhere(val,Qₕ)
   φ_field = AlgoimCallLevelSetFunction(_φ,∇(_φ))
-  sm = Gridap.CellData.KDTreeSearch(num_nearest_vertices=3)
+
+  sm = Gridap.CellData.KDTreeSearch(num_nearest_vertices=5)
   iφ_field = map(local_views(φ_field.values)) do iφ
     Gridap.CellData.Interpolable(iφ,searchmethod=sm)
   end |> GridapDistributed.DistributedInterpolable
+  
   cpps = compute_closest_point_projections(
     fmodel,Vₕ,φ_field,order,max_refinement_level)
-  writevtk(Ω,"Ω",cellfields=["cpp"=>cpps,"phi"=>iφ_field∘cpps],nsubcells=3)
+  dists = compute_distance_fe_function(
+    fmodel,Qₕ,Vₕ,φ_field,order,max_refinement_level)
+
+  writevtk(Ω,"Ω_field",cellfields=["cpp"=>cpps,"dist"=>dists,"phi"=>iφ_field∘cpps],nsubcells=1)
 
-  # [TODO]: Compute distance function
-  # [TODO]: Check visualisation of cpp and distance function
+  # [TODO]: Visualization for order >= 3 and num_subcells > 1 is strange.
+  #         This might be due to oscillations of the polynomial, but not sure.
 
   true
 end