From 464df19cbbd3fc3387a5c912b69fcd4b8463d2fa Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 21 Feb 2024 01:32:26 -0800
Subject: [PATCH] add integrals for (1) and (3)

---
 docs/photonics/examples/effective_area.py | 32 +++++++++++++++++++++--
 femwell/maxwell/waveguide.py              |  9 +++----
 2 files changed, 33 insertions(+), 8 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index ddafb02a..13349303 100644
--- 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
index d6e750e8..7d8548cb 100644
--- 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]]