diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index a26d7dd2..b30eeffc 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -23,7 +23,7 @@ jobs:
           python-version: "3.12"
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.9.2'
+          version: '1.10.1'
       - uses: julia-actions/cache@v1
       - name: Install JuliaFormatter and format
         run: |
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 1735af7e..110c3e45 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -48,7 +48,7 @@ jobs:
       
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.9.3'
+          version: '1.10.1'
       - uses: julia-actions/cache@v1
         with:
           cache-name: ${{ runner.os }}-${{ steps.get-date.outputs.date }}
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:
diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl
index b51087ba..0010a40a 100644
--- a/docs/julia/heater_3d.jl
+++ b/docs/julia/heater_3d.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/lithium_niobate_phase_shifter.jl b/docs/julia/lithium_niobate_phase_shifter.jl
index c78ecd43..366b7ff0 100644
--- a/docs/julia/lithium_niobate_phase_shifter.jl
+++ b/docs/julia/lithium_niobate_phase_shifter.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/microstrip.jl b/docs/julia/microstrip.jl
index 27ffe145..4d534aac 100644
--- a/docs/julia/microstrip.jl
+++ b/docs/julia/microstrip.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/thermal_fill.jl b/docs/julia/thermal_fill.jl
new file mode 100644
index 00000000..581c36bd
--- /dev/null
+++ b/docs/julia/thermal_fill.jl
@@ -0,0 +1,212 @@
+# ---
+# 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.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")
+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()
+
+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
+
+# %% [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]
+
+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]
+# 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 ∘ τ,
+    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(
+    "Effective thermal conductivity: ",
+    (power / temperature_difference^2) /
+    ((unitcell_thickness * unitcell_width) / unitcell_length),
+)
diff --git a/docs/julia/thermal_simple.jl b/docs/julia/thermal_simple.jl
index f03a7302..83389206 100644
--- a/docs/julia/thermal_simple.jl
+++ b/docs/julia/thermal_simple.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/waveguide.jl b/docs/julia/waveguide.jl
index 9e80f3ce..2f7684d8 100644
--- a/docs/julia/waveguide.jl
+++ b/docs/julia/waveguide.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/waveguide_anisotropic.jl b/docs/julia/waveguide_anisotropic.jl
index 4bf62e8b..39596207 100644
--- a/docs/julia/waveguide_anisotropic.jl
+++ b/docs/julia/waveguide_anisotropic.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/waveguide_bent.jl b/docs/julia/waveguide_bent.jl
index f171918a..beb82e09 100644
--- a/docs/julia/waveguide_bent.jl
+++ b/docs/julia/waveguide_bent.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/waveguide_bent_coupling.jl b/docs/julia/waveguide_bent_coupling.jl
index a0271536..9a3e22f6 100644
--- a/docs/julia/waveguide_bent_coupling.jl
+++ b/docs/julia/waveguide_bent_coupling.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]
diff --git a/docs/julia/waveguide_overlap_integral.jl b/docs/julia/waveguide_overlap_integral.jl
index 6600d927..ed8c5b21 100644
--- a/docs/julia/waveguide_overlap_integral.jl
+++ b/docs/julia/waveguide_overlap_integral.jl
@@ -11,7 +11,7 @@
 #   kernelspec:
 #     display_name: base
 #     language: julia
-#     name: julia-1.9
+#     name: julia-1.10
 # ---
 
 # %% [markdown]