Fem: Add electric charge density object

src/Mod/Fem/femsolver/elmer/equations/electrostatic_writer.py

@@ -30,6 +30,8 @@ __url__ = "https://www.freecad.org"
 ## \addtogroup FEM
 #  @{
 
+from FreeCAD import Units
+
 from .. import sifio
 
 
@@ -133,8 +135,8 @@ class ESwriter:
                     elif obj.BoundaryCondition == "Neumann":
                         self.write.boundary(
                             name,
-                            "Surface Charge Density",
-                            obj.SurfaceChargeDensity.getValueAs("C/m^2").Value,
+                            "Electric Flux",
+                            obj.ElectricFluxDensity.getValueAs("C/m^2").Value,
                         )
                     if obj.PotentialConstant:
                         self.write.boundary(name, "Potential Constant", True)
@@ -146,5 +148,60 @@ class ESwriter:
                         self.write.boundary(name, "Capacitance Body", obj.CapacitanceBody)
                 self.write.handled(obj)
 
+        for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
+            if obj.Mode not in ["Interface", "Total Interface"]:
+                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.handled(obj)
+
+    def handleElectrostaticBodyForces(self):
+        for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
+            if obj.Mode not in ["Source", "Total Source"]:
+                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.handled(obj)
+
 
 ##  @}

src/Mod/Fem/femsolver/elmer/writer.py

@@ -471,6 +471,7 @@ class Writer:
         if activeIn:
             ESW.handleElectrostaticConstants()
             ESW.handleElectrostaticBndConditions()
+            ESW.handleElectrostaticBodyForces()
             ESW.handleElectrostaticMaterial(activeIn)
 
     # -------------------------------------------------------------------------------------------