From 5341fdebc841a99a96174a1b268d81732ab31c6a Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Tue, 19 May 2026 09:16:01 +1000 Subject: [PATCH 1/7] Started implementing driver --- test/hdiv_stabilisation.jl | 94 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 test/hdiv_stabilisation.jl diff --git a/test/hdiv_stabilisation.jl b/test/hdiv_stabilisation.jl new file mode 100644 index 0000000..f9f43ee --- /dev/null +++ b/test/hdiv_stabilisation.jl @@ -0,0 +1,94 @@ + +using Gridap +using GridapEmbedded +using Gridap.Geometry, Gridap.Arrays +using Gridap.FESpaces, Gridap.CellData + +function get_aggregates(cutgeo) + bgmodel = get_background_model(cutgeo) + geo = get_geometry(cutgeo) + strategy = AggregateAllCutCells() + aggregates = aggregate(strategy,cutgeo,geo) + colors = color_aggregates(aggregates,bgmodel) + return inverse_table(colors) +end + +function projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + mass(u,v) = ∫(u⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, mass, mass + ) + return P +end + +R = 1.4 +geo = disk(R;x0=Point(1.0,1.7)) + +n = 8 +order = 1 +qdegree = 2*order + +bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) +cutgeo = cut(bgmodel,geo) + +Ω_ac = Triangulation(cutgeo,ACTIVE) +V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order)) +Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) +X = MultiFieldFESpace([V,Q]) + +ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) +Ωp = PatchTriangulation(bgmodel, ptopo) +dΩp = Measure(Ωp, qdegree) + +Wd = FESpaces.PatchFESpace(bgmodel,Ωp,VectorValue{2,Float64},order) +Πd = projection_operator(ptopo, V, W, Ωp, dΩp) +W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order) +Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) + +Ω = Triangulation(cutgeo,PHYSICAL) +Ω_cut = Triangulation(cutgeo,CUT_IN) +Γ = EmbeddedBoundary(cutgeo) +dΩ = Measure(Ω, qdegree) +dΩ_cut = Measure(Ω_cut, qdegree) +dΓ = Measure(Γ, 2*qdegree) + +η = 1.0 +u_exact(x) = VectorValue(x[1], -x[2]) +p_exact(x) = x[1]*x[2] +f(x) = η*u_exact(x) + ∇(p_exact)(x) +g(x) = -(∇⋅u_exact)(x) + +γ = 10.0 +hΓ = CellField(γ ./ get_array(∫(1.0)dΓ),Γ) +n_Γ = get_normal_vector(Γ) +a(u,v) = ∫(u⋅v)dΩ + ∫(hΓ*(u⋅n_Γ)*(v⋅n_Γ))dΓ +b(u,q) = ∫((∇⋅u)*q)dΩ +btilde(u,q) = b(u,q) + ∫((u⋅n_Γ)*q)dΓ +l((v,q)) = ∫(v⋅f + g⋅q)dΩ - ∫((v⋅n_Γ)*p_exact)dΓ + ∫(hΓ*(v⋅n_Γ)⋅(u_exact⋅n_Γ))dΓ + +τ = 1.0 +sd(u,v) = ∫(τ*(u⋅v))dΩ_cut +s0(p,q) = ∫(τ*p*q)dΩ_cut + +function weakform(x,y) + u, p = x + v, q = y + Πu, Πv = Π(u), Π(v) + Xp = FESpaces.PatchFESpace(X,ptopo) + data = FESpaces.collect_and_merge_cell_matrix_and_vector( + (X, X , a(u,v) + btilde(v,p) + b(u,q) + s(u,v), l(y), zero(X)), + (X, Xp , -1*s(u,Πv), DomainContribution(), zero(X)), + (Xp, X , -1*s(Πu,v), DomainContribution(), zero(Xp)), + (Xp, Xp, s(Πu,Πv), DomainContribution(), zero(Xp)), + ) + assem = SparseMatrixAssembler(X,X) + A, B = assemble_matrix_and_vector(assem,data) + return AffineFEOperator(X,X,A,B) +end + +op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) +uh, ph = solve(op) + +l2_err_u = sqrt(sum(∫( (uh-u_exact)⋅(uh-u_exact) )dΩ)) +l2_err_p = sqrt(sum(∫( (ph-p_exact)*(ph-p_exact) )dΩ)) From a28af055b50ea99f0e56cff2609db30995785d00 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Sun, 24 May 2026 23:23:26 +1000 Subject: [PATCH 2/7] Finished driver --- test/hdiv_stabilisation.jl | 68 ++++++++++++++++++++++++++++---------- 1 file changed, 51 insertions(+), 17 deletions(-) diff --git a/test/hdiv_stabilisation.jl b/test/hdiv_stabilisation.jl index f9f43ee..a21e229 100644 --- a/test/hdiv_stabilisation.jl +++ b/test/hdiv_stabilisation.jl @@ -22,12 +22,22 @@ function projection_operator(ptopo, V, W, Ω, dΩ) return P end -R = 1.4 -geo = disk(R;x0=Point(1.0,1.7)) +function div_projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + lhs(u,v) = ∫(u⋅Π(v,Ω))dΩ + rhs(u,v) = ∫((∇⋅u)⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, lhs, rhs + ) + return P +end + +R = 0.4 +geo = disk(R;x0=Point(0.5,0.5)) -n = 8 +n = 12 order = 1 -qdegree = 2*order +qdegree = 2*order+2 bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) cutgeo = cut(bgmodel,geo) @@ -41,10 +51,11 @@ ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) Ωp = PatchTriangulation(bgmodel, ptopo) dΩp = Measure(Ωp, qdegree) -Wd = FESpaces.PatchFESpace(bgmodel,Ωp,VectorValue{2,Float64},order) -Πd = projection_operator(ptopo, V, W, Ωp, dΩp) +Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) +Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order) Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) +Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) Ω = Triangulation(cutgeo,PHYSICAL) Ω_cut = Triangulation(cutgeo,CUT_IN) @@ -55,17 +66,17 @@ dΓ = Measure(Γ, 2*qdegree) η = 1.0 u_exact(x) = VectorValue(x[1], -x[2]) -p_exact(x) = x[1]*x[2] +p_exact(x) = x[1] + x[2] f(x) = η*u_exact(x) + ∇(p_exact)(x) g(x) = -(∇⋅u_exact)(x) γ = 10.0 hΓ = CellField(γ ./ get_array(∫(1.0)dΓ),Γ) n_Γ = get_normal_vector(Γ) -a(u,v) = ∫(u⋅v)dΩ + ∫(hΓ*(u⋅n_Γ)*(v⋅n_Γ))dΓ -b(u,q) = ∫((∇⋅u)*q)dΩ +a(u,v) = ∫(η*(u⋅v))dΩ + ∫(hΓ*(u⋅n_Γ)*(v⋅n_Γ))dΓ +b(u,q) = ∫(-(∇⋅u)*q)dΩ btilde(u,q) = b(u,q) + ∫((u⋅n_Γ)*q)dΓ -l((v,q)) = ∫(v⋅f + g⋅q)dΩ - ∫((v⋅n_Γ)*p_exact)dΓ + ∫(hΓ*(v⋅n_Γ)⋅(u_exact⋅n_Γ))dΓ +l((v,q)) = ∫(v⋅f + g⋅q)dΩ + ∫(hΓ*(v⋅n_Γ)⋅(u_exact⋅n_Γ))dΓ #- ∫((v⋅n_Γ)*p_exact)dΓ τ = 1.0 sd(u,v) = ∫(τ*(u⋅v))dΩ_cut @@ -74,13 +85,16 @@ s0(p,q) = ∫(τ*p*q)dΩ_cut function weakform(x,y) u, p = x v, q = y - Πu, Πv = Π(u), Π(v) + Πu, Πv = Πd(u), Πd(v) + Πp, Πq = Π0(p), Π0(q) + divu, divv = ∇⋅u, ∇⋅v + Πdivu, Πdivv = Π0div(u), Π0div(v) Xp = FESpaces.PatchFESpace(X,ptopo) data = FESpaces.collect_and_merge_cell_matrix_and_vector( - (X, X , a(u,v) + btilde(v,p) + b(u,q) + s(u,v), l(y), zero(X)), - (X, Xp , -1*s(u,Πv), DomainContribution(), zero(X)), - (Xp, X , -1*s(Πu,v), DomainContribution(), zero(Xp)), - (Xp, Xp, s(Πu,Πv), DomainContribution(), zero(Xp)), + (X, X , a(u,v) + btilde(v,p) + b(u,q) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y), zero(X)), + (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution(), zero(X)), + (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution(), zero(Xp)), + (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution(), zero(Xp)), ) assem = SparseMatrixAssembler(X,X) A, B = assemble_matrix_and_vector(assem,data) @@ -90,5 +104,25 @@ end op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) uh, ph = solve(op) -l2_err_u = sqrt(sum(∫( (uh-u_exact)⋅(uh-u_exact) )dΩ)) -l2_err_p = sqrt(sum(∫( (ph-p_exact)*(ph-p_exact) )dΩ)) +eu = uh-u_exact +ep = ph-p_exact +l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) +l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) + +writevtk( + Ω, "hdiv_cut"; append=false, + cellfields=[ + "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact,"eu"=>eu,"ep"=>ep + ], +) + +aggregates = get_aggregates(cutgeo) +cell_to_agg = flatten_partition(aggregates,num_cells(bgmodel)) +writevtk( + Triangulation(bgmodel), "aggregates"; append=false, + celldata = ["agg"=>cell_to_agg] +) + +writevtk( + Γ, "boundary"; append=false, +) From 22a4b70e5cef7fb5e4180f14da65a477bd3e5c3b Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 25 May 2026 21:47:00 +1000 Subject: [PATCH 3/7] Completed drivers --- test/darcy_BGP_flux.jl | 141 ++++++++++++++++++ ...iv_stabilisation.jl => darcy_BGP_mixed.jl} | 68 ++++++--- 2 files changed, 186 insertions(+), 23 deletions(-) create mode 100644 test/darcy_BGP_flux.jl rename test/{hdiv_stabilisation.jl => darcy_BGP_mixed.jl} (65%) diff --git a/test/darcy_BGP_flux.jl b/test/darcy_BGP_flux.jl new file mode 100644 index 0000000..bf27f15 --- /dev/null +++ b/test/darcy_BGP_flux.jl @@ -0,0 +1,141 @@ + +using Gridap +using GridapEmbedded +using Gridap.Geometry, Gridap.Arrays +using Gridap.FESpaces, Gridap.CellData + +function get_aggregates(cutgeo) + geo = get_geometry(cutgeo) + strategy = AggregateAllCutCells() + aggregates = aggregate(strategy,cutgeo,geo) + return Table(filter(agg -> length(agg) > 1, inverse_table(aggregates))) +end + +function projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + mass(u,v) = ∫(u⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, mass, mass + ) + return P +end + +function div_projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + lhs(u,v) = ∫(u⋅Π(v,Ω))dΩ + rhs(u,v) = ∫((∇⋅u)⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, lhs, rhs + ) + return P +end + +function generate_mesh(n,R=0.4,x0=0.835) + geo1 = disk(R;x0=Point(0.5,0.5),name="disk") + geo2 = plane(x0=Point(x0,0.5),v=Point(1.0,0.0),name="plane") + geo = intersect(geo1,geo2,name="domain") + + bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) + cutgeo = cut(bgmodel,geo) + return cutgeo +end + +n = 12 +cutgeo = generate_mesh(n) +bgmodel = get_background_model(cutgeo) + +order = 1 +qdegree = 2*(order+1) + +Ω_ac = Triangulation(cutgeo,ACTIVE) +V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) +Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) +Λ = ConstantFESpace(Ω_ac) +X = MultiFieldFESpace([V,Q,Λ]) + +ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) +Ωp = PatchTriangulation(bgmodel, ptopo) +dΩp = Measure(Ωp, qdegree) + +Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) +Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) +W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) +Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) +Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) + +Ω = Triangulation(cutgeo,PHYSICAL) +Ω_cut = Triangulation(cutgeo,CUT_IN) +Γu = EmbeddedBoundary(cutgeo) + +dΩ = Measure(Ω, qdegree) +dΩ_cut = Measure(Ω_cut, qdegree) +dΓu = Measure(Γu, 2*qdegree) + +η = 1.0 +u_exact(x) = VectorValue(x[1], -x[2]) +p_exact(x) = x[1] + x[2] +f(x) = η*u_exact(x) + ∇(p_exact)(x) +g(x) = -(∇⋅u_exact)(x) + +γ = 10.0 * (1/n) +n_Γu = get_normal_vector(Γu) +a(u,v) = ∫(η*(u⋅v))dΩ + ∫(γ*(u⋅n_Γu)*(v⋅n_Γu))dΓu +b(u,q) = ∫(-(∇⋅u)*q)dΩ +btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu +c(λ,q) = ∫(λ*q)dΩ +l((v,q,μ)) = ∫(v⋅f + g⋅q)dΩ + ∫(γ*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu + c(μ,p_exact) + +τ = 1.0 +sd(u,v) = ∫(τ*(u⋅v))dΩp +s0(p,q) = ∫(τ*(p*q))dΩp + +function weakform(x,y) + u, p, λ = x + v, q, μ = y + Πu, Πv = Πd(u), Πd(v) + Πp, Πq = Π0(p), Π0(q) + divu, divv = ∇⋅u, ∇⋅v + Πdivu, Πdivv = Π0div(u), Π0div(v) + Xp = FESpaces.PatchFESpace(X,ptopo) + data = FESpaces.collect_and_merge_cell_matrix_and_vector( + (X, X , a(u,v) + btilde(v,p) + b(u,q) + c(λ,q) + c(μ,p) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), + (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), + (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), + (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), + ) + assem = SparseMatrixAssembler(X,X) + A, B = assemble_matrix_and_vector(assem,data) + return AffineFEOperator(X,X,A,B) +end + +op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) +uh, ph = solve(op) + +eu = uh-u_exact +ep = ph-p_exact +l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) +l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) + +writevtk( + Ω, "hdiv_cut"; append=false, + cellfields=[ + "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact,"eu"=>eu,"ep"=>ep + ], +) + +aggregates = get_aggregates(cutgeo) +cell_to_agg = flatten_partition(aggregates,num_cells(bgmodel)) +writevtk( + Triangulation(bgmodel), "aggregates"; append=false, + celldata = ["agg"=>cell_to_agg] +) + +function normal_at_centroid(Γ) + n = get_normal_vector(Γ) + x = CellPoint(map(mean,get_cell_coordinates(Γ)),Γ,PhysicalDomain()) + return n(x) +end + +writevtk( + Γu, "boundary_u"; append=false, celldata=["n"=>normal_at_centroid(Γu)] +) diff --git a/test/hdiv_stabilisation.jl b/test/darcy_BGP_mixed.jl similarity index 65% rename from test/hdiv_stabilisation.jl rename to test/darcy_BGP_mixed.jl index a21e229..1965d67 100644 --- a/test/hdiv_stabilisation.jl +++ b/test/darcy_BGP_mixed.jl @@ -5,12 +5,10 @@ using Gridap.Geometry, Gridap.Arrays using Gridap.FESpaces, Gridap.CellData function get_aggregates(cutgeo) - bgmodel = get_background_model(cutgeo) geo = get_geometry(cutgeo) strategy = AggregateAllCutCells() aggregates = aggregate(strategy,cutgeo,geo) - colors = color_aggregates(aggregates,bgmodel) - return inverse_table(colors) + return Table(filter(agg -> length(agg) > 1, inverse_table(aggregates))) end function projection_operator(ptopo, V, W, Ω, dΩ) @@ -32,18 +30,25 @@ function div_projection_operator(ptopo, V, W, Ω, dΩ) return P end -R = 0.4 -geo = disk(R;x0=Point(0.5,0.5)) +function generate_mesh(n,R=0.4,x0=0.835) + geo1 = disk(R;x0=Point(0.5,0.5),name="disk") + geo2 = plane(x0=Point(x0,0.5),v=Point(1.0,0.0),name="plane") + geo = intersect(geo1,geo2,name="domain") + + bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) + cutgeo = cut(bgmodel,geo) + return cutgeo +end n = 12 -order = 1 -qdegree = 2*order+2 +cutgeo = generate_mesh(n) +bgmodel = get_background_model(cutgeo) -bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) -cutgeo = cut(bgmodel,geo) +order = 1 +qdegree = 2*(order+1) Ω_ac = Triangulation(cutgeo,ACTIVE) -V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order)) +V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) X = MultiFieldFESpace([V,Q]) @@ -53,16 +58,21 @@ dΩp = Measure(Ωp, qdegree) Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) -W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order) +W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) Ω = Triangulation(cutgeo,PHYSICAL) Ω_cut = Triangulation(cutgeo,CUT_IN) Γ = EmbeddedBoundary(cutgeo) +Γu = EmbeddedBoundary(cutgeo, "disk","domain") +Γp = EmbeddedBoundary(cutgeo, "plane","domain") + dΩ = Measure(Ω, qdegree) dΩ_cut = Measure(Ω_cut, qdegree) dΓ = Measure(Γ, 2*qdegree) +dΓu = Measure(Γu, 2*qdegree) +dΓp = Measure(Γp, 2*qdegree) η = 1.0 u_exact(x) = VectorValue(x[1], -x[2]) @@ -70,17 +80,17 @@ p_exact(x) = x[1] + x[2] f(x) = η*u_exact(x) + ∇(p_exact)(x) g(x) = -(∇⋅u_exact)(x) -γ = 10.0 -hΓ = CellField(γ ./ get_array(∫(1.0)dΓ),Γ) -n_Γ = get_normal_vector(Γ) -a(u,v) = ∫(η*(u⋅v))dΩ + ∫(hΓ*(u⋅n_Γ)*(v⋅n_Γ))dΓ +γ = 10.0 * (1/n) +n_Γu = get_normal_vector(Γu) +n_Γp = get_normal_vector(Γp) +a(u,v) = ∫(η*(u⋅v))dΩ + ∫(γ*(u⋅n_Γu)*(v⋅n_Γu))dΓu b(u,q) = ∫(-(∇⋅u)*q)dΩ -btilde(u,q) = b(u,q) + ∫((u⋅n_Γ)*q)dΓ -l((v,q)) = ∫(v⋅f + g⋅q)dΩ + ∫(hΓ*(v⋅n_Γ)⋅(u_exact⋅n_Γ))dΓ #- ∫((v⋅n_Γ)*p_exact)dΓ +btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu +l((v,q)) = ∫(v⋅f + g⋅q)dΩ + ∫(γ*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu - ∫((v⋅n_Γp)*p_exact)dΓp τ = 1.0 -sd(u,v) = ∫(τ*(u⋅v))dΩ_cut -s0(p,q) = ∫(τ*p*q)dΩ_cut +sd(u,v) = ∫(τ*(u⋅v))dΩp +s0(p,q) = ∫(τ*p*q)dΩp function weakform(x,y) u, p = x @@ -91,10 +101,10 @@ function weakform(x,y) Πdivu, Πdivv = Π0div(u), Π0div(v) Xp = FESpaces.PatchFESpace(X,ptopo) data = FESpaces.collect_and_merge_cell_matrix_and_vector( - (X, X , a(u,v) + btilde(v,p) + b(u,q) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y), zero(X)), - (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution(), zero(X)), - (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution(), zero(Xp)), - (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution(), zero(Xp)), + (X, X , a(u,v) + btilde(v,p) + b(u,q) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), + (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), + (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), + (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), ) assem = SparseMatrixAssembler(X,X) A, B = assemble_matrix_and_vector(assem,data) @@ -123,6 +133,18 @@ writevtk( celldata = ["agg"=>cell_to_agg] ) +function normal_at_centroid(Γ) + n = get_normal_vector(Γ) + x = CellPoint(map(mean,get_cell_coordinates(Γ)),Γ,PhysicalDomain()) + return n(x) +end + writevtk( Γ, "boundary"; append=false, ) +writevtk( + Γu, "boundary_u"; append=false, celldata=["n"=>normal_at_centroid(Γu)] +) +writevtk( + Γp, "boundary_p"; append=false, celldata=["n"=>normal_at_centroid(Γp)] +) From 76e67e5210bd40480b4228c8804935b1b4d80d7b Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 25 May 2026 22:05:25 +1000 Subject: [PATCH 4/7] Minor --- test/darcy_BGP_mixed.jl | 175 ++++++++++++++++++++++------------------ 1 file changed, 95 insertions(+), 80 deletions(-) diff --git a/test/darcy_BGP_mixed.jl b/test/darcy_BGP_mixed.jl index 1965d67..2ce703a 100644 --- a/test/darcy_BGP_mixed.jl +++ b/test/darcy_BGP_mixed.jl @@ -1,4 +1,5 @@ +using Test using Gridap using GridapEmbedded using Gridap.Geometry, Gridap.Arrays @@ -40,94 +41,109 @@ function generate_mesh(n,R=0.4,x0=0.835) return cutgeo end -n = 12 -cutgeo = generate_mesh(n) -bgmodel = get_background_model(cutgeo) - -order = 1 -qdegree = 2*(order+1) - -Ω_ac = Triangulation(cutgeo,ACTIVE) -V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) -Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) -X = MultiFieldFESpace([V,Q]) - -ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) -Ωp = PatchTriangulation(bgmodel, ptopo) -dΩp = Measure(Ωp, qdegree) - -Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) -Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) -W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) -Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) -Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) - -Ω = Triangulation(cutgeo,PHYSICAL) -Ω_cut = Triangulation(cutgeo,CUT_IN) -Γ = EmbeddedBoundary(cutgeo) -Γu = EmbeddedBoundary(cutgeo, "disk","domain") -Γp = EmbeddedBoundary(cutgeo, "plane","domain") +function driver(n,order,u_exact,p_exact,η=1.0,γ=10.0,τ=10.0) + cutgeo = generate_mesh(n) + bgmodel = get_background_model(cutgeo) + qdegree = 2*(order+1) + + Ω_ac = Triangulation(cutgeo,ACTIVE) + V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) + Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) + X = MultiFieldFESpace([V,Q]) + + ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) + Ωp = PatchTriangulation(bgmodel, ptopo) + + Ω = Triangulation(cutgeo,PHYSICAL) + Γu = EmbeddedBoundary(cutgeo, "disk","domain") + Γp = EmbeddedBoundary(cutgeo, "plane","domain") + + dΩ = Measure(Ω, qdegree) + dΩp = Measure(Ωp, qdegree) + dΓu = Measure(Γu, 2*qdegree) + dΓp = Measure(Γp, 2*qdegree) + + Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) + Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) + W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) + Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) + Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) + + f(x) = η*u_exact(x) + ∇(p_exact)(x) + g(x) = -(∇⋅u_exact)(x) + + h = γ/n + n_Γu = get_normal_vector(Γu) + n_Γp = get_normal_vector(Γp) + a(u,v) = ∫(η*(u⋅v))dΩ + ∫(h*(u⋅n_Γu)*(v⋅n_Γu))dΓu + b(u,q) = ∫(-(∇⋅u)*q)dΩ + btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu + l((v,q)) = ∫(v⋅f + g⋅q)dΩ + ∫(h*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu - ∫((v⋅n_Γp)*p_exact)dΓp + + sd(u,v) = ∫(τ*(u⋅v))dΩp + s0(p,q) = ∫(τ*p*q)dΩp + + function weakform(x,y) + u, p = x + v, q = y + Πu, Πv = Πd(u), Πd(v) + Πp, Πq = Π0(p), Π0(q) + divu, divv = ∇⋅u, ∇⋅v + Πdivu, Πdivv = Π0div(u), Π0div(v) + Xp = FESpaces.PatchFESpace(X,ptopo) + data = FESpaces.collect_and_merge_cell_matrix_and_vector( + (X, X , a(u,v) + btilde(v,p) + b(u,q) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), + (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), + (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), + (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), + ) + assem = SparseMatrixAssembler(X,X) + A, B = assemble_matrix_and_vector(assem,data) + return AffineFEOperator(X,X,A,B) + end + + op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) + uh, ph = solve(op) + + eu = uh-u_exact + ep = ph-p_exact + l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) + l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) + + return cutgeo, uh, ph, l2_err_u, l2_err_p +end -dΩ = Measure(Ω, qdegree) -dΩ_cut = Measure(Ω_cut, qdegree) -dΓ = Measure(Γ, 2*qdegree) -dΓu = Measure(Γu, 2*qdegree) -dΓp = Measure(Γp, 2*qdegree) +function convergence(ns,order,u_exact,p_exact) + l2_u = zeros(Float64,length(ns)) + l2_p = zeros(Float64,length(ns)) + for (i,n) in enumerate(ns) + _, _, _, l2_u[i], l2_p[i] = driver(n,order,u_exact,p_exact) + end + slope_u = log.(l2_u[1:end-1]./l2_u[2:end]) ./ log.(ns[2:end]./ns[1:end-1]) + slope_p = log.(l2_p[1:end-1]./l2_p[2:end]) ./ log.(ns[2:end]./ns[1:end-1]) + return l2_u, l2_p, slope_u, slope_p +end -η = 1.0 +# Manufactured solution u_exact(x) = VectorValue(x[1], -x[2]) p_exact(x) = x[1] + x[2] -f(x) = η*u_exact(x) + ∇(p_exact)(x) -g(x) = -(∇⋅u_exact)(x) - -γ = 10.0 * (1/n) -n_Γu = get_normal_vector(Γu) -n_Γp = get_normal_vector(Γp) -a(u,v) = ∫(η*(u⋅v))dΩ + ∫(γ*(u⋅n_Γu)*(v⋅n_Γu))dΓu -b(u,q) = ∫(-(∇⋅u)*q)dΩ -btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu -l((v,q)) = ∫(v⋅f + g⋅q)dΩ + ∫(γ*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu - ∫((v⋅n_Γp)*p_exact)dΓp - -τ = 1.0 -sd(u,v) = ∫(τ*(u⋅v))dΩp -s0(p,q) = ∫(τ*p*q)dΩp - -function weakform(x,y) - u, p = x - v, q = y - Πu, Πv = Πd(u), Πd(v) - Πp, Πq = Π0(p), Π0(q) - divu, divv = ∇⋅u, ∇⋅v - Πdivu, Πdivv = Π0div(u), Π0div(v) - Xp = FESpaces.PatchFESpace(X,ptopo) - data = FESpaces.collect_and_merge_cell_matrix_and_vector( - (X, X , a(u,v) + btilde(v,p) + b(u,q) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), - (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), - (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), - (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), - ) - assem = SparseMatrixAssembler(X,X) - A, B = assemble_matrix_and_vector(assem,data) - return AffineFEOperator(X,X,A,B) -end - -op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) -uh, ph = solve(op) +cutgeo, uh, ph, l2_err_u, l2_err_p = driver(8,1,u_exact,p_exact) +@test l2_err_u < 1.e-10 +@test l2_err_p < 1.e-10 -eu = uh-u_exact -ep = ph-p_exact -l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) -l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) +u_conv(x) = VectorValue(x[1] + sin(π*x[2]), -x[2] + sin(π*x[1])) +p_conv(x) = sin(π*x[1]) - sin(π*x[2]) +l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) writevtk( - Ω, "hdiv_cut"; append=false, + Ω, "darcy_BGP_mixed"; append=false, cellfields=[ - "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact,"eu"=>eu,"ep"=>ep + "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact, + "eu"=>uh-u_exact,"ep"=>ph-p_exact ], ) -aggregates = get_aggregates(cutgeo) -cell_to_agg = flatten_partition(aggregates,num_cells(bgmodel)) +cell_to_agg = flatten_partition(get_aggregates(cutgeo),num_cells(bgmodel)) writevtk( Triangulation(bgmodel), "aggregates"; append=false, celldata = ["agg"=>cell_to_agg] @@ -139,9 +155,8 @@ function normal_at_centroid(Γ) return n(x) end -writevtk( - Γ, "boundary"; append=false, -) +Γu = EmbeddedBoundary(cutgeo, "disk","domain") +Γp = EmbeddedBoundary(cutgeo, "plane","domain") writevtk( Γu, "boundary_u"; append=false, celldata=["n"=>normal_at_centroid(Γu)] ) From 5353c7d1cac6dd6b4b510997c8e230b20843b0eb Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 25 May 2026 22:14:49 +1000 Subject: [PATCH 5/7] Final drivers --- test/GridapEmbeddedTests/DarcyBGPflux.jl | 169 ++++++++++++++++++ .../DarcyBGPmixed.jl} | 9 + test/darcy_BGP_flux.jl | 141 --------------- 3 files changed, 178 insertions(+), 141 deletions(-) create mode 100644 test/GridapEmbeddedTests/DarcyBGPflux.jl rename test/{darcy_BGP_mixed.jl => GridapEmbeddedTests/DarcyBGPmixed.jl} (93%) delete mode 100644 test/darcy_BGP_flux.jl diff --git a/test/GridapEmbeddedTests/DarcyBGPflux.jl b/test/GridapEmbeddedTests/DarcyBGPflux.jl new file mode 100644 index 0000000..04c2cf9 --- /dev/null +++ b/test/GridapEmbeddedTests/DarcyBGPflux.jl @@ -0,0 +1,169 @@ +module DarcyBGPflux +# Darcy problem: BGP stabilisation with flux boundary conditions (pure pressure trace). +# A Lagrange multiplier fixes the mean pressure to close the system. +# +# Ref: Badia, Boschman, Martín, Nilsson, Ruiz-Baier, Zahedi (2026) +# "Divergence-free unfitted finite element discretisations for the Darcy problem" +# arXiv:2603.26212 + +using Test +using Gridap +using GridapEmbedded +using Gridap.Geometry, Gridap.Arrays +using Gridap.FESpaces, Gridap.CellData + +function get_aggregates(cutgeo) + geo = get_geometry(cutgeo) + strategy = AggregateAllCutCells() + aggregates = aggregate(strategy,cutgeo,geo) + return Table(filter(agg -> length(agg) > 1, inverse_table(aggregates))) +end + +function projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + mass(u,v) = ∫(u⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, mass, mass + ) + return P +end + +function div_projection_operator(ptopo, V, W, Ω, dΩ) + Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) + lhs(u,v) = ∫(u⋅Π(v,Ω))dΩ + rhs(u,v) = ∫((∇⋅u)⋅Π(v,Ω))dΩ + P = LocalOperator( + LocalSolveMap(), ptopo, W, V, lhs, rhs + ) + return P +end + +function generate_mesh(n,R=0.4,x0=0.835) + geo1 = disk(R;x0=Point(0.5,0.5),name="disk") + geo2 = plane(x0=Point(x0,0.5),v=Point(1.0,0.0),name="plane") + geo = intersect(geo1,geo2,name="domain") + + bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) + cutgeo = cut(bgmodel,geo) + return cutgeo +end + +function driver(n,order,u_exact,p_exact,η=1.0,γ=10.0,τ=1.0) + cutgeo = generate_mesh(n) + bgmodel = get_background_model(cutgeo) + qdegree = 2*(order+1) + + Ω_ac = Triangulation(cutgeo,ACTIVE) + V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) + Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) + Λ = ConstantFESpace(Ω_ac) + X = MultiFieldFESpace([V,Q,Λ]) + + ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) + Ωp = PatchTriangulation(bgmodel, ptopo) + + Ω = Triangulation(cutgeo,PHYSICAL) + Γu = EmbeddedBoundary(cutgeo) + + dΩ = Measure(Ω, qdegree) + dΩp = Measure(Ωp, qdegree) + dΓu = Measure(Γu, 2*qdegree) + + Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) + Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) + W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) + Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) + Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) + + f(x) = η*u_exact(x) + ∇(p_exact)(x) + g(x) = -(∇⋅u_exact)(x) + + h = γ/n + n_Γu = get_normal_vector(Γu) + a(u,v) = ∫(η*(u⋅v))dΩ + ∫(h*(u⋅n_Γu)*(v⋅n_Γu))dΓu + b(u,q) = ∫(-(∇⋅u)*q)dΩ + btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu + c(λ,q) = ∫(λ*q)dΩ + l((v,q,μ)) = ∫(v⋅f + g⋅q)dΩ + ∫(h*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu + c(μ,p_exact) + + sd(u,v) = ∫(τ*(u⋅v))dΩp + s0(p,q) = ∫(τ*(p*q))dΩp + + function weakform(x,y) + u, p, λ = x + v, q, μ = y + Πu, Πv = Πd(u), Πd(v) + Πp, Πq = Π0(p), Π0(q) + divu, divv = ∇⋅u, ∇⋅v + Πdivu, Πdivv = Π0div(u), Π0div(v) + Xp = FESpaces.PatchFESpace(X,ptopo) + data = FESpaces.collect_and_merge_cell_matrix_and_vector( + (X, X , a(u,v) + btilde(v,p) + b(u,q) + c(λ,q) + c(μ,p) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), + (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), + (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), + (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), + ) + assem = SparseMatrixAssembler(X,X) + A, B = assemble_matrix_and_vector(assem,data) + return AffineFEOperator(X,X,A,B) + end + + op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) + uh, ph, λh = solve(op) + + eu = uh-u_exact + ep = ph-p_exact + l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) + l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) + + return cutgeo, uh, ph, l2_err_u, l2_err_p +end + +function convergence(ns,order,u_exact,p_exact) + l2_u = zeros(Float64,length(ns)) + l2_p = zeros(Float64,length(ns)) + for (i,n) in enumerate(ns) + _, _, _, l2_u[i], l2_p[i] = driver(n,order,u_exact,p_exact) + end + slope_u = log.(l2_u[1:end-1]./l2_u[2:end]) ./ log.(ns[2:end]./ns[1:end-1]) + slope_p = log.(l2_p[1:end-1]./l2_p[2:end]) ./ log.(ns[2:end]./ns[1:end-1]) + return l2_u, l2_p, slope_u, slope_p +end + +# Manufactured solution +u_exact(x) = VectorValue(x[1], -x[2]) +p_exact(x) = x[1] + x[2] +cutgeo, uh, ph, l2_err_u, l2_err_p = driver(8,1,u_exact,p_exact) +@test l2_err_u < 1.e-10 +@test l2_err_p < 1.e-10 + +u_conv(x) = VectorValue(x[1] + sin(π*x[2]), -x[2] + sin(π*x[1])) +p_conv(x) = sin(π*x[1]) - sin(π*x[2]) +l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) + +writevtk( + Ω, "darcy_BGP_flux"; append=false, + cellfields=[ + "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact, + "eu"=>uh-u_exact,"ep"=>ph-p_exact + ], +) + +cell_to_agg = flatten_partition(get_aggregates(cutgeo),num_cells(bgmodel)) +writevtk( + Triangulation(bgmodel), "aggregates"; append=false, + celldata = ["agg"=>cell_to_agg] +) + +function normal_at_centroid(Γ) + n = get_normal_vector(Γ) + x = CellPoint(map(mean,get_cell_coordinates(Γ)),Γ,PhysicalDomain()) + return n(x) +end + +Γu = EmbeddedBoundary(cutgeo) +writevtk( + Γu, "boundary_u"; append=false, celldata=["n"=>normal_at_centroid(Γu)] +) + +end # module \ No newline at end of file diff --git a/test/darcy_BGP_mixed.jl b/test/GridapEmbeddedTests/DarcyBGPmixed.jl similarity index 93% rename from test/darcy_BGP_mixed.jl rename to test/GridapEmbeddedTests/DarcyBGPmixed.jl index 2ce703a..552a610 100644 --- a/test/darcy_BGP_mixed.jl +++ b/test/GridapEmbeddedTests/DarcyBGPmixed.jl @@ -1,3 +1,10 @@ +module DarcyBGPmixed +# Darcy problem: BGP stabilisation with mixed boundary conditions (flux + pressure traces +# imposed on different parts of the unfitted boundary). +# +# Ref: Badia, Boschman, Martín, Nilsson, Ruiz-Baier, Zahedi (2026) +# "Divergence-free unfitted finite element discretisations for the Darcy problem" +# arXiv:2603.26212 using Test using Gridap @@ -163,3 +170,5 @@ writevtk( writevtk( Γp, "boundary_p"; append=false, celldata=["n"=>normal_at_centroid(Γp)] ) + +end # module \ No newline at end of file diff --git a/test/darcy_BGP_flux.jl b/test/darcy_BGP_flux.jl deleted file mode 100644 index bf27f15..0000000 --- a/test/darcy_BGP_flux.jl +++ /dev/null @@ -1,141 +0,0 @@ - -using Gridap -using GridapEmbedded -using Gridap.Geometry, Gridap.Arrays -using Gridap.FESpaces, Gridap.CellData - -function get_aggregates(cutgeo) - geo = get_geometry(cutgeo) - strategy = AggregateAllCutCells() - aggregates = aggregate(strategy,cutgeo,geo) - return Table(filter(agg -> length(agg) > 1, inverse_table(aggregates))) -end - -function projection_operator(ptopo, V, W, Ω, dΩ) - Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) - mass(u,v) = ∫(u⋅Π(v,Ω))dΩ - P = LocalOperator( - LocalSolveMap(), ptopo, W, V, mass, mass - ) - return P -end - -function div_projection_operator(ptopo, V, W, Ω, dΩ) - Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) - lhs(u,v) = ∫(u⋅Π(v,Ω))dΩ - rhs(u,v) = ∫((∇⋅u)⋅Π(v,Ω))dΩ - P = LocalOperator( - LocalSolveMap(), ptopo, W, V, lhs, rhs - ) - return P -end - -function generate_mesh(n,R=0.4,x0=0.835) - geo1 = disk(R;x0=Point(0.5,0.5),name="disk") - geo2 = plane(x0=Point(x0,0.5),v=Point(1.0,0.0),name="plane") - geo = intersect(geo1,geo2,name="domain") - - bgmodel = CartesianDiscreteModel((0,1,0,1),(n,n)) - cutgeo = cut(bgmodel,geo) - return cutgeo -end - -n = 12 -cutgeo = generate_mesh(n) -bgmodel = get_background_model(cutgeo) - -order = 1 -qdegree = 2*(order+1) - -Ω_ac = Triangulation(cutgeo,ACTIVE) -V = FESpace(Ω_ac,ReferenceFE(raviart_thomas,Float64,order); conformity=:Hdiv) -Q = FESpace(Ω_ac,ReferenceFE(lagrangian,Float64,order); conformity=:L2) -Λ = ConstantFESpace(Ω_ac) -X = MultiFieldFESpace([V,Q,Λ]) - -ptopo = PatchTopology(get_grid_topology(bgmodel), get_aggregates(cutgeo)) -Ωp = PatchTriangulation(bgmodel, ptopo) -dΩp = Measure(Ωp, qdegree) - -Wd = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:RT) -Πd = projection_operator(ptopo, V, Wd, Ωp, dΩp) -W0 = FESpaces.PatchFESpace(bgmodel,Ωp,Float64,order;space=:Q) -Π0 = projection_operator(ptopo, Q, W0, Ωp, dΩp) -Π0div = div_projection_operator(ptopo, V, W0, Ωp, dΩp) - -Ω = Triangulation(cutgeo,PHYSICAL) -Ω_cut = Triangulation(cutgeo,CUT_IN) -Γu = EmbeddedBoundary(cutgeo) - -dΩ = Measure(Ω, qdegree) -dΩ_cut = Measure(Ω_cut, qdegree) -dΓu = Measure(Γu, 2*qdegree) - -η = 1.0 -u_exact(x) = VectorValue(x[1], -x[2]) -p_exact(x) = x[1] + x[2] -f(x) = η*u_exact(x) + ∇(p_exact)(x) -g(x) = -(∇⋅u_exact)(x) - -γ = 10.0 * (1/n) -n_Γu = get_normal_vector(Γu) -a(u,v) = ∫(η*(u⋅v))dΩ + ∫(γ*(u⋅n_Γu)*(v⋅n_Γu))dΓu -b(u,q) = ∫(-(∇⋅u)*q)dΩ -btilde(u,q) = b(u,q) + ∫((u⋅n_Γu)*q)dΓu -c(λ,q) = ∫(λ*q)dΩ -l((v,q,μ)) = ∫(v⋅f + g⋅q)dΩ + ∫(γ*(v⋅n_Γu)⋅(u_exact⋅n_Γu))dΓu + c(μ,p_exact) - -τ = 1.0 -sd(u,v) = ∫(τ*(u⋅v))dΩp -s0(p,q) = ∫(τ*(p*q))dΩp - -function weakform(x,y) - u, p, λ = x - v, q, μ = y - Πu, Πv = Πd(u), Πd(v) - Πp, Πq = Π0(p), Π0(q) - divu, divv = ∇⋅u, ∇⋅v - Πdivu, Πdivv = Π0div(u), Π0div(v) - Xp = FESpaces.PatchFESpace(X,ptopo) - data = FESpaces.collect_and_merge_cell_matrix_and_vector( - (X, X , a(u,v) + btilde(v,p) + b(u,q) + c(λ,q) + c(μ,p) + sd(u,v) - s0(divu,q) - s0(divv,p), l(y)), - (X, Xp , s0(divu,Πq) + s0(Πdivv,p) - sd(u,Πv), DomainContribution()), - (Xp, X , s0(Πdivu,q) + s0(divv,Πp) - sd(Πu,v), DomainContribution()), - (Xp, Xp, sd(Πu,Πv) - s0(Πdivu,Πq) - s0(Πdivv,Πp), DomainContribution()), - ) - assem = SparseMatrixAssembler(X,X) - A, B = assemble_matrix_and_vector(assem,data) - return AffineFEOperator(X,X,A,B) -end - -op = weakform(get_trial_fe_basis(X),get_fe_basis(X)) -uh, ph = solve(op) - -eu = uh-u_exact -ep = ph-p_exact -l2_err_u = sqrt(sum(∫( eu ⋅ eu )dΩ)) -l2_err_p = sqrt(sum(∫( ep ⋅ ep )dΩ)) - -writevtk( - Ω, "hdiv_cut"; append=false, - cellfields=[ - "uh"=>uh,"ph"=>ph,"u_exact"=>u_exact,"p_exact"=>p_exact,"eu"=>eu,"ep"=>ep - ], -) - -aggregates = get_aggregates(cutgeo) -cell_to_agg = flatten_partition(aggregates,num_cells(bgmodel)) -writevtk( - Triangulation(bgmodel), "aggregates"; append=false, - celldata = ["agg"=>cell_to_agg] -) - -function normal_at_centroid(Γ) - n = get_normal_vector(Γ) - x = CellPoint(map(mean,get_cell_coordinates(Γ)),Γ,PhysicalDomain()) - return n(x) -end - -writevtk( - Γu, "boundary_u"; append=false, celldata=["n"=>normal_at_centroid(Γu)] -) From 19735ccb054f69d36069837a95bcf938667bf8ce Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 25 May 2026 22:15:21 +1000 Subject: [PATCH 6/7] Minor --- test/GridapEmbeddedTests/DarcyBGPflux.jl | 2 +- test/GridapEmbeddedTests/DarcyBGPmixed.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/GridapEmbeddedTests/DarcyBGPflux.jl b/test/GridapEmbeddedTests/DarcyBGPflux.jl index 04c2cf9..e271f80 100644 --- a/test/GridapEmbeddedTests/DarcyBGPflux.jl +++ b/test/GridapEmbeddedTests/DarcyBGPflux.jl @@ -139,7 +139,7 @@ cutgeo, uh, ph, l2_err_u, l2_err_p = driver(8,1,u_exact,p_exact) u_conv(x) = VectorValue(x[1] + sin(π*x[2]), -x[2] + sin(π*x[1])) p_conv(x) = sin(π*x[1]) - sin(π*x[2]) -l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) +# l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) writevtk( Ω, "darcy_BGP_flux"; append=false, diff --git a/test/GridapEmbeddedTests/DarcyBGPmixed.jl b/test/GridapEmbeddedTests/DarcyBGPmixed.jl index 552a610..8cc3767 100644 --- a/test/GridapEmbeddedTests/DarcyBGPmixed.jl +++ b/test/GridapEmbeddedTests/DarcyBGPmixed.jl @@ -140,7 +140,7 @@ cutgeo, uh, ph, l2_err_u, l2_err_p = driver(8,1,u_exact,p_exact) u_conv(x) = VectorValue(x[1] + sin(π*x[2]), -x[2] + sin(π*x[1])) p_conv(x) = sin(π*x[1]) - sin(π*x[2]) -l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) +# l2_u, l2_p, su, sp = convergence([8,16,32,64],1,u_conv,p_conv) writevtk( Ω, "darcy_BGP_mixed"; append=false, From dab697bd140c70ddd485684232d4f93e94860a39 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 3 Jun 2026 17:16:43 +1000 Subject: [PATCH 7/7] Actiavate HDiv tests --- NEWS.md | 4 ++++ test/GridapEmbeddedTests/runtests.jl | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/NEWS.md b/NEWS.md index de74695..e1ab35b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Added drivers for Darcy flow using the BGP method. Since PR[#122](https://github.com/gridap/GridapEmbedded.jl/pull/122). + ### Fixed - Fix for issue #115, which is a bug in the logic of the `ACTIVE` mask and `GhostSkeleton` construction. Since PR [#123](https://github.com/gridap/GridapEmbedded.jl/pull/123). diff --git a/test/GridapEmbeddedTests/runtests.jl b/test/GridapEmbeddedTests/runtests.jl index f6d52df..16d9736 100644 --- a/test/GridapEmbeddedTests/runtests.jl +++ b/test/GridapEmbeddedTests/runtests.jl @@ -20,4 +20,9 @@ using Test @time @testset "TraceFEM" begin include("TraceFEMTests.jl") end +@time @testset "DarcyBGP" begin + include("DarcyBGPflux.jl") + include("DarcyBGPmixed.jl") +end + end