Added amr driver

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

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")

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