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/ccd2887443cefb54124f7436a2ec4a6c90d7fe9f2d42ce4a64e3d9642fcd3411.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_9174/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/27694f92b73794d94eb615ac301026774de12769cf73bf7f8c3e51c634552ae7.png
Named boundaries invalidated by a call to Mesh.refined()
Error: 4.2371e-05+0j
../../_images/b2f4791ba5267d2140ac8ce20b0715d4fa27c167fb5e84d5bfa025053acef4a6.png ../../_images/e9e75edc5756fb3d1eda8d8fb9aa17a66470cda1b4c30fa756a4d98a4cd70dd6.png
Error: -0.00011698+0j
../../_images/8b6ade8d8e55591f7a08cdb7af3d38964f27f18e72a1e1ad4bb1188f99faba98.png ../../_images/9855a60fae781a9b563998736402c1e0eb95cd751d863c5ee89752b18920f10c.png
Error: -0.00014164+0j
../../_images/800d32c100d674a6de7b7fed860da9c3e5e8c56fe939d53f1b2727cd2da367b8.png ../../_images/9a98aefa58cf640b595197775cfc9675db3a234fc55dc20746739b1df61ba98f.png
Error: -0.00014494+0j
../../_images/f53692f175fc5a00881a5537fc9a3356243485cf5501fbc0f1f52dd0244e6703.png ../../_images/067ac52cb1799fbd812a33e0d061a7f2a8688643aad16d3fbd65bb2f79ec4ecb.png
Error: -0.00014891+0j
../../_images/46a6f219340794d6debc8d6243cff365a3092f1ba62cabebb56be73546e4f8b3.png ../../_images/29cebdd872f6826b503e595d747d56fbebf20bd0a3fe0006f320938b353d1470.png
Error: -0.00013566+0j
../../_images/4660883b1fa6a4b9aae22a925f7a52a2aa5ddd065edd8210f3f2721370cd8398.png ../../_images/139254eba07ea162ccc7146272c3cb8c5cff5a04328b58798d2cb0ba9541c007.png
Error: -0.00013255+0j
../../_images/153c20d66e99ccbcc635da1184f9d705b6470bac8e587bfbce8eb333f869b068.png ../../_images/120b45ae84627197f7d19a2ba1cd8d0c1d7ba345b5f22dc8995433a2091a5d3b.png
Error: -9.8289e-05+0j
../../_images/2b2ea96b19036376e7d3957f54b1a41607271f84afb877e6c2b7c06ceeab1dad.png ../../_images/16a5a4fbe9c43a242c146e2dedb0c10f4faa05781cf0f3c07a3e4f2ead8ef916.png
Error: -5.0584e-05+0j
../../_images/e0775c22a73786bee9c8497d07b0c8f37ebacec1e28f2f649c044d518e32c872.png ../../_images/19f217bdf6d8e7e980dcec0701f0f682e032762adeafe811477fa568ca178991.png
Error: -4.1234e-05+0j
../../_images/af2c09f26844c1e6ffbb09cb1ac747d450ea78cfea3de836b9e432a2c031980e.png ../../_images/2d4bea551e349bd4ecfd1de850df71c850a8a9804a439d49450d85df4cf85470.png
Error: -3.4453e-05+0j
../../_images/ae273054ab7ffe620a021a6ef5903996268138f7807f527155fe1147616142e5.png ../../_images/a193ebfb36a4445f03efab8a02694d73b02e60dc777a4dbace59beecd43c8637.png
Error: -2.3801e-05+0j
../../_images/f6aff68d1fbb0a12ba26cb01be00998d37a231b598668c09ebed65f6df3e7f64.png ../../_images/72e19c81bcb66e66c26b4174dcea2f6861033e3bc5eeb14061f2c90d53567a35.png
Error: -1.5849e-05+0j
../../_images/d97f1ebc1eee211d653badbd5793cc0e9ba25d2640722836b11da6ad5ded7bd9.png ../../_images/2d07701efb253ba5aa7a168263918906d65ba5b6ab520af2feb7c9ed18bbe8dc.png
Error: -1.139e-05+0j
../../_images/ba7c5afdf06f2bd1d5a8555ed6452402cc3818d705b1b90d2f03cbad5c03537b.png ../../_images/ea0918e145f6182029450fa4853101ceaf9c3a54220eea7f33e05742d941403e.png
Error: -9.823e-06+0j
../../_images/2251fa71f03ff43729ae39e48093b682dfa42324190fe5537d9f624a587765cd.png ../../_images/bfeb7262afaa65ae649cfef6572c572438b357a0c43c7d87ae80c02851938cd3.png
Error: -7.9428e-06+0j
../../_images/549bda0775a1eca8fc496ad80d3f5d1b0f9d6f37807d6cb09cf32186c3d04614.png