Added more tests

test/AdaptivityTests/macro_poisson.jl

@@ -0,0 +1,64 @@
+using Gridap
+import Gridap.ReferenceFEs: LagrangianRefFE, Quadrature
+import Gridap.Adaptivity: RedRefinementRule, MacroReferenceFE, CompositeQuadrature
+import Gridap.Adaptivity: get_polytope, num_subcells
+import FillArrays: Fill
+
+
+order, ncell = 2, 64
+bc_tags = Dict(:dirichlet_tags => "boundary")
+
+u((x, y)) = sin(3.2x * (x - y))cos(x + 4.3y) + sin(4.6 * (x + 2y))cos(2.6(y - 2x))
+f(x) = -Δ(u)(x)
+
+model = CartesianDiscreteModel((0, 1, 0, 1), (ncell, ncell))
+u_reffe = ReferenceFE(lagrangian, Float64, order)
+U = TrialFESpace(FESpace(model, u_reffe; bc_tags...), u)
+
+rrule = RedRefinementRule(QUAD)
+poly = get_polytope(rrule)
+reffe = LagrangianRefFE(Float64, poly, 1)
+sub_reffes = Fill(reffe, num_subcells(rrule))
+v_reffe = MacroReferenceFE(rrule, sub_reffes)
+V = FESpace(model, v_reffe; bc_tags...)
+
+Ω = Triangulation(model)
+macro_quad = Quadrature(poly, CompositeQuadrature(), rrule, 2order)
+dΩ = Measure(Ω, macro_quad)
+dΩ⁺ = Measure(Ω, 10)
+
+a(u, v) = ∫(∇(v) ⋅ ∇(u))dΩ
+l(v) = ∫(f * v)dΩ
+op = AffineFEOperator(a, l, U, V)
+eh = solve(op) - u
+fem_l2err = sqrt(∑(∫(eh * eh)dΩ))    # 0.5745, too large!
+
+a1(u, v) = ∫(u * v)dΩ
+l1(v) = ∫(u * v)dΩ
+op1 = AffineFEOperator(a1, l1, U, V)
+eh1 = solve(op1) - u
+fem_l2err1 = sqrt(∑(∫(eh1 * eh1)dΩ))  # 1.52e-5, reasonable
+
+############################################################################################
+
+l2_norm(a) = sum(∫(a⋅a)*dΩ)
+l2_error(a,b) = l2_norm(a-b)
+
+sol(x) = sum(x)
+V = FESpace(model, v_reffe)
+
+uh = interpolate(sol, V)
+l2_error(uh, sol)
+
+∇uh = gradient(uh)
+l2_error(∇uh, ∇(sol))
+
+Ωrr = Triangulation(rrule.ref_grid)
+
+cmaps = get_cell_map(rrule)
+cjacs = lazy_map(∇,cmaps)
+cinvjacs = lazy_map(inv,cjacs)
+
+pts = Gridap.CellData.get_data(get_cell_points(Ωrr))
+
+lazy_map(evaluate,cinvjacs,pts)

test/AdaptivityTests/macro_stokes.jl

@@ -0,0 +1,60 @@
+
+using Gridap
+using Gridap.Adaptivity, Gridap.Geometry, Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Arrays
+using Gridap.ReferenceFEs
+
+using FillArrays
+
+# Barycentric working for D = 2,3
+# PowellSabin working for D = 2 
+
+# Parameters
+Dc = 2
+reftype = 2
+
+u_sol(x) = (Dc == 2) ? VectorValue(x[1],-x[2]) : VectorValue(x[1],-x[2],0.0)
+p_sol(x) = (x[1] - 1.0/2.0)
+
+domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1)
+nc = (Dc == 2) ? (4,4) : (2,2,2)
+model = simplexify(CartesianDiscreteModel(domain,nc))
+
+min_order = (reftype == 1) ? Dc : Dc-1
+order = max(2,min_order)
+poly  = (Dc == 2) ? TRI : TET
+rrule = (reftype == 1) ? Adaptivity.BarycentricRefinementRule(poly) : Adaptivity.PowellSabinRefinementRule(poly)
+
+reffes = Fill(LagrangianRefFE(VectorValue{Dc,Float64},poly,order),Adaptivity.num_subcells(rrule))
+reffe_u = Adaptivity.MacroReferenceFE(rrule,reffes)
+reffe_p = LagrangianRefFE(Float64,poly,order-1)
+
+qdegree = 2*order
+quad  = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,qdegree)
+
+V = FESpace(model,reffe_u,dirichlet_tags=["boundary"])
+Q = FESpace(model,reffe_p,conformity=:L2,constraint=:zeromean)
+U = TrialFESpace(V,u_sol)
+
+Ω = Triangulation(model)
+dΩ = Measure(Ω,quad)
+@assert abs(sum(∫(p_sol)dΩ)) < 1.e-15
+@assert abs(sum(∫(divergence(u_sol))dΩ)) < 1.e-15
+
+X = MultiFieldFESpace([U,Q])
+Y = MultiFieldFESpace([V,Q])
+
+f(x) = -Δ(u_sol)(x) + ∇(p_sol)(x)
+lap(u,v) = ∫(∇(u)⊙∇(v))dΩ
+a((u,p),(v,q)) = lap(u,v) + ∫(divergence(u)*q - divergence(v)*p)dΩ
+l((v,q)) = ∫(f⋅v)dΩ
+
+op = AffineFEOperator(a,l,X,Y)
+xh = solve(op)
+uh, ph = xh
+sum(∫(uh⋅uh)dΩ)
+sum(∫(ph)dΩ)
+
+eh_u = uh - u_sol
+eh_p = ph - p_sol
+sum(∫(eh_u⋅eh_u)dΩ)
+sum(∫(eh_p*eh_p)dΩ)