add integrals for (1) and (3)

HelgeGehring · Feb 21, 2024 · 464df19 · 464df19
1 parent bdcd730
commit 464df19
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 8 deletions.
diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
@@ -22,7 +22,7 @@
 import shapely.affinity
 from matplotlib.ticker import MultipleLocator
 from shapely.ops import clip_by_rect
-from skfem import Basis, ElementTriP0
+from skfem import Basis, ElementTriP0, Functional
 from skfem.io.meshio import from_meshio
 
 from femwell.maxwell.waveguide import compute_modes
@@ -97,7 +97,35 @@
                 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)
+                p_list.append(mode.Sz)
+
+                @Functional
+                def I(w):
+                    return 1
+
+                @Functional
+                def Sz(w):
+                    return w["Sz"]
+
+                @Functional
+                def Sz2(w):
+                    return w["Sz"] ** 2
+
+                Sz_basis, Sz_vec = mode.Sz
+
+                print(
+                    "int(Sz)", Sz.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec))
+                )  # 1 as it's normalized
+                print("int_core(1)", I.assemble(Sz_basis.with_elements("core")))  # area of core
+                print(
+                    "int_core(Sz)",
+                    Sz.assemble(
+                        Sz_basis.with_elements("core"),
+                        Sz=Sz_basis.with_elements("core").interpolate(Sz_vec),
+                    ),
+                )
+                print("int(Sz^2)", Sz2.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec)))
+
                 break
         else:
             print(f"no TM mode found for {width}")

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
@@ -92,18 +92,15 @@ def poynting(self):
 
         return poynting_basis, P_proj
 
-    @cached_property
-    def Px(self):
+    def Sx(self):
         basis, _P = self.poynting
         return basis.split_bases()[0], _P[basis.split_indices()[0]]
 
-    @cached_property
-    def Py(self):
+    def Sy(self):
         basis, _P = self.poynting
         return basis.split_bases()[1], _P[basis.split_indices()[1]]
 
-    @cached_property
-    def Pz(self):
+    def Sz(self):
         basis, _P = self.poynting
         return basis.split_bases()[2], _P[basis.split_indices()[2]]