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/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py
index 2d18f57f99..9a7ecbae2e 100644
--- a/src/Mod/Fem/femmesh/meshsetsgetter.py
+++ b/src/Mod/Fem/femmesh/meshsetsgetter.py
@@ -533,25 +533,24 @@ 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)
-        #     ]
+        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_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"]
+            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["HeatFluxFaces"] = result
 
     # ********************************************************************************************
     # ********************************************************************************************
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")