diff --git a/src/Mod/Fem/femobjects/constraint_electricchargedensity.py b/src/Mod/Fem/femobjects/constraint_electricchargedensity.py
index fb2b9a6b16..ed4d0f7634 100644
--- a/src/Mod/Fem/femobjects/constraint_electricchargedensity.py
+++ b/src/Mod/Fem/femobjects/constraint_electricchargedensity.py
@@ -29,6 +29,9 @@ __url__ = "https://www.freecad.org"
 #  \ingroup FEM
 #  \brief constraint electric charge density object
 
+
+from FreeCAD import Units
+
 from . import base_fempythonobject
 
 _PropHelper = base_fempythonobject._PropHelper
@@ -85,3 +88,43 @@ class ConstraintElectricChargeDensity(base_fempythonobject.BaseFemPythonObject):
         )
 
         return prop
+
+    def get_total_source_density(self, obj):
+        """
+        Calculate density for `Total Source` mode.
+        """
+        size = 0
+        items = []
+        for feat, sub_elem in obj.References:
+            for name in sub_elem:
+                sub = feat.getSubObject(name)
+                if sub.ShapeType == "Solid":
+                    size += sub.Volume
+                    items.append(name)
+                elif sub.ShapeType == "Face":
+                    size += sub.Area
+                    items.append(name)
+
+        if items:
+            vol = Units.Quantity(f"{size} mm^3")
+            return obj.TotalCharge / vol
+
+    def get_total_interface_density(self, obj):
+        """
+        Calculate density for `Total Interface` mode.
+        """
+        size = 0
+        items = []
+        for feat, sub_elem in obj.References:
+            for name in sub_elem:
+                sub = feat.getSubObject(name)
+                if sub.ShapeType == "Face":
+                    size += sub.Area
+                    items.append(name)
+                elif sub.ShapeType == "Edge":
+                    size += sub.Length
+                    items.append(name)
+
+        if items:
+            area = Units.Quantity(f"{size} mm^2")
+            return obj.TotalCharge / area
diff --git a/src/Mod/Fem/femsolver/elmer/equations/electrostatic_writer.py b/src/Mod/Fem/femsolver/elmer/equations/electrostatic_writer.py
index 0d1e93903c..752704860f 100644
--- a/src/Mod/Fem/femsolver/elmer/equations/electrostatic_writer.py
+++ b/src/Mod/Fem/femsolver/elmer/equations/electrostatic_writer.py
@@ -152,58 +152,38 @@ class ESwriter:
                 self.write.handled(obj)
 
         for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
-            if obj.Mode not in ["Interface", "Total Interface"]:
-                continue
+            match obj.Mode:
+                case "Interface":
+                    density = obj.InterfaceChargeDensity
+                case "Total Interface":
+                    density = obj.Proxy.get_total_interface_density(obj)
+                case _:
+                    continue
 
-            size = 0
-            items = []
             for feat, sub_elem in obj.References:
                 for name in sub_elem:
-                    sub = feat.getSubObject(name)
-                    if sub.ShapeType == "Face":
-                        size += sub.Area
-                        items.append(name)
-                    elif sub.ShapeType == "Edge":
-                        size += sub.Length
-                        items.append(name)
-
-            if items:
-                if obj.Mode == "Interface":
-                    density = obj.InterfaceChargeDensity.getValueAs("C/m^2").Value
-                elif obj.Mode == "Total Interface":
-                    area = Units.Quantity(f"{size} mm^2")
-                    density = (obj.TotalCharge / area).getValueAs("C/m^2").Value
-                for name in items:
                     self.write.boundary(name, "! FreeCAD Name", obj.Label)
-                    self.write.boundary(name, "Surface Charge Density", round(density, 6))
+                    self.write.boundary(
+                        name, "Surface Charge Density", round(density.getValueAs("C/m^2").Value, 6)
+                    )
                     self.write.handled(obj)
 
     def handleElectrostaticBodyForces(self):
         for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
-            if obj.Mode not in ["Source", "Total Source"]:
-                continue
+            match obj.Mode:
+                case "Source":
+                    density = obj.SourceChargeDensity
+                case "Total Source":
+                    density = obj.Proxy.get_total_source_density(obj)
+                case _:
+                    continue
 
-            size = 0
-            items = []
             for feat, sub_elem in obj.References:
                 for name in sub_elem:
-                    sub = feat.getSubObject(name)
-                    if sub.ShapeType == "Solid":
-                        size += sub.Volume
-                        items.append(name)
-                    elif sub.ShapeType == "Face":
-                        size += sub.Area
-                        items.append(name)
-
-            if items:
-                if obj.Mode == "Source":
-                    density = obj.SourceChargeDensity.getValueAs("C/m^3").Value
-                elif obj.Mode == "Total Source":
-                    vol = Units.Quantity(f"{size} mm^3")
-                    density = (obj.TotalCharge / vol).getValueAs("C/m^3").Value
-                for name in items:
                     self.write.bodyForce(name, "! FreeCAD Name", obj.Label)
-                    self.write.bodyForce(name, "Charge Density", round(density, 6))
+                    self.write.bodyForce(
+                        name, "Charge Density", round(density.getValueAs("C/m^3").Value, 6)
+                    )
                     self.write.handled(obj)