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 9a7ecbae2e..e063218a6d 100644
--- a/src/Mod/Fem/femmesh/meshsetsgetter.py
+++ b/src/Mod/Fem/femmesh/meshsetsgetter.py
@@ -342,21 +342,6 @@ class MeshSetsGetter:
     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"])
-        """
-
         if not self.femnodes_mesh:
             self.femnodes_mesh = self.femmesh.Nodes
         if not self.femelement_table:
@@ -368,17 +353,14 @@ class MeshSetsGetter:
 
         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"]))
+            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))
+
+            femobj["PressureFaces"] = result
 
     def get_constraints_electrostatic_faces(self):
         if not self.member.cons_electrostatic:
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")