From e29f5920b66fb1263df1715b63e09c84e2e9b70a Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 18 Dec 2023 17:57:55 +0100 Subject: [PATCH 1/8] Adds Algoim utils + tests --- Project.toml | 3 + src/AlgoimUtils/AlgoimUtils.jl | 479 ++++++++++++++++++ src/Exports.jl | 14 + src/GridapEmbedded.jl | 2 + test/AlgoimUtilsTests/AlgoimInterfaceTests.jl | 108 ++++ test/AlgoimUtilsTests/ClosestPointTests.jl | 44 ++ test/AlgoimUtilsTests/DualQuadraturesTests.jl | 99 ++++ test/AlgoimUtilsTests/PoissonAlgoimTests.jl | 64 +++ .../AlgoimUtilsTests/QuadratureDegreeTests.jl | 167 ++++++ test/AlgoimUtilsTests/VisualizationTests.jl | 43 ++ .../VolumeConservationTests.jl | 69 +++ test/AlgoimUtilsTests/runtests.jl | 19 + test/runtests.jl | 2 + 13 files changed, 1113 insertions(+) create mode 100644 src/AlgoimUtils/AlgoimUtils.jl create mode 100644 test/AlgoimUtilsTests/AlgoimInterfaceTests.jl create mode 100644 test/AlgoimUtilsTests/ClosestPointTests.jl create mode 100644 test/AlgoimUtilsTests/DualQuadraturesTests.jl create mode 100644 test/AlgoimUtilsTests/PoissonAlgoimTests.jl create mode 100644 test/AlgoimUtilsTests/QuadratureDegreeTests.jl create mode 100644 test/AlgoimUtilsTests/VisualizationTests.jl create mode 100644 test/AlgoimUtilsTests/VolumeConservationTests.jl create mode 100644 test/AlgoimUtilsTests/runtests.jl diff --git a/Project.toml b/Project.toml index 7bb24c5..b9ce2fc 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,9 @@ version = "0.8.2" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +Algoim = "0eb9048c-21de-4c7a-bfac-056de1940b74" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" @@ -13,6 +15,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +algoimWrapper_jll = "3c43aa7b-5398-51f3-8d75-8f051e6faa4d" [compat] AbstractTrees = "0.3.3, 0.4" diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl new file mode 100644 index 0000000..2ae817a --- /dev/null +++ b/src/AlgoimUtils/AlgoimUtils.jl @@ -0,0 +1,479 @@ +module AlgoimUtils + +using CxxWrap +using Algoim +using Gridap + +import Base.prod +@inline prod(x::Point{N,T}) where {N,T} = prod(x.data) + +import Algoim: to_array +@inline to_array(x::Point{N,T}) where {N,T} = collect(x.data) + +using Gridap.Helpers +using Gridap.ReferenceFEs +using Gridap.Arrays: collect1d, CompressedArray, Table +using Gridap.Adaptivity +using Gridap.Fields: testitem +using Gridap.CellData +using Gridap.CellData: GenericCellField, get_data +using Gridap.CellData: _point_to_cell_cache, _point_to_cell! +using Gridap.Geometry + +using GridapEmbedded.Interfaces +using GridapEmbedded.Interfaces: Simplex +using GridapEmbedded.AgFEM: _aggregate_by_threshold_barrier + +using MiniQhull +using FillArrays + +import Algoim: CachedLevelSetValue +import Algoim: CachedLevelSetGradient +import Algoim: AlgoimCallLevelSetFunction +import Algoim: normal +import Gridap.ReferenceFEs: Quadrature +import GridapEmbedded.AgFEM: aggregate + +export TriangulationAndMeasure +export algoim +export Quadrature +export is_cell_active +export restrict_measure +export aggregate_narrow_band +export init_bboxes +export fill_cpp_data +export fill_cpp_data_raw +export compute_closest_point_projections +export compute_normal_displacement +export compute_distance_fe_function +export delaunaytrian +export convexhull + +struct Algoim <: QuadratureName end +const algoim = Algoim() + +(φ::CachedLevelSetValue{<:CellField})(p,i::Float32) = begin + _p = Point(p) + (carr,cgetindex,ceval) = φ.c + evaluate!(ceval,getindex!(cgetindex,carr,Int(i)),_p) +end + +(φ::CachedLevelSetGradient{<:CellField})(p,i::Float32) = begin + _p = Point(p) + (carr,cgetindex,ceval) = φ.c + _val = evaluate!(ceval,getindex!(cgetindex,carr,Int(i)),_p) + ConstCxxRef(to_uvector(to_array(_val))) +end + +function AlgoimCallLevelSetFunction(φ::CellField,∇φ::CellField) + φtrian = get_triangulation(φ) + ∇φtrian = get_triangulation(∇φ) + xφ = testitem(testitem(get_cell_coordinates(φtrian))) + x∇φ = testitem(testitem(get_cell_coordinates(∇φtrian))) + cellφ = get_array(φ) + cell∇φ = get_array(∇φ) + cache_φ = (cellφ,array_cache(cellφ),return_cache(testitem(cellφ),xφ)) + cache_∇φ = (cell∇φ,array_cache(cell∇φ),return_cache(testitem(cell∇φ),x∇φ)) + AlgoimCallLevelSetFunction{typeof(φ),typeof(∇φ),typeof(cache_φ),typeof(cache_∇φ)}(φ,∇φ,cache_φ,cache_∇φ) +end + +function normal(phi::AlgoimCallLevelSetFunction,x::AbstractVector{<:Point},cell_id::Int=1) + map(xi->normal(phi,xi,cell_id),x) +end + +function normal(phi::AlgoimCallLevelSetFunction,trian::Triangulation) + f = [ x -> normal(phi,x,i) for i in 1:num_cells(trian) ] + GenericCellField(f,trian,PhysicalDomain()) +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) + gx/norm(gx) +end + +function Quadrature(trian::Grid,::Algoim,phi::LevelSetFunction,degree::Int;kwargs...) + ctype_polytope = map(get_polytope,get_reffes(trian)) + @notimplementedif !all(map(is_n_cube,ctype_polytope)) + cell_to_coords = get_cell_coordinates(trian) + cell_to_bboxes = collect1d(lazy_map(a->(a[1],a[end]),cell_to_coords)) + jls = JuliaFunctionLevelSet(phi,Val{num_dims(trian)}()) + cell_to_quad = map(enumerate(cell_to_bboxes)) do (cell_id,bbox) + bbmin, bbmax = bbox + Quadrature(cell_id,bbmin,bbmax,jls,phi,degree;kwargs...) + end + CompressedArray(cell_to_quad,1:length(cell_to_quad)) +end + +function Quadrature(trian::Grid,::Algoim, + phi1::LevelSetFunction,phi2::LevelSetFunction, + degree::Int;kwargs...) + ctype_polytope = map(get_polytope,get_reffes(trian)) + @notimplementedif !all(map(is_n_cube,ctype_polytope)) + cell_to_coords = get_cell_coordinates(trian) + cell_to_bboxes = collect1d(lazy_map(a->(a[1],a[end]),cell_to_coords)) + jls1 = JuliaFunctionLevelSet(phi1,Val{num_dims(trian)}()) + jls2 = JuliaFunctionLevelSet(phi2,Val{num_dims(trian)}()) + cell_to_quad = map(enumerate(cell_to_bboxes)) do (cell_id,bbox) + bbmin, bbmax = bbox + Quadrature(cell_id,bbmin,bbmax,jls1,jls2,phi1,phi2,degree;kwargs...) + end + CompressedArray(cell_to_quad,1:length(cell_to_quad)) +end + +function Quadrature(cell_id::Int, + xmin::Point{N,T}, + xmax::Point{N,T}, + jls::LevelSetFunction, + phi::LevelSetFunction, + degree::Int; + phase::Int=CUT) where {N,T} + coords, weights = fill_quad_data(jls,phi,xmin,xmax,phase,degree,cell_id) + GenericQuadrature(coords,weights,"Algoim quadrature of degree $degree") +end + +function Quadrature(cell_id::Int, + xmin::Point{N,T}, + xmax::Point{N,T}, + jls1::LevelSetFunction, + jls2::LevelSetFunction, + phi1::LevelSetFunction, + phi2::LevelSetFunction, + degree::Int; + phase1::Int=IN, + phase2::Int=IN) where {N,T} + coords, weights = fill_quad_data(jls1,jls2,phi1,phi2,xmin,xmax, + phase1,phase2,degree,cell_id) + GenericQuadrature(coords,weights,"Algoim quadrature of degree $degree") +end + +function is_cell_active(meas::Measure) + has_non_empty_quad(x) = num_points(x) > 0 + lazy_map(has_non_empty_quad,get_data(meas.quad)) +end + +function restrict_measure(meas::Measure,trian::Triangulation) + ocell_quad = meas.quad + cell_quad = lazy_map(Reindex(ocell_quad.cell_quad),trian.tface_to_mface) + cell_point = lazy_map(Reindex(ocell_quad.cell_point),trian.tface_to_mface) + cell_weight = lazy_map(Reindex(ocell_quad.cell_weight),trian.tface_to_mface) + dds = ocell_quad.data_domain_style + ids = ocell_quad.integration_domain_style + Measure(CellQuadrature(cell_quad,cell_point,cell_weight,trian,dds,ids)) +end + +function TriangulationAndMeasure(Ωbg::Triangulation,quad::Tuple) + msg = "TriangulationAndMeasure can only receive the background triangulation" + @notimplementedif num_cells(get_background_model(Ωbg)) != num_cells(Ωbg) msg + dΩbg = Measure(Ωbg,quad,data_domain_style=PhysicalDomain()) + # RMK: This is a hack, but algoim interface does not let you + # know if a (given) cell intersects the interior of the level + # set. The hack consists in inferring this from the size of + # each quadrature. I do not expect this hack implies a lot of + # extra operations with regards to the proper way to do it. + cell_to_is_active = is_cell_active(dΩbg) + Ωᵃ = Triangulation(Ωbg,cell_to_is_active) + dΩᵃ = restrict_measure(dΩbg,Ωᵃ) + Ωᵃ,dΩᵃ,cell_to_is_active +end + +function aggregate(bgtrian,cell_to_is_active,cell_to_is_cut,in_or_out) + n_cells = length(cell_to_is_active) + @assert n_cells == length(cell_to_is_cut) + + cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + !isa ? 0.0 : (isc ? 0.0 : 1.0) + end + + cell_to_inoutcut = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + !isa ? OUT : (isc ? CUT : IN) + end + + cell_to_coords = get_cell_coordinates(bgtrian) + model = get_background_model(bgtrian) + topo = get_grid_topology(model) + D = num_cell_dims(model) + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + # A hack follows to avoid constructing the actual facet_to_inoutcut array + facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) + + threshold = 1.0 + _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) +end + +function aggregate_narrow_band(bgtrian,cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut,in_or_out) + n_cells = length(cell_to_is_in_narrow) + @assert n_cells == length(cell_to_is_active) + @assert n_cells == length(cell_to_is_cut) + + cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + ( isa & !isc ) ? 1.0 : 0.0 + end + + cell_to_inoutcut = lazy_map(cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut) do isn, isa, isc + !isn ? OUT : ( ( isa & !isc ) ? IN : CUT ) + end + + cell_to_coords = get_cell_coordinates(bgtrian) + model = get_background_model(bgtrian) + topo = get_grid_topology(model) + D = num_cell_dims(model) + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + # A hack follows to avoid constructing the actual facet_to_inoutcut array + facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) + + threshold = 1.0 + _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) +end + +using Gridap.Geometry: get_cell_to_parent_cell +using Gridap.CellData: get_cell_quadrature +using GridapEmbedded.AgFEM: compute_subcell_bbox +import GridapEmbedded.AgFEM: init_bboxes + +function init_bboxes(cell_to_coords,cut_measure::Measure) + bgcell_to_cbboxes = init_bboxes(cell_to_coords) + quad = get_cell_quadrature(cut_measure) + trian = get_triangulation(quad) + model = get_active_model(trian) + ccell_to_bgcell = get_cell_to_parent_cell(model) + for (cc,bc) in enumerate(ccell_to_bgcell) + bgcell_to_cbboxes[bc] = compute_subcell_bbox(quad.cell_point[cc]) + end + bgcell_to_cbboxes +end + +function compute_closest_point_projections(Ω::Triangulation,φ; + cppdegree::Int=2,trim::Bool=false,limitstol::Float64=1.0e-8) + compute_closest_point_projections(get_active_model(Ω),φ, + cppdegree=cppdegree,trim=trim,limitstol=limitstol) +end + +function compute_closest_point_projections(model::CartesianDiscreteModel, + φ::AlgoimCallLevelSetFunction; + cppdegree::Int=2, + trim::Bool=false, + limitstol::Float64=1.0e-8) + cdesc = get_cartesian_descriptor(model) + partition = Int32[cdesc.partition...] + xmin = cdesc.origin + xmax = xmin + Point(cdesc.sizes .* partition) + fill_cpp_data(φ,partition,xmin,xmax,cppdegree,trim,limitstol) +end + +function compute_closest_point_projections(fespace::FESpace, + φ::AlgoimCallLevelSetFunction, + order::Int; + cppdegree::Int=2, + trim::Bool=false, + limitstol::Float64=1.0e-8) + trian = get_triangulation(fespace) + model = get_active_model(trian) + # TO-ANSWER: Do I need the rmodel in this scope or not? + rmodel = refine(model,order) + cps = compute_closest_point_projections( + rmodel,φ,order,cppdegree=cppdegree,trim=trim,limitstol=limitstol) + msg = "Is the FE space order the same as the input order?" + @assert length(cps) == num_free_dofs(fespace) msg + cps = node_to_dof_order(cps,fespace,rmodel,order) +end + +function compute_closest_point_projections(model::AdaptedDiscreteModel, + φ::AlgoimCallLevelSetFunction, + order::Int; + cppdegree::Int=2, + trim::Bool=false, + limitstol::Float64=1.0e-8) + reffe = ReferenceFE(lagrangian,Float64,order) + rfespace = TestFESpace(model,reffe) + _rφ = interpolate_everywhere(φ.φ,rfespace) + rφ = AlgoimCallLevelSetFunction(_rφ,∇(_rφ)) + cdesc = get_cartesian_descriptor(get_model(model)) + partition = Int32[cdesc.partition...] + xmin = cdesc.origin + xmax = xmin + Point(cdesc.sizes .* partition) + fill_cpp_data(rφ,partition,xmin,xmax,cppdegree,trim,limitstol) +end + +function node_to_dof_order(cps, + fespace::FESpace, + rmodel::AdaptedDiscreteModel, + order::Int) + D = num_dims(rmodel) + cdesc = get_cartesian_descriptor(get_model(rmodel)) + partition = cdesc.partition + orders = tfill(order,Val{D}()) + ones = tfill(1,Val{D}()) + range = CartesianIndices(orders.+1) .- CartesianIndex(ones) # 0-based + ldof_to_lnode = get_ldof_to_lnode(orders,D) + o2n_faces_map = rmodel.glue.o2n_faces_map + node_partition = partition .+ 1 + ncells = num_cells(rmodel.parent) + cell_node_ids = lazy_map(1:ncells) do cellid + anchor_node = o2n_faces_map[cellid][1] + node_ijk = CartesianIndices(partition)[anchor_node] + range_cis = node_ijk .+ range + node_ids = LinearIndices(node_partition)[range_cis] + node_ids[ldof_to_lnode] + end + cell_dofs_ids = fespace.cell_dofs_ids + c1 = array_cache(fespace.cell_dofs_ids) + c2 = array_cache(cell_node_ids) + ncps = similar(cps) + for cellid in 1:ncells + dof_ids = getindex!(c1,cell_dofs_ids,cellid) + node_ids = getindex!(c2,cell_node_ids,cellid) + ncps[dof_ids] = cps[node_ids] + end + ncps +end + +function get_ldof_to_lnode(orders,D) + + # Generate indices of n-faces and order s.t. + # (1) dimension-increasing (2) lexicographic + bin_rang_nfaces = tfill(0:1,Val{D}()) + bin_ids_nfaces = vec(collect(Iterators.product(bin_rang_nfaces...))) + sum_bin_ids_nfaces = sum.(bin_ids_nfaces) + bin_ids_nfaces = permute!(bin_ids_nfaces,sortperm(sum_bin_ids_nfaces)) + + # Generate LIs of basis funs s.t. order by n-faces + lids_b = LinearIndices(Tuple([orders[i]+1 for i=1:D])) + + eet = eltype(eltype(bin_ids_nfaces)) + f(x) = Tuple( x[i] == one(eet) ? (0:0) : (1:orders[i]:(orders[i]+1)) for i in 1:length(x) ) + g(x) = Tuple( x[i] == one(eet) ? (2:orders[i]) : (0:0) for i in 1:length(x) ) + rang_nfaces = map(f,bin_ids_nfaces) + rang_own_dofs = map(g,bin_ids_nfaces) + + perm = Int64[] + for i = 1:length(bin_ids_nfaces) + cis_nfaces = CartesianIndices(rang_nfaces[i]) + cis_own_dofs = CartesianIndices(rang_own_dofs[i]) + for ci in cis_nfaces + ci = ci .+ cis_own_dofs + perm = vcat(perm,reshape(lids_b[ci],length(ci))) + end + end + + perm +end + +function compute_normal_displacement( + cps::AbstractVector{<:Point}, + phi::AlgoimCallLevelSetFunction, + fun, + dt::Float64, + Ω::Triangulation) + # Note that cps must be (scalar) DoF-numbered, not lexicographic-numbered + 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(Ω)) + 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 +end + +signed_distance(φ::Function,x,y) = sign(φ(y))*norm(x-y) + +signed_distance(φ::T,x,y) where {T<:Number} = sign(φ)*norm(x-y) + +function _compute_signed_distance( + φ::AlgoimCallLevelSetFunction{<:Function,<:Function}, + cps::Vector{<:Point},cos::Matrix{<:Point}) + _dist(x,y) = signed_distance(φ.φ,x,y) + lazy_map(_dist,cps,cos) +end + +function _compute_signed_distance( + φ::AlgoimCallLevelSetFunction{<:CellField,<:CellField}, + cps::Vector{<:Point},cos::Matrix{<:Point}) + φs = get_free_dof_values(φ.φ) + lazy_map(signed_distance,φs,cps,cos) +end + +function compute_distance_fe_function( + bgmodel::CartesianDiscreteModel, + fespace::FESpace, + φ::AlgoimCallLevelSetFunction, + order::Int; + cppdegree::Int=2) + cps = compute_closest_point_projections( + fespace,φ,order,cppdegree=cppdegree) + rmodel = refine(bgmodel,order) + cos = get_node_coordinates(rmodel) + cos = node_to_dof_order(cos,fespace,rmodel,order) + dists = _compute_signed_distance(φ,cps,cos) + FEFunction(fespace,dists) +end + +abstract type QhullType end + +struct DelaunayTrian <: QhullType end +const delaunaytrian = DelaunayTrian() + +struct ConvexHull <: QhullType end +const convexhull = ConvexHull() + +get_flags(::QhullType) = @abstractmethod +get_flags(::DelaunayTrian) = "qhull d Qt Qbb Qc Qz" +get_flags(::ConvexHull) = "qhull Qt Qc" + +get_dimension(::QhullType,dim) = @abstractmethod +get_dimension(::DelaunayTrian,dim) = dim +get_dimension(::ConvexHull,dim) = dim-1 + +import Gridap.Visualization: visualization_data + +function visualization_data(meas::Measure,filename;cellfields=Dict(),qhulltype=DelaunayTrian()) + node_coordinates = collect(Iterators.flatten(meas.quad.cell_point.values)) + grid = _to_grid(node_coordinates,qhulltype) + ndata = Dict() + for (k,v) in cellfields + ndata[k] = lazy_map(v,node_coordinates) + end + visualization_data(grid,filename,nodaldata=ndata) +end + +function visualization_data(meas::Vector{<:Measure},filename;cellfields=Dict(),qhulltype=DelaunayTrian()) + node_coordinates = vcat(map(m->collect(Iterators.flatten(m.quad.cell_point.values)),meas)...) + grid = _to_grid(node_coordinates,qhulltype) + ndata = Dict() + for (k,v) in cellfields + ndata[k] = lazy_map(v,node_coordinates) + end + visualization_data(grid,filename,nodaldata=ndata) +end + +function _to_grid(node_coordinates::Vector{<:Point{Dp,Tp}},qhulltype) where {Dp,Tp} + d = get_dimension(qhulltype,Dp) + connectivity = delaunay(reinterpret(node_coordinates),get_flags(qhulltype))[1:(d+1),:] + cell_node_ids = Table(collect(eachcol(connectivity))) + reffes = [LagrangianRefFE(Float64,Simplex(Val{d}()),1)] + cell_types = collect(Fill(Int8(1),length(cell_node_ids))) + UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types) +end + +end # module \ No newline at end of file diff --git a/src/Exports.jl b/src/Exports.jl index 8a6f2a5..cb8ea2b 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -49,4 +49,18 @@ end @publish AgFEM compute_cell_bboxes @publish AgFEM compute_cell_to_dface_bboxes +@publish AlgoimUtils algoim +@publish AlgoimUtils fill_quad_data +@publish AlgoimUtils AlgoimCallLevelSetFunction +@publish AlgoimUtils TriangulationAndMeasure +@publish AlgoimUtils normal +@publish AlgoimUtils init_bboxes +@publish AlgoimUtils fill_cpp_data +@publish AlgoimUtils fill_cpp_data_raw +@publish AlgoimUtils compute_closest_point_projections +@publish AlgoimUtils compute_normal_displacement +@publish AlgoimUtils compute_distance_fe_function +@publish AlgoimUtils delaunaytrian +@publish AlgoimUtils convexhull + @publish MomentFittedQuadratures momentfitted \ No newline at end of file diff --git a/src/GridapEmbedded.jl b/src/GridapEmbedded.jl index b84e57e..175cee0 100644 --- a/src/GridapEmbedded.jl +++ b/src/GridapEmbedded.jl @@ -10,6 +10,8 @@ include("AgFEM/AgFEM.jl") include("MomentFittedQuadratures/MomentFittedQuadratures.jl") +include("AlgoimUtils/AlgoimUtils.jl") + include("Exports.jl") end # module diff --git a/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl b/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl new file mode 100644 index 0000000..d1f64ec --- /dev/null +++ b/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl @@ -0,0 +1,108 @@ +module AlgoimInterfaceTests + +using Test +using Gridap +using Gridap.ReferenceFEs +using GridapEmbedded +using GridapEmbedded.Interfaces + +function run_case_raw_point(degree,phase) + + u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)) - 1.0 + gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5)) + + phi = AlgoimCallLevelSetFunction(u,gradu) + + xmin = Point(-1.1,-1.1) + xmax = Point(1.1,1.1) + + _, weights = fill_quad_data(phi,xmin,xmax,phase,degree) + + s = ∑(weights) + if phase == IN + @test s ≈ 1.6055602335693198 + elseif phase == CUT + @test s ≈ 5.04463087731472 + else + error() + end + +end + +function run_case_function(cells,degree,phase) + + domain = (-1.1,1.1,-1.1,1.1) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)) - 1.0 + gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5)) + + phi = AlgoimCallLevelSetFunction(u,gradu) + quad = Quadrature(algoim,phi,degree,phase=phase) + dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain()) + + s = ∫(1)dΩ + if phase == IN + @test ∑(s) ≈ 1.570796326817732 + elseif phase == CUT + @test ∑(s) ≈ 4.844224109684756 + else + error() + end + +end + +function run_case_fe_function(cells,order,degree,phase) + + domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(Ω,reffe) + U = TrialFESpace(V) + + u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.33*0.33)) - 1.0 + gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.33*0.33)) + + uₕ = interpolate(u,U) + ∇uₕ = ∇(uₕ) + + phi = AlgoimCallLevelSetFunction(uₕ,∇uₕ) + quad = Quadrature(algoim,phi,degree,phase=phase) + dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain()) + + s = ∫(1)dΩ + if phase == IN + @test ∑(s) ≈ 0.6911505144863921 + elseif phase == CUT + @test ∑(s) ≈ 4.382658365204242 + else + error() + end + +end + +f(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.33*0.33)) - 1.0 +gradf(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.33*0.33)) + +phi = AlgoimCallLevelSetFunction(f,gradf) +@test phi(VectorValue(1.0,0.5,0.33)) == 2.0 +@test all(phi.∇φ(VectorValue(1.0,0.5,0.33)) .≈ VectorValue(2.0, 4.0, 6.0606060606060606)) + +run_case_raw_point(3,IN) +run_case_raw_point(3,CUT) + +run_case_function((64,64),3,IN) +run_case_function((64,66),3,IN) + +run_case_function((64,64),3,CUT) +run_case_function((64,66),3,CUT) + +run_case_fe_function((32,32,32),2,3,IN) +run_case_fe_function((32,32,32),2,3,CUT) + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/ClosestPointTests.jl b/test/AlgoimUtilsTests/ClosestPointTests.jl new file mode 100644 index 0000000..0779288 --- /dev/null +++ b/test/AlgoimUtilsTests/ClosestPointTests.jl @@ -0,0 +1,44 @@ +module ClosestPointTests + +using Test +using CxxWrap +using Gridap +using GridapEmbedded + +const IN = -1 +const OUT = 1 +const CUT = 0 + +function run_case(degree) + + xmin = Point(-1.1,-1.1,-1.1) + xmax = Point(1.1,1.1,1.1) + partition = Int32[16,16,16] + + domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1) + bgmodel = CartesianDiscreteModel(domain,partition) + Ω = Triangulation(bgmodel) + + reffeᵠ = ReferenceFE(lagrangian,Float64,1) + V = TestFESpace(Ω,reffeᵠ,conformity=:L2) + + f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.5*0.5)) - 1.0 + function gradf(x::V) where {V} + V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.5*0.5)]) + end + + fₕ = interpolate_everywhere(f,V) + phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ)) + + coords = fill_cpp_data(phi,partition,xmin,xmax,degree) + @test maximum(f.(coords)) < 1.0e-4 + +end + +run_case(2) +run_case(3) +run_case(4) +run_case(5) +run_case(-1) + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/DualQuadraturesTests.jl b/test/AlgoimUtilsTests/DualQuadraturesTests.jl new file mode 100644 index 0000000..13f2fad --- /dev/null +++ b/test/AlgoimUtilsTests/DualQuadraturesTests.jl @@ -0,0 +1,99 @@ +module DualQuadraturesTests + +using Test +using Gridap +using Gridap.ReferenceFEs +using GridapEmbedded +using GridapEmbedded.Interfaces + +# 3D funs + +u¹(x) = ( (x[1]-0.25) * (x[1]-0.25) + + (x[2]-0.25) * (x[2]-0.25) + + (x[3]-0.25) * (x[3]-0.25) ) / ( 0.20 * 0.20 ) - 1.0 + +u²(x) = ( (x[1]-0.75) * (x[1]-0.75) + + (x[2]-0.75) * (x[2]-0.75) + + (x[3]-0.75) * (x[3]-0.75) ) / ( 0.22 * 0.22 ) - 1.0 + +u³(x) = ( (x[1]-0.50) * (x[1]-0.50) + + (x[2]-0.50) * (x[2]-0.50) + + (x[3]-0.50) * (x[3]-0.50) ) / ( 0.40 * 0.40 ) - 1.0 + +# 2D funs + +v¹(x) = ( (x[1]-0.25) * (x[1]-0.25) + + (x[2]-0.25) * (x[2]-0.25) ) / ( 0.20 * 0.20 ) - 1.0 + +v²(x) = ( (x[1]-0.75) * (x[1]-0.75) + + (x[2]-0.75) * (x[2]-0.75) ) / ( 0.22 * 0.22 ) - 1.0 + +v³(x) = ( (x[1]-0.50) * (x[1]-0.50) + + (x[2]-0.50) * (x[2]-0.50) ) / ( 0.40 * 0.40 ) - 1.0 + +function run_case_fe_function(domain,cells,order,degree,phase¹,phase²,u¹,u²) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(Ω,reffe) + U = TrialFESpace(V) + + uₕ¹ = interpolate(u¹,U) + ∇uₕ¹ = ∇(uₕ¹) + + uₕ² = interpolate(u²,U) + ∇uₕ² = ∇(uₕ²) + + phi¹ = AlgoimCallLevelSetFunction(uₕ¹,∇uₕ¹) + phi² = AlgoimCallLevelSetFunction(uₕ²,∇uₕ²) + + quad = Quadrature(algoim,phi¹,phi²,degree,phase1=phase¹,phase2=phase²) + dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain()) + + s = ∑(∫(1)dΩ) + if ( phase¹ == IN ) && ( phase² == OUT ) + @test s ≈ 0.03351032164371535 + elseif ( phase¹ == OUT ) && ( phase² == IN ) + @test s ≈ 0.04460223810211001 + elseif ( phase¹ == IN ) && ( phase² == CUT ) + @test_skip s ≈ 0.11301756191551385 # Different precision in GitHub's CI workflow + elseif ( phase¹ == OUT ) && ( phase² == CUT ) + @test_skip s ≈ 1.8976021283682756 # Different precision in GitHub's CI workflow + end + +end + +domain3D = (0.0,1.1,0.0,1.1,0.0,1.1) + +run_case_fe_function(domain3D,(24,24,24),2,5,IN,OUT,u¹,u²) +# run_case_fe_function(domain3D,(24,24,24),2,5,OUT,IN,u¹,u²) + +run_case_fe_function(domain3D,(24,24,24),2,5,IN,CUT,u¹,u³) +# run_case_fe_function(domain3D,(24,24,24),2,5,OUT,CUT,u¹,u³) + +# domain2D = (0.0,1.1,0.0,1.1) +# domain3D = (0.0,1.1,0.0,1.1,0.0,1.1) +# order = 2 +# phaseA = IN +# phaseB = OUT +# u = u¹ +# v = u³ + +# T = pi*16/25 # intersection 3D +# # T = pi*4/5 # intersection 2D +# # T = pi*(0.2^2+0.22^2) # volumes 2D +# # T = (4/3)*pi*(0.2^3+0.22^3) # volumes 3D +# # T = 4*pi*(0.2^2+0.22^2) # non-intersecting surfaces 3D + +# for deg in [3,4,5] +# @info "Case: degree = $deg" +# for n in [6,12,24,48] +# A = run_case_fe_function(domain3D,(n,n,n),2,deg,phaseA,CUT,u¹,u³) +# B = run_case_fe_function(domain3D,(n,n,n),2,deg,phaseB,CUT,u¹,u³) +# @show A+B-T +# end +# end + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/PoissonAlgoimTests.jl b/test/AlgoimUtilsTests/PoissonAlgoimTests.jl new file mode 100644 index 0000000..e359b91 --- /dev/null +++ b/test/AlgoimUtilsTests/PoissonAlgoimTests.jl @@ -0,0 +1,64 @@ +module PoissonAlgoimTests + +using Test +using Gridap +using Gridap.ReferenceFEs +using GridapEmbedded +using GridapEmbedded.Interfaces + +function run_poisson(domain,cells,order,solution_degree) + + φ(x) = x[1]-x[2] + ∇φ(x) = VectorValue(1.0,-1.0,0.0) + phi = AlgoimCallLevelSetFunction(φ,∇φ) + + u(x) = sum(x)^solution_degree + f(x) = -Δ(u)(x) + ud(x) = u(x) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + degree = order == 1 ? 3 : 2*order + vquad = Quadrature(algoim,phi,degree,phase=IN) + Ωᵃ,dΩᵃ = TriangulationAndMeasure(Ω,vquad) + squad = Quadrature(algoim,phi,degree) + _,dΓ = TriangulationAndMeasure(Ω,squad) + n_Γ = normal(phi,Ω) + + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(Ωᵃ,reffe,conformity=:H1,dirichlet_tags="boundary") + U = TrialFESpace(V,ud) + + # Nitsche method + γᵈ = 2.0*order^2 + h = (domain[2]-domain[1])/cells[1] + + a(u,v) = + ∫( ∇(v)⋅∇(u) )dΩᵃ + + ∫( (γᵈ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )dΓ + + l(v) = + ∫( v*f )dΩᵃ + + ∫( (γᵈ/h)*v*ud - (n_Γ⋅∇(v))*ud )dΓ + + # FE problem + op = AffineFEOperator(a,l,U,V) + uₕ = solve(op) + + eₕ = u - uₕ + + l2(u) = √(∑( ∫( u*u )dΩᵃ )) + h1(u) = √(∑( ∫( u*u + ∇(u)⋅∇(u) )dΩᵃ )) + + el2 = l2(eₕ) + eh1 = h1(eₕ) + + @test el2 < 1.0e-014 + @test eh1 < 1.0e-013 + +end + +run_poisson((-1.1,1.1,-1.1,1.1,-1.1,1.1),(8,8,8),1,1) + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/QuadratureDegreeTests.jl b/test/AlgoimUtilsTests/QuadratureDegreeTests.jl new file mode 100644 index 0000000..530a792 --- /dev/null +++ b/test/AlgoimUtilsTests/QuadratureDegreeTests.jl @@ -0,0 +1,167 @@ +module QuadratureDegreeTests + +using Test +using Gridap +using Gridap.ReferenceFEs +using GridapEmbedded +using GridapEmbedded.Interfaces + +φ(x) = x[1]-x[2] +∇φ(x) = VectorValue(1.0,-1.0,0.0) +phi = AlgoimCallLevelSetFunction(φ,∇φ) + +function run_diffusion_reaction_bulk(domain,cells,order,degree) + + u(x) = sum(x)^order + f(x) = u(x)-Δ(u)(x) + ud(x) = u(x) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + vquad = Quadrature(algoim,phi,degree,phase=IN) + Ωᵃ,dΩᵃ,is_a = TriangulationAndMeasure(Ω,vquad) + squad = Quadrature(algoim,phi,degree) + _,dΓ,is_c = TriangulationAndMeasure(Ω,squad) + n_Γ = normal(phi,Ω) + + cell_to_cellin = aggregate(Ω,is_a,is_c,IN) + + reffe = ReferenceFE(lagrangian,Float64,order) + Vstd = TestFESpace(Ωᵃ,reffe,conformity=:H1,dirichlet_tags="boundary") + + V = AgFEMSpace(Vstd,cell_to_cellin) + U = TrialFESpace(V,ud) + + # Nitsche method + γᵈ = 5.0*order^2 + h = (domain[2]-domain[1])/cells[1] + + a(u,v) = + ∫( v⋅u )dΩᵃ + + ∫( ∇(v)⋅∇(u) )dΩᵃ + + ∫( (γᵈ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )dΓ + + l(v) = + ∫( v*f )dΩᵃ + + ∫( (γᵈ/h)*v*ud - (n_Γ⋅∇(v))*ud )dΓ + + # FE problem + op = AffineFEOperator(a,l,U,V) + uₕ = solve(op) + + eₕ = u - uₕ + + l2(u) = √(∑( ∫( u*u )dΩᵃ )) + h1(u) = √(∑( ∫( u*u + ∇(u)⋅∇(u) )dΩᵃ )) + + el2 = l2(eₕ) + eh1 = h1(eₕ) + + el2, eh1 + +end + +left_project(u,n) = u - n⊗(n⋅u) +∇ᵈ(u,n) = left_project(∇(u),n) + +function run_diffusion_reaction_surface(domain,cells,order,degree) + + u(x) = sum(x)^order + f(x) = u(x)-Δ(u)(x) + ud(x) = u(x) + + model = CartesianDiscreteModel(domain,cells) + Ω = Triangulation(model) + + squad = Quadrature(algoim,phi,degree) + Ωc,dΓ = TriangulationAndMeasure(Ω,squad) + + n_Γ = normal(phi,Ωc) + dΩc = Measure(Ωc,2*order) + + reffe = ReferenceFE(lagrangian,Float64,order) + + V = TestFESpace(Ωc,reffe,dirichlet_tags="boundary") + U = TrialFESpace(V,ud) + + m(u,v) = ∫( v⋅u )dΓ + a(u,v) = ∫( ∇ᵈ(v,n_Γ)⋅∇ᵈ(u,n_Γ) )dΓ + + # Γ-orthogonal derivative volume stabilisation + h = (domain[2]-domain[1])/cells[1] + γˢ = 10.0*h + s(u,v) = ∫( γˢ*((n_Γ⋅∇(v))⊙(n_Γ⋅∇(u))) )dΩc + + aₛ(u,v) = m(u,v) + a(u,v) + s(u,v) + l(v) = ∫( v⊙f )dΓ + + op = AffineFEOperator(aₛ,l,U,V) + uₕ = solve(op) + eₕ = uₕ - u + + l2(v) = √(∑(∫(v⋅v)dΓ)) + l2(eₕ),l2(∇ᵈ(eₕ,n_Γ)) + +end + +domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1) +partition = (8,8,8) + +# BULK tests + +order = 1 +degree = 2 +el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree) +@test el2 < 1.0e-015 +@test eh1 < 1.0e-014 + +order = 2 +degree = 3 +el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree) +@test el2 < 1.0e-013 +@test eh1 < 1.0e-012 + +# For order > 2, expect accuracy loss with Lagrangian AgFEM + +# order = 3 +# degree = 5 +# el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree) +# @test el2 < 1.0e-011 +# @test eh1 < 1.0e-009 + +# order = 4 +# degree = 6 +# el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree) +# @test el2 < 1.0e-007 +# @test eh1 < 1.0e-005 + +# SURFACE tests + +order = 1 +degree = 2 +el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree) +@test el2 < 1.0e-015 +@test eh1 < 1.0e-014 + +order = 2 +degree = 3 +el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree) +@test el2 < 1.0e-013 +@test eh1 < 1.0e-012 + +# For order > 2, expect accuracy loss with Lagrangian AgFEM + +# order = 3 +# degree = 5 +# el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree) +# @test el2 < 1.0e-009 +# @test eh1 < 1.0e-008 + +# order = 4 +# degree = 6 +# el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree) +# @test el2 < 1.0e-007 +# @test eh1 < 1.0e-005 + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/VisualizationTests.jl b/test/AlgoimUtilsTests/VisualizationTests.jl new file mode 100644 index 0000000..81cd7ea --- /dev/null +++ b/test/AlgoimUtilsTests/VisualizationTests.jl @@ -0,0 +1,43 @@ +module VisualizationTests + + using CxxWrap + using Gridap + using GridapEmbedded + using Gridap.ReferenceFEs + + const IN = -1 + const OUT = 1 + const CUT = 0 + + order = 1 + domain = (-1.1,1.1,-1.1,1.1) + + f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)) - 1.0 + function gradf(x::V) where {V} + V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5)]) + end + + n = 12 + degree = 3 + cppdegree = -1 + + partition = Int32[n,n] + bgmodel = CartesianDiscreteModel(domain,partition) + Ω = Triangulation(bgmodel) + + reffeᵠ = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(Ω,reffeᵠ) + + fₕ = interpolate_everywhere(f,V) + phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ)) + + squad = Quadrature(algoim,phi,degree,phase=CUT) + _,dΓ = TriangulationAndMeasure(Ω,squad) + + vquad = Quadrature(algoim,phi,degree,phase=IN) + _,dΩ = TriangulationAndMeasure(Ω,vquad) + + writevtk(dΓ,"res_sur",cellfields=["f"=>fₕ],qhulltype=convexhull) + writevtk([dΩ,dΓ],"res_vol",cellfields=["f"=>fₕ]) + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/VolumeConservationTests.jl b/test/AlgoimUtilsTests/VolumeConservationTests.jl new file mode 100644 index 0000000..374d26d --- /dev/null +++ b/test/AlgoimUtilsTests/VolumeConservationTests.jl @@ -0,0 +1,69 @@ +module VolumeConservationTests + +using Test +using CxxWrap +using Gridap +using GridapEmbedded + +using Gridap.ReferenceFEs + +const IN = -1 +const OUT = 1 +const CUT = 0 + +domain = (-1.1,1.1,-1.1,1.1) + +f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)) - 1.0 +function gradf(x::V) where {V} + V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5)]) +end + +function run_case(n::Int,order::Int,cppdegree::Int) + + partition = Int32[n,n] + bgmodel = CartesianDiscreteModel(domain,partition) + Ω = Triangulation(bgmodel) + + reffeᵠ = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(Ω,reffeᵠ) + + fₕ = interpolate_everywhere(f,V) + phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ)) + + degree = order == 1 ? 3 : 2*order + vquad = Quadrature(algoim,phi,degree,phase=IN) + dΩ = Measure(Ω,vquad,data_domain_style=PhysicalDomain()) + vol_1 = ∑(∫(1)dΩ) + # @show vol_1 + + cps = compute_closest_point_projections(V,phi,order,cppdegree=cppdegree) + + reffeᶜ = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + W = TestFESpace(Ω,reffeᶜ) + + g(x) = VectorValue(1.0,0.0) + gₕ = interpolate_everywhere(g,W) + + dt = 0.2 + _phi₂ = compute_normal_displacement(cps,phi,gₕ,dt,Ω) + _phi₂ = get_free_dof_values(fₕ) - _phi₂ + _phi₂ = FEFunction(V,_phi₂) + phi₂ = AlgoimCallLevelSetFunction(_phi₂,∇(_phi₂)) + + vquad = Quadrature(algoim,phi₂,degree,phase=IN) + dΩ₂ = Measure(Ω,vquad,data_domain_style=PhysicalDomain()) + # writevtk(dΩ₂,"res_vol") + vol_2 = ∑(∫(1)dΩ₂) + # @show vol_2 + + abs(vol_1-vol_2)/abs(vol_1) + +end + +order = 2 # Affects CP accuracy _and_ volume conservation +cppdegree = 2 # Affects CP accuracy, but not volume conservation + +@test run_case(12,order,cppdegree) < 1.0e-5 +# @test run_case(24,order,cppdegree) < 1.0e-6 + +end # module \ No newline at end of file diff --git a/test/AlgoimUtilsTests/runtests.jl b/test/AlgoimUtilsTests/runtests.jl new file mode 100644 index 0000000..f775609 --- /dev/null +++ b/test/AlgoimUtilsTests/runtests.jl @@ -0,0 +1,19 @@ +module AlgoimTests + +using Test + +@testset "AlgoimInterface" begin include("AlgoimInterfaceTests.jl") end + +@testset "QuadratureDegree" begin include("QuadratureDegreeTests.jl") end + +@testset "DualQuadratures" begin include("DualQuadraturesTests.jl") end + +@testset "PoissonAlgoim" begin include("PoissonAlgoimTests.jl") end + +@testset "ClosestPoint" begin include("ClosestPointTests.jl") end + +@testset "VolumeConservation" begin include("VolumeConservationTests.jl") end + +@testset "Visualization" begin include("VisualizationTests.jl") end + +end # module \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 8b1363f..83cf4d2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,8 @@ if QHULL_LOADED @time @testset "MomentFittedQuadratures" begin include("MomentFittedQuadraturesTests/runtests.jl") end + @time @testset "AlgoimUtils" begin include("AlgoimUtilsTests/runtests.jl") end + include(joinpath(@__DIR__,"..","examples","runexamples.jl")) else From 26a8e8a0856a91f0d92c72c157f0b1379cc36c69 Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 18 Dec 2023 18:05:03 +0100 Subject: [PATCH 2/8] Relocated some procedures --- src/AgFEM/AgFEM.jl | 1 + src/AgFEM/AggregateBoundingBoxes.jl | 12 +++++ src/AgFEM/CellAggregation.jl | 57 ++++++++++++++++++++++ src/AlgoimUtils/AlgoimUtils.jl | 73 ----------------------------- src/Exports.jl | 2 +- 5 files changed, 71 insertions(+), 74 deletions(-) diff --git a/src/AgFEM/AgFEM.jl b/src/AgFEM/AgFEM.jl index 2bf2684..72fed37 100644 --- a/src/AgFEM/AgFEM.jl +++ b/src/AgFEM/AgFEM.jl @@ -18,6 +18,7 @@ using GridapEmbedded.Interfaces: compute_subcell_to_inout export aggregate export color_aggregates +export aggregate_narrow_band export AggregateAllCutCells export compute_cell_bboxes export compute_cell_to_dface_bboxes diff --git a/src/AgFEM/AggregateBoundingBoxes.jl b/src/AgFEM/AggregateBoundingBoxes.jl index 4604144..80f26b0 100644 --- a/src/AgFEM/AggregateBoundingBoxes.jl +++ b/src/AgFEM/AggregateBoundingBoxes.jl @@ -19,6 +19,18 @@ function init_bboxes(cell_to_coords,cut::EmbeddedDiscretization;in_or_out=IN) bgcell_to_cbboxes end +function init_bboxes(cell_to_coords,cut_measure::Measure) + bgcell_to_cbboxes = init_bboxes(cell_to_coords) + quad = get_cell_quadrature(cut_measure) + trian = get_triangulation(quad) + model = get_active_model(trian) + ccell_to_bgcell = get_cell_to_parent_cell(model) + for (cc,bc) in enumerate(ccell_to_bgcell) + bgcell_to_cbboxes[bc] = compute_subcell_bbox(quad.cell_point[cc]) + end + bgcell_to_cbboxes +end + function init_cut_bboxes(cut,ccell_to_bgcell,in_or_out) subcell_to_inout = compute_subcell_to_inout(cut,cut.geo) subcell_is_in = lazy_map(i->i==in_or_out,subcell_to_inout) diff --git a/src/AgFEM/CellAggregation.jl b/src/AgFEM/CellAggregation.jl index c41df61..90df2fb 100644 --- a/src/AgFEM/CellAggregation.jl +++ b/src/AgFEM/CellAggregation.jl @@ -303,3 +303,60 @@ function _color_aggregates_barrier(cell_to_cellin,cell_to_faces,face_to_cells) cell_to_color end + +# Specialised methods for Algoim quadratures + +function aggregate(bgtrian,cell_to_is_active,cell_to_is_cut,in_or_out) + n_cells = length(cell_to_is_active) + @assert n_cells == length(cell_to_is_cut) + + cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + !isa ? 0.0 : (isc ? 0.0 : 1.0) + end + + cell_to_inoutcut = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + !isa ? OUT : (isc ? CUT : IN) + end + + cell_to_coords = get_cell_coordinates(bgtrian) + model = get_background_model(bgtrian) + topo = get_grid_topology(model) + D = num_cell_dims(model) + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + # A hack follows to avoid constructing the actual facet_to_inoutcut array + facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) + + threshold = 1.0 + _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) +end + +function aggregate_narrow_band(bgtrian,cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut,in_or_out) + n_cells = length(cell_to_is_in_narrow) + @assert n_cells == length(cell_to_is_active) + @assert n_cells == length(cell_to_is_cut) + + cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc + ( isa & !isc ) ? 1.0 : 0.0 + end + + cell_to_inoutcut = lazy_map(cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut) do isn, isa, isc + !isn ? OUT : ( ( isa & !isc ) ? IN : CUT ) + end + + cell_to_coords = get_cell_coordinates(bgtrian) + model = get_background_model(bgtrian) + topo = get_grid_topology(model) + D = num_cell_dims(model) + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + # A hack follows to avoid constructing the actual facet_to_inoutcut array + facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) + + threshold = 1.0 + _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) +end \ No newline at end of file diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl index 2ae817a..9f057c3 100644 --- a/src/AlgoimUtils/AlgoimUtils.jl +++ b/src/AlgoimUtils/AlgoimUtils.jl @@ -22,7 +22,6 @@ using Gridap.Geometry using GridapEmbedded.Interfaces using GridapEmbedded.Interfaces: Simplex -using GridapEmbedded.AgFEM: _aggregate_by_threshold_barrier using MiniQhull using FillArrays @@ -32,15 +31,12 @@ import Algoim: CachedLevelSetGradient import Algoim: AlgoimCallLevelSetFunction import Algoim: normal import Gridap.ReferenceFEs: Quadrature -import GridapEmbedded.AgFEM: aggregate export TriangulationAndMeasure export algoim export Quadrature export is_cell_active export restrict_measure -export aggregate_narrow_band -export init_bboxes export fill_cpp_data export fill_cpp_data_raw export compute_closest_point_projections @@ -177,77 +173,8 @@ function TriangulationAndMeasure(Ωbg::Triangulation,quad::Tuple) Ωᵃ,dΩᵃ,cell_to_is_active end -function aggregate(bgtrian,cell_to_is_active,cell_to_is_cut,in_or_out) - n_cells = length(cell_to_is_active) - @assert n_cells == length(cell_to_is_cut) - - cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc - !isa ? 0.0 : (isc ? 0.0 : 1.0) - end - - cell_to_inoutcut = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc - !isa ? OUT : (isc ? CUT : IN) - end - - cell_to_coords = get_cell_coordinates(bgtrian) - model = get_background_model(bgtrian) - topo = get_grid_topology(model) - D = num_cell_dims(model) - cell_to_faces = get_faces(topo,D,D-1) - face_to_cells = get_faces(topo,D-1,D) - # A hack follows to avoid constructing the actual facet_to_inoutcut array - facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) - - threshold = 1.0 - _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) -end - -function aggregate_narrow_band(bgtrian,cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut,in_or_out) - n_cells = length(cell_to_is_in_narrow) - @assert n_cells == length(cell_to_is_active) - @assert n_cells == length(cell_to_is_cut) - - cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc - ( isa & !isc ) ? 1.0 : 0.0 - end - - cell_to_inoutcut = lazy_map(cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut) do isn, isa, isc - !isn ? OUT : ( ( isa & !isc ) ? IN : CUT ) - end - - cell_to_coords = get_cell_coordinates(bgtrian) - model = get_background_model(bgtrian) - topo = get_grid_topology(model) - D = num_cell_dims(model) - cell_to_faces = get_faces(topo,D,D-1) - face_to_cells = get_faces(topo,D-1,D) - # A hack follows to avoid constructing the actual facet_to_inoutcut array - facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) - - threshold = 1.0 - _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) -end - using Gridap.Geometry: get_cell_to_parent_cell using Gridap.CellData: get_cell_quadrature -using GridapEmbedded.AgFEM: compute_subcell_bbox -import GridapEmbedded.AgFEM: init_bboxes - -function init_bboxes(cell_to_coords,cut_measure::Measure) - bgcell_to_cbboxes = init_bboxes(cell_to_coords) - quad = get_cell_quadrature(cut_measure) - trian = get_triangulation(quad) - model = get_active_model(trian) - ccell_to_bgcell = get_cell_to_parent_cell(model) - for (cc,bc) in enumerate(ccell_to_bgcell) - bgcell_to_cbboxes[bc] = compute_subcell_bbox(quad.cell_point[cc]) - end - bgcell_to_cbboxes -end function compute_closest_point_projections(Ω::Triangulation,φ; cppdegree::Int=2,trim::Bool=false,limitstol::Float64=1.0e-8) diff --git a/src/Exports.jl b/src/Exports.jl index cb8ea2b..66974a1 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -44,6 +44,7 @@ end @publish AgFEM AgFEMSpace @publish AgFEM aggregate @publish AgFEM color_aggregates +@publish AgFEM aggregate_narrow_band @publish AgFEM AggregateCutCellsByThreshold @publish AgFEM AggregateAllCutCells @publish AgFEM compute_cell_bboxes @@ -54,7 +55,6 @@ end @publish AlgoimUtils AlgoimCallLevelSetFunction @publish AlgoimUtils TriangulationAndMeasure @publish AlgoimUtils normal -@publish AlgoimUtils init_bboxes @publish AlgoimUtils fill_cpp_data @publish AlgoimUtils fill_cpp_data_raw @publish AlgoimUtils compute_closest_point_projections From 74c93c6d7e60867c7a0177aedccef29bcdc32519 Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Tue, 19 Dec 2023 12:27:54 +0100 Subject: [PATCH 3/8] Perturbing VisualizationTest for testing engine --- test/AlgoimUtilsTests/VisualizationTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/AlgoimUtilsTests/VisualizationTests.jl b/test/AlgoimUtilsTests/VisualizationTests.jl index 81cd7ea..9118178 100644 --- a/test/AlgoimUtilsTests/VisualizationTests.jl +++ b/test/AlgoimUtilsTests/VisualizationTests.jl @@ -17,7 +17,7 @@ module VisualizationTests V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5)]) end - n = 12 + n = 14 degree = 3 cppdegree = -1 From 8ab1ebeb64a1373ddd3281ed2d5fc3885d7be65a Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Thu, 4 Jan 2024 18:50:22 +0100 Subject: [PATCH 4/8] Update NEWS.md --- NEWS.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index eeb9b68..e149230 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,7 +4,12 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.8.3] - 2024-08-02 +## Unreleased + +### Added +- Interface to Algoim v0.2, which provides algoim's high-order quadrature algorithms for domains implicitly-defined by multivariate polynomials and high-order accurate algorithms for computing closest points on implicitly-defined surfaces. Since PR [#76](https://github.com/gridap/GridapEmbedded.jl/pull/76). + +## [0.8.3] - 2024-01-02 ### Added - Support for FillArrays v1. Since PR [#75](https://github.com/gridap/GridapEmbedded.jl/pull/75). From e5d54765d3b26e060f5ce92d00dfab33e742726e Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Fri, 5 Jan 2024 19:30:11 +0100 Subject: [PATCH 5/8] Run tests for AlgoimUtils only in x86 --- .github/workflows/ci.yml | 35 +++++++++++++++++++++++++++++++ test/AlgoimUtilsTests/runtests.jl | 14 ++++++------- test/runtests.jl | 2 -- 3 files changed, 42 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b047a8d..f316fd1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,6 +70,41 @@ jobs: - run: | julia --color=yes --project=. --check-bounds=yes --depwarn=error -e ' (1,) .== 1; include("test/GridapEmbeddedTests/runtests.jl")' + algoimutils: + name: AlgoimUtils ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.8' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - run: | + julia --color=yes --project=. --check-bounds=yes --depwarn=error -e ' + using Pkg; Pkg.instantiate()' + - run: | + julia --color=yes --project=. --check-bounds=yes --depwarn=error -e ' + (1,) .== 1; include("test/AlgoimUtilsTests/runtests.jl")' docs: name: Documentation runs-on: ubuntu-latest diff --git a/test/AlgoimUtilsTests/runtests.jl b/test/AlgoimUtilsTests/runtests.jl index f775609..f24d881 100644 --- a/test/AlgoimUtilsTests/runtests.jl +++ b/test/AlgoimUtilsTests/runtests.jl @@ -2,18 +2,18 @@ module AlgoimTests using Test -@testset "AlgoimInterface" begin include("AlgoimInterfaceTests.jl") end +@time @testset "AlgoimInterface" begin include("AlgoimInterfaceTests.jl") end -@testset "QuadratureDegree" begin include("QuadratureDegreeTests.jl") end +@time @testset "QuadratureDegree" begin include("QuadratureDegreeTests.jl") end -@testset "DualQuadratures" begin include("DualQuadraturesTests.jl") end +@time @testset "DualQuadratures" begin include("DualQuadraturesTests.jl") end -@testset "PoissonAlgoim" begin include("PoissonAlgoimTests.jl") end +@time @testset "PoissonAlgoim" begin include("PoissonAlgoimTests.jl") end -@testset "ClosestPoint" begin include("ClosestPointTests.jl") end +@time @testset "ClosestPoint" begin include("ClosestPointTests.jl") end -@testset "VolumeConservation" begin include("VolumeConservationTests.jl") end +@time @testset "VolumeConservation" begin include("VolumeConservationTests.jl") end -@testset "Visualization" begin include("VisualizationTests.jl") end +@time @testset "Visualization" begin include("VisualizationTests.jl") end end # module \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 83cf4d2..8b1363f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,8 +21,6 @@ if QHULL_LOADED @time @testset "MomentFittedQuadratures" begin include("MomentFittedQuadraturesTests/runtests.jl") end - @time @testset "AlgoimUtils" begin include("AlgoimUtilsTests/runtests.jl") end - include(joinpath(@__DIR__,"..","examples","runexamples.jl")) else From 67a271077b059dd4eb6a0b27c35eef5ff7aef459 Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 15 Jan 2024 13:35:55 +0100 Subject: [PATCH 6/8] Optimised Algoim visualisation --- src/AlgoimUtils/AlgoimUtils.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl index 9f057c3..c4d36d1 100644 --- a/src/AlgoimUtils/AlgoimUtils.jl +++ b/src/AlgoimUtils/AlgoimUtils.jl @@ -379,7 +379,9 @@ function visualization_data(meas::Measure,filename;cellfields=Dict(),qhulltype=D grid = _to_grid(node_coordinates,qhulltype) ndata = Dict() for (k,v) in cellfields - ndata[k] = lazy_map(v,node_coordinates) + pts = get_cell_points(meas) + eval = evaluate(v,pts) + ndata[k] = collect(Iterators.flatten(eval)) end visualization_data(grid,filename,nodaldata=ndata) end @@ -389,7 +391,9 @@ function visualization_data(meas::Vector{<:Measure},filename;cellfields=Dict(),q grid = _to_grid(node_coordinates,qhulltype) ndata = Dict() for (k,v) in cellfields - ndata[k] = lazy_map(v,node_coordinates) + pts = map(m->get_cell_points(m),meas) + eval = map(p->evaluate(v,p),pts) + ndata[k] = vcat(map(e->collect(Iterators.flatten(e)),eval)...) end visualization_data(grid,filename,nodaldata=ndata) end From bf8f9da0d53e3c76706af482c423c8a5cad1c431 Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 22 Jan 2024 11:47:50 +0100 Subject: [PATCH 7/8] Updated authors and news --- LICENSE.md | 38 +++++++++++++++++++++++++++----------- NEWS.md | 2 +- Project.toml | 4 ++-- docs/make.jl | 2 +- 4 files changed, 31 insertions(+), 15 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index bc8fbce..41792ce 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,7 +1,7 @@ GridapEmbedded.jl Copyright and License == -Copyright (c) 2020 by [Francesc Verdugo](mailto:fverdugo@cimne.upc.edu) and [Santiago Badia](mailto:santiago.badia@monash.edu). +Copyright (c) 2020 by [Francesc Verdugo](mailto:f.verdugo.rojano@vu.nl), [Eric Neiva](mailto:eric.neiva@college-de-france.fr) and [Santiago Badia](mailto:santiago.badia@monash.edu). GridapEmbedded.jl is licensed under the MIT Expat License @@ -28,19 +28,35 @@ SOFTWARE. Citing `Gridap` == -In order to give credit to the `Gridap` contributors, we simply ask you to cite the refence below in any publication in which you have made use of `Gridap` packages. +In order to give credit to the `Gridap` contributors, we simply ask you to cite the references below in any publication in which you have made use of `Gridap` packages. F. Verdugo and S. Badia. *A user-guide to Gridap -- grid-based approximation of partial differential equations in Julia*, 2019. [arXiv:1910.01412](https://arxiv.org/abs/1910.01412) ``` -@article{gridap_guide_2019, - author={Francesc Verdugo and Santiago Badia}, - journal = {{arXiv}}, - title = {{A user-guide to Gridap -- grid-based approximation of partial differential equations in Julia}}, - year = {2019}, - eprint={1910.01412}, - archivePrefix={arXiv}, - primaryClass={cs.MS}, +@article{Badia2020, + doi = {10.21105/joss.02520}, + url = {https://doi.org/10.21105/joss.02520}, + year = {2020}, + publisher = {The Open Journal}, + volume = {5}, + number = {52}, + pages = {2520}, + author = {Santiago Badia and Francesc Verdugo}, + title = {Gridap: An extensible Finite Element toolbox in Julia}, + journal = {Journal of Open Source Software} +} + +@article{Verdugo2022, + doi = {10.1016/j.cpc.2022.108341}, + url = {https://doi.org/10.1016/j.cpc.2022.108341}, + year = {2022}, + month = jul, + publisher = {Elsevier {BV}}, + volume = {276}, + pages = {108341}, + author = {Francesc Verdugo and Santiago Badia}, + title = {The software design of Gridap: A Finite Element package based on the Julia {JIT} compiler}, + journal = {Computer Physics Communications} } ``` @@ -92,4 +108,4 @@ By making a contribution to this project, I certify that: Contact == -Please, contact the project administrators, [Santiago Badia](mailto:santiago.badia@monash.edu) and [Francesc Verdugo](mailto:fverdugo@cimne.upc.edu), for further questions about licenses and terms of use. +Please, contact the project administrators, [Francesc Verdugo](mailto:f.verdugo.rojano@vu.nl), [Eric Neiva](mailto:eric.neiva@college-de-france.fr) and [Santiago Badia](mailto:santiago.badia@monash.edu), for further questions about licenses and terms of use. diff --git a/NEWS.md b/NEWS.md index e149230..82c962d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## Unreleased +## [0.9.0] - 2024-01-22 ### Added - Interface to Algoim v0.2, which provides algoim's high-order quadrature algorithms for domains implicitly-defined by multivariate polynomials and high-order accurate algorithms for computing closest points on implicitly-defined surfaces. Since PR [#76](https://github.com/gridap/GridapEmbedded.jl/pull/76). diff --git a/Project.toml b/Project.toml index 3d20492..2a82563 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapEmbedded" uuid = "8838a6a3-0006-4405-b874-385995508d5d" -authors = ["Francesc Verdugo "] -version = "0.8.3" +authors = ["Francesc Verdugo ", "Eric Neiva ", "Santiago Badia "] +version = "0.9.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/docs/make.jl b/docs/make.jl index c09398c..3bad67a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,7 +8,7 @@ makedocs(; ], repo="https://github.com/gridap/GridapEmbedded.jl/blob/{commit}{path}#L{line}", sitename="GridapEmbedded.jl", - authors="Francesc Verdugo ", + authors="Francesc Verdugo , Eric Neiva and Santiago Badia ", ) deploydocs(; From cea454fe5564060414c06b717e9a4a84f4018906 Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 22 Jan 2024 11:49:09 +0100 Subject: [PATCH 8/8] Small amendment to citing --- LICENSE.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index 41792ce..652d476 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -30,8 +30,6 @@ Citing `Gridap` In order to give credit to the `Gridap` contributors, we simply ask you to cite the references below in any publication in which you have made use of `Gridap` packages. -F. Verdugo and S. Badia. *A user-guide to Gridap -- grid-based approximation of partial differential equations in Julia*, 2019. [arXiv:1910.01412](https://arxiv.org/abs/1910.01412) - ``` @article{Badia2020, doi = {10.21105/joss.02520},