diff --git a/src/Mod/Fem/femsolver/calculix/writer.py b/src/Mod/Fem/femsolver/calculix/writer.py
index a2378c748f..cc4193ec67 100644
--- a/src/Mod/Fem/femsolver/calculix/writer.py
+++ b/src/Mod/Fem/femsolver/calculix/writer.py
@@ -895,40 +895,33 @@ class FemInputWriterCcx(writerbase.FemInputWriter):
             caller_method_name=sys._getframe().f_code.co_name,
         )
 
-    # TODO move code parts from this method to base writer module
-    # into get_constraints_heatflux_faces method
     def write_faceheatflux_constraints_heatflux(self, f, femobj, heatflux_obj):
         if heatflux_obj.ConstraintType == "Convection":
-            f.write("*FILM\n")
-            for o, elem_tup in heatflux_obj.References:
-                for elem in elem_tup:
-                    ho = o.Shape.getElement(elem)
-                    if ho.ShapeType == "Face":
-                        v = self.mesh_object.FemMesh.getccxVolumesByFace(ho)
-                        f.write("** Heat flux on face {}\n".format(elem))
-                        for i in v:
-                            # SvdW: add factor to force heatflux to units system of t/mm/s/K
-                            # OvG: Only write out the VolumeIDs linked to a particular face
-                            f.write("{},F{},{},{}\n".format(
-                                i[0],
-                                i[1],
-                                heatflux_obj.AmbientTemp,
-                                heatflux_obj.FilmCoef * 0.001
-                            ))
+            heatflux_key_word = "FILM"
+            heatflux_facetype = "F"
+            # SvdW: add factor to force heatflux to units system of t/mm/s/K
+            heatflux_values = "{},{}".format(
+                heatflux_obj.AmbientTemp,
+                heatflux_obj.FilmCoef * 0.001
+            )
         elif heatflux_obj.ConstraintType == "DFlux":
-            f.write("*DFLUX\n")
-            for o, elem_tup in heatflux_obj.References:
-                for elem in elem_tup:
-                    ho = o.Shape.getElement(elem)
-                    if ho.ShapeType == "Face":
-                        v = self.mesh_object.FemMesh.getccxVolumesByFace(ho)
-                        f.write("** Heat flux on face {}\n".format(elem))
-                        for i in v:
-                            f.write("{},S{},{}\n".format(
-                                i[0],
-                                i[1],
-                                heatflux_obj.DFlux * 0.001
-                            ))
+            heatflux_key_word = "DFLUX"
+            heatflux_facetype = "S"
+            heatflux_values = "{}".format(heatflux_obj.DFlux * 0.001)
+
+        f.write("*{}\n".format(heatflux_key_word))
+        for ref_shape in femobj["HeatFluxFaceTable"]:
+            elem_string = ref_shape[0]
+            face_table = ref_shape[1]
+            f.write("** Heat flux on face {}\n".format(elem_string))
+            for i in face_table:
+                # OvG: Only write out the VolumeIDs linked to a particular face
+                f.write("{},{}{},{}\n".format(
+                    i[0],
+                    heatflux_facetype,
+                    i[1],
+                    heatflux_values
+                ))
 
     # ********************************************************************************************
     # constraints fluidsection
diff --git a/src/Mod/Fem/femsolver/writerbase.py b/src/Mod/Fem/femsolver/writerbase.py
index 538c69662f..1c7bb88fe6 100644
--- a/src/Mod/Fem/femsolver/writerbase.py
+++ b/src/Mod/Fem/femsolver/writerbase.py
@@ -436,8 +436,25 @@ class FemInputWriter():
                         )
 
     def get_constraints_heatflux_faces(self):
-        pass
-        # TODO implement
+        # TODO: use meshtools to get the surfaces (or move to mesh tools)
+        # see constraint contact or constrint 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.heatflux_objects:
+            # 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 = "{}:{}".format(o.Name, elem)
+                        face_table = self.mesh_object.FemMesh.getccxVolumesByFace(ho)
+                        femobj["HeatFluxFaceTable"].append((elem_info, face_table))
 
     # ********************************************************************************************
     # ********************************************************************************************
diff --git a/src/Mod/Fem/femtest/data/calculix/thermomech_spine.inp b/src/Mod/Fem/femtest/data/calculix/thermomech_spine.inp
index 9e03315ebe..c3510dee50 100644
--- a/src/Mod/Fem/femtest/data/calculix/thermomech_spine.inp
+++ b/src/Mod/Fem/femtest/data/calculix/thermomech_spine.inp
@@ -169,22 +169,22 @@ FemConstraintTemperature,11,11,310.93
 ** written by write_constraints_heatflux function
 ** FemConstraintHeatflux
 *FILM
-** Heat flux on face Face3
+** Heat flux on face Box:Face3
 3,F2,255.3722,0.005678
 7,F3,255.3722,0.005678
 11,F2,255.3722,0.005678
 12,F2,255.3722,0.005678
-** Heat flux on face Face4
+** Heat flux on face Box:Face4
 1,F2,255.3722,0.005678
 2,F4,255.3722,0.005678
 4,F4,255.3722,0.005678
 9,F4,255.3722,0.005678
-** Heat flux on face Face5
+** Heat flux on face Box:Face5
 1,F4,255.3722,0.005678
 3,F4,255.3722,0.005678
 4,F3,255.3722,0.005678
 10,F3,255.3722,0.005678
-** Heat flux on face Face6
+** Heat flux on face Box:Face6
 6,F4,255.3722,0.005678
 7,F4,255.3722,0.005678
 8,F4,255.3722,0.005678