using Gridap
using GridapGmsh
using GridapMakie
using CairoMakie
Solving the 2D Poisson equation - seeking the potential two coaxial electrodes with applied potentials. This works as expected...
write("square_hole.geo", """SetFactory("OpenCASCADE");
Rectangle(1) = {.0, 0, 0, 1., 1.};
Rectangle(2) = {.4, 0.4, 0, .2, .2};
bnd1() = Boundary{ Surface{1}; };
bnd2() = Boundary{ Surface{2}; };
BooleanDifference(3) = { Surface{1}; Delete; }{ Surface{2}; Delete; };
Physical Surface("Vacuum", 1) = {3};
Physical Curve("C0", 10) = {bnd1()};
Physical Curve("C1", 20) = {bnd2()};
""")
run(`gmsh -1 -2 -v 0 square_hole.geo`);
model = GmshDiscreteModel("square_hole.msh");
Info : Reading 'square_hole.msh'... Info : 17 entities Info : 112 nodes Info : 224 elements Info : Done reading 'square_hole.msh'
# set boundary conditions
u0(x) = 1.
u1(x) = 0.
f(x) = 0.
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags=["C0", "C1"])
U = TrialFESpace(V0,[u0,u1])
TrialFESpace()
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ
b(v) = ∫( v*f )*dΩ
op = AffineFEOperator(a,b,U,V0)
uh = solve(op)
SingleFieldFEFunction(): num_cells: 184 DomainStyle: ReferenceDomain() Triangulation: BodyFittedTriangulation() Triangulation id: 5697224318339329232
fig = Figure()
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = "y",
title = "grid", aspect=DataAspect())
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = "y",
title = "uh", aspect=DataAspect())
plot!(ax1, Ω)
wireframe!(ax1, Ω, color=:black, linewidth=1)
ax1.aspect = DataAspect()
plt = plot!(ax2, Ω, uh, colorrange=(0, 1.))
Colorbar(fig[1,3], plt)
fig
This appears to be broken. The potential of the small inner electrode is not accounted for.
write("square_hole_small.geo", """SetFactory("OpenCASCADE");
Rectangle(1) = {.0, 0, 0, 1., 1.};
Rectangle(2) = {.4, 0.4, 0, .05, .05};
bnd1() = Boundary{ Surface{1}; };
bnd2() = Boundary{ Surface{2}; };
BooleanDifference(3) = { Surface{1}; Delete; }{ Surface{2}; Delete; };
Physical Surface("Vacuum", 1) = {3};
Physical Curve("C0", 10) = {bnd1()};
Physical Curve("C1", 20) = {bnd2()};
""")
run(`gmsh -1 -2 -v 0 square_hole_small.geo`);
model = GmshDiscreteModel("square_hole_small.msh");
Info : Reading 'square_hole_small.msh'... Info : 17 entities Info : 148 nodes Info : 296 elements Info : Done reading 'square_hole_small.msh'
u0(x) = 1.
u1(x) = 0.
f(x) = 0.
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags=["C0", "C1"])
U = TrialFESpace(V0,[u0,u1])
TrialFESpace()
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ
b(v) = ∫( v*f )*dΩ
op = AffineFEOperator(a,b,U,V0)
uh = solve(op)
SingleFieldFEFunction(): num_cells: 260 DomainStyle: ReferenceDomain() Triangulation: BodyFittedTriangulation() Triangulation id: 8393209368221062039
fig = Figure()
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = "y",
title = "grid", aspect=DataAspect())
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = "y",
title = "uh", aspect=DataAspect())
plot!(ax1, Ω)
wireframe!(ax1, Ω, color=:black, linewidth=1)
ax1.aspect = DataAspect()
plt = plot!(ax2, Ω, uh, colorrange=(0, 1.))
Colorbar(fig[1,3], plt)
fig