diff --git a/src/Mod/Fem/App/FemConstraintInitialTemperature.cpp b/src/Mod/Fem/App/FemConstraintInitialTemperature.cpp
index 0a6420a0a8..f4e4b7ad0b 100644
--- a/src/Mod/Fem/App/FemConstraintInitialTemperature.cpp
+++ b/src/Mod/Fem/App/FemConstraintInitialTemperature.cpp
@@ -35,9 +35,6 @@ PROPERTY_SOURCE(Fem::ConstraintInitialTemperature, Fem::Constraint)
 ConstraintInitialTemperature::ConstraintInitialTemperature()
 {
     ADD_PROPERTY(initialTemperature, (300.0));
-
-    References.setStatus(App::Property::ReadOnly, true);
-    References.setStatus(App::Property::Hidden, true);
 }
 
 App::DocumentObjectExecReturn* ConstraintInitialTemperature::execute()
diff --git a/src/Mod/Fem/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py
index 1b267eb540..e0127662ca 100644
--- a/src/Mod/Fem/femmesh/meshsetsgetter.py
+++ b/src/Mod/Fem/femmesh/meshsetsgetter.py
@@ -144,6 +144,7 @@ class MeshSetsGetter:
         self.get_constraints_sectionprint_faces()
         self.get_constraints_transform_nodes()
         self.get_constraints_temperature_nodes()
+        self.get_constraints_initialtemperature_nodes()
         self.get_constraints_electrostatic_nodes()
 
         # constraints sets with constraint data
@@ -245,6 +246,15 @@ class MeshSetsGetter:
             print_obj_info(femobj["Object"])
             femobj["Nodes"] = meshtools.get_femnodes_by_femobj_with_references(self.femmesh, femobj)
 
+    def get_constraints_initialtemperature_nodes(self):
+        if not self.member.cons_initialtemperature:
+            return
+        # get nodes
+        for femobj in self.member.cons_initialtemperature:
+            # femobj --> dict, FreeCAD document object is femobj["Object"]
+            print_obj_info(femobj["Object"])
+            femobj["Nodes"] = meshtools.get_femnodes_by_femobj_with_references(self.femmesh, femobj)
+
     def get_constraints_fluidsection_nodes(self):
         if not self.member.geos_fluidsection:
             return
diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_initialtemperature.py b/src/Mod/Fem/femsolver/calculix/write_constraint_initialtemperature.py
index c3bf3e5182..5a5b02add0 100644
--- a/src/Mod/Fem/femsolver/calculix/write_constraint_initialtemperature.py
+++ b/src/Mod/Fem/femsolver/calculix/write_constraint_initialtemperature.py
@@ -32,10 +32,31 @@ def get_analysis_types():
     return ["thermomech"]
 
 
+def get_sets_name():
+    return "constraints_initial_temperature_node_sets"
+
+
 def get_constraint_title():
     return "Initial temperature constraint"
 
 
+def write_meshdata_constraint(f, femobj, inittemp_obj, ccxwriter):
+    if inittemp_obj.References and len(inittemp_obj.References) > 0:
+        f.write(f"*NSET,NSET={inittemp_obj.Name}\n")
+        for n in femobj["Nodes"]:
+            f.write(f"{n},\n")
+    else:
+        return
+
+
+def get_before_write_meshdata_constraint():
+    return ""
+
+
+def get_after_write_meshdata_constraint():
+    return ""
+
+
 def get_before_write_constraint():
     return "*INITIAL CONDITIONS,TYPE=TEMPERATURE\n"
 
@@ -48,11 +69,9 @@ def write_constraint(f, femobj, inittemp_obj, ccxwriter):
 
     # floats read from ccx should use {:.13G}, see comment in writer module
 
-    f.write(
-        "{},{}\n".format(
-            ccxwriter.ccx_nall, Units.Quantity(inittemp_obj.initialTemperature.getValueAs("K"))
-        )
-    )
+    initialtemp = inittemp_obj.initialTemperature.getValueAs("K")
 
-
-# Should only be one object in the analysis
+    if inittemp_obj.References and len(inittemp_obj.References) > 0:
+        f.write(f"{inittemp_obj.Name},{initialtemp}\n")
+    else:
+        f.write(f"{ccxwriter.ccx_nall},{initialtemp}\n")
diff --git a/src/Mod/Fem/femsolver/calculix/writer.py b/src/Mod/Fem/femsolver/calculix/writer.py
index 8f830bf159..5d1e477885 100644
--- a/src/Mod/Fem/femsolver/calculix/writer.py
+++ b/src/Mod/Fem/femsolver/calculix/writer.py
@@ -159,6 +159,7 @@ class FemInputWriterCcx(writerbase.FemInputWriter):
         self.write_constraints_meshsets(inpfile, self.member.cons_planerotation, con_planerotation)
         self.write_constraints_meshsets(inpfile, self.member.cons_transform, con_transform)
         self.write_constraints_meshsets(inpfile, self.member.cons_temperature, con_temperature)
+        self.write_constraints_meshsets(inpfile, self.member.cons_initialtemperature, con_itemp)
         self.write_constraints_meshsets(
             inpfile, self.member.cons_electricchargedensity, con_electricchargedensity
         )
diff --git a/src/Mod/Fem/femtest/data/calculix/thermomech_bimetal.inp b/src/Mod/Fem/femtest/data/calculix/thermomech_bimetal.inp
index 1267d1aa0b..370b8d3712 100644
--- a/src/Mod/Fem/femtest/data/calculix/thermomech_bimetal.inp
+++ b/src/Mod/Fem/femtest/data/calculix/thermomech_bimetal.inp
@@ -7040,6 +7040,10 @@ Evolumes
 1520,
 1521,
 
+***********************************************************
+** constraints initial temperature node sets
+** ConstraintInitialTemperature
+
 ** Physical constants for SI(mm) unit system with Kelvins
 *PHYSICAL CONSTANTS, ABSOLUTE ZERO=0, STEFAN BOLTZMANN=5.670374419e-11
 
diff --git a/src/Mod/Fem/femtools/checksanalysis.py b/src/Mod/Fem/femtools/checksanalysis.py
index b6f67f5c1a..ea457348a7 100644
--- a/src/Mod/Fem/femtools/checksanalysis.py
+++ b/src/Mod/Fem/femtools/checksanalysis.py
@@ -226,8 +226,6 @@ def check_member_for_solver_calculix(analysis, solver, mesh, member):
         if not member.cons_initialtemperature:
             if not member.geos_fluidsection:
                 message += "Thermomechanical analysis: No initial temperature defined.\n"
-        if len(member.cons_initialtemperature) > 1:
-            message += "Thermomechanical analysis: Only one initial temperature is allowed.\n"
 
     # constraints
     # fixed