Add basic mesh for consistency checks

docs/julia/fresnel.jl

@@ -0,0 +1,115 @@
+using Gridap
+using Gridap.Geometry
+using GridapGmsh
+using GridapMakie, GLMakie
+using PyCall
+
+import MUMPS
+using MPI, SparseArrays
+
+MPI.Init()
+
+GLMakie.inline!(true)
+
+solver = LinearFESolver(BackslashSolver())
+solver = LinearFESolver(LUSolver())
+
+model = GmshDiscreteModel("mesh.msh")
+Ω = Triangulation(model)
+dΩ = Measure(Ω, 4)
+Γ_in = BoundaryTriangulation(model, tags = "PML_top___volume1")
+dΓ_in = Measure(Γ_in, 4)
+Γ_out = BoundaryTriangulation(model, tags = "PML_bottom___volume1")
+dΓ_out = Measure(Γ_out, 4)
+#fig = plot(Ω)
+#fig.axis.aspect = DataAspect()
+#wireframe!(Ω, color = :black, linewidth = 1)
+#display(fig)
+
+mu_r = 1
+k0 = 1e4 # 1 THz
+
+labels = get_face_labeling(model)
+
+n1 = 1
+n2 = 1.2
+epsilons = ["volume1" => n1^2, "PML_top" => n2^2, "PML_bottom" => n2^2]
+ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]
+τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
+
+ε_total = (ε ∘ τ) + (x -> 10im * max(0, x[3] - 1)^2) + (x -> 10im * min(0, x[3] + 1)^2)
+
+function lhs_maxwell(u, v)
+    0.5 * ∫((1 / ε_total * curl(u) ⋅ curl(v)) - k0^2 * u ⋅ v)dΩ
+end
+
+input_form_rhs(v) = ∫(-2 * ((x -> VectorValue(1, 0, 0)) ⋅ v))dΓ_in
+
+V = TestFESpace(
+    model,
+    ReferenceFE(nedelec, Float64, 1),
+    conformity = :HCurl,
+    vector_type = Vector{ComplexF64},
+    #dirichlet_tags = ["cladding___boundary", "substrate___boundary"],
+    #dirichlet_masks=[(true,false,false),(true,false,false)]
+)
+#U = TrialFESpace(V, x -> VectorValue(0, 0, 0))
+U = TrialFESpace(V)
+println(num_free_dofs(V))
+
+#op = AffineFEOperator(lhs_maxwell, input_form_rhs, U, V)
+
+#sol = solve(solver, op)
+
+assem = Gridap.FESpaces.SparseMatrixAssembler(U, V)
+A = assemble_matrix(lhs_maxwell, assem, U, V)
+B = assemble_vector(input_form_rhs, assem, U)
+x = MUMPS.solve(A, B)[:, 1]
+sol = FEFunction(U, x)
+
+Base.real(x::VectorValue{3,ComplexF64}) = VectorValue(real.(x.data))
+Base.imag(x::VectorValue{3,ComplexF64}) = VectorValue(imag.(x.data))
+writevtk(
+    Ω,
+    "fresnel_results",
+    cellfields = [
+        "epsilon_real" => ε ∘ τ,
+        "epsilon_imag" => imag(ε_total),
+        "field_real" => real(sol),
+        "field_imag" => imag(sol),
+    ],
+)
+
+#fig = plot(Ω, real(sol))
+#display(fig)
+#fig = plot(Ω, abs(sol ⋅ sol))
+#display(fig)
+
+#sum_in = abs(sum(∫(norm ∘ sol)dΓ_in))
+#sum_out = abs(sum(∫(norm ∘ sol)dΓ_out))
+#println(sum_in)
+#println(sum_out)
+
+
+#display(
+#    plot([imag(sol(Gridap.Point(0, 0, x)) ⋅ VectorValue(1, 0, 0)) for x = -0.25:0.01:0.25]),
+#)
+#
+#println(
+#    (
+#        minimum([abs((norm ∘ sol)(Gridap.Point(0.5, 0.5, x))) for x = -1:0.01:1]) + maximum([abs((norm ∘ sol)(Gridap.Point(0.5, 0.5, x))) for x = -1:0.01:1])
+#    ) / 2,
+#)
+#println(minimum([abs((norm ∘ sol)(Gridap.Point(0.5, 0.5, x))) for x = -1:0.01:1]))
+#println(maximum([abs((norm ∘ sol)(Gridap.Point(0.5, 0.5, x))) for x = -1:0.01:1]))
+#println(maximum([abs((norm ∘ sol)(Gridap.Point(0.5, 0.5, x))) for x = -1:0.01:1]))
+#
+#transmission =
+#    maximum([imag(sol(Gridap.Point(0.5, 0.5, x)) ⋅ VectorValue(1, 0, 0)) for x = -1:0.01:0])
+#in_plus_reflection =
+#    maximum([imag(sol(Gridap.Point(0.5, 0.5, x)) ⋅ VectorValue(1, 0, 0)) for x = 0:0.01:1])
+#
+#println(transmission / in_plus_reflection)
+#
+#reflection = abs((n1 - n2) / (n1 + n2))
+#println((1 - reflection) / (1 + reflection))

docs/julia/fresnel.py

@@ -0,0 +1,130 @@
+from meshwell.model import Model
+from meshwell.prism import Prism
+from shapely import Polygon, box
+
+if __name__ == "__main__":
+    model = Model()
+
+"""
+DEFINE ENTITIES
+"""
+
+simvolx = 1
+simvoly = 1
+simvolz = 0.5
+pml_z = 0.25
+boundary_thickness = 0.25
+
+vol1_polygon = box(xmin=-simvolx / 2, ymin=-simvoly / 2, xmax=simvolx / 2, ymax=simvoly / 2)
+
+vol1 = Prism(
+    polygons=vol1_polygon,
+    buffers={
+        -simvolz / 2: 0.0,
+        simvolz / 2: 0.0,
+    },
+    model=model,
+    physical_name="volume1",
+    mesh_order=3,
+    resolution={"resolution": 0.1, "SizeMax": 1.0, "DistMax": 1.0},
+)
+
+PML_bot = Prism(
+    polygons=vol1_polygon,
+    buffers={
+        -simvolz / 2 - pml_z: 0.0,
+        -simvolz / 2: 0.0,
+    },
+    model=model,
+    physical_name="PML_bottom",
+    mesh_order=2,
+)
+
+PML_top = Prism(
+    polygons=vol1_polygon,
+    buffers={
+        simvolz / 2: 0.0,
+        simvolz / 2 + pml_z: 0.0,
+    },
+    model=model,
+    physical_name="PML_top",
+    mesh_order=1,
+)
+
+"""
+BOUNDARIES
+"""
+
+right_polygon = box(
+    xmin=simvolx / 2, ymin=-simvoly / 2, xmax=simvolx / 2 + boundary_thickness, ymax=simvoly / 2
+)
+left_polygon = box(
+    xmin=-simvolx / 2 - boundary_thickness,
+    ymin=-simvoly / 2,
+    xmax=-simvolx / 2,
+    ymax=simvoly / 2,
+)
+front_polygon = box(
+    xmin=-simvolx / 2, ymin=simvoly / 2, xmax=simvolx / 2, ymax=simvoly / 2 + boundary_thickness
+)
+back_polygon = box(
+    xmin=-simvolx / 2,
+    ymin=-simvoly / 2 - boundary_thickness,
+    xmax=simvolx / 2,
+    ymax=-simvoly / 2,
+)
+
+boundary_buffers = {-simvolz / 2 - pml_z: 0.0, simvolz / 2 + pml_z: 0.0}
+
+right = Prism(
+    polygons=right_polygon,
+    buffers=boundary_buffers,
+    model=model,
+    physical_name="right",
+    mesh_bool=False,
+    mesh_order=0,
+)
+
+left = Prism(
+    polygons=left_polygon,
+    buffers=boundary_buffers,
+    model=model,
+    physical_name="left",
+    mesh_bool=False,
+    mesh_order=0,
+)
+
+up = Prism(
+    polygons=front_polygon,
+    buffers=boundary_buffers,
+    model=model,
+    physical_name="up",
+    mesh_bool=False,
+    mesh_order=0,
+)
+
+down = Prism(
+    polygons=back_polygon,
+    buffers=boundary_buffers,
+    model=model,
+    physical_name="down",
+    mesh_bool=False,
+    mesh_order=0,
+)
+
+"""
+ASSEMBLE AND NAME ENTITIES
+"""
+entities = [vol1, PML_bot, PML_top, up, down, left, right]
+
+mesh = model.mesh(
+    entities_list=entities,
+    verbosity=0,
+    global_scaling=60e-6,
+    filename="mesh.msh",
+    periodic_entities=[
+        (x + "___" + s1, x + "___" + s2)
+        for x in ("volume1", "PML_bot", "PML_top")
+        for (s1, s2) in (("left", "right"), ("up", "down"))
+    ],
+)