From a81a688e415175650b7fac35a06a3de2015f74cf Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Tue, 24 Sep 2024 17:53:22 +0200 Subject: [PATCH] Compute distance FE function for Algoim LSs with DistributedCellFields --- src/AlgoimUtils/AlgoimUtils.jl | 123 +++++++++++++++------ test/AlgoimUtilsTests/ClosestPointTests.jl | 35 ++++-- 2 files changed, 115 insertions(+), 43 deletions(-) diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl index 8b994c7..21941a4 100644 --- a/src/AlgoimUtils/AlgoimUtils.jl +++ b/src/AlgoimUtils/AlgoimUtils.jl @@ -87,6 +87,8 @@ function AlgoimCallLevelSetFunction(φ::CellField,∇φ::CellField) end struct DistributedAlgoimCallLevelSetFunction{A<:AbstractArray{<:AlgoimCallLevelSetFunction}} <: GridapType + values ::DistributedCellField + gradients::DistributedCellField levelsets::A end @@ -94,7 +96,7 @@ local_views(a::DistributedAlgoimCallLevelSetFunction) = a.levelsets function AlgoimCallLevelSetFunction(φ::DistributedCellField,∇φ::DistributedCellField) levelsets = map((v,g)->AlgoimCallLevelSetFunction(v,g),local_views(φ),local_views(∇φ)) - DistributedAlgoimCallLevelSetFunction(levelsets) + DistributedAlgoimCallLevelSetFunction(φ,∇φ,levelsets) end function normal(phi::AlgoimCallLevelSetFunction,x::AbstractVector{<:Point},cell_id::Int=1) @@ -112,6 +114,11 @@ function normal(phi::AlgoimCallLevelSetFunction,trian::DistributedTriangulation) DistributedCellField(normals,trian) end +function normal(phi::DistributedAlgoimCallLevelSetFunction,trian::DistributedTriangulation) + normals = map((φ,t)->normal(φ,t),local_views(phi),local_views(trian)) + DistributedCellField(normals,trian) +end + function normal(ls::AlgoimCallLevelSetFunction{<:CellField,<:CellField},x::Point,cell_id::Int=1) (carr,cgetindex,ceval) = ls.cache_∇φ gx = evaluate!(ceval,getindex!(cgetindex,carr,cell_id),x) @@ -339,8 +346,7 @@ function compute_closest_point_projections(model::DistributedCartesianDiscreteMo cmin = Int32[cmin...] cpartition = Int32[cdesc.partition...] cmax = cmin + cpartition - fill_cpp_data(φ,gpartition, - xmin,xmax, + fill_cpp_data(φ,gpartition,xmin,xmax, cppdegree,trim,limitstol, rmin=cmin,rmax=cmax) end @@ -354,32 +360,29 @@ function compute_closest_point_projections(model::DistributedCartesianDiscreteMo cppdegree::Int=2, trim::Bool=false, limitstol::Float64=1.0e-8) - - @notimplemented - # gdesc = model.metadata.descriptor - # xmin = gdesc.origin - # gpartition = Int32[gdesc.partition...] - # xmax = xmin + Point(gdesc.sizes .* gpartition) - # cpps = map(local_views(model),local_views(φ)) do m,f - # cdesc = get_cartesian_descriptor(m) - # cmin = ( cdesc.origin - gdesc.origin ) - # cmin = round.(cmin.data ./ gdesc.sizes) .+ 1 - # cmin = Int32[cmin...] - # cpartition = Int32[cdesc.partition...] - # cmax = cmin + cpartition - # fill_cpp_data(f,gpartition,xmin,xmax,cppdegree,trim,limitstol) - # end - # node_gids = get_face_gids(model,0) - # PVector(cpps,partition(node_gids)) |> consistent! |> wait - # cpps + gdesc = model.metadata.descriptor + xmin = gdesc.origin + gpartition = Int32[gdesc.partition...] + xmax = xmin + Point(gdesc.sizes .* gpartition) + cpps = map(local_views(model),local_views(φ)) do m,f + cdesc = get_cartesian_descriptor(m) + cmin = ( cdesc.origin - gdesc.origin ) + cmin = round.(cmin.data ./ gdesc.sizes) + cmin = Int32[cmin...] + cpartition = Int32[cdesc.partition...] + cmax = cmin + cpartition + fill_cpp_data(f,gpartition,xmin,xmax, + cppdegree,trim,limitstol, + rmin=cmin,rmax=cmax) + end + node_gids = get_face_gids(model,0) + PVector(cpps,partition(node_gids)) |> consistent! |> wait + cpps end function compute_closest_point_projections(fespace::FESpace, - φ::AlgoimCallLevelSetFunction, - order::Int; - cppdegree::Int=2, - trim::Bool=false, - limitstol::Float64=1.0e-8) + φ::Union{AlgoimCallLevelSetFunction,DistributedAlgoimCallLevelSetFunction}, + order::Int;cppdegree::Int=2,trim::Bool=false,limitstol::Float64=1.0e-8) trian = get_triangulation(fespace) model = get_background_model(trian) cps = compute_closest_point_projections( @@ -397,11 +400,11 @@ function compute_closest_point_projections(model::CartesianDiscreteModel, xmin = cdesc.origin xmax = xmin + Point(cdesc.sizes .* cdesc.partition) partition = Int32[cdesc.partition...] .* Int32(order) - fill_cpp_data(φ,partition,xmin,xmax,cppdegree,trim,limitstol) + fill_cpp_data(φ,partition,xmin,xmax,cppdegree,trim,limitstol,order=order) end function compute_closest_point_projections(model::DistributedCartesianDiscreteModel, - φ::AlgoimCallLevelSetFunction, + φ::AlgoimCallLevelSetFunction{<:Function,<:Function}, order::Int; cppdegree::Int=2, trim::Bool=false, @@ -416,9 +419,41 @@ function compute_closest_point_projections(model::DistributedCartesianDiscreteMo cmin = Int.(round.(cmin.data ./ gdesc.sizes) .* order) cpartition = cdesc.partition .* order cmax = cmin .+ cpartition - fill_cpp_data(φ,gpartition, - xmin,xmax, - cppdegree,trim,limitstol, + fill_cpp_data(φ,gpartition,xmin,xmax,cppdegree,trim,limitstol, + order=order,rmin=Int32[cmin...],rmax=Int32[cmax...]) + end +end + +function compute_closest_point_projections(model::DistributedCartesianDiscreteModel, + φ::AlgoimCallLevelSetFunction{<:CellField,<:CellField}, + order::Int; + cppdegree::Int=2, + trim::Bool=false, + limitstol::Float64=1.0e-8) + @notimplemented "The Algoim LS function must be defined from a DistributedCellField" +end + +function compute_closest_point_projections(model::DistributedCartesianDiscreteModel, + φ::DistributedAlgoimCallLevelSetFunction, + order::Int; + cppdegree::Int=2, + trim::Bool=false, + limitstol::Float64=1.0e-8) + gdesc = model.metadata.descriptor + xmin = gdesc.origin + xmax = xmin + Point(gdesc.sizes .* gdesc.partition) + gpartition = gdesc.partition .* order + gsizes = gdesc.sizes ./ order + ggrid = vec(map(i->xmin+Point((Tuple(i).-1).*gsizes),CartesianIndices(gpartition.+1))) + gvals = evaluate(φ.values,ggrid) + gpartition = Int32[gpartition...] + map(local_views(model)) do m + cdesc = get_cartesian_descriptor(m) + cmin = ( cdesc.origin - gdesc.origin ) + cmin = Int.(round.(cmin.data ./ gdesc.sizes) .* order) + cpartition = cdesc.partition .* order + cmax = cmin .+ cpartition + fill_cpp_data(gvals,gpartition,xmin,xmax,cppdegree,trim,limitstol, rmin=Int32[cmin...],rmax=Int32[cmax...]) end end @@ -518,7 +553,7 @@ end function compute_normal_displacement( cps::AbstractArray, - phi::AlgoimCallLevelSetFunction, + phi::Union{AlgoimCallLevelSetFunction,DistributedAlgoimCallLevelSetFunction}, fun::DistributedCellField, dt::Float64, Ω::DistributedTriangulation) @@ -526,7 +561,7 @@ function compute_normal_displacement( map(cps) do cp fun_xs = evaluate(fun,cp) nΓ_xs = evaluate(normal(phi,Ω),cp) - lazy_map(fun_xs,nΓ_xs) do fun_x,nΓ_x + map(fun_xs,nΓ_xs) do fun_x,nΓ_x dt * ( fun_x ⋅ nΓ_x ) end end @@ -680,6 +715,28 @@ function compute_distance_fe_function( FEFunction(fespace,dists) end +function compute_distance_fe_function( + bgmodel::DistributedCartesianDiscreteModel, + fespace::FESpace, + φ::DistributedAlgoimCallLevelSetFunction, + order::Int; + cppdegree::Int=2) + cps = compute_closest_point_projections( + fespace,φ,order,cppdegree=cppdegree) + _dists = map(cps, + local_views(fespace), + local_views(bgmodel), + local_views(φ)) do cp,fs,bg,φl + rm = refine(bg,order) + cos = get_node_coordinates(rm) + cos = node_to_dof_order(cos,fs,bg,order) + _compute_signed_distance(φl,cp,cos) + end + dists = PVector(_dists,partition(fespace.gids)) + consistent!(dists) |> wait + FEFunction(fespace,dists) +end + function compute_distance_fe_function( fespace_scalar_type::FESpace, fespace_vector_type::FESpace, diff --git a/test/AlgoimUtilsTests/ClosestPointTests.jl b/test/AlgoimUtilsTests/ClosestPointTests.jl index 8e1027e..ef491e7 100644 --- a/test/AlgoimUtilsTests/ClosestPointTests.jl +++ b/test/AlgoimUtilsTests/ClosestPointTests.jl @@ -16,7 +16,7 @@ const CUT = 0 function run_case(distribute,parts,degree) ranks = distribute(LinearIndices((prod(parts),))) - cells = Int32[4,4] + cells = Int32[16,16] domain = (-1.1,1.1,-1.1,1.1) bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells) @@ -32,19 +32,34 @@ function run_case(distribute,parts,degree) phi = AlgoimCallLevelSetFunction(f,∇f) fₕ = compute_distance_fe_function(bgmodel,V,phi,order,cppdegree=degree) - writevtk(Ω,"distance",cellfields=["f"=>fₕ]) - + # id = parts[1] + # writevtk(Ω,"distance",cellfields=["f"=>fₕ]) φₕ = AlgoimCallLevelSetFunction(fₕ,∇(fₕ)) - # cp₋ = compute_closest_point_projections(V,φₕ,order) - cpps = compute_closest_point_projections(V,phi,order) - nΓ = normal(phi,Ω) + gₕ = compute_distance_fe_function(bgmodel,V,φₕ,order,cppdegree=degree) + # id = parts[1] + # writevtk(Ω,"distance_$id",cellfields=["g"=>gₕ]) + + # cpps = compute_closest_point_projections(bgmodel,φₕ) + cpps = compute_closest_point_projections(V,φₕ,order) + # cpps = compute_closest_point_projections(V,phi,order) + # map(cpps) do cpp + # @show cpp + # end + nΓ = normal(φₕ,Ω) dt = 0.1 - dₕ = compute_normal_displacement(cpps,phi,nΓ,dt,Ω) - map(dₕ) do d - @show d - end + _dₕ = compute_normal_displacement(cpps,φₕ,nΓ,dt,Ω) + dₕ_fv = PVector(_dₕ,partition(V.gids)) + consistent!(dₕ_fv) |> wait + dₕ = FEFunction(V,dₕ_fv) + id = parts[1] + writevtk(Ω,"disps_$id",cellfields=["d"=>dₕ]) + + # dₕ = compute_normal_displacement(cpps,phi,nΓ,dt,Ω) + # map(dₕ) do d + # @show d + # end end