From 890e65be576099209c374a4c6112436193a00b1a Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Thu, 28 Dec 2023 00:42:36 -0800
Subject: [PATCH 01/19] feat: add calculate_poynting to Mode

---
 femwell/maxwell/waveguide.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 68bffc52..3d766e6e 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -175,6 +175,26 @@ def E4(w):
     def calculate_propagation_loss(self, distance):
         return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance
 
+    def calculate_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 = 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 = new_basis.project(Px, dtype=np.complex64)
+        Py_proj = new_basis.project(Py, dtype=np.complex64)
+        Pz_proj = new_basis.project(Pz, dtype=np.complex64)
+
+        return new_basis, np.array([Px_proj, Py_proj, Pz_proj])
+
     def calculate_power(self, elements=None):
         if not elements:
             basis = self.basis

From d5443de80acab1b901a9b7ae9c59258831bed51d Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Thu, 28 Dec 2023 01:29:03 -0800
Subject: [PATCH 02/19] doc: implement Poynting plot for effective_area example

---
 docs/photonics/examples/effective_area.py | 45 +++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 74ac9924..ca026c95 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/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.calculate_poynting())
                 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,44 @@
 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, P_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, P_fig_f = p_dict[str(h_fig_f)][idx_fig_f]
+
+fig, ax = plt.subplots(1, 2)
+
+basis_fig_h.plot(np.abs(P_fig_h[2]), 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(P_fig_f[2]), 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()

From e7615e09303bf44e8c247a7b45fcb7ea8f8348ac Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Thu, 28 Dec 2023 01:54:02 -0800
Subject: [PATCH 03/19] fix: black refactor

---
 docs/photonics/examples/effective_area.py | 24 +++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index ca026c95..fc99173b 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -208,24 +208,32 @@
 fig, ax = plt.subplots(1, 2)
 
 basis_fig_h.plot(np.abs(P_fig_h[2]), 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_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')
+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_h.mesh.restrict(subdomain).draw(
+        ax=ax[0], boundaries_only=True, color="k", linewidth=1.0
+    )
 
 basis_fig_f.plot(np.abs(P_fig_f[2]), 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_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')
+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)
+    basis_fig_f.mesh.restrict(subdomain).draw(
+        ax=ax[1], boundaries_only=True, color="k", linewidth=1.0
+    )
 
 fig.tight_layout()
 plt.show()

From bb320eed2e1949187590fa04be102c790a21960a Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Fri, 29 Dec 2023 00:09:43 -0800
Subject: [PATCH 04/19] feat: poynting calc with use the same element as solved
 Ez

---
 femwell/maxwell/waveguide.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 3d766e6e..b02d0662 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -181,6 +181,9 @@ def calculate_poynting(self):
         # 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)
         new_basis = self.basis.with_element(ElementDG(ElementTriP1()))
 
         # Calculation of the Poynting vector
@@ -350,7 +353,9 @@ 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.")
 

From dc4eb9361d2fb6f7cb61a522efadab2708c685a8 Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Fri, 29 Dec 2023 00:16:06 -0800
Subject: [PATCH 05/19] feat: add poynting as a Mode (cached) property

---
 docs/photonics/examples/effective_area.py |  2 +-
 femwell/maxwell/waveguide.py              | 46 ++++++++++++-----------
 2 files changed, 25 insertions(+), 23 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index fc99173b..226cc23d 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -97,7 +97,7 @@
                 neff_list.append(np.real(mode.n_eff))
                 aeff_list.append(mode.calculate_effective_area())
                 tm_list.append(mode.transversality)
-                p_list.append(mode.calculate_poynting())
+                p_list.append(mode.poynting)
                 break
         else:
             print(f"no TM mode found for {width}")
diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index b02d0662..b771baee 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -69,6 +69,30 @@ 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)
+        new_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 = new_basis.project(Px, dtype=np.complex64)
+        Py_proj = new_basis.project(Py, dtype=np.complex64)
+        Pz_proj = new_basis.project(Pz, dtype=np.complex64)
+
+        return new_basis, np.array([Px_proj, Py_proj, Pz_proj])
+
     @cached_property
     def te_fraction(self):
         """TE-fraction of the mode"""
@@ -175,28 +199,6 @@ def E4(w):
     def calculate_propagation_loss(self, distance):
         return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance
 
-    def calculate_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)
-        new_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 = new_basis.project(Px, dtype=np.complex64)
-        Py_proj = new_basis.project(Py, dtype=np.complex64)
-        Pz_proj = new_basis.project(Pz, dtype=np.complex64)
-
-        return new_basis, np.array([Px_proj, Py_proj, Pz_proj])
 
     def calculate_power(self, elements=None):
         if not elements:

