diff --git a/src/Mod/Fem/Gui/TaskFemConstraintContact.cpp b/src/Mod/Fem/Gui/TaskFemConstraintContact.cpp
index 65019b1533..ddd0e47358 100644
--- a/src/Mod/Fem/Gui/TaskFemConstraintContact.cpp
+++ b/src/Mod/Fem/Gui/TaskFemConstraintContact.cpp
@@ -227,8 +227,10 @@ void TaskFemConstraintContact::addToSelectionSlave()
         }
         for (const auto& subName : subNames) {  // for every selected sub element
             bool addMe = true;
-            if (subName.substr(0, 4) != "Face") {
-                QMessageBox::warning(this, tr("Selection error"), tr("Only faces can be picked"));
+            if ((subName.substr(0, 4) != "Face") && (subName.substr(0, 4) != "Edge")) {
+                QMessageBox::warning(this,
+                                     tr("Selection error"),
+                                     tr("Only faces can be picked (edges in 2D models)"));
                 return;
             }
             for (auto itr = std::ranges::find(SubElements, subName); itr != SubElements.end();
@@ -352,8 +354,10 @@ void TaskFemConstraintContact::addToSelectionMaster()
         }
         for (const auto& subName : subNames) {  // for every selected sub element
             bool addMe = true;
-            if (subName.substr(0, 4) != "Face") {
-                QMessageBox::warning(this, tr("Selection error"), tr("Only faces can be picked"));
+            if ((subName.substr(0, 4) != "Face") && (subName.substr(0, 4) != "Edge")) {
+                QMessageBox::warning(this,
+                                     tr("Selection error"),
+                                     tr("Only faces can be picked (edges in 2D models)"));
                 return;
             }
             for (auto itr = std::ranges::find(SubElements.begin(), SubElements.end(), subName);
diff --git a/src/Mod/Fem/Gui/TaskFemConstraintHeatflux.cpp b/src/Mod/Fem/Gui/TaskFemConstraintHeatflux.cpp
index 52dc84abcb..9978a32a3d 100644
--- a/src/Mod/Fem/Gui/TaskFemConstraintHeatflux.cpp
+++ b/src/Mod/Fem/Gui/TaskFemConstraintHeatflux.cpp
@@ -281,10 +281,11 @@ void TaskFemConstraintHeatflux::addToSelection()
 
         if (!subNames.empty()) {
             for (const auto& subName : subNames) {
-                if (subName.substr(0, 4) != "Face") {
-                    QMessageBox::warning(this,
-                                         tr("Selection error"),
-                                         tr("Selection must only consist of faces!"));
+                if ((subName.substr(0, 4) != "Face") && (subName.substr(0, 4) != "Edge")) {
+                    QMessageBox::warning(
+                        this,
+                        tr("Selection error"),
+                        tr("Selection must only consist of faces! (edges in 2D models)"));
                     return;
                 }
             }
@@ -345,10 +346,11 @@ void TaskFemConstraintHeatflux::removeFromSelection()
 
         if (!subNames.empty()) {
             for (const auto& subName : subNames) {
-                if (subName.substr(0, 4) != "Face") {
-                    QMessageBox::warning(this,
-                                         tr("Selection error"),
-                                         tr("Selection must only consist of faces!"));
+                if ((subName.substr(0, 4) != "Face") && (subName.substr(0, 4) != "Edge")) {
+                    QMessageBox::warning(
+                        this,
+                        tr("Selection error"),
+                        tr("Selection must only consist of faces! (edges in 2D models)"));
                     return;
                 }
             }
diff --git a/src/Mod/Fem/Gui/TaskFemConstraintPressure.cpp b/src/Mod/Fem/Gui/TaskFemConstraintPressure.cpp
index 02c9878444..86694f4bfc 100644
--- a/src/Mod/Fem/Gui/TaskFemConstraintPressure.cpp
+++ b/src/Mod/Fem/Gui/TaskFemConstraintPressure.cpp
@@ -146,8 +146,10 @@ void TaskFemConstraintPressure::addToSelection()
 
         for (const auto& subName : subNames) {  // for every selected sub element
             bool addMe = true;
-            if (subName.substr(0, 4) != "Face") {
-                QMessageBox::warning(this, tr("Selection error"), tr("Only faces can be picked"));
+            if ((subName.substr(0, 4) != "Face") && (subName.substr(0, 4) != "Edge")) {
+                QMessageBox::warning(this,
+                                     tr("Selection error"),
+                                     tr("Only faces (edges in 2D models) can be picked"));
                 return;
             }
             for (auto itr = std::ranges::find(SubElements, subName); itr != SubElements.end();
diff --git a/src/Mod/Fem/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py
index e0127662ca..ce88857418 100644
--- a/src/Mod/Fem/femmesh/meshsetsgetter.py
+++ b/src/Mod/Fem/femmesh/meshsetsgetter.py
@@ -101,6 +101,19 @@ class MeshSetsGetter:
         self.femelement_count_test = True
         self.mat_geo_sets = []
 
+        self._load_tables()
+
+    def _load_tables(self):
+        if self.mesh_object:
+            if not self.femnodes_mesh:
+                self.femnodes_mesh = self.femmesh.Nodes
+            if not self.femelement_table:
+                self.femelement_table = meshtools.get_femelement_table(self.femmesh)
+            if not self.femnodes_ele_table:
+                self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
+                    self.femnodes_mesh, self.femelement_table
+                )
+
     # ********************************************************************************************
     # ********************************************************************************************
     # use set for node sets to be sure all nodes are unique
@@ -338,130 +351,48 @@ class MeshSetsGetter:
 
     # ********************************************************************************************
     # ********************************************************************************************
+    def _get_ccx_elements(self, obj):
+        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))
+
+        return result
+
     # faces sets
     def get_constraints_pressure_faces(self):
-        if not self.member.cons_pressure:
-            return
-        # TODO see comments in get_constraints_force_nodeloads()
-        # it applies here too. Mhh it applies to all constraints ...
-
-        """
-        # deprecated version
-        # get the faces and face numbers
         for femobj in self.member.cons_pressure:
-            # femobj --> dict, FreeCAD document object is femobj["Object"]
-            femobj["PressureFaces"] = meshtools.get_pressure_obj_faces_depreciated(
-                self.femmesh,
-                femobj
-            )
-            # print(femobj["PressureFaces"])
-        """
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
 
-        if not self.femnodes_mesh:
-            self.femnodes_mesh = self.femmesh.Nodes
-        if not self.femelement_table:
-            self.femelement_table = meshtools.get_femelement_table(self.femmesh)
-        if not self.femnodes_ele_table:
-            self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
-                self.femnodes_mesh, self.femelement_table
-            )
-
-        for femobj in self.member.cons_pressure:
-            # femobj --> dict, FreeCAD document object is femobj["Object"]
-            print_obj_info(femobj["Object"])
-            pressure_faces = meshtools.get_pressure_obj_faces(
-                self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-            )
-            # the data model is for compatibility reason with deprecated version
-            # get_pressure_obj_faces_depreciated returns the face ids in a tuple per ref_shape
-            # some_string was the reference_shape_element_string in deprecated method
-            # [(some_string, [ele_id, ele_face_id], [ele_id, ele_face_id], ...])]
-            some_string = "{}: face load".format(femobj["Object"].Name)
-            femobj["PressureFaces"] = [(some_string, pressure_faces)]
-            FreeCAD.Console.PrintLog("{}\n".format(femobj["PressureFaces"]))
+            femobj["PressureFaces"] = result
 
     def get_constraints_electrostatic_faces(self):
-        if not self.member.cons_electrostatic:
-            return
-        if not self.femnodes_mesh:
-            self.femnodes_mesh = self.femmesh.Nodes
-        if not self.femelement_table:
-            self.femelement_table = meshtools.get_femelement_table(self.femmesh)
-        if not self.femnodes_ele_table:
-            self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
-                self.femnodes_mesh, self.femelement_table
-            )
-
         for femobj in self.member.cons_electrostatic:
-            if femobj["Object"].BoundaryCondition == "Neumann":
-                print_obj_info(femobj["Object"])
+            obj = femobj["Object"]
+            if obj.BoundaryCondition == "Neumann":
+                result = self._get_ccx_elements(obj)
 
-                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:
-            return
-        if not self.femnodes_mesh:
-            self.femnodes_mesh = self.femmesh.Nodes
-        if not self.femelement_table:
-            self.femelement_table = meshtools.get_femelement_table(self.femmesh)
-        if not self.femnodes_ele_table:
-            self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
-                self.femnodes_mesh, self.femelement_table
-            )
-
         for femobj in self.member.cons_electricchargedensity:
-            if femobj["Object"].Mode in ["Interface", "Total Interface"]:
-                print_obj_info(femobj["Object"])
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
 
-                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:
-            return
-        if not self.femnodes_mesh:
-            self.femnodes_mesh = self.femmesh.Nodes
-        if not self.femelement_table:
-            self.femelement_table = meshtools.get_femelement_table(self.femmesh)
-        if not self.femnodes_ele_table:
-            self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
-                self.femnodes_mesh, self.femelement_table
-            )
-
         for femobj in self.member.cons_contact:
-            # femobj --> dict, FreeCAD document object is femobj["Object"]
-            print_obj_info(femobj["Object"])
-            contact_slave_faces, contact_master_faces = meshtools.get_contact_obj_faces(
-                self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-            )
-            # [ele_id, ele_face_id], [ele_id, ele_face_id], ...]
-            # whereas the ele_face_id might be ccx specific
-            femobj["ContactSlaveFaces"] = contact_slave_faces
-            femobj["ContactMasterFaces"] = contact_master_faces
-            # FreeCAD.Console.PrintLog("{}\n".format(femobj["ContactSlaveFaces"]))
-            # FreeCAD.Console.PrintLog("{}\n".format(femobj["ContactMasterFaces"]))
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
+
+            femobj["ContactSlaveFaces"] = result[:-1]
+            femobj["ContactMasterFaces"] = result[-1:]
 
     # information in the regard of element faces constraints
     # forum post: https://forum.freecad.org/viewtopic.php?f=18&t=42783&p=370286#p366723
@@ -471,29 +402,12 @@ class MeshSetsGetter:
     #                from one side of the geometric face are needed
 
     def get_constraints_tie_faces(self):
-        if not self.member.cons_tie:
-            return
-        if not self.femnodes_mesh:
-            self.femnodes_mesh = self.femmesh.Nodes
-        if not self.femelement_table:
-            self.femelement_table = meshtools.get_femelement_table(self.femmesh)
-        if not self.femnodes_ele_table:
-            self.femnodes_ele_table = meshtools.get_femnodes_ele_table(
-                self.femnodes_mesh, self.femelement_table
-            )
-
         for femobj in self.member.cons_tie:
-            # femobj --> dict, FreeCAD document object is femobj["Object"]
-            print_obj_info(femobj["Object"])
-            slave_faces, master_faces = meshtools.get_tie_obj_faces(
-                self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
-            )
-            # [ele_id, ele_face_id], [ele_id, ele_face_id], ...]
-            # whereas the ele_face_id might be ccx specific
-            femobj["TieSlaveFaces"] = slave_faces
-            femobj["TieMasterFaces"] = master_faces
-            # FreeCAD.Console.PrintLog("{}\n".format(femobj["ContactSlaveFaces"]))
-            # FreeCAD.Console.PrintLog("{}\n".format(femobj["ContactMasterFaces"]))
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
+
+            femobj["TieSlaveFaces"] = result[:-1]
+            femobj["TieMasterFaces"] = result[-1:]
 
     def get_constraints_sectionprint_faces(self):
         if not self.member.cons_sectionprint:
@@ -543,27 +457,11 @@ class MeshSetsGetter:
                         )
 
     def get_constraints_heatflux_faces(self):
-        if not self.member.cons_heatflux:
-            return
-        # TODO: use meshtools to get the surfaces (or move to mesh tools)
-        # see constraint contact or constraint tie and constraint force
-        # heatflux_obj_face_table: see force_obj_node_load_table
-        #     [
-        #         ("refshape_name:elemname", face_table),
-        #         ...,
-        #         ("refshape_name:elemname", face_table)
-        #     ]
         for femobj in self.member.cons_heatflux:
-            # femobj --> dict, FreeCAD document object is femobj["Object"]
-            heatflux_obj = femobj["Object"]
-            femobj["HeatFluxFaceTable"] = []
-            for o, elem_tup in heatflux_obj.References:
-                for elem in elem_tup:
-                    ho = o.Shape.getElement(elem)
-                    if ho.ShapeType == "Face":
-                        elem_info = f"{o.Name}:{elem}"
-                        face_table = self.mesh_object.FemMesh.getccxVolumesByFace(ho)
-                        femobj["HeatFluxFaceTable"].append((elem_info, face_table))
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
+
+            femobj["HeatFluxFaces"] = result
 
     # ********************************************************************************************
     # ********************************************************************************************
@@ -581,16 +479,11 @@ class MeshSetsGetter:
             self.get_solid_element_sets(self.member.cons_centrif)
 
     def get_constraints_bodyheatsource_elements(self):
-        # get element ids and write them into the femobj
-        if not self.member.cons_bodyheatsource:
-            return
-        if (
-            len(self.member.cons_bodyheatsource) == 1
-            and not self.member.cons_bodyheatsource[0]["Object"].References
-        ):
-            self.member.cons_bodyheatsource[0]["FEMElements"] = self.ccx_evolumes
-        else:
-            self.get_solid_element_sets(self.member.cons_bodyheatsource)
+        for femobj in self.member.cons_bodyheatsource:
+            obj = femobj["Object"]
+            result = self._get_ccx_elements(obj)
+
+            femobj["BodyHeatSourceElements"] = result
 
     # ********************************************************************************************
     # ********************************************************************************************
diff --git a/src/Mod/Fem/femmesh/meshtools.py b/src/Mod/Fem/femmesh/meshtools.py
index 0c7919e972..e5b1075194 100644
--- a/src/Mod/Fem/femmesh/meshtools.py
+++ b/src/Mod/Fem/femmesh/meshtools.py
@@ -291,7 +291,9 @@ def get_ccxelement_faces_elements_from_binary_search(bit_pattern_dict):
     return faces
 
 
-def get_ccxelement_edges_from_binary_search(bit_pattern_dict):
+def get_ccxelement_edges_from_binary_search(bit_pattern_dict, sets_getter):
+    shell_mode = sets_getter.solver_obj.ModelSpace == "3D"
+    offset = 2 if shell_mode else 0
     tria3_mask = {0b011: 1, 0b110: 2, 0b101: 3}
     tria6_mask = {0b001011: 1, 0b010110: 2, 0b100101: 3}
     quad4_mask = {0b0011: 1, 0b0110: 2, 0b1100: 3, 0b1001: 4}
@@ -307,7 +309,7 @@ def get_ccxelement_edges_from_binary_search(bit_pattern_dict):
         mask_dict = vol_dict[bit_pattern_dict[ele][0]]
         for key in mask_dict:
             if (key & bit_pattern_dict[ele][1]) == key:
-                faces.append([ele, mask_dict[key]])
+                faces.append([ele, mask_dict[key] + offset])
     # print("EDGES:", faces)
     FreeCAD.Console.PrintMessage(f"found Edges: {len(faces)}\n")
 
@@ -1502,9 +1504,66 @@ def get_ref_shape_node_sum_geom_table(node_geom_table):
 # ************************************************************************************************
 # ***** methods for retrieving element face sets *************************************************
 # ***** charged faces ****************************************************************************
-def get_charge_density_obj_elements(femmesh, femelement_table, femnodes_ele_table, femobj):
+def pair_obj_reference(obj_ref):
+    pairs = []
+    for feat, ref in obj_ref:
+        for sub_ref in ref:
+            sub = (feat, sub_ref)
+            pairs.append(sub)
+
+    return pairs
+
+
+def get_ccx_elements(sets_getter, ref_pair):
+    ref_obj, sub_ref = ref_pair
+    geom_type = ref_obj.getSubObject(sub_ref).ShapeType
+    elem = []
+    is_sub_element = False
+    model_dim = 0
+    if is_solid_femmesh(sets_getter.femmesh):
+        model_dim = 3
+    elif is_face_femmesh(sets_getter.femmesh):
+        model_dim = 2
+    elif is_edge_femmesh(sets_getter.femmesh):
+        model_dim = 1
+
+    match model_dim:
+        case 3:
+            match geom_type:
+                case "Solid":
+                    elem = get_ccx_elements_by_references(sets_getter, ref_pair)
+                    is_sub_element = False
+                case "Face" | "Edge" | "Vertex":
+                    elem = get_ccx_subelements_by_references(sets_getter, ref_pair)
+                    is_sub_element = True
+        case 2:
+            match geom_type:
+                case "Face":
+                    elem = get_ccx_elements_by_references(sets_getter, ref_pair)
+                    is_sub_element = False
+                case "Edge" | "Vertex":
+                    elem = get_ccx_subelements_by_references(sets_getter, ref_pair)
+                    is_sub_element = True
+        case 1:
+            match geom_type:
+                case "Edge":
+                    is_sub_element = False
+                    elem = get_ccx_elements_by_references(sets_getter, ref_pair)
+                case "Vertex":
+                    elem = get_ccx_subelements_by_references(sets_getter, ref_pair)
+                    is_sub_element = True
+        case 0:
+            match geom_type:
+                case "Vertex":
+                    elem = get_ccx_elements_by_references(sets_getter, ref_pair)
+                    is_sub_element = False
+
+    return (*elem, is_sub_element)
+
+
+def get_ccx_elements_by_references(sets_getter, femobj_ref):
     node_set = []
-    res = []
+    result = []
     # TODO get elements from mesh groups
     # if femmesh.GroupCount:
     #     node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, "Node")
@@ -1515,36 +1574,33 @@ def get_charge_density_obj_elements(femmesh, femelement_table, femnodes_ele_tabl
     #             "from existent finite element mesh group data.\n"
     #         )
     if not node_set:
+        elem = []
         FreeCAD.Console.PrintLog(
             "    Finite element mesh nodes will be retrieved "
             "by searching the appropriate nodes in the finite element mesh.\n"
         )
-        for feat, ref in femobj["Object"].References:
-            for sub_ref in ref:
-                sub = (feat, (sub_ref,))
-                node_set = get_femnodes_by_references(femmesh, [sub])
-                charged_volume_node_set = sorted(set(node_set))
+        feat, sub_ref = femobj_ref
+        sub = (feat, (sub_ref,))
+        node_set = get_femnodes_by_references(sets_getter.femmesh, [sub])
+        charged_volume_node_set = sorted(set(node_set))
 
-                bit_pattern_dict = get_bit_pattern_dict(
-                    femelement_table, femnodes_ele_table, charged_volume_node_set
-                )
-                sh = feat.getSubObject(sub_ref)
-                if sh.ShapeType == "Solid":
-                    charged_elem = get_ccxelement_volumes_elements_from_binary_search(
-                        bit_pattern_dict
-                    )
-                elif sh.ShapeType == "Face":
-                    charged_elem = get_ccxelement_faces_elements_from_binary_search(
-                        bit_pattern_dict
-                    )
-                res.append((sub, charged_elem))
+        bit_pattern_dict = get_bit_pattern_dict(
+            sets_getter.femelement_table, sets_getter.femnodes_ele_table, charged_volume_node_set
+        )
+        sh = feat.getSubObject(sub_ref)
+        if sh.ShapeType == "Solid":
+            elem = get_ccxelement_volumes_elements_from_binary_search(bit_pattern_dict)
+        elif sh.ShapeType == "Face":
+            elem = get_ccxelement_faces_elements_from_binary_search(bit_pattern_dict)
 
-    return res
+        result = (sub, elem)
+
+    return result
 
 
-def get_charge_density_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj):
+def get_ccx_subelements_by_references(sets_getter, femobj_ref):
     node_set = []
-    res = []
+    result = []
     # TODO get elements from mesh groups
     # if femmesh.GroupCount:
     #     node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, "Node")
@@ -1555,28 +1611,28 @@ def get_charge_density_obj_faces(femmesh, femelement_table, femnodes_ele_table,
     #             "from existent finite element mesh group data.\n"
     #         )
     if not node_set:
+        sub_elem = []
         FreeCAD.Console.PrintLog(
             "    Finite element mesh nodes will be retrieved "
             "by searching the appropriate nodes in the finite element mesh.\n"
         )
-        for feat, ref in femobj["Object"].References:
-            for sub_ref in ref:
-                sub = (feat, (sub_ref,))
-                node_set = get_femnodes_by_references(femmesh, [sub])
-                charged_face_node_set = sorted(set(node_set))
+        feat, sub_ref = femobj_ref
+        sub = (feat, (sub_ref,))
+        node_set = get_femnodes_by_references(sets_getter.femmesh, [sub])
+        charged_face_node_set = sorted(set(node_set))
 
-                bit_pattern_dict = get_bit_pattern_dict(
-                    femelement_table, femnodes_ele_table, charged_face_node_set
-                )
-                sh = feat.getSubObject(sub_ref)
-                if sh.ShapeType == "Face":
-                    charged_faces = get_ccxelement_faces_from_binary_search(bit_pattern_dict)
-                elif sh.ShapeType == "Edge":
-                    charged_faces = get_ccxelement_edges_from_binary_search(bit_pattern_dict)
+        bit_pattern_dict = get_bit_pattern_dict(
+            sets_getter.femelement_table, sets_getter.femnodes_ele_table, charged_face_node_set
+        )
+        sh = feat.getSubObject(sub_ref)
+        if sh.ShapeType == "Face":
+            sub_elem = get_ccxelement_faces_from_binary_search(bit_pattern_dict)
+        elif sh.ShapeType == "Edge":
+            sub_elem = get_ccxelement_edges_from_binary_search(bit_pattern_dict, sets_getter)
 
-                res.append((sub, charged_faces))
+        result = (sub, sub_elem)
 
-    return res
+    return result
 
 
 # ***** pressure faces ***************************************************************************
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_bodyheatsource.py b/src/Mod/Fem/femsolver/calculix/write_constraint_bodyheatsource.py
index 65d48a7c24..35fde847d1 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_bodyheatsource.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_bodyheatsource.py
@@ -58,12 +58,12 @@ def get_after_write_constraint():
 
 
 def write_meshdata_constraint(f, femobj, bodyheatsource_obj, ccxwriter):
+
     f.write(f"*ELSET,ELSET={bodyheatsource_obj.Name}\n")
-    if isinstance(femobj["FEMElements"], str):
-        f.write("{}\n".format(femobj["FEMElements"]))
-    else:
-        for e in femobj["FEMElements"]:
-            f.write(f"{e},\n")
+    for refs, surf, is_sub_el in femobj["BodyHeatSourceElements"]:
+        if not is_sub_el:
+            for elem in surf:
+                f.write(f"{elem},\n")
 
 
 def write_constraint(f, femobj, bodyheatsource_obj, ccxwriter):
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_contact.py b/src/Mod/Fem/femsolver/calculix/write_constraint_contact.py
index f602af6ba6..c46873c351 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_contact.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_contact.py
@@ -57,12 +57,23 @@ def get_after_write_constraint():
 def write_meshdata_constraint(f, femobj, contact_obj, ccxwriter):
     # slave DEP
     f.write(f"*SURFACE, NAME=DEP{contact_obj.Name}\n")
-    for i in femobj["ContactSlaveFaces"]:
-        f.write(f"{i[0]},S{i[1]}\n")
+    for refs, surf, is_sub_el in femobj["ContactSlaveFaces"]:
+        if is_sub_el:
+            for elem, fno in surf:
+                f.write(f"{elem},S{fno}\n")
+        else:
+            for elem in surf:
+                f.write(f"{elem},S2\n")
+
     # master IND
     f.write(f"*SURFACE, NAME=IND{contact_obj.Name}\n")
-    for i in femobj["ContactMasterFaces"]:
-        f.write(f"{i[0]},S{i[1]}\n")
+    for refs, surf, is_sub_el in femobj["ContactMasterFaces"]:
+        if is_sub_el:
+            for elem, fno in surf:
+                f.write(f"{elem},S{fno}\n")
+        else:
+            for elem in surf:
+                f.write(f"{elem},S2\n")
 
 
 def write_constraint(f, femobj, contact_obj, ccxwriter):
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")
 
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_heatflux.py b/src/Mod/Fem/femsolver/calculix/write_constraint_heatflux.py
index 057b71c1e6..3925b0cc89 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_heatflux.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_heatflux.py
@@ -82,10 +82,7 @@ def write_meshdata_constraint(f, femobj, heatflux_obj, ccxwriter):
         heatflux_amplitude = ""
 
     f.write(f"*{heatflux_key_word}{heatflux_amplitude}\n")
-    for ref_shape in femobj["HeatFluxFaceTable"]:
-        elem_string = ref_shape[0]
-        face_table = ref_shape[1]
-        f.write(f"** Heat flux on face {elem_string}\n")
-        for i in face_table:
-            # OvG: Only write out the VolumeIDs linked to a particular face
-            f.write(f"{i[0]},{heatflux_facetype}{i[1]}{heatflux_facesubtype},{heatflux_values}\n")
+    for feat, surf, is_sub_el in femobj["HeatFluxFaces"]:
+        f.write("** {0.Name}.{1[0]}\n".format(*feat))
+        for face, fno in surf:
+            f.write(f"{face},{heatflux_facetype}{fno}{heatflux_facesubtype},{heatflux_values}\n")
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_pressure.py b/src/Mod/Fem/femsolver/calculix/write_constraint_pressure.py
index d389da8099..80ab6187c8 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_pressure.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_pressure.py
@@ -54,21 +54,13 @@ def write_meshdata_constraint(f, femobj, prs_obj, ccxwriter):
         f.write("*DLOAD\n")
     rev = -1 if prs_obj.Reversed else 1
     # the pressure has to be output in MPa
-    pressure_quantity = FreeCAD.Units.Quantity(prs_obj.Pressure.getValueAs("MPa"))
-    press_rev = rev * pressure_quantity
-    for ref_shape in femobj["PressureFaces"]:
-        # the loop is needed for compatibility reason
-        # in deprecated method get_pressure_obj_faces_depreciated
-        # the face ids where per ref_shape
-        f.write("** " + ref_shape[0] + "\n")
-        for face, fno in ref_shape[1]:
-            if fno > 0:  # solid mesh face
-                f.write(f"{face},P{fno},{press_rev}\n")
-            # on shell mesh face: fno == 0
-            # normal of element face == face normal
-            elif fno == 0:
-                f.write(f"{face},P,{press_rev}\n")
-            # on shell mesh face: fno == -1
-            # normal of element face opposite direction face normal
-            elif fno == -1:
-                f.write(f"{face},P,{-1 * press_rev}\n")
+    pressure = prs_obj.Pressure.getValueAs("MPa").Value
+    pressure *= rev
+    for feat, surf, is_sub_el in femobj["PressureFaces"]:
+        f.write("** {0.Name}.{1[0]}\n".format(*feat))
+        if is_sub_el:
+            for elem, fno in surf:
+                f.write(f"{elem},P{fno},{pressure}\n")
+        else:
+            for elem in surf:
+                f.write(f"{elem},P,{-1*pressure}\n")
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_tie.py b/src/Mod/Fem/femsolver/calculix/write_constraint_tie.py
index 4f96280f0e..70d7e73070 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_tie.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_tie.py
@@ -60,12 +60,23 @@ def get_after_write_constraint():
 def write_meshdata_constraint(f, femobj, tie_obj, ccxwriter):
     # slave DEP
     f.write(f"*SURFACE, NAME=TIE_DEP{tie_obj.Name}\n")
-    for i in femobj["TieSlaveFaces"]:
-        f.write(f"{i[0]},S{i[1]}\n")
+    for refs, surf, is_sub_el in femobj["TieSlaveFaces"]:
+        if is_sub_el:
+            for elem, fno in surf:
+                f.write(f"{elem},S{fno}\n")
+        else:
+            for elem in surf:
+                f.write(f"{elem},S2\n")
+
     # master IND
     f.write(f"*SURFACE, NAME=TIE_IND{tie_obj.Name}\n")
-    for i in femobj["TieMasterFaces"]:
-        f.write(f"{i[0]},S{i[1]}\n")
+    for refs, surf, is_sub_el in femobj["TieMasterFaces"]:
+        if is_sub_el:
+            for elem, fno in surf:
+                f.write(f"{elem},S{fno}\n")
+        else:
+            for elem in surf:
+                f.write(f"{elem},S2\n")
 
 
 def write_constraint(f, femobj, tie_obj, ccxwriter):
diff --git a/src/Mod/Fem/femtaskpanels/task_constraint_tie.py b/src/Mod/Fem/femtaskpanels/task_constraint_tie.py
index 3b1cdfbdb0..5f3e633dd7 100644
--- a/src/Mod/Fem/femtaskpanels/task_constraint_tie.py
+++ b/src/Mod/Fem/femtaskpanels/task_constraint_tie.py
@@ -63,7 +63,7 @@ class _TaskPanel(base_femtaskpanel._BaseTaskPanel):
 
         # geometry selection widget
         self.selectionWidget = selection_widgets.GeometryElementsSelection(
-            obj.References, ["Face"], False, False
+            obj.References, ["Edge", "Face"], False, False
         )
 
         # form made from param and selection widget
diff --git a/src/Mod/Fem/femtest/data/calculix/box_static.inp b/src/Mod/Fem/femtest/data/calculix/box_static.inp
index 26f6972664..b6e81fab79 100644
--- a/src/Mod/Fem/femtest/data/calculix/box_static.inp
+++ b/src/Mod/Fem/femtest/data/calculix/box_static.inp
@@ -558,7 +558,7 @@ FemConstraintFixed,3
 ** constraints pressure element face loads
 ** FemConstraintPressure
 *DLOAD
-** FemConstraintPressure: face load
+** Box.Face2
 3,P1,1000.0
 6,P3,1000.0
 9,P2,1000.0
diff --git a/src/Mod/Fem/femtest/data/calculix/constraint_sectionprint.inp b/src/Mod/Fem/femtest/data/calculix/constraint_sectionprint.inp
index 4c195149ef..72e170fe3a 100644
--- a/src/Mod/Fem/femtest/data/calculix/constraint_sectionprint.inp
+++ b/src/Mod/Fem/femtest/data/calculix/constraint_sectionprint.inp
@@ -3427,7 +3427,7 @@ SOF, SOM, SOAREA
 ** constraints pressure element face loads
 ** ConstraintPressure
 *DLOAD
-** ConstraintPressure: face load
+** BooleanFragments.Face1
 846,P3,100.0
 861,P1,100.0
 875,P4,100.0
diff --git a/src/Mod/Fem/femtest/data/calculix/constraint_transform_beam_hinged.inp b/src/Mod/Fem/femtest/data/calculix/constraint_transform_beam_hinged.inp
index d8a738e1b5..2621f8463e 100644
--- a/src/Mod/Fem/femtest/data/calculix/constraint_transform_beam_hinged.inp
+++ b/src/Mod/Fem/femtest/data/calculix/constraint_transform_beam_hinged.inp
@@ -3657,7 +3657,7 @@ FemConstraintDisplacement,1,1,0.0
 ** constraints pressure element face loads
 ** FemConstraintPressure
 *DLOAD
-** FemConstraintPressure: face load
+** CompoundFilter.Face8
 1384,P3,10.0
 1394,P2,10.0
 1397,P2,10.0
diff --git a/src/Mod/Fem/femtest/data/calculix/material_multiple_tensionrod_twoboxes.inp b/src/Mod/Fem/femtest/data/calculix/material_multiple_tensionrod_twoboxes.inp
index 17882d7c3a..757454bd73 100644
--- a/src/Mod/Fem/femtest/data/calculix/material_multiple_tensionrod_twoboxes.inp
+++ b/src/Mod/Fem/femtest/data/calculix/material_multiple_tensionrod_twoboxes.inp
@@ -1251,7 +1251,7 @@ ConstraintFixed,3
 ** constraints pressure element face loads
 ** ConstraintPressure
 *DLOAD
-** ConstraintPressure: face load
+** MultiMatCompSolid.Face11
 61,P3,1000.0
 74,P3,1000.0
 75,P3,1000.0
diff --git a/src/Mod/Fem/femtest/data/calculix/material_nonlinear.inp b/src/Mod/Fem/femtest/data/calculix/material_nonlinear.inp
index 9ee5e70a93..c68f03634c 100644
--- a/src/Mod/Fem/femtest/data/calculix/material_nonlinear.inp
+++ b/src/Mod/Fem/femtest/data/calculix/material_nonlinear.inp
@@ -20024,7 +20024,7 @@ ConstraintFixed,3
 ** constraints pressure element face loads
 ** ConstraintPressure
 *DLOAD
-** ConstraintPressure: face load
+** Hole_Plate.Face2
 4796,P3,-130.0
 4799,P4,-130.0
 4872,P4,-130.0