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
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 = 1

mesh = from_meshio(
    mesh_from_OrderedDict(polygons_paper, {}, default_resolution_max=0.08, filename="mesh.msh")
)
mesh.draw().show()
../../_images/71ce420133c5ab3326399eb82a5b1856446b873872d6b553cdcb30aad2fc4d3a.png
errors = []
nelements = []

for i in range(20):
    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, order=1, metallic_boundaries=boundaries[num]
    )
    modes[0].plot(modes[0].E.real, direction="x")
    plt.show()
    error = modes[0].n_eff.real - neff_values_paper[num]
    errors.append(error)
    nelements.append(mesh.nelements)
    elements_to_refine = adaptive_theta(modes[0].eval_error_estimator(), theta=0.5)
    print(
        f"Error {i:2d}: {error: .7f},\tElements: {mesh.nelements:10d}\tSelected: {len(elements_to_refine)}"
    )

    fig, ax = plt.subplots(1, 2)
    for a in ax:
        a.set_aspect("equal")
    mesh.draw(ax=ax[0])
    mesh_restricted = mesh.restrict(elements=elements_to_refine)
    mesh_restricted.draw(color="red", linewidth=2, ax=ax[0])
    mesh = mesh.refined(elements_to_refine)
    mesh.draw(ax=ax[1])
    mesh_restricted.draw(color="red", linewidth=2, ax=ax[1]).show()

fig, axs = plt.subplots(2, 1)
axs[0].plot(errors, label="Error")
axs[0].axhline(y=0, color="r", linestyle="-")
axs[0].set_ylabel("Error")
axs[1].plot(nelements)
axs[1].semilogy()
axs[1].set_xlabel("Iteration")
axs[1].set_ylabel("Number of elements")
plt.show()
../../_images/92045d881bcc847342e1055ba8765079305a51d0ca3801235380d75d5fd20fc6.png
Named boundaries invalidated by a call to Mesh.refined()
Error  0: -0.0009402,	Elements:        484	Selected: 8
../../_images/a49ec3798e4f97c76de0fb498a35f53f4ec80dfb695396c4867dd6c76b8cbdf1.png ../../_images/13f0d7a2f2c9505e6e751912d67115f7b459c748cba6713e0e45ef61757cfe3c.png
Error  1: -0.0009728,	Elements:        534	Selected: 26
../../_images/18160fe7f9f9b65b08cfdbb27264417da8d498e74060b97937944001f939b6e7.png ../../_images/911a93c4b56cc7e661800ad1d2b8bae7c12253e55dc60b92d41a1a6b593ab24f.png
Error  2: -0.0009282,	Elements:        697	Selected: 11
../../_images/2e1b30e4d058f386aaeba46798f7d46e214a493487b2f71038197bd6f84bfade.png ../../_images/57359dbbe8c740896b64fe9ac3650440c3345266c23b7aef0866ce615ec2ecb5.png
Error  3: -0.0008647,	Elements:        778	Selected: 72
../../_images/58de2995750128d8a7f31c29aca36a343a38a8b7cc8fbc683397e348b45169fa.png ../../_images/7228083db16a276ba29cf9abc5aa101229a64d0b26a3d93e5d69c6daae8fb49d.png
Error  4: -0.0005305,	Elements:       1174	Selected: 44
../../_images/a2ebf1381ecdc18a8c86913b1ca2cefd13e1edb58ff5eef88cf2f7bc176b854e.png ../../_images/4d844b82b58141a8b1e142132ce2424930eed62ae8e3d53d69eea74fda140144.png
Error  5: -0.0004795,	Elements:       1442	Selected: 120
../../_images/ef6b4fcd4699753ef73a66bdc3267e2881c2f9856feb59881dce75ed8dc49d18.png ../../_images/58028ad6d958b7db9f8498c8f614f2b32a7166667ce4d19999d4fa2c598b35f9.png
Error  6: -0.0003474,	Elements:       2147	Selected: 100
../../_images/f0d5c7d356175a4df7b7b75d16c50fccefb9e66c7f2fde7f76c42a3e0285fbd6.png ../../_images/86419fab264c1862e97f0726f71d1810cf78b07a72378cc696dad43b29d3ded5.png
Error  7: -0.0002198,	Elements:       2732	Selected: 248
../../_images/ab656afd34f36a9476161cd4bc8745896d2029b54e898661340ef87cfda99906.png ../../_images/0a83c705d1693be5e8dbd6c8956494c54005cbf65ebb46688d2017b1442e3cf4.png
Error  8: -0.0001822,	Elements:       4046	Selected: 247
../../_images/ae762b809dd7dde9f6ae67605b0e7e4d304eb0af3889464df878260dfd8493db.png ../../_images/1399ae04988d3ecb5288358131930157fe991a508618ccf989bef48a32398998.png
Error  9: -0.0001337,	Elements:       5551	Selected: 482
../../_images/ad846891bfaf46d447893c9a0fd2d5b0f360e4ad9fd7d621f2a9da190822ef85.png ../../_images/f19cca5b623bdb0c40720aeeda0b2bc6a3e422e1bc57a4aa87cb656f2efbc302.png
Error 10: -0.0000891,	Elements:       8213	Selected: 488
../../_images/8d4c5dcfb2c33126384c07c6af1c85df26d9d99f357cb9e8834b331c3ebbacfd.png ../../_images/c67e991995ca90bd46a513a0a44c20aab8b579e85e1ef7ca7a058ed82085bc09.png
Error 11: -0.0000643,	Elements:      11156	Selected: 994
../../_images/fc90ee43d50ce989c01a531c18e07ef0da123c4365807d6f8bf7f52b08af175f.png ../../_images/05240e6a811e78dd7553790936fd6b11cc16c90dc0672568475fe507232e3946.png
Error 12: -0.0000486,	Elements:      16390	Selected: 853
../../_images/de579bb21ffbbb0a251c8555c26850926bf12387119c157e79bbd6593d03478a.png ../../_images/d1110d4b79d1f3777044344b6e1bf4e40a854714e427ebf7166cf404c280ea19.png
Error 13: -0.0000372,	Elements:      21480	Selected: 1874
../../_images/12cfc773e92720d6bb49917f867fbfd7d60a64564b0a8ed3d9cc528956354096.png ../../_images/24e5e3d16bd8dca7289ba5252e9417f0891e9dbd85c988ce496824091fee612e.png
Error 14: -0.0000250,	Elements:      31266	Selected: 1417
../../_images/2cc9746b920a0bed6c28b86f271f0bd74320b238e8c3459c707b452e300f89e4.png ../../_images/c7a43fbc10c7bd2996bbd016b16e4941d8030e66a2baefb34c52e0f1b619a64f.png
Error 15: -0.0000200,	Elements:      40060	Selected: 2756
../../_images/625727730a32df27a93ccf480d30fa97bd11a57d01bafaa1bfde518c65dcddb5.png ../../_images/f15773e67d1bc7b5e3953e434724166ae4ef56115373dcb9f9e296f39263d806.png
Error 16: -0.0000140,	Elements:      54708	Selected: 4203
../../_images/141cc4231495fc4d23102e77edbf75fa5e7e0e772a524a143abbdff671878edc.png ../../_images/887a0c775d849a107cd1c6ce73aa28106ef21aab9db936f0d78f00a60377be95.png
Error 17: -0.0000108,	Elements:      78269	Selected: 3754
../../_images/289ba143b831a3e06b7df6b65472fff154bed0a52013d665dcf4c3c999b7bacb.png ../../_images/4b7935d9f6a939f2b1365f1ab2752dc38d2f34665cb3c9eab81225596937fd09.png
Error 18: -0.0000079,	Elements:      99379	Selected: 6361
../../_images/53cfa88af4d2ddb2b45024fe281dc60a1b2e5e96ea941621e3cdd47281273bec.png ../../_images/069ae46c397dee5a595af9fe5b8950d9ee68c51df4307d1e20be0c68aca6347d.png
Error 19: -0.0000062,	Elements:     133490	Selected: 7457
../../_images/af09146804617a14090ef0b29a8d04b99a0d1c174db4e53bb317ff7a69450c69.png ../../_images/6eee8acb8399eff9e14c5420cf402c8f1e471e926a9f9894b0d7b7a821c372c9.png