From b12affefa8dbbdc53505b904e6a6d7755eb06575 Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Fri, 29 Dec 2023 00:16:36 -0800
Subject: [PATCH 06/19] refactor: black compliant

---
 femwell/maxwell/waveguide.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index b771baee..3d8d975c 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -199,7 +199,6 @@ def E4(w):
     def calculate_propagation_loss(self, distance):
         return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance
 
-
     def calculate_power(self, elements=None):
         if not elements:
             basis = self.basis
@@ -357,7 +356,8 @@ def compute_modes(
     elif order == 2:
         # element = ElementTriN2() * ElementTriP2()
         from skfem.element import ElementQuad2
-        element = ElementTriN2() * ElementQuad2#ElementTriP2()
+
+        element = ElementTriN2() * ElementQuad2  # ElementTriP2()
     else:
         raise AssertionError("Only order 1 and 2 implemented by now.")
 

From 560f62745a2209ceca5cf92a8b3b5e18ecf797ec Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Fri, 29 Dec 2023 00:21:14 -0800
Subject: [PATCH 07/19] fix: change basis name to make it more clear

---
 femwell/maxwell/waveguide.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 3d8d975c..2625e7f6 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -79,7 +79,7 @@ def poynting(self):
 
         # New basis will be the discontinuous variant used by the solved Ez-field
         (Et, Et_basis), (En, En_basis) = self.basis.split(self.E)
-        new_basis = self.basis.with_element(ElementDG(ElementTriP1()))
+        poynting_basis = self.basis.with_element(ElementDG(ElementTriP1()))
 
         # Calculation of the Poynting vector
         Px = Ey * Hz - Ez * Hy
@@ -87,11 +87,11 @@ def poynting(self):
         Pz = Ex * Hy - Ey * Hx
 
         # Projection of the Poynting vector on the new basis
-        Px_proj = new_basis.project(Px, dtype=np.complex64)
-        Py_proj = new_basis.project(Py, dtype=np.complex64)
-        Pz_proj = new_basis.project(Pz, dtype=np.complex64)
+        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 new_basis, np.array([Px_proj, Py_proj, Pz_proj])
+        return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj])
 
     @cached_property
     def te_fraction(self):

From 5ab12bdf67bda30f20498df763d607c9babf16f9 Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Fri, 29 Dec 2023 00:26:52 -0800
Subject: [PATCH 08/19] feat: add poynting component as Mode properties

---
 docs/photonics/examples/effective_area.py | 10 +++++-----
 femwell/maxwell/waveguide.py              | 15 +++++++++++++++
 2 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 226cc23d..ddafb02a 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -97,7 +97,7 @@
                 neff_list.append(np.real(mode.n_eff))
                 aeff_list.append(mode.calculate_effective_area())
                 tm_list.append(mode.transversality)
-                p_list.append(mode.poynting)
+                p_list.append(mode.Pz)
                 break
         else:
             print(f"no TM mode found for {width}")
@@ -197,17 +197,17 @@
 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, P_fig_h = p_dict[str(h_fig_h)][idx_fig_h]
+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, P_fig_f = p_dict[str(h_fig_f)][idx_fig_f]
+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(P_fig_h[2]), ax=ax[0], aspect="equal")
+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)"
 )
@@ -221,7 +221,7 @@
         ax=ax[0], boundaries_only=True, color="k", linewidth=1.0
     )
 
-basis_fig_f.plot(np.abs(P_fig_f[2]), ax=ax[1], aspect="equal")
+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)"
 )
diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 2625e7f6..7b4e0298 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -93,6 +93,21 @@ def poynting(self):
 
         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"""

From eaa5a36a0a136f4d4ae5230af53a08a910caf26e Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Sat, 30 Dec 2023 16:25:09 -0800
Subject: [PATCH 09/19] fix: restore ElementTriP2 instead of Quad2 for order=2

---
 femwell/maxwell/waveguide.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 7b4e0298..10e383d4 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -369,10 +369,7 @@ def compute_modes(
     if order == 1:
         element = ElementTriN1() * ElementTriP1()
     elif order == 2:
-        # element = ElementTriN2() * ElementTriP2()
-        from skfem.element import ElementQuad2
-
-        element = ElementTriN2() * ElementQuad2  # ElementTriP2()
+        element = ElementTriN2() * ElementTriP2()
     else:
         raise AssertionError("Only order 1 and 2 implemented by now.")
 

From 98cb85359c28d5cf17411b1c0d6874e5c3761cce Mon Sep 17 00:00:00 2001
From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com>
Date: Sat, 30 Dec 2023 16:54:57 -0800
Subject: [PATCH 10/19] feat: for poynting, auto extraction of Ez elements

---
 femwell/maxwell/waveguide.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 10e383d4..5a720831 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -78,8 +78,7 @@ def poynting(self):
         (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()))
+        poynting_basis = self.basis.with_element(ElementDG(self.basis.elem.elems[1]))
 
         # Calculation of the Poynting vector
         Px = Ey * Hz - Ez * Hy

From 9f0397a6b90e936d5da526daa389336e1c26cb32 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 10 Jan 2024 11:01:57 -0800
Subject: [PATCH 11/19] use vectorbasis instead of three individual bases for
 the poynting vector

---
 femwell/maxwell/waveguide.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 5a720831..d6e750e8 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -78,7 +78,9 @@ def poynting(self):
         (Hx, Hy), Hz = self.basis.interpolate(self.H)
 
         # New basis will be the discontinuous variant used by the solved Ez-field
-        poynting_basis = self.basis.with_element(ElementDG(self.basis.elem.elems[1]))
+        poynting_basis = self.basis.with_element(
+            ElementVector(ElementDG(self.basis.elem.elems[1]), 3)
+        )
 
         # Calculation of the Poynting vector
         Px = Ey * Hz - Ez * Hy
@@ -86,26 +88,24 @@ def poynting(self):
         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)
+        P_proj = poynting_basis.project(np.stack([Px, Py, Pz], axis=0), dtype=np.complex64)
 
-        return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj])
+        return poynting_basis, P_proj
 
     @cached_property
     def Px(self):
         basis, _P = self.poynting
-        return basis, _P[0]
+        return basis.split_bases()[0], _P[basis.split_indices()[0]]
 
     @cached_property
     def Py(self):
         basis, _P = self.poynting
-        return basis, _P[1]
+        return basis.split_bases()[1], _P[basis.split_indices()[1]]
 
     @cached_property
     def Pz(self):
         basis, _P = self.poynting
-        return basis, _P[2]
+        return basis.split_bases()[2], _P[basis.split_indices()[2]]
 
     @cached_property
     def te_fraction(self):

From 9f8c5b7c8367657c65dd3dddbd60c56470e8f2a8 Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Thu, 18 Jan 2024 05:00:24 +0000
Subject: [PATCH 12/19] Bump actions/cache from 3 to 4

Bumps [actions/cache](https://github.com/actions/cache) from 3 to 4.
- [Release notes](https://github.com/actions/cache/releases)
- [Changelog](https://github.com/actions/cache/blob/main/RELEASES.md)
- [Commits](https://github.com/actions/cache/compare/v3...v4)

---
updated-dependencies:
- dependency-name: actions/cache
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
---
 .github/workflows/docs.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 172929f9..30d6db07 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -60,7 +60,7 @@ jobs:
         run: |
           pip install juliacall
           julia -e 'using Pkg; Pkg.add("PythonCall")'
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         with:
           path: |
             docs/_build/.jupyter_cache

From b49765902469faf37b758693d385319fd8f30dac Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Tue, 30 Jan 2024 00:48:56 -0800
Subject: [PATCH 13/19] try not pinning sphinx/docutils

---
 .github/workflows/docs.yml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 30d6db07..1735af7e 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -43,8 +43,6 @@ jobs:
           mamba install julia
           pip install jupyter-book
           pip install --upgrade jupytext
-          pip install sphinx=="5.3.0"
-          pip install docutils=="0.17.1"
           pip install gdsfactory gplugins
           pip list
       

From c9aabb8716e7d721de9596d6b73350bfd1dbc92a Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Thu, 8 Feb 2024 04:30:56 +0000
Subject: [PATCH 14/19] Bump pre-commit/action from 3.0.0 to 3.0.1

Bumps [pre-commit/action](https://github.com/pre-commit/action) from 3.0.0 to 3.0.1.
- [Release notes](https://github.com/pre-commit/action/releases)
- [Commits](https://github.com/pre-commit/action/compare/v3.0.0...v3.0.1)

---
updated-dependencies:
- dependency-name: pre-commit/action
  dependency-type: direct:production
  update-type: version-update:semver-patch
...

Signed-off-by: dependabot[bot] <support@github.com>
---
 .github/workflows/build.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 1d930728..a26d7dd2 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -28,7 +28,7 @@ jobs:
       - name: Install JuliaFormatter and format
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
-      - uses: pre-commit/action@v3.0.0
+      - uses: pre-commit/action@v3.0.1
 
   flake8-lint:
     runs-on: ubuntu-latest

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 15/19] 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]]
 

From bed87732e537cdde950acd5f4b9db8ce49146da4 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 21 Feb 2024 17:09:50 -0800
Subject: [PATCH 16/19] start example on effective thermal conductivity

---
 docs/julia/thermal_fill.jl | 190 +++++++++++++++++++++++++++++++++++++
 1 file changed, 190 insertions(+)
 create mode 100644 docs/julia/thermal_fill.jl

diff --git a/docs/julia/thermal_fill.jl b/docs/julia/thermal_fill.jl
new file mode 100644
index 00000000..b4214893
--- /dev/null
+++ b/docs/julia/thermal_fill.jl
@@ -0,0 +1,190 @@
+# ---
+# jupyter:
+#   jupytext:
+#     custom_cell_magics: kql
+#     formats: jl:percent,ipynb
+#     text_representation:
+#       extension: .jl
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: base
+#     language: julia
+#     name: julia-1.9
+# ---
+
+# %% [markdown]
+# # Effective thermal conductivity
+
+# %% tags=["hide-input", "thebe-init", "remove-stderr"]
+using PyCall
+np = pyimport("numpy")
+shapely = pyimport("shapely")
+shapely.affinity = pyimport("shapely.affinity")
+OrderedDict = pyimport("collections").OrderedDict
+meshwell = pyimport("meshwell")
+Prism = pyimport("meshwell.prism").Prism
+Model = pyimport("meshwell.model").Model
+
+model = Model()
+
+unitcell_width = 2
+unitcell_length = 2
+unitcell_thickness = 2
+
+fill_width = 2
+fill_length = 1
+fill_thickness = 2
+
+fill_polygon = shapely.box(
+    xmin = unitcell_length / 2 - fill_length / 2,
+    ymin = unitcell_width / 2 - fill_width / 2,
+    xmax = unitcell_length / 2 + fill_length / 2,
+    ymax = unitcell_width / 2 + fill_width / 2,
+)
+
+fill = Prism(
+    polygons = fill_polygon,
+    buffers = Dict(0 => 0.0, fill_thickness => 0.0),
+    model = model,
+    physical_name = "fill",
+    mesh_order = 1,
+    resolution = Dict("resolution" => 0.2, "SizeMax" => 1.0, "DistMax" => 1.0),
+)
+
+unitcell_polygon =
+    shapely.box(xmin = 0, ymin = 0, xmax = unitcell_length, ymax = unitcell_width)
+unitcell_buffer = Dict(0 => 0.0, unitcell_thickness => 0.0)
+
+unitcell = Prism(
+    polygons = unitcell_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "unitcell",
+    mesh_order = 2,
+)
+
+"""
+BOUNDARIES
+"""
+
+right_polygon = shapely.box(
+    xmin = unitcell_length,
+    ymin = 0,
+    xmax = unitcell_length + 1,
+    ymax = unitcell_width,
+)
+right = Prism(
+    polygons = right_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "right",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+left_polygon = shapely.box(xmin = -1, ymin = 0, xmax = 0, ymax = unitcell_width)
+left = Prism(
+    polygons = left_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "left",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+up_polygon = shapely.box(
+    xmin = 0,
+    ymin = unitcell_width,
+    xmax = unitcell_length,
+    ymax = unitcell_width + 1,
+)
+up = Prism(
+    polygons = up_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "up",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+down_polygon = shapely.box(xmin = 0, ymin = -1, xmax = unitcell_length, ymax = 0)
+down = Prism(
+    polygons = down_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "down",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+"""
+ASSEMBLE AND NAME ENTITIES
+"""
+entities = [unitcell, fill, up, down, left, right]
+
+mesh = model.mesh(
+    entities_list = entities,
+    verbosity = 0,
+    filename = "mesh.msh",
+    default_characteristic_length = 0.2,
+)
+
+
+
+# %% tags=["hide-input", "thebe-init", "remove-stderr"]
+using Gridap
+using GridapGmsh
+using Gridap.Geometry
+using GridapMakie, CairoMakie
+
+using Femwell.Maxwell.Electrostatic
+using Femwell.Thermal
+
+
+# %%
+thermal_conductivities = ["unitcell" => 1, "fill" => 100.0]
+# %% [markdown]
+# First we load the mesh, get the labels, define a CellField to evaluate the constants on
+# and create functions which work on the tags
+
+# %% tags=["remove-stderr", "hide-output"]
+model = GmshDiscreteModel("mesh.msh")
+Ω = Triangulation(model)
+dΩ = Measure(Ω, 1)
+
+labels = get_face_labeling(model)
+tags = get_face_tag(labels, num_cell_dims(model))
+τ = CellField(tags, Ω)
+
+thermal_conductivities =
+    Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_conductivities)
+ϵ_conductivities(tag) = thermal_conductivities[tag]
+
+# %% [markdown]
+# The next step is to define the boundary conditions, this can be done simply via julia-dicts:
+
+# %% tags=["remove-stderr", "hide-output"]
+boundary_temperatures = Dict("left___unitcell" => 100.0, "right___unitcell" => 0.0)
+
+# %% tags=["remove-stderr"]
+T0 = calculate_temperature(
+    ϵ_conductivities ∘ τ,
+    CellField(x -> 0, Ω),
+    boundary_temperatures,
+    order = 2,
+)
+temperature_difference = abs(sum(values(boundary_temperatures) .* [-1, 1]))
+power = abs(
+    sum(
+        ∫(
+            (ϵ_conductivities ∘ τ) *
+            (norm ∘ gradient(temperature(T0))) *
+            (norm ∘ gradient(temperature(T0))),
+        )dΩ,
+    ),
+)
+#println(∑(∫(1)dΩ))
+(power / temperature_difference^2) /
+((unitcell_thickness * unitcell_width) / unitcell_length)

From c87e426ac952a083464e6e1c80a9b752ff3e225a Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 21 Feb 2024 17:28:28 -0800
Subject: [PATCH 17/19] add text

---
 docs/julia/thermal_fill.jl | 56 ++++++++++++++++++++++++++------------
 1 file changed, 39 insertions(+), 17 deletions(-)

diff --git a/docs/julia/thermal_fill.jl b/docs/julia/thermal_fill.jl
index b4214893..581c36bd 100644
--- a/docs/julia/thermal_fill.jl
+++ b/docs/julia/thermal_fill.jl
@@ -11,12 +11,32 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
 # # Effective thermal conductivity
 
+# Usually, after the design of the circuitry, the empty space on chips gets filled with equally spaced small rectangles in each layer.
+# Optically, these structures are supposed to be far enough away that their influence on the structures can be neglected.
+# But for thermal considerations, those fill structures can have an impact on the temperature distribution on the chip and thus e.g. on the crosstalk between thermal phase shifters.
+# As it's computationally challenging to include all the small cuboids in the model (which is especially for the meshing a major challenge),
+# a preferable approach is to consider the filled area as a homogenous area of higher thermal conductivity.
+# For this, we calculate the effective thermal conductivity of the filled area by examining a single unit cell.
+
+# To have an intuitively understandable problem, we consider half of the unit cell to be filled with a highly thermally conductive material (metal/silicon) surrounded by a material with low thermal conductance (e.g. silicon dioxide)
+# Let's start with the mesh!
+
+# %% 
+
+unitcell_width = 2
+unitcell_length = 2
+unitcell_thickness = 2
+
+fill_width = 2
+fill_length = 1
+fill_thickness = 2
+
 # %% tags=["hide-input", "thebe-init", "remove-stderr"]
 using PyCall
 np = pyimport("numpy")
@@ -29,14 +49,6 @@ Model = pyimport("meshwell.model").Model
 
 model = Model()
 
-unitcell_width = 2
-unitcell_length = 2
-unitcell_thickness = 2
-
-fill_width = 2
-fill_length = 1
-fill_thickness = 2
-
 fill_polygon = shapely.box(
     xmin = unitcell_length / 2 - fill_length / 2,
     ymin = unitcell_width / 2 - fill_width / 2,
@@ -142,14 +154,12 @@ using GridapMakie, CairoMakie
 using Femwell.Maxwell.Electrostatic
 using Femwell.Thermal
 
+# %% [markdown]
+# For simplicity, we define the thermal conductivity of the unitcell to be 1 and the thermal conductivity of the fill to be 100, which is almost negligible in comparison.
 
 # %%
 thermal_conductivities = ["unitcell" => 1, "fill" => 100.0]
-# %% [markdown]
-# First we load the mesh, get the labels, define a CellField to evaluate the constants on
-# and create functions which work on the tags
 
-# %% tags=["remove-stderr", "hide-output"]
 model = GmshDiscreteModel("mesh.msh")
 Ω = Triangulation(model)
 dΩ = Measure(Ω, 1)
@@ -163,11 +173,20 @@ thermal_conductivities =
 ϵ_conductivities(tag) = thermal_conductivities[tag]
 
 # %% [markdown]
-# The next step is to define the boundary conditions, this can be done simply via julia-dicts:
+# We define the temperatures at both sides of the unitcell to define the direction in which we want to estimate the thermal conductivity.
+# In other directions the boundaries are considered to be insulating/mirroring.
 
 # %% tags=["remove-stderr", "hide-output"]
 boundary_temperatures = Dict("left___unitcell" => 100.0, "right___unitcell" => 0.0)
 
+
+# %% [markdown]
+# We start with calculating the temperature distribution.
+# From this, we calculate, analog to how we calculate resistances for electrical simulations first the integral
+# $\int σ \left|\frac{\mathrm{d}T}{\mathrm{d}\vec{x}}\right|^2 dA which gives an equivalent of the power
+# and calculate from this the effective thermal conductivity.
+# We expect the thermal conductivity almost twice as high as the thermal conductivity of the material with lower thermal conductivity.
+
 # %% tags=["remove-stderr"]
 T0 = calculate_temperature(
     ϵ_conductivities ∘ τ,
@@ -185,6 +204,9 @@ power = abs(
         )dΩ,
     ),
 )
-#println(∑(∫(1)dΩ))
-(power / temperature_difference^2) /
-((unitcell_thickness * unitcell_width) / unitcell_length)
+
+println(
+    "Effective thermal conductivity: ",
+    (power / temperature_difference^2) /
+    ((unitcell_thickness * unitcell_width) / unitcell_length),
+)

From ab4b6025f4fd0b66c65ebf163309a2c3d59bd0ef Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 21 Feb 2024 17:30:37 -0800
Subject: [PATCH 18/19] add example to toc

---
 docs/_toc.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/_toc.yml b/docs/_toc.yml
index 00ea7c7b..3c27db04 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -52,6 +52,7 @@ parts:
     - file: julia/lithium_niobate_phase_shifter.jl
     - file: julia/thermal_simple.jl
     - file: julia/heater_3d.jl
+    - file: julia/thermal_fill.jl
 
   - caption: Benchmarks
     chapters:

From 17a1e82c67a497cba538a0f4ae510b4d52a7ce36 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Wed, 21 Feb 2024 17:35:19 -0800
Subject: [PATCH 19/19] fix

---
 docs/julia/thermal_fill.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/julia/thermal_fill.jl b/docs/julia/thermal_fill.jl
index 581c36bd..af329b21 100644
--- a/docs/julia/thermal_fill.jl
+++ b/docs/julia/thermal_fill.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.10
+#     name: julia-1.9
 # ---
 
 # %% [markdown]