cleanup

Author: HelgeGehring

docs/photonics/examples/refinement.py

@@ -27,7 +27,7 @@
 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.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
 
 # %%
@@ -71,17 +71,16 @@
         epsilon[basis0.get_dofs(elements=subdomain)] = e
 
     modes = compute_modes(
-        basis0, epsilon, wavelength=1.5, num_modes=1, order=1, metallic_boundaries=boundaries[num]
+        basis0, epsilon, wavelength=1.5, order=1, metallic_boundaries=boundaries[num]
     )
-    modes[0].show(modes[0].E.real, direction="x")
+    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)
     print(f"Error {i:2d}: {error: .7f}, Elements: {mesh.nelements:10d}")
 
-    elements_to_refine = adaptive_theta(
-        eval_error_estimator(modes[0].basis, modes[0].E, epsilon), theta=0.5
-    )
+    elements_to_refine = adaptive_theta(modes[0].eval_error_estimator(), theta=0.5)
 
     fig, ax = plt.subplots()
     mesh.restrict(elements=elements_to_refine).draw(color="red", linewidth=2, ax=ax)

femwell/maxwell/waveguide.py

@@ -463,6 +463,37 @@ def plot_intensity(
 
         return fig, ax
 
+    def eval_error_estimator(self):
+        # facet jump
+        fbasis = [InteriorFacetBasis(self.basis.mesh, self.basis.elem, side=i) for i in [0, 1]]
+        fbasis_epsilon = [
+            InteriorFacetBasis(
+                self.basis.mesh, self.basis_epsilon_r.elem, side=i, quadrature=fbasis[0].quadrature
+            )
+            for i in [0, 1]
+        ]
+        w = {f"u{str(i + 1)}": fbasis[i].interpolate(self.E) for i in [0, 1]}
+        w2 = {f"epsilon{str(i + 1)}": fbasis_epsilon[i].interpolate(self.epsilon_r) for i in [0, 1]}
+
+        # norm_0 = np.linalg.norm(w["u1"][0]*w2["epsilon1"])
+        # norm_1 = np.linalg.norm(grad(w["u1"][1]))*10
+        norm_0 = norm_1 = 1
+
+        @Functional
+        def edge_jump(w):
+            return w.h * (
+                np.abs(dot(grad(w["u1"][1]) - grad(w["u2"][1]), w.n)) ** 2 / norm_1**2
+                + np.abs(dot(w["u1"][0] * w["epsilon1"][0] - w["u2"][0] * w["epsilon2"][0], w.n))
+                ** 2
+                / norm_0**2
+            )
+
+        tmp = np.zeros(self.basis.mesh.facets.shape[1])
+        tmp[fbasis[0].find] = edge_jump.elemental(fbasis[0], **w, **w2)
+        eta_E = np.sum(0.5 * tmp[self.basis.mesh.t2f], axis=0)
+
+        return eta_E
+
 
 @dataclass(frozen=True)
 class Modes:
@@ -800,44 +831,6 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi
     return fig, axs
 
 
-def eval_error_estimator(basis, u, epsilon):
-    # @Functional
-    # def interior_residual(w):
-    #     h = w.h
-    #     x, y = w.x
-    #     return h**2  # * load_func(x, y) ** 2
-
-    # eta_K = interior_residual.elemental(basis, w=basis.interpolate(u))
-
-    # facet jump
-    fbasis = [InteriorFacetBasis(basis.mesh, basis.elem, side=i) for i in [0, 1]]
-    fbasis_epsilon = [
-        InteriorFacetBasis(basis.mesh, ElementTriP0(), side=i, quadrature=fbasis[0].quadrature)
-        for i in [0, 1]
-    ]
-    w = {f"u{str(i + 1)}": fbasis[i].interpolate(u) for i in [0, 1]}
-    w2 = {f"epsilon{str(i + 1)}": fbasis_epsilon[i].interpolate(epsilon) for i in [0, 1]}
-
-    # norm_0 = np.linalg.norm(w["u1"][0]*w2["epsilon1"])
-    # norm_1 = np.linalg.norm(grad(w["u1"][1]))*10
-    norm_0 = norm_1 = 1
-
-    @Functional
-    def edge_jump(w):
-        return w.h * (
-            np.abs(dot(grad(w["u1"][1]) - grad(w["u2"][1]), w.n)) ** 2 / norm_1**2
-            + np.abs(dot(w["u1"][0] * w["epsilon1"][0] - w["u2"][0] * w["epsilon2"][0], w.n)) ** 2
-            / norm_0**2
-        )
-
-    tmp = np.zeros(basis.mesh.facets.shape[1])
-    tmp[fbasis[0].find] = edge_jump.elemental(fbasis[0], **w, **w2)
-    eta_E = np.sum(0.5 * tmp[basis.mesh.t2f], axis=0)
-
-    return eta_E
-    # return eta_K + eta_E
-
-
 if __name__ == "__main__":
     from collections import OrderedDict