Merge pull request #78 from gridap/optimise_algoimutils

Optimise computation of normal displacement in dynamic simulations

NEWS.md

@@ -6,9 +6,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## Unreleased
 
+### Changed
+
+- Optimise computation of normal displacement in dynamic simulations. Since PR [#78](https://github.com/gridap/GridapEmbedded.jl/pull/78).
+
 ### Fixed
 
-- Compute signed distance FE function in 3D [#77](https://github.com/gridap/GridapEmbedded.jl/pull/77).
+- Compute signed distance FE function in 3D. Since PR [#77](https://github.com/gridap/GridapEmbedded.jl/pull/77).
 
 ## [0.9.0] - 2024-01-22
 

src/AlgoimUtils/AlgoimUtils.jl

@@ -41,6 +41,7 @@ export fill_cpp_data
 export fill_cpp_data_raw
 export compute_closest_point_projections
 export compute_normal_displacement
+export compute_normal_displacement!
 export compute_distance_fe_function
 export delaunaytrian
 export convexhull
@@ -323,6 +324,40 @@ function compute_normal_displacement(
   disps
 end
 
+function compute_normal_displacement!(
+    cell_to_points,
+    cps::AbstractVector{<:Point},
+    phi::AlgoimCallLevelSetFunction,
+    fun,
+    dt::Float64,
+    Ω::Triangulation)
+  # Note that cps must be (scalar) DoF-numbered, not lexicographic-numbered
+  if isnothing(cell_to_points)
+    searchmethod = KDTreeSearch()
+    cache1 = _point_to_cell_cache(searchmethod,Ω)
+    x_to_cell(x) = _point_to_cell!(cache1, x)
+    point_to_cell = lazy_map(x_to_cell, cps)
+    cell_to_points, _ = make_inverse_table(point_to_cell, num_cells(Ω))
+  end
+  cell_to_xs = lazy_map(Broadcasting(Reindex(cps)), cell_to_points)
+  cell_point_xs = CellPoint(cell_to_xs, Ω, PhysicalDomain())
+  fun_xs = evaluate(fun,cell_point_xs)
+  nΓ_xs = evaluate(normal(phi,Ω),cell_point_xs)
+  cell_point_disp = lazy_map(Broadcasting(⋅),fun_xs,nΓ_xs)
+  cache_vals = array_cache(cell_point_disp)
+  cache_ctop = array_cache(cell_to_points)
+  disps = zeros(Float64,length(cps))
+  for cell in 1:length(cell_to_points)
+    pts = getindex!(cache_ctop,cell_to_points,cell)
+    vals = getindex!(cache_vals,cell_point_disp,cell)
+    for (i,pt) in enumerate(pts)
+      val = vals[i]
+      disps[pt] = dt * val
+    end
+  end
+  disps, cell_to_points
+end
+
 signed_distance(φ::Function,x,y) = sign(φ(y))*norm(x-y)
 
 signed_distance(φ::T,x,y) where {T<:Number} = sign(φ)*norm(x-y)

src/Exports.jl

@@ -59,6 +59,7 @@ end
 @publish AlgoimUtils fill_cpp_data_raw
 @publish AlgoimUtils compute_closest_point_projections
 @publish AlgoimUtils compute_normal_displacement
+@publish AlgoimUtils compute_normal_displacement!
 @publish AlgoimUtils compute_distance_fe_function
 @publish AlgoimUtils delaunaytrian
 @publish AlgoimUtils convexhull