diff --git a/src/Mod/Fem/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py
index e0127662ca..053ea28a97 100644
--- a/src/Mod/Fem/femmesh/meshsetsgetter.py
+++ b/src/Mod/Fem/femmesh/meshsetsgetter.py
@@ -393,15 +393,15 @@ class MeshSetsGetter:
             )
 
         for femobj in self.member.cons_electrostatic:
-            if femobj["Object"].BoundaryCondition == "Neumann":
-                print_obj_info(femobj["Object"])
+            obj = femobj["Object"]
+            if obj.BoundaryCondition == "Neumann":
+                print_obj_info(obj)
+                result = []
+                ref_data = meshtools.pair_obj_reference(obj.References)
+                for ref_pair in ref_data:
+                    result.append(meshtools.get_ccx_elements(self, ref_pair))
 
-                charged_faces = meshtools.get_charge_density_obj_faces(
-                    self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-                )
-                some_string = "{}: face electric flux".format(femobj["Object"].Name)
-                femobj["ElectricFluxFaces"] = [(some_string, charged_faces)]
-                FreeCAD.Console.PrintLog("{}\n".format(femobj["ElectricFluxFaces"]))
+                femobj["ElectricFluxFaces"] = result
 
     def get_constraints_electricchargedensity_faces(self):
         if not self.member.cons_electricchargedensity:
@@ -416,27 +416,17 @@ class MeshSetsGetter:
             )
 
         for femobj in self.member.cons_electricchargedensity:
-            if femobj["Object"].Mode in ["Interface", "Total Interface"]:
-                print_obj_info(femobj["Object"])
+            obj = femobj["Object"]
+            print_obj_info(obj)
+            result = []
+            ref_data = meshtools.pair_obj_reference(obj.References)
+            for ref_pair in ref_data:
+                result.append(meshtools.get_ccx_elements(self, ref_pair))
 
-                charged_faces = meshtools.get_charge_density_obj_faces(
-                    self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-                )
-                some_string = "{}: face electric charge density".format(femobj["Object"].Name)
-                femobj["ChargeDensityFaces"] = [(some_string, charged_faces)]
-                FreeCAD.Console.PrintLog("{}\n".format(femobj["ChargeDensityFaces"]))
-
-            elif femobj["Object"].Mode in ["Source", "Total Source"]:
-                print_obj_info(femobj["Object"])
-
-                charged_volumes = meshtools.get_charge_density_obj_elements(
-                    self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-                )
-                some_string = "{}: Elements with electric charge density".format(
-                    femobj["Object"].Name
-                )
-                femobj["ChargeDensityElements"] = (some_string, charged_volumes)
-                FreeCAD.Console.PrintLog("{}\n".format(femobj["ChargeDensityElements"]))
+            if obj.Mode in ["Interface", "Total Interface"]:
+                femobj["ChargeDensityFaces"] = result
+            elif obj.Mode in ["Source", "Total Source"]:
+                femobj["ChargeDensityElements"] = result
 
     def get_constraints_contact_faces(self):
         if not self.member.cons_contact:
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py b/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py
index 28903d558e..0b903a70ae 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py
@@ -31,7 +31,7 @@ def get_analysis_types():
 
 
 def get_sets_name():
-    return "constraints_electricchargedensity_node_sets"
+    return "constraints_electricchargedensity_element_sets"
 
 
 def get_constraint_title():
@@ -40,30 +40,15 @@ def get_constraint_title():
 
 def write_meshdata_constraint(f, femobj, den_obj, ccxwriter):
 
-    # print("MESHDATA: ", femobj)
     if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic":
         return
 
     if den_obj.Mode in ["Source", "Total Source"]:
-        msg, data = femobj["ChargeDensityElements"]
-        f.write(f"** {msg}\n")
         f.write(f"*ELSET,ELSET={den_obj.Name}\n")
-        for ref, elem in data:
-            for e in elem:
-                f.write(f"{e},\n")
-
-
-#        if isinstance(femobj["FEMElements"], str):
-#            f.write("{}\n".format(femobj["FEMElements"]))
-#        else:
-#            for e in femobj["FEMElements"]:
-#                f.write(f"{e},\n")
-
-
-#    if femobj["Object"].BoundaryCondition == "Dirichlet":
-#        f.write(f"*NSET,NSET={den_obj.Name}\n")
-#        for n in femobj["Nodes"]:
-#            f.write(f"{n},\n")
+        for refs, surf, is_sub_el in femobj["ChargeDensityElements"]:
+            if not is_sub_el:
+                for elem in surf:
+                    f.write(f"{elem},\n")
 
 
 def get_before_write_meshdata_constraint():
@@ -85,10 +70,6 @@ def get_after_write_constraint():
 def write_constraint(f, femobj, den_obj, ccxwriter):
 
     # floats read from ccx should use {:.13G}, see comment in writer module
-    #    if den_obj.BoundaryCondition == "Dirichlet":
-    #        f.write("*BOUNDARY\n")
-    #        f.write("{},11,11,{:.13G}\n".format(den_obj.Name, den_obj.Potential.getValueAs("mV").Value))
-    #        f.write("\n")
     if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic":
         return
 
@@ -108,18 +89,17 @@ def write_constraint(f, femobj, den_obj, ccxwriter):
         f.write("\n")
 
     elif den_obj.Mode in ["Interface", "Total Interface"]:
-        # check internal interface
         density = density.getValueAs("C/mm^2").Value
+        # check internal interface
         internal = _check_shared_interface(den_obj)
-        for feat, refs in femobj["ChargeDensityFaces"]:
-            f.write("** " + feat + "\n")
+        for feat, surf, is_sub_el in femobj["ChargeDensityFaces"]:
+            f.write("** {0.Name}.{1[0]}\n".format(*feat))
             f.write("*DFLUX\n")
-            for ref in refs:
-                d = density
-                if ref[0] in internal:
-                    d = density / 2
-                for face, fno in ref[1]:
-                    f.write("{},S{},{:.13G}\n".format(face, fno, d))
+            d = density
+            if feat in internal:
+                d = density / 2
+            for face, fno in surf:
+                f.write("{},S{},{:.13G}\n".format(face, fno, d))
 
         f.write("\n")
 
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py b/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py
index 1b5b9ef65f..1f86835dfe 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py
@@ -80,15 +80,14 @@ def write_constraint(f, femobj, pot_obj, ccxwriter):
         density = pot_obj.ElectricFluxDensity.getValueAs("C/mm^2").Value
         # check internal interface
         internal = _check_shared_interface(pot_obj)
-        for feat, refs in femobj["ElectricFluxFaces"]:
-            f.write("** " + feat + "\n")
+        for feat, surf, is_sub_el in femobj["ElectricFluxFaces"]:
+            f.write("** {0.Name}.{1[0]}\n".format(*feat))
             f.write("*DFLUX\n")
-            for ref in refs:
-                d = density
-                if ref[0] in internal:
-                    d = density / 2
-                for face, fno in ref[1]:
-                    f.write("{},S{},{:.13G}\n".format(face, fno, d))
+            d = density
+            if feat in internal:
+                d = density / 2
+            for face, fno in surf:
+                f.write("{},S{},{:.13G}\n".format(face, fno, d))
 
     f.write("\n")