Implement the Poynting vector calculation

docs/photonics/examples/effective_area.py

@@ -48,10 +48,12 @@
 neff_dict = dict()
 aeff_dict = dict()
 tm_dict = dict()
+p_dict = dict()
 for h in h_list:
     neff_list = []
     aeff_list = []
     tm_list = []
+    p_list = []
     for width in w_list:
         width = width * 1e-3
         nc = shapely.geometry.box(
@@ -95,12 +97,14 @@
                 neff_list.append(np.real(mode.n_eff))
                 aeff_list.append(mode.calculate_effective_area())
                 tm_list.append(mode.transversality)
+                p_list.append(mode.Pz)
                 break
         else:
             print(f"no TM mode found for {width}")
     neff_dict[str(h)] = neff_list
     aeff_dict[str(h)] = aeff_list
     tm_dict[str(h)] = tm_list
+    p_dict[str(h)] = p_list
 # %% [markdown]
 # Plot the result
 
@@ -184,3 +188,52 @@
 ax4.legend()
 
 plt.show()
+
+# %% [markdown]
+# We can also plot the modes to compare them with the reference article
+
+# %%
+# Data figure h
+w_fig_h = 500
+idx_fig_h = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_h))
+h_fig_h = 0.7
+basis_fig_h, Pz_fig_h = p_dict[str(h_fig_h)][idx_fig_h]
+
+# Data figure f
+w_fig_f = 600
+idx_fig_f = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_f))
+h_fig_f = 0.5
+basis_fig_f, Pz_fig_f = p_dict[str(h_fig_f)][idx_fig_f]
+
+fig, ax = plt.subplots(1, 2)
+
+basis_fig_h.plot(np.abs(Pz_fig_h), ax=ax[0], aspect="equal")
+ax[0].set_title(
+    f"Poynting vector $S_z$\nfor h = {h_fig_h}μm & w = {w_fig_h}nm\n(Reproduction of Fig.1.h)"
+)
+ax[0].set_xlim(-w_fig_h * 1e-3 / 2 - 0.1, w_fig_h * 1e-3 / 2 + 0.1)
+ax[0].set_ylim(capital_h - 0.1, capital_h + h_fig_h + 0.1)
+# Turn off the axis wrt to the article figure
+ax[0].axis("off")
+# Add the contour
+for subdomain in basis_fig_h.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
+    basis_fig_h.mesh.restrict(subdomain).draw(
+        ax=ax[0], boundaries_only=True, color="k", linewidth=1.0
+    )
+
+basis_fig_f.plot(np.abs(Pz_fig_f), ax=ax[1], aspect="equal")
+ax[1].set_title(
+    f"Poynting vector $S_z$\nfor h = {h_fig_f}μm & w = {w_fig_f}nm\n(Reproduction of Fig.1.f)"
+)
+ax[1].set_xlim(-w_list[idx_fig_f] * 1e-3 / 2 - 0.1, w_list[idx_fig_f] * 1e-3 / 2 + 0.1)
+ax[1].set_ylim(capital_h - 0.1, capital_h + h_fig_f + 0.1)
+# Turn off the axis wrt to the article figure
+ax[1].axis("off")
+# Add the contour
+for subdomain in basis_fig_f.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
+    basis_fig_f.mesh.restrict(subdomain).draw(
+        ax=ax[1], boundaries_only=True, color="k", linewidth=1.0
+    )
+
+fig.tight_layout()
+plt.show()

femwell/maxwell/waveguide.py

@@ -69,6 +69,45 @@ def n_eff(self):
         """Effective refractive index of the mode"""
         return self.k / self.k0
 
+    @cached_property
+    def poynting(self):
+        """Poynting vector of the mode"""
+
+        # Extraction of the fields
+        (Ex, Ey), Ez = self.basis.interpolate(self.E)
+        (Hx, Hy), Hz = self.basis.interpolate(self.H)
+
+        # New basis will be the discontinuous variant used by the solved Ez-field
+        (Et, Et_basis), (En, En_basis) = self.basis.split(self.E)
+        poynting_basis = self.basis.with_element(ElementDG(ElementTriP1()))
+
+        # Calculation of the Poynting vector
+        Px = Ey * Hz - Ez * Hy
+        Py = Ez * Hx - Ex * Hz
+        Pz = Ex * Hy - Ey * Hx
+
+        # Projection of the Poynting vector on the new basis
+        Px_proj = poynting_basis.project(Px, dtype=np.complex64)
+        Py_proj = poynting_basis.project(Py, dtype=np.complex64)
+        Pz_proj = poynting_basis.project(Pz, dtype=np.complex64)
+
+        return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj])
+
+    @cached_property
+    def Px(self):
+        basis, _P = self.poynting
+        return basis, _P[0]
+
+    @cached_property
+    def Py(self):
+        basis, _P = self.poynting
+        return basis, _P[1]
+
+    @cached_property
+    def Pz(self):
+        basis, _P = self.poynting
+        return basis, _P[2]
+
     @cached_property
     def te_fraction(self):
         """TE-fraction of the mode"""
@@ -330,7 +369,10 @@ def compute_modes(
     if order == 1:
         element = ElementTriN1() * ElementTriP1()
     elif order == 2:
-        element = ElementTriN2() * ElementTriP2()
+        # element = ElementTriN2() * ElementTriP2()
+        from skfem.element import ElementQuad2
+
+        element = ElementTriN2() * ElementQuad2  # ElementTriP2()
     else:
         raise AssertionError("Only order 1 and 2 implemented by now.")