From 8920d76d9393f283e5ca43f30914bc1ffa57413f Mon Sep 17 00:00:00 2001
From: marioalexis <mario.passaglia@gmail.com>
Date: Wed, 15 May 2024 01:20:59 -0300
Subject: [PATCH] Fem: Load heat flux from CalculiX results - fixes #12117

---
 src/Mod/Fem/App/FemVTKTools.cpp               |  1 +
 src/Mod/Fem/feminout/importCcxFrdResults.py   | 21 +++++++++++++++++++
 src/Mod/Fem/feminout/importToolsFem.py        |  5 +++++
 src/Mod/Fem/femobjects/result_mechanical.py   | 10 +++++++++
 .../femsolver/calculix/write_step_output.py   |  9 +++++---
 .../data/calculix/thermomech_bimetall.inp     |  2 +-
 6 files changed, 44 insertions(+), 4 deletions(-)

diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp
index fd5df42bc9..fbb6175913 100644
--- a/src/Mod/Fem/App/FemVTKTools.cpp
+++ b/src/Mod/Fem/App/FemVTKTools.cpp
@@ -659,6 +659,7 @@ std::map<std::string, std::string> _getFreeCADMechResultVectorProperties()
     resFCVecProp["PS1Vector"] = "Major Principal Stress Vector";
     resFCVecProp["PS2Vector"] = "Intermediate Principal Stress Vector";
     resFCVecProp["PS3Vector"] = "Minor Principal Stress Vector";
+    resFCVecProp["HeatFlux"] = "Heat Flux";
 
     return resFCVecProp;
 }
diff --git a/src/Mod/Fem/feminout/importCcxFrdResults.py b/src/Mod/Fem/feminout/importCcxFrdResults.py
index e2d4f63699..effc473f87 100644
--- a/src/Mod/Fem/feminout/importCcxFrdResults.py
+++ b/src/Mod/Fem/feminout/importCcxFrdResults.py
@@ -311,6 +311,7 @@ def read_frd_result(
     mode_strain = {}
     mode_peeq = {}
     mode_temp = {}
+    mode_heatflux = {}
     mode_massflow = {}
     mode_networkpressure = {}
 
@@ -322,6 +323,7 @@ def read_frd_result(
     mode_strain_found = False
     mode_peeq_found = False
     mode_temp_found = False
+    mode_heatflux_found = False
     mode_massflow_found = False
     mode_networkpressure_found = False
     end_of_section_found = False
@@ -641,6 +643,19 @@ def read_frd_result(
             temperature = float(line[13:25])
             mode_temp[elem] = (temperature)
 
+        # Check if we found heat flux section
+        if line[5:9] == "FLUX":
+            mode_heatflux_found = True
+        if mode_heatflux_found and (line[1:3] == "-1"):
+            # we found a heat_flux line
+            elem = int(line[4:13])
+            mode_heatflux_x = float(line[13:25])
+            mode_heatflux_y = float(line[25:37])
+            mode_heatflux_z = float(line[37:49])
+            mode_heatflux[elem] = FreeCAD.Vector(mode_heatflux_x, mode_heatflux_y, mode_heatflux_z)
+
+
+
         # Check if we found a mass flow section
         if line[5:11] == "MAFLOW":
             mode_massflow_found = True
@@ -713,6 +728,12 @@ def read_frd_result(
                 mode_temp_found = False
                 node_element_section = False
 
+            if mode_heatflux_found:
+                mode_results["heatflux"] = mode_heatflux
+                mode_heatflux = {}
+                mode_heatflux_found = False
+                node_element_section = False
+
             if mode_massflow_found:
                 mode_results["mflow"] = mode_massflow
                 mode_massflow = {}
diff --git a/src/Mod/Fem/feminout/importToolsFem.py b/src/Mod/Fem/feminout/importToolsFem.py
index 09d39ca32d..2c0aa9c671 100644
--- a/src/Mod/Fem/feminout/importToolsFem.py
+++ b/src/Mod/Fem/feminout/importToolsFem.py
@@ -455,6 +455,11 @@ def fill_femresult_mechanical(
                 res_obj.Temperature = list(map((lambda x: x), Temperature.values()))
             res_obj.Time = step_time
 
+    if "heatflux" in result_set:
+        HeatFlux = result_set["heatflux"]
+        if HeatFlux:
+            res_obj.HeatFlux = list(map((lambda x: x), HeatFlux.values()))
+
     # fill res_obj.MassFlow
     if "mflow" in result_set:
         MassFlow = result_set["mflow"]
diff --git a/src/Mod/Fem/femobjects/result_mechanical.py b/src/Mod/Fem/femobjects/result_mechanical.py
index d23cf72d87..6887b09679 100644
--- a/src/Mod/Fem/femobjects/result_mechanical.py
+++ b/src/Mod/Fem/femobjects/result_mechanical.py
@@ -230,6 +230,16 @@ class ResultMechanical(base_fempythonobject.BaseFemPythonObject):
             "Temperature field",
             True
         )
+        obj.addProperty(
+            "App::PropertyVectorList",
+            "HeatFlux",
+            "NodeData",
+            "List of heat flux vectors",
+            True
+        )
+        obj.setPropertyStatus("HeatFlux", "LockDynamic")
+
+
         obj.setPropertyStatus("Temperature", "LockDynamic")
         obj.addProperty(
             "App::PropertyFloatList",
diff --git a/src/Mod/Fem/femsolver/calculix/write_step_output.py b/src/Mod/Fem/femsolver/calculix/write_step_output.py
index e03bb9293c..2ad190e3cc 100644
--- a/src/Mod/Fem/femsolver/calculix/write_step_output.py
+++ b/src/Mod/Fem/femsolver/calculix/write_step_output.py
@@ -51,10 +51,13 @@ def write_step_output(f, ccxwriter):
         f.write("U\n")
     if not ccxwriter.member.geos_fluidsection:
         f.write("*EL FILE\n")
+        variables = "S, E"
+        if ccxwriter.analysis_type == "thermomech":
+            variables += ", HFL"
         if ccxwriter.solver_obj.MaterialNonlinearity == "nonlinear":
-            f.write("S, E, PEEQ\n")
-        else:
-            f.write("S, E\n")
+            variables += ", PEEQ"
+
+        f.write(variables + "\n")
 
         # dat file
         # reaction forces: freecad.org/tracker/view.php?id=2934
diff --git a/src/Mod/Fem/femtest/data/calculix/thermomech_bimetall.inp b/src/Mod/Fem/femtest/data/calculix/thermomech_bimetall.inp
index 5954c13b35..6f638b535b 100644
--- a/src/Mod/Fem/femtest/data/calculix/thermomech_bimetall.inp
+++ b/src/Mod/Fem/femtest/data/calculix/thermomech_bimetall.inp
@@ -7111,7 +7111,7 @@ ConstraintTemperatureNormal,11,11,273.0
 *NODE FILE
 U, NT
 *EL FILE
-S, E
+S, E, HFL
 ** outputs --> dat file
 ** reaction forces for Constraint fixed
 *NODE PRINT, NSET=ConstraintFixed, TOTALS=ONLY