From c7eeb9477667f3f5a88538c57aff58e5802c9e30 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 3 Dec 2024 09:45:29 +1100 Subject: [PATCH 01/10] Some optimisations to NVB --- src/Adaptivity/EdgeBasedRefinement.jl | 69 +++++++++++++++------------ test/AdaptivityTests/estimators.jl | 18 +++++++ 2 files changed, 56 insertions(+), 31 deletions(-) create mode 100644 test/AdaptivityTests/estimators.jl diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl index 724f74fd4..fc50d45ea 100644 --- a/src/Adaptivity/EdgeBasedRefinement.jl +++ b/src/Adaptivity/EdgeBasedRefinement.jl @@ -78,7 +78,6 @@ function refine_edge_based_topology( c2n_map_new = get_refined_cell_to_vertex_map(topo,rrules,faces_list) polys_new, cell_type_new = _get_cell_polytopes(rrules) orientation = NonOriented() - return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,orientation) end @@ -109,17 +108,28 @@ The new vertices are ordered by parent dimension and face id (in that order). function get_new_coordinates_from_faces(p::Union{Polytope{D},GridTopology{D}},faces_list::Tuple) where {D} @check length(faces_list) == D+1 - nN_new = sum(x->length(x),faces_list) + nN_new = sum(length,faces_list) coords_old = get_vertex_coordinates(p) coords_new = Vector{eltype(coords_old)}(undef,nN_new) n = 1 - for (d,dfaces) in enumerate(faces_list) - if length(dfaces) > 0 - nf = length(dfaces) - d2n_map = get_faces(p,d-1,0) - coords_new[n:n+nf-1] .= map(f -> sum(coords_old[d2n_map[f]])/length(d2n_map[f]), dfaces) - n += nf + # Nodes + if !isempty(faces_list[1]) + for node in faces_list[1] + coords_new[n] = coords_old[node] + n += 1 + end + end + # Faces (d > 0) + for (d,dfaces) in enumerate(faces_list[2:end]) + if !isempty(dfaces) + d2n_map = get_faces(p,d,0) + cache = array_cache(d2n_map) + for face in dfaces + face_nodes = getindex!(cache,d2n_map,face) + coords_new[n] = sum(coords_old[face_nodes])/length(face_nodes) + n += 1 + end end end @@ -255,12 +265,13 @@ function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTo setup_edge_based_rrules(method, topo, collect(1:num_faces(topo,Dc))) end -function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc +function setup_edge_based_rrules( + method::NVBRefinement, topo::UnstructuredGridTopology{Dc}, cells_to_refine::AbstractArray{<:Integer} +) where Dc nE = num_faces(topo,1) - c2e_map = get_faces(topo,Dc,1) - c2e_map_cache = array_cache(c2e_map) - e2c_map = get_faces(topo,1,Dc) - polys = topo.polytopes + c2e_map = get_faces(topo,Dc,1) + e2c_map = get_faces(topo,1,Dc) + polys = topo.polytopes cell_types = topo.cell_type cell_color = copy(cell_types) # WHITE # Hardcoded for TRI @@ -274,43 +285,39 @@ function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTo c_to_longest_edge_lid = method.cell_to_longest_edge_lid is_refined = falses(nE) # Loop over cells and mark edges to refine i.e. is_refined - # The reason to not loop directly on c is that we need to change c within - # a given iteration of the for loop - for i in 1:length(cells_to_refine) - c = cells_to_refine[i] - e_longest = c_to_longest_edge_gid[c] + for i in eachindex(cells_to_refine) # Has to terminate because an edge is marked each iteration or we skip an # iteration due to a boundary cell - while !is_refined[e_longest] - is_refined[e_longest] = true - c_nbor_lid = findfirst(c′ -> c′ != c, e2c_map[e_longest]) - if isnothing(c_nbor_lid) # We've reach the boundary + c = cells_to_refine[i] + e = c_to_longest_edge_gid[c] + while !is_refined[e] + is_refined[e] = true + e_cells = view(e2c_map,e) + if length(e_cells) == 1 # We've reach the boundary continue else - # Get the longest edge of the neighbor - c_nbor_gid = e2c_map[e_longest][c_nbor_lid] - e_longest = c_to_longest_edge_gid[c_nbor_gid] - # Set the current cell gid to that of the neighbor - c = c_nbor_gid + # Propagate to neighboring cell + c = ifelse(e_cells[1] == c, e_cells[2], e_cells[1]) + e = c_to_longest_edge_gid[c] end end end # Loop over cells and refine based on marked edges for c in 1:length(c2e_map) - c_edges = getindex!(c2e_map_cache, c2e_map, c) + c_edges = view(c2e_map,c) refined_edge_lids = findall(is_refined[c_edges]) - # GREEN refinement because only one edge should be bisected if length(refined_edge_lids) == 1 + # GREEN refinement because only one edge should be bisected ref_edge = refined_edge_lids[1] cell_color[c] = GREEN + Int8(ref_edge-1) - # BLUE refinement: two bisected elseif length(refined_edge_lids) == 2 + # BLUE refinement: two bisected long_ref_edge_lid = c_to_longest_edge_lid[c] short_ref_edge_lid = setdiff(refined_edge_lids, long_ref_edge_lid)[1] blue_idx = BLUE_dict[(long_ref_edge_lid, short_ref_edge_lid)] cell_color[c] = BLUE + Int8(blue_idx - 1) - # DOUBLE BLUE refinement: three bisected edges (somewhat rare) elseif length(refined_edge_lids) == 3 + # DOUBLE BLUE refinement: three bisected edges (somewhat rare) long_ref_edge_lid = c_to_longest_edge_lid[c] cell_color[c] = BLUE_DOUBLE + Int(long_ref_edge_lid - 1) end diff --git a/test/AdaptivityTests/estimators.jl b/test/AdaptivityTests/estimators.jl new file mode 100644 index 000000000..f15a605d3 --- /dev/null +++ b/test/AdaptivityTests/estimators.jl @@ -0,0 +1,18 @@ + +using Gridap, Gridap.Geometry, Gridap.Adaptivity + +function LShapedModel(n) + model = CartesianDiscreteModel((0,1,0,1),(n,n)) + cell_coords = map(mean,get_cell_coordinates(model)) + l_shape_filter(x) = (x[1] < 0.5) || (x[2] < 0.5) + mask = map(l_shape_filter,cell_coords) + return simplexify(DiscreteModelPortion(model,mask)) +end + +model = LShapedModel(10) + +method = Adaptivity.NVBRefinement(model) +cells_to_refine = [collect(1:10)...,collect(20:30)...] +fmodel = refine(method,model;cells_to_refine) + +writevtk(Triangulation(fmodel),"tmp/fmodel";append=false) From c73e24f0549ff53220dbc97a04346d9ccb51ad46 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 16:15:28 +1100 Subject: [PATCH 02/10] Added DorflerMarking --- src/Adaptivity/AdaptiveMeshRefinement.jl | 97 +++++++++++++++++++ src/Adaptivity/Adaptivity.jl | 3 + ...tors.jl => AdaptiveMeshRefinementTests.jl} | 25 ++++- 3 files changed, 124 insertions(+), 1 deletion(-) create mode 100644 src/Adaptivity/AdaptiveMeshRefinement.jl rename test/AdaptivityTests/{estimators.jl => AdaptiveMeshRefinementTests.jl} (52%) diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl new file mode 100644 index 000000000..5b57d24d8 --- /dev/null +++ b/src/Adaptivity/AdaptiveMeshRefinement.jl @@ -0,0 +1,97 @@ + +struct DorflerMarking + θ :: Float64 + ν :: Float64 + strategy :: Symbol + function DorflerMarking( + θ::Float64; + ν::Float64 = 0.5, + strategy::Symbol = :sort + ) + @assert 0 < θ < 1 + @assert strategy ∈ (:sort,:binsort,:quickmark) "Strategy not recognized. Available values are (:sort)" + new(θ,ν,strategy) + end +end + +mark(m::DorflerMarking, η::Vector{<:Real}) = mark(Val(m.strategy), m, η) + +function mark(::Val{:sort}, m::DorflerMarking, η::Vector{<:Real}) + target = m.θ * sum(η) + perm = sortperm(η, rev=true, alg=QuickSort) + s = zero(eltype(η)) + k = 0 + while s < target + k += 1 + s += η[perm[k]] + end + return perm[1:k] +end + +function mark(::Val{:binsort}, m::DorflerMarking, η::Vector{<:Real}) + target = m.θ * sum(η) + M = maximum(η) + N = length(η) + + # Find minimal K such that + # νᴷ⁺¹ M ≤ (1 - θ) / (θ * N) * target + K = 0 + a = M*m.ν + b = (1 - m.θ) / (m.θ * N * M) * target + while a > b + K += 1 + a *= m.ν + end + + # Sort into bins Bk such that ηi ∈ Bk if + # νᴷ⁺¹ M ≤ ηi < νᴷ M + bins = zeros(Int,N) + sums = zeros(Float64,K+1) + lbs = zeros(Float64,K+1) + lbs[1] = M * m.ν + for i in 2:K + lbs[i] = lbs[i-1] * m.ν + end + lbs[end] = 0.0 + + for (i,ηi) in enumerate(η) + k = 1 + while ηi < lbs[k] + k += 1 + end + bins[i] = k + sums[k] += ηi + end + + # Find minimal set of bins that gets over target + k = 0 + s = 0.0 + while s < target + k += 1 + s += sums[k] + end + + return findall(i -> i <= k, bins) +end + +function mark(::Val{:quickmark}, m::DorflerMarking, η::Vector{<:Real}) + function quickmark!(η, perm, l, u, target) :: Int + m = (u - l) ÷ 2 + sort!(view(perm,l:u), by=i->η[i], rev=true, alg=PartialQuickSort(m)) + + p = l + m + t = η[perm[p]] + σ = sum(η[perm[l:p-1]]) + + (σ >= target) && return quickmark!(η, perm, l, p, target) + (σ + t >= target) && return p + return quickmark!(η, perm, p + 1, u, target - σ - t) + end + + N = length(η) + l, u = 1, N + perm = collect(1:N) + target = m.θ * sum(η) + m = quickmark!(η, perm, l, u, target) + return perm[1:m] +end diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl index d80eb0c3e..5b85e0460 100644 --- a/src/Adaptivity/Adaptivity.jl +++ b/src/Adaptivity/Adaptivity.jl @@ -40,6 +40,8 @@ export AdaptedTriangulation export Triangulation, is_change_possible, best_target, get_adapted_model export change_domain, move_contributions +export DorflerMarking, mark + include("RefinementRules.jl") include("FineToCoarseFields.jl") include("OldToNewFields.jl") @@ -51,5 +53,6 @@ include("MacroFEs.jl") include("CompositeQuadratures.jl") include("EdgeBasedRefinement.jl") include("SimplexifyRefinement.jl") +include("AdaptiveMeshRefinement.jl") end # module diff --git a/test/AdaptivityTests/estimators.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl similarity index 52% rename from test/AdaptivityTests/estimators.jl rename to test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index f15a605d3..da240a6a7 100644 --- a/test/AdaptivityTests/estimators.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -1,5 +1,27 @@ - +using Test using Gridap, Gridap.Geometry, Gridap.Adaptivity +using DataStructures + +# Marking tests + +function test_marking() + for strategy in (:sort,:binsort,:quickmark) + for n in (1000,10000) + for θ in (0.3,0.5) + η = abs.(randn(n)) + m = DorflerMarking(θ;strategy) + I = Adaptivity.mark(m,η) + @test sum(η[I]) > θ * sum(η) + println("Strategy: $strategy, θ: $θ, n: $n") + println(" > N marked = $(length(I)), val marked = $(sum(η[I]) / sum(η))") + end + end + end +end + +test_marking() + +# AMR tests function LShapedModel(n) model = CartesianDiscreteModel((0,1,0,1),(n,n)) @@ -16,3 +38,4 @@ cells_to_refine = [collect(1:10)...,collect(20:30)...] fmodel = refine(method,model;cells_to_refine) writevtk(Triangulation(fmodel),"tmp/fmodel";append=false) + From 31175b04805faa78faff176c68b1da380ab9ec7a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 16:30:47 +1100 Subject: [PATCH 03/10] Added documentation --- src/Adaptivity/AdaptiveMeshRefinement.jl | 44 ++++++++++++++++++- .../AdaptiveMeshRefinementTests.jl | 4 +- 2 files changed, 45 insertions(+), 3 deletions(-) diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl index 5b57d24d8..2f5d83f34 100644 --- a/src/Adaptivity/AdaptiveMeshRefinement.jl +++ b/src/Adaptivity/AdaptiveMeshRefinement.jl @@ -1,4 +1,46 @@ +""" + struct DorflerMarking + θ :: Float64 + ν :: Float64 + strategy :: Symbol + end + + DorflerMarking(θ::Float64; ν::Float64 = 0.5, strategy::Symbol = :quickmark) + +Implements the Dorfler marking strategy. Given a vector `η` of real positive numbers, +the marking strategy find a subset of indices `I` such that + + sum(η[I]) > θ * sum(η) + +where `0 < θ < 1` is a threshold parameter. + +For more details, see the following reference: + +`Dörfler marking with minimal cardinality is a linear complexity problem`, Pfeiler et al. (2020) + +The marking algorithm is controlled by the `strategy` parameter, which can take +the following values: + +- `:sort`: Optimal cardinality, O(N log N) complexity. See Algorithm 2 in the reference. +- `:binsort`: Quasi-optimal cardinality, O(N) complexity. See Algorithm 7 in the reference. +- `:quickmark`: Optimal cardinality, O(N) complexity. See Algorithm 10 in the reference. + +# Arguments + +- `θ::Float64`: The threshold parameter. Between 0 and 1. +- `ν::Float64`: Extra parameter for `:binsort`. Default is 0.5. +- `strategy::Symbol`: The marking strategy. Default is `:quickmark`. + +# Usage + +```julia +η = abs.(randn(1000)) +m = DorflerMarking(0.5) +I = mark(m,η) +``` + +""" struct DorflerMarking θ :: Float64 ν :: Float64 @@ -6,7 +48,7 @@ struct DorflerMarking function DorflerMarking( θ::Float64; ν::Float64 = 0.5, - strategy::Symbol = :sort + strategy::Symbol = :quickmark ) @assert 0 < θ < 1 @assert strategy ∈ (:sort,:binsort,:quickmark) "Strategy not recognized. Available values are (:sort)" diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index da240a6a7..40776de45 100644 --- a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -4,7 +4,7 @@ using DataStructures # Marking tests -function test_marking() +function test_dorfler_marking() for strategy in (:sort,:binsort,:quickmark) for n in (1000,10000) for θ in (0.3,0.5) @@ -19,7 +19,7 @@ function test_marking() end end -test_marking() +@testset "Dorfler marking" test_dorfler_marking() # AMR tests From 60be5128409aebbdf5ced02e8e3c714068ae4f12 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 17:27:19 +1100 Subject: [PATCH 04/10] Added amr driver --- src/Adaptivity/AdaptiveMeshRefinement.jl | 42 +++++++++++++++ src/Adaptivity/Adaptivity.jl | 2 +- .../AdaptiveMeshRefinementTests.jl | 51 +++++++++++++++++-- 3 files changed, 89 insertions(+), 6 deletions(-) diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl index 2f5d83f34..0dfc9173d 100644 --- a/src/Adaptivity/AdaptiveMeshRefinement.jl +++ b/src/Adaptivity/AdaptiveMeshRefinement.jl @@ -1,4 +1,6 @@ +# Marking strategies + """ struct DorflerMarking θ :: Float64 @@ -56,6 +58,12 @@ struct DorflerMarking end end +""" + mark(m::DorflerMarking, η::Vector{<:Real}) -> Vector{Int} + +Given a vector `η` of real positive numbers, returns a subset of indices `I` such that +satisfying the Dorfler marking condition. +""" mark(m::DorflerMarking, η::Vector{<:Real}) = mark(Val(m.strategy), m, η) function mark(::Val{:sort}, m::DorflerMarking, η::Vector{<:Real}) @@ -137,3 +145,37 @@ function mark(::Val{:quickmark}, m::DorflerMarking, η::Vector{<:Real}) m = quickmark!(η, perm, l, u, target) return perm[1:m] end + +# Estimators + +function estimate(f::Function, uh) + collect_estimator(f(uh)) +end + +function collect_estimator(c::DomainContribution) + trians = get_domains(c) + bgmodel = get_background_model(first(trians)) + msg = "Estimator not implemented for mixed background models" + @notimplementedif !all([bgmodel == get_background_model(trian) for trian in trians]) msg + + Dc = num_cell_dims(bgmodel) + η = zeros(Float64,num_cells(bgmodel)) + for trian in trians + glue = get_glue(trian,Val(Dc)) + collect_estimator!(η,glue,get_contribution(c,trian)) + end + + return η +end + +function collect_estimator!(η, glue::FaceToFaceGlue, c) + cache = array_cache(c) + for (face,bgcell) in enumerate(glue.tface_to_mface) + η[bgcell] += getindex!(cache,c,face) + end +end + +function collect_estimator!(η, glue::SkeletonPair, c) + collect_estimator!(η,glue.plus,c) + collect_estimator!(η,glue.minus,c) +end diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl index 5b85e0460..32a5d1594 100644 --- a/src/Adaptivity/Adaptivity.jl +++ b/src/Adaptivity/Adaptivity.jl @@ -40,7 +40,7 @@ export AdaptedTriangulation export Triangulation, is_change_possible, best_target, get_adapted_model export change_domain, move_contributions -export DorflerMarking, mark +export DorflerMarking, mark, estimate include("RefinementRules.jl") include("FineToCoarseFields.jl") diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index 40776de45..6a0165941 100644 --- a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -31,11 +31,52 @@ function LShapedModel(n) return simplexify(DiscreteModelPortion(model,mask)) end -model = LShapedModel(10) +function amr_step(model) + order = 1 + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe) + + Ω = Triangulation(model) + Γ = Boundary(model) + Λ = Skeleton(model) + + dΩ = Measure(Ω,2*order) + dΓ = Measure(Γ,2*order) + dΛ = Measure(Λ,2*order) + + hK = CellField(sqrt.(collect(get_array(∫(1)dΩ))),Ω) + + f(x) = 1.0 / ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2) + a(u,v) = ∫(∇(u)⋅∇(v))dΩ + l(v) = ∫(f*v)dΩ + ηh(u) = ∫(hK*f)dΩ + ∫(hK*∇(u)⋅∇(u))dΓ + ∫(hK*jump(∇(u))⋅jump(∇(u)))dΛ + + op = AffineFEOperator(a,l,V,V) + uh = solve(op) + η = estimate(ηh,uh) + + m = DorflerMarking(0.5) + I = Adaptivity.mark(m,η) + + method = Adaptivity.NVBRefinement(model) + fmodel = Adaptivity.get_model(refine(method,model;cells_to_refine=I)) -method = Adaptivity.NVBRefinement(model) -cells_to_refine = [collect(1:10)...,collect(20:30)...] -fmodel = refine(method,model;cells_to_refine) + return fmodel, uh, η, I +end -writevtk(Triangulation(fmodel),"tmp/fmodel";append=false) +model = LShapedModel(10) +for i in 1:10 + fmodel, uh, η, I = amr_step(model) + is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model)) + Ω = Triangulation(model) + writevtk( + Ω,"tmp/model_$(i-1)",append=false, + cellfields = [ + "uh" => uh, + "η" => CellField(η,Ω), + "is_refined" => CellField(is_refined,Ω) + ], + ) + model = fmodel +end From 66d44a326b9c3b9fc00d5000268dc390de4b8397 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 17:30:59 +1100 Subject: [PATCH 05/10] Minor --- test/AdaptivityTests/AdaptiveMeshRefinementTests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index 6a0165941..81fcc3dcc 100644 --- a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -64,9 +64,10 @@ function amr_step(model) return fmodel, uh, η, I end +nsteps = 10 model = LShapedModel(10) -for i in 1:10 +for i in 1:nsteps fmodel, uh, η, I = amr_step(model) is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model)) Ω = Triangulation(model) From 2ba6583774c9faeb5c8f9a4a31d3afd364766dc2 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 17:43:46 +1100 Subject: [PATCH 06/10] More docs --- docs/src/Adaptivity.md | 12 ++++++++++++ src/Adaptivity/AdaptiveMeshRefinement.jl | 7 +++++++ 2 files changed, 19 insertions(+) diff --git a/docs/src/Adaptivity.md b/docs/src/Adaptivity.md index 377cd5d23..5a698528a 100644 --- a/docs/src/Adaptivity.md +++ b/docs/src/Adaptivity.md @@ -80,6 +80,18 @@ The API is given by the following methods: MacroReferenceFE ``` +## Adaptive Mesh Refinement + +One of the main uses of mesh refinement is Adaptive Mesh Refinement, where the mesh is refined only in regions of interest. + +The typical AMR workflow is the so-called `solve-estimate-mark-refine` loop. Since estimators will generally be problem-dependent, we only aim to provide some generic tools that can be combined by the user: + +```@docs + DorflerMarking + mark + estimate +``` + ## Notes for users Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon). diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl index 0dfc9173d..78ca189f0 100644 --- a/src/Adaptivity/AdaptiveMeshRefinement.jl +++ b/src/Adaptivity/AdaptiveMeshRefinement.jl @@ -148,6 +148,13 @@ end # Estimators +""" + estimate(f::Function, uh::Function) -> Vector{Float64} + +Given a functional `f` and a function `uh`, such that `f(uh)` produces a +scalar-valued `DomainContribution`, collects the estimator values for +each cell in the background model. +""" function estimate(f::Function, uh) collect_estimator(f(uh)) end From ce22f58316751f01d62e1804f8b26bc7ae183535 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 6 Dec 2024 17:49:38 +1100 Subject: [PATCH 07/10] Minor --- src/Adaptivity/AdaptiveMeshRefinement.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl index 78ca189f0..40f9eea42 100644 --- a/src/Adaptivity/AdaptiveMeshRefinement.jl +++ b/src/Adaptivity/AdaptiveMeshRefinement.jl @@ -2,13 +2,13 @@ # Marking strategies """ - struct DorflerMarking - θ :: Float64 - ν :: Float64 - strategy :: Symbol - end + struct DorflerMarking + θ :: Float64 + ν :: Float64 + strategy :: Symbol + end - DorflerMarking(θ::Float64; ν::Float64 = 0.5, strategy::Symbol = :quickmark) + DorflerMarking(θ::Float64; ν::Float64 = 0.5, strategy::Symbol = :quickmark) Implements the Dorfler marking strategy. Given a vector `η` of real positive numbers, the marking strategy find a subset of indices `I` such that @@ -19,7 +19,7 @@ where `0 < θ < 1` is a threshold parameter. For more details, see the following reference: -`Dörfler marking with minimal cardinality is a linear complexity problem`, Pfeiler et al. (2020) +"Dörfler marking with minimal cardinality is a linear complexity problem", Pfeiler et al. (2020) The marking algorithm is controlled by the `strategy` parameter, which can take the following values: From 9cf8f080ee01e29127575c185ce6663353f4d265 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 9 Dec 2024 12:55:16 +1100 Subject: [PATCH 08/10] Added tests --- .../AdaptiveMeshRefinementTests.jl | 83 +++++++++++++------ test/AdaptivityTests/runtests.jl | 4 + 2 files changed, 60 insertions(+), 27 deletions(-) diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index 81fcc3dcc..c512e11ae 100644 --- a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -1,3 +1,5 @@ +module AdaptiveMeshRefinementTests + using Test using Gridap, Gridap.Geometry, Gridap.Adaptivity using DataStructures @@ -19,8 +21,6 @@ function test_dorfler_marking() end end -@testset "Dorfler marking" test_dorfler_marking() - # AMR tests function LShapedModel(n) @@ -31,53 +31,82 @@ function LShapedModel(n) return simplexify(DiscreteModelPortion(model,mask)) end -function amr_step(model) - order = 1 +l2_norm(he,xh,dΩ) = ∫(he*(xh*xh))*dΩ +l2_norm(xh,dΩ) = ∫(xh*xh)*dΩ + +function amr_step(model,u_exact;order=1) reffe = ReferenceFE(lagrangian,Float64,order) - V = TestFESpace(model,reffe) + V = TestFESpace(model,reffe;dirichlet_tags=["boundary"]) + U = TrialFESpace(V,u_exact) Ω = Triangulation(model) Γ = Boundary(model) Λ = Skeleton(model) - dΩ = Measure(Ω,2*order) + dΩ = Measure(Ω,4*order) dΓ = Measure(Γ,2*order) dΛ = Measure(Λ,2*order) hK = CellField(sqrt.(collect(get_array(∫(1)dΩ))),Ω) - - f(x) = 1.0 / ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2) + + nΓ = get_normal_vector(Γ) + nΛ = get_normal_vector(Λ) + + ∇u(x) = ∇(u_exact)(x) + f(x) = -Δ(u_exact)(x) a(u,v) = ∫(∇(u)⋅∇(v))dΩ l(v) = ∫(f*v)dΩ - ηh(u) = ∫(hK*f)dΩ + ∫(hK*∇(u)⋅∇(u))dΓ + ∫(hK*jump(∇(u))⋅jump(∇(u)))dΛ + ηh(u) = l2_norm(hK,(f + Δ(u)),dΩ) + l2_norm(hK,(∇(u) - ∇u)⋅nΓ,dΓ) + l2_norm(hK,jump(∇(u)⋅nΛ),dΛ) - op = AffineFEOperator(a,l,V,V) + op = AffineFEOperator(a,l,U,V) uh = solve(op) η = estimate(ηh,uh) - m = DorflerMarking(0.5) + m = DorflerMarking(0.8) I = Adaptivity.mark(m,η) method = Adaptivity.NVBRefinement(model) fmodel = Adaptivity.get_model(refine(method,model;cells_to_refine=I)) - return fmodel, uh, η, I + error = sum(l2_norm(uh - u_exact,dΩ)) + return fmodel, uh, η, I, error end -nsteps = 10 -model = LShapedModel(10) +function test_amr(nsteps,order) + model = LShapedModel(10) -for i in 1:nsteps - fmodel, uh, η, I = amr_step(model) - is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model)) - Ω = Triangulation(model) - writevtk( - Ω,"tmp/model_$(i-1)",append=false, - cellfields = [ - "uh" => uh, - "η" => CellField(η,Ω), - "is_refined" => CellField(is_refined,Ω) - ], - ) - model = fmodel + ϵ = 1e-2 + r(x) = ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2) + u_exact(x) = 1.0 / (ϵ + r(x)) + + vtk = false + last_error = Inf + for i in 1:nsteps + fmodel, uh, η, I, error = amr_step(model,u_exact;order) + if vtk + is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model)) + Ω = Triangulation(model) + writevtk( + Ω,"tmp/model_$(i-1)",append=false, + cellfields = [ + "uh" => uh, + "η" => CellField(η,Ω), + "is_refined" => CellField(is_refined,Ω), + "u_exact" => CellField(u_exact,Ω), + ], + ) + end + + println("Error: $error, Error η: $(sum(η))") + @test (i < 3) || (error < last_error) + last_error = error + model = fmodel + end end + +############################################################################################ + +@testset "Dorfler marking" test_dorfler_marking() +@testset "AMR - Poisson" test_amr(20,2) + +end # module \ No newline at end of file diff --git a/test/AdaptivityTests/runtests.jl b/test/AdaptivityTests/runtests.jl index 87437afc3..2d2a569da 100644 --- a/test/AdaptivityTests/runtests.jl +++ b/test/AdaptivityTests/runtests.jl @@ -23,4 +23,8 @@ end include("MacroFEStokesTests.jl") end +@testset "AMR" begin + include("AdaptiveMeshRefinementTests.jl") +end + end # module \ No newline at end of file From c5969723c5125eda6fe4bec050c412e6b6285652 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 9 Dec 2024 15:23:22 +1100 Subject: [PATCH 09/10] Minor --- test/AdaptivityTests/AdaptiveMeshRefinementTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl index c512e11ae..6f0915d0e 100644 --- a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl +++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl @@ -56,7 +56,7 @@ function amr_step(model,u_exact;order=1) f(x) = -Δ(u_exact)(x) a(u,v) = ∫(∇(u)⋅∇(v))dΩ l(v) = ∫(f*v)dΩ - ηh(u) = l2_norm(hK,(f + Δ(u)),dΩ) + l2_norm(hK,(∇(u) - ∇u)⋅nΓ,dΓ) + l2_norm(hK,jump(∇(u)⋅nΛ),dΛ) + ηh(u) = l2_norm(hK*(f + Δ(u)),dΩ) + l2_norm(hK*(∇(u) - ∇u)⋅nΓ,dΓ) + l2_norm(jump(hK*∇(u)⋅nΛ),dΛ) op = AffineFEOperator(a,l,U,V) uh = solve(op) From bcfba317cb0494482672e17420c00d942345122a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 9 Dec 2024 16:57:18 +1100 Subject: [PATCH 10/10] Updated NEWS --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index e9287dbdc..524a03951 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,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). +## [Unreleased] + +### Added + +- Added AMR-related methods `mark` and `estimate` to `Adaptivity` module. Implemented Dorfler marking strategy. Since PR[#1063](https://github.com/gridap/Gridap.jl/pull/1063). + ## [0.18.8] - 2024-12-2 ### Added