From 835a76a40cbe237d8ec6f0c841f8372d29d7d18e Mon Sep 17 00:00:00 2001 From: Sanath Keshav Date: Wed, 30 Oct 2024 12:55:29 +0100 Subject: [PATCH] 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 | 2 +- README.md | 15 ++- include/material_models/LinearThermal.h | 115 ++++++++++++++++++ .../material_models/LinearThermalIsotropic.h | 47 ------- include/setup.h | 4 +- 5 files changed, 130 insertions(+), 53 deletions(-) create mode 100644 include/material_models/LinearThermal.h delete mode 100644 include/material_models/LinearThermalIsotropic.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d80ea8a..5a7abb0 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/README.md b/README.md index faf8c09..b32f81d 100644 --- a/README.md +++ b/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 diff --git a/include/material_models/LinearThermal.h b/include/material_models/LinearThermal.h new file mode 100644 index 0000000..bb8fc47 --- /dev/null +++ b/include/material_models/LinearThermal.h @@ -0,0 +1,115 @@ +#ifndef LINEARTHERMAL_H +#define LINEARTHERMAL_H + +#include "matmodel.h" +#include // For Eigen's aligned_allocator + +class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> { + public: + LinearThermalIsotropic(vector l_e, json materialProperties) + : ThermalModel(l_e) + { + try { + conductivity = materialProperties["conductivity"].get>(); + } 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[n_mat]; + + for (size_t i = 0; i < n_mat; ++i) { + phase_stiffness[i] = Matrix::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 conductivity; +}; + +class LinearThermalTriclinic : public ThermalModel, public LinearModel<1> { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Ensure proper alignment for Eigen structures + + LinearThermalTriclinic(vector l_e, json materialProperties) + : ThermalModel(l_e) + { + vector K_keys = { + "K_11", "K_12", "K_13", + "K_22", "K_23", + "K_33"}; + + try { + n_mat = materialProperties.at("K_11").get>().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>(); + if (values.size() != n_mat) { + throw std::runtime_error("Inconsistent size for material property: " + K_keys[k]); + } + K_constants.row(k) = Eigen::Map(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[n_mat]; + for (size_t i = 0; i < n_mat; ++i) { + phase_stiffness[i] = Matrix::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> K_mats; + MatrixXd K_constants; +}; + +#endif // LINEARTHERMAL_H diff --git a/include/material_models/LinearThermalIsotropic.h b/include/material_models/LinearThermalIsotropic.h deleted file mode 100644 index 35446b8..0000000 --- a/include/material_models/LinearThermalIsotropic.h +++ /dev/null @@ -1,47 +0,0 @@ -#ifndef LINEARTHERMALISOTROPIC_H -#define LINEARTHERMALISOTROPIC_H - -#include "matmodel.h" - -class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> { - public: - LinearThermalIsotropic(vector l_e, json materialProperties) - : ThermalModel(l_e) - { - try { - conductivity = materialProperties["conductivity"].get>(); - } 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[n_mat]; - - for (size_t i = 0; i < n_mat; ++i) { - phase_stiffness[i] = Matrix::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 conductivity; -}; - -#endif // LINEARTHERMALISOTROPIC_H diff --git a/include/setup.h b/include/setup.h index beecd1c..a9fec2f 100644 --- a/include/setup.h +++ b/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"); }