added linear thermal triclinic model (#32)

* added linear thermal triclinic model similar to linear elastic triclinic model

* fixed readme.md indent causing linting error

CMakeLists.txt

@@ -133,7 +133,7 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/solver.h
         include/setup.h
 
-        include/material_models/LinearThermalIsotropic.h
+        include/material_models/LinearThermal.h
         include/material_models/LinearElastic.h
         include/material_models/PseudoPlastic.h
         include/material_models/J2Plasticity.h

README.md

@@ -178,13 +178,20 @@ FANS requires a JSON input file specifying the problem parameters. Example input
 
 ```json
 "method": "cg",
-"TOL": 1e-10,
-"n_it": 100
+"error_parameters":{
+    "measure": "Linfinity",
+    "type": "absolute",
+    "tolerance": 1e-10
+},
+"n_it": 100,
 ```
 
 - `method`: This indicates the numerical method to be used for solving the system of equations. `cg` stands for the Conjugate Gradient method, and `fp` stands for the Fixed Point method.
-- `TOL`: This sets the tolerance level for the solver. It defines the convergence criterion which is based on the L-infinity norm of the nodal finite element residual; the solver iterates until the solution meets this tolerance.
-- `n_it`: This specifies the maximum number of iterations allowed for the FANS solver.
+- `error_parameters`: This section defines the error parameters for the solver. Error control is applied on the finite element nodal residual of the problem.
+  - `measure`: Specifies the norm used to measure the error. Options include `Linfinity`, `L1`, or `L2`.
+  - `type`: Defines the type of error measurement. Options are `absolute` or `relative`.
+  - `tolerance`: Sets the tolerance level for the solver, defining the convergence criterion based on the chosen error measure. The solver iterates until the solution meets this tolerance.
+- `n_it`: Specifies the maximum number of iterations allowed for the FANS solver.
 
 ### Macroscale Loading Conditions
 

include/material_models/LinearThermal.h

@@ -0,0 +1,115 @@
+#ifndef LINEARTHERMAL_H
+#define LINEARTHERMAL_H
+
+#include "matmodel.h"
+#include <Eigen/StdVector> // For Eigen's aligned_allocator
+
+class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> {
+  public:
+    LinearThermalIsotropic(vector<double> l_e, json materialProperties)
+        : ThermalModel(l_e)
+    {
+        try {
+            conductivity = materialProperties["conductivity"].get<vector<double>>();
+        } catch (const std::exception &e) {
+            throw std::runtime_error("Missing material properties for the requested material model.");
+        }
+        n_mat = conductivity.size();
+
+        double kappa_ref = (*max_element(conductivity.begin(), conductivity.end()) +
+                            *min_element(conductivity.begin(), conductivity.end())) /
+                           2;
+        kapparef_mat = kappa_ref * Matrix3d::Identity();
+
+        Matrix3d phase_kappa;
+        phase_stiffness = new Matrix<double, 8, 8>[n_mat];
+
+        for (size_t i = 0; i < n_mat; ++i) {
+            phase_stiffness[i] = Matrix<double, 8, 8>::Zero();
+            phase_kappa        = conductivity[i] * Matrix3d::Identity();
+
+            for (int p = 0; p < 8; ++p) {
+                phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e * 0.1250;
+            }
+        }
+    }
+
+    void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override
+    {
+        sigma(i + 0, 0) = conductivity[mat_index] * eps(i + 0, 0);
+        sigma(i + 1, 0) = conductivity[mat_index] * eps(i + 1, 0);
+        sigma(i + 2, 0) = conductivity[mat_index] * eps(i + 2, 0);
+    }
+
+  private:
+    vector<double> conductivity;
+};
+
+class LinearThermalTriclinic : public ThermalModel, public LinearModel<1> {
+  public:
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Ensure proper alignment for Eigen structures
+
+    LinearThermalTriclinic(vector<double> l_e, json materialProperties)
+        : ThermalModel(l_e)
+    {
+        vector<string> K_keys = {
+            "K_11", "K_12", "K_13",
+            "K_22", "K_23",
+            "K_33"};
+
+        try {
+            n_mat                = materialProperties.at("K_11").get<vector<double>>().size();
+            size_t num_constants = K_keys.size();
+
+            // Initialize matrix to hold all constants
+            K_constants.resize(num_constants, n_mat);
+
+            // Load material constants into matrix
+            for (size_t k = 0; k < num_constants; ++k) {
+                const auto &values = materialProperties.at(K_keys[k]).get<vector<double>>();
+                if (values.size() != n_mat) {
+                    throw std::runtime_error("Inconsistent size for material property: " + K_keys[k]);
+                }
+                K_constants.row(k) = Eigen::Map<const RowVectorXd>(values.data(), values.size());
+            }
+        } catch (const std::exception &e) {
+            throw std::runtime_error("Missing or inconsistent material properties for the requested material model.");
+        }
+
+        // Assemble conductivity matrices for each material
+        K_mats.resize(n_mat);
+        kapparef_mat = Matrix3d::Zero();
+
+        for (size_t i = 0; i < n_mat; ++i) {
+            Matrix3d K_i;
+            K_i << K_constants(0, i), K_constants(1, i), K_constants(2, i),
+                K_constants(1, i), K_constants(3, i), K_constants(4, i),
+                K_constants(2, i), K_constants(4, i), K_constants(5, i);
+
+            K_mats[i] = K_i;
+            kapparef_mat += K_i;
+        }
+
+        kapparef_mat /= n_mat;
+
+        // Compute phase stiffness matrices
+        phase_stiffness = new Matrix<double, 8, 8>[n_mat];
+        for (size_t i = 0; i < n_mat; ++i) {
+            phase_stiffness[i] = Matrix<double, 8, 8>::Zero();
+            for (int p = 0; p < 8; ++p) {
+                phase_stiffness[i] += B_int[p].transpose() * K_mats[i] * B_int[p] * v_e * 0.1250;
+            }
+        }
+    }
+
+    void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override
+    {
+        sigma.segment<3>(i) = K_mats[mat_index] * eps.segment<3>(i);
+    }
+
+  private:
+    std::vector<Matrix3d, Eigen::aligned_allocator<Matrix3d>> K_mats;
+    MatrixXd                                                  K_constants;
+};
+
+#endif // LINEARTHERMAL_H

include/setup.h

@@ -2,7 +2,7 @@
 #include "solverFP.h"
 
 // Thermal models
-#include "material_models/LinearThermalIsotropic.h"
+#include "material_models/LinearThermal.h"
 
 // Mechanical models
 #include "material_models/LinearElastic.h"
@@ -17,6 +17,8 @@ Matmodel<1> *createMatmodel(const Reader &reader)
 {
     if (reader.matmodel == "LinearThermalIsotropic") {
         return new LinearThermalIsotropic(reader.l_e, reader.materialProperties);
+    } else if (reader.matmodel == "LinearThermalTriclinic") {
+        return new LinearThermalTriclinic(reader.l_e, reader.materialProperties);
     } else {
         throw std::invalid_argument(reader.matmodel + " is not a valid matmodel for thermal problem");
     }