Mesh refinement

Mesh refinement#

Hide code cell source
from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import shapely
import shapely.affinity
from skfem import Basis, ElementTriP0, adaptive_theta
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes, eval_error_estimator
from femwell.mesh import mesh_from_OrderedDict
polygons_paper = OrderedDict(
    left=shapely.LineString(((0, 0), (0, 1))),
    right=shapely.LineString(((1, 0), (1, 1))),
    top=shapely.LineString(((0, 1), (1, 1))),
    core=shapely.box(0, 0, 0.5, 0.5),
    clad=shapely.box(0, 0, 1, 1),
)

boundaries = [
    [lambda x: x[0] == np.min(x[0]), lambda x: x[0] == np.max(x[0])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[0] == np.max(x[0])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[1] == np.max(x[1])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[1] == np.max(x[1])],
]

epsilons_paper = [
    {"core": 2.25, "clad": 1},
    {"core": 8, "clad": 1},
    {"core": 1, "clad": 2.25},
    {"core": 1, "clad": 8},
]

neff_values_paper = [1.27627404, 2.65679692, 1.387926425, 2.761465320]
num = 3

mesh = from_meshio(
    mesh_from_OrderedDict(polygons_paper, {}, default_resolution_max=0.1, filename="mesh.msh")
)
mesh.draw().show()
../../_images/827f13aeb4339aedfbc51034ee8160b371a35ace9a82a66675f76512589a2b9a.png
for i in range(16):
    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros()
    for subdomain, e in epsilons_paper[num].items():
        epsilon[basis0.get_dofs(elements=subdomain)] = e

    modes = compute_modes(
        basis0, epsilon, wavelength=1.5, num_modes=1, order=1, metallic_boundaries=boundaries[num]
    )
    modes[0].show(modes[0].E.real, direction="x")
    print(f"Error: {modes[0].n_eff - neff_values_paper[num]:.5g}")

    mesh = mesh.refined(adaptive_theta(eval_error_estimator(modes[0].basis, modes[0].E), theta=0.5))
    mesh.draw().show()
/tmp/ipykernel_9218/275219315.py:10: DeprecationWarning: The behavior of passing an array directly to `show` is deprecated and will be removed in the future. Use `plot` instead.
  modes[0].show(modes[0].E.real, direction="x")
../../_images/0b63eade47d9abbb4c862819892404071613e5102a0b218417596836d571e6b7.png
Named boundaries invalidated by a call to Mesh.refined()
Error: 4.2371e-05+0j
../../_images/d504a3b67e6af4b9d9585cc436c61130fd1e59a06726ee974c2f825ef1f6b6c4.png ../../_images/3d8c1bbd297b4e63635a164e0dca2792f5aa306319f54d48f9d50da5c9a4551c.png
Error: -0.00011698+0j
../../_images/78fb8704f729d61a19832da532b0edbdafcb09fbe1fa69f78818a07e21909f16.png ../../_images/7c242c1a0c2745172640c08342ca6999eb6d052c41869d740fb5384c5de71be3.png
Error: -0.00014164+0j
../../_images/154c6a86d9e325969f5dae6238e037cc9d0fa38581b0085320faf3e029b9cad8.png ../../_images/94e505f13fd0a94c21dbf0dddcbf40ac72f521c1fae86ce0b0b6b4abc53b76ec.png
Error: -0.00014494+0j
../../_images/10d67735df38d6b617f004a119493769943ed80397a656b0ae75415a55931fa4.png ../../_images/823173421fe59f912ab188923f84c47fafd187ed12d9c7d697986146a8afbba5.png
Error: -0.00014891+0j
../../_images/e48877dcbdfd16feb30377d81efe116b37efeb93c5d4b72ac9d098b82403023f.png ../../_images/2b23d5a1a32980a79a939e573b784b8a4310c9bd6bcccbb3890cff9e4cb471f5.png
Error: -0.00013566+0j
../../_images/ade3ec2bf4c5aa53067fcff59355c921191774dd697a4cd2dc7c06ee56f72804.png ../../_images/7813352d2b88d2e7e5d44e81622b4d706a344e41e34ffd2e4fe4d61a0eefbecb.png
Error: -0.00013255+0j
../../_images/7a255a08e61095eb002d57f8166cd6dace42961607078a3fe66d312cc5879e15.png ../../_images/8072b327bd7bf2dcf64c61e0187ef4be526176361bbf2ee258704369196895fd.png
Error: -9.8289e-05+0j
../../_images/ee645845ab654571f3951e51da136b975cbe13b65071109a0e1e7cc6d477a781.png ../../_images/f34952b8c3615779d4ef058ea4f93afecbe417d9ddec8fe9fc41bc3eae72e32e.png
Error: -5.0584e-05+0j
../../_images/12ecd09733349857c8c0d3134f868a22a1391a578e3c4a12c17654f74e4ca2c9.png ../../_images/403975708d3f5c128fc64f31e6fc0ceae633e0b516f12871a4d3c402f9ed0a3f.png
Error: -4.1234e-05+0j
../../_images/e4a0723a2701d005c145761afa0cc23c856efcb1355ffbfb5fe3ef66389c639f.png ../../_images/e90af069cc826aefa3b2c44ae545e43fff29f3158391b8a67e7474435b507a82.png
Error: -3.4453e-05+0j
../../_images/b0f27245dc1e0245998a5e4cd71bd369d207fd3987a1f2665c25142b068393cc.png ../../_images/1c324b2f30c9cec6e2ceb03f0e13cbaf3feaae2234cd27251ec626513b025bb2.png
Error: -2.3801e-05+0j
../../_images/737af4962d6347b4d6d52ac81c86778f9d155e7a57f11ad1017d50ec0cc0986c.png ../../_images/a5cf79851ccbfe5bddc27d85afd92f2bf7e5dcda71a75b48ba0a861abb58c2e6.png
Error: -1.5849e-05+0j
../../_images/928dc9e7f07beccfb846d183ec24ef135a1242702206516019d6c7efe01f9dfb.png ../../_images/134b843b7c48feab02261fcbb1e2cd23bda5b7c4458b57dcd52950ca11fb5576.png
Error: -1.139e-05+0j
../../_images/b5deb5ec1cd1167981f49b688f06c07d471845bd5acdf02a4d13cafc1cec616a.png ../../_images/8d9492f2e8de8ca26d7d5926ff3683e6f52298376959237aca3b3667388cff2d.png
Error: -9.823e-06+0j
../../_images/4f5ce1b98ada6bf40e06684775feea8ba4c949fe9ca5a8fb1d8366c8cb7c59be.png ../../_images/0cd1e23294e400e437630a2d80864de4b5a3eee1245bc9a178ce4d1ed51a6742.png
Error: -7.9428e-06+0j
../../_images/9d0bb86f75698b112be9b3195e44ad6024466db35536b84c02b4d44c7fbcd030.png