|
| 1 | +module DarcyTests |
| 2 | + |
| 3 | +using Test |
| 4 | +using Gridap |
| 5 | +import Gridap: ∇, divergence |
| 6 | + |
| 7 | +u(x) = VectorValue(2*x[1],x[1]+x[2]) |
| 8 | + |
| 9 | +divergence(::typeof(u)) = (x) -> 3 |
| 10 | + |
| 11 | +p(x) = x[1]-x[2] |
| 12 | + |
| 13 | +∇p(x) = VectorValue(1,-1) |
| 14 | + |
| 15 | +∇(::typeof(p)) = ∇p |
| 16 | + |
| 17 | +f(x) = u(x) + ∇p(x) |
| 18 | + |
| 19 | +domain = (0,1,0,1) |
| 20 | +partition = (4,4) |
| 21 | +order = 2 |
| 22 | +model = CartesianDiscreteModel(domain,partition) |
| 23 | + |
| 24 | +V = FESpace( |
| 25 | + reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64}, |
| 26 | + conformity=:Hdiv, model=model, dirichlet_tags=[5,6]) |
| 27 | + |
| 28 | +Q = FESpace( |
| 29 | + reffe=:QLagrangian, order=order-1, valuetype=Float64, |
| 30 | + conformity=:L2, model=model) |
| 31 | + |
| 32 | +U = TrialFESpace(V,u) |
| 33 | +P = TrialFESpace(Q) |
| 34 | + |
| 35 | +Y = MultiFieldFESpace([V, Q]) |
| 36 | +X = MultiFieldFESpace([U, P]) |
| 37 | + |
| 38 | +trian = Triangulation(model) |
| 39 | +degree = 2 |
| 40 | +quad = CellQuadrature(trian,degree) |
| 41 | + |
| 42 | +neumanntags = [7,8] |
| 43 | +btrian = BoundaryTriangulation(model,neumanntags) |
| 44 | +degree = 2*order |
| 45 | +bquad = CellQuadrature(btrian,degree) |
| 46 | +nb = get_normal_vector(btrian) |
| 47 | + |
| 48 | +function a(y,x) |
| 49 | + v, q = y |
| 50 | + u, p = x |
| 51 | + v*u - (∇*v)*p + q*(∇*u) |
| 52 | +end |
| 53 | + |
| 54 | +function l(y) |
| 55 | + v, q = y |
| 56 | + v*f + q*(∇*u) |
| 57 | +end |
| 58 | + |
| 59 | +function l_Γ(y) |
| 60 | + v, q = y |
| 61 | + -(v*nb)*p |
| 62 | +end |
| 63 | + |
| 64 | +t_Ω = AffineFETerm(a,l,trian,quad) |
| 65 | +t_Γ = FESource(l_Γ,btrian,bquad) |
| 66 | +op = AffineFEOperator(Y,X,t_Ω,t_Γ) |
| 67 | +xh = solve(op) |
| 68 | +uh, ph = xh |
| 69 | + |
| 70 | +eu = u - uh |
| 71 | +ep = p - ph |
| 72 | + |
| 73 | +l2(v) = v*v |
| 74 | +h1(v) = v*v + ∇(v)*∇(v) |
| 75 | + |
| 76 | +eu_l2 = sum(integrate(l2(eu),trian,quad)) |
| 77 | +ep_l2 = sum(integrate(l2(ep),trian,quad)) |
| 78 | +ep_h1 = sum(integrate(h1(ep),trian,quad)) |
| 79 | + |
| 80 | +tol = 1.0e-9 |
| 81 | +@test eu_l2 < tol |
| 82 | +@test ep_l2 < tol |
| 83 | +@test ep_h1 < tol |
| 84 | + |
| 85 | +end # module |
0 commit comments