using Gridap
using GridapGmsh
using GridapMakie
using CairoMakie

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

# 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])

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)

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

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

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

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)

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