diff --git a/CMakeLists.txt b/CMakeLists.txt index 1918e68..d80ea8a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -134,7 +134,7 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER include/setup.h include/material_models/LinearThermalIsotropic.h - include/material_models/LinearElasticIsotropic.h + include/material_models/LinearElastic.h include/material_models/PseudoPlastic.h include/material_models/J2Plasticity.h ) diff --git a/include/material_models/LinearElastic.h b/include/material_models/LinearElastic.h new file mode 100644 index 0000000..39f54c8 --- /dev/null +++ b/include/material_models/LinearElastic.h @@ -0,0 +1,145 @@ +#ifndef LINEARELASTIC_H +#define LINEARELASTIC_H + +#include "matmodel.h" +#include // For Eigen's aligned_allocator + +class LinearElasticIsotropic : public MechModel, public LinearModel<3> { + public: + LinearElasticIsotropic(vector l_e, map> materialProperties) + : MechModel(l_e) + { + try { + bulk_modulus = materialProperties["bulk_modulus"]; + mu = materialProperties["shear_modulus"]; + } catch (const std::exception &e) { + throw std::runtime_error("Missing material properties for the requested material model."); + } + + n_mat = bulk_modulus.size(); + lambda.resize(n_mat); + mu.resize(n_mat); + + for (int i = 0; i < n_mat; ++i) { + lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i]; + } + + double lambda_ref = (*max_element(lambda.begin(), lambda.end()) + + *min_element(lambda.begin(), lambda.end())) / + 2; + double mu_ref = (*max_element(mu.begin(), mu.end()) + + *min_element(mu.begin(), mu.end())) / + 2; + + kapparef_mat = Matrix::Zero(); + kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref); + kapparef_mat += 2 * mu_ref * Matrix::Identity(); + + phase_stiffness = new Matrix[n_mat]; + Matrix phase_kappa; + + for (int i = 0; i < n_mat; i++) { + phase_kappa.setZero(); + phase_stiffness[i] = Matrix::Zero(); + + phase_kappa.topLeftCorner(3, 3).setConstant(lambda[i]); + phase_kappa += 2 * mu[i] * Matrix::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 + { + double buf1 = lambda[mat_index] * (eps(i, 0) + eps(i + 1, 0) + eps(i + 2, 0)); + double buf2 = 2 * mu[mat_index]; + sigma(i + 0, 0) = buf1 + buf2 * eps(i + 0, 0); + sigma(i + 1, 0) = buf1 + buf2 * eps(i + 1, 0); + sigma(i + 2, 0) = buf1 + buf2 * eps(i + 2, 0); + sigma(i + 3, 0) = buf2 * eps(i + 3, 0); + sigma(i + 4, 0) = buf2 * eps(i + 4, 0); + sigma(i + 5, 0) = buf2 * eps(i + 5, 0); + } + + private: + vector bulk_modulus; + vector lambda; + vector mu; +}; + +class LinearElasticTriclinic : public MechModel, public LinearModel<3> { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Ensure proper alignment for Eigen structures + + LinearElasticTriclinic(const vector &l_e, const map> &materialProperties) + : MechModel(l_e) + { + vector C_keys = { + "C_11", "C_12", "C_13", "C_14", "C_15", "C_16", + "C_22", "C_23", "C_24", "C_25", "C_26", + "C_33", "C_34", "C_35", "C_36", + "C_44", "C_45", "C_46", + "C_55", "C_56", + "C_66"}; + + n_mat = materialProperties.at("C_11").size(); + size_t num_constants = C_keys.size(); + + // Initialize matrix to hold all constants + MatrixXd C_constants(num_constants, n_mat); + + // Load material constants into matrix + for (size_t k = 0; k < num_constants; ++k) { + const auto &values = materialProperties.at(C_keys[k]); + if (values.size() != n_mat) { + throw std::runtime_error("Inconsistent size for material property: " + C_keys[k]); + } + C_constants.row(k) = Eigen::Map(values.data(), values.size()); + } + + // Assemble stiffness matrices for each material + C_mats.resize(n_mat); + kapparef_mat = Matrix::Zero(); + + for (size_t i = 0; i < n_mat; ++i) { + Matrix C_i = Matrix::Zero(); + int k = 0; // Index for C_constants + + // Assign constants to the upper triangular part + for (int row = 0; row < 6; ++row) { + for (int col = row; col < 6; ++col) { + C_i(row, col) = C_constants(k++, i); + } + } + + // Symmetrize the matrix + C_i = C_i.selfadjointView(); + + C_mats[i] = C_i; + kapparef_mat += C_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() * C_mats[i] * B_int[p] * v_e * 0.1250; + } + } + } + + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override + { + sigma.segment<6>(i) = C_mats[mat_index] * eps.segment<6>(i); + } + + private: + std::vector, Eigen::aligned_allocator>> C_mats; +}; + +#endif // LINEARELASTIC_H diff --git a/include/material_models/LinearElasticIsotropic.h b/include/material_models/LinearElasticIsotropic.h deleted file mode 100644 index 684ba79..0000000 --- a/include/material_models/LinearElasticIsotropic.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef LINEARELASTICISOTROPIC_H -#define LINEARELASTICISOTROPIC_H - -#include "matmodel.h" - -class LinearElasticIsotropic : public MechModel, public LinearModel<3> { - public: - LinearElasticIsotropic(vector l_e, map> materialProperties) - : MechModel(l_e) - { - try { - bulk_modulus = materialProperties["bulk_modulus"]; - mu = materialProperties["shear_modulus"]; - } catch (const std::exception &e) { - throw std::runtime_error("Missing material properties for the requested material model."); - } - - n_mat = bulk_modulus.size(); - lambda.resize(n_mat); - mu.resize(n_mat); - - for (int i = 0; i < n_mat; ++i) { - lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i]; - } - - double lambda_ref = (*max_element(lambda.begin(), lambda.end()) + - *min_element(lambda.begin(), lambda.end())) / - 2; - double mu_ref = (*max_element(mu.begin(), mu.end()) + - *min_element(mu.begin(), mu.end())) / - 2; - - kapparef_mat = Matrix::Zero(); - kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref); - kapparef_mat += 2 * mu_ref * Matrix::Identity(); - - phase_stiffness = new Matrix[n_mat]; - Matrix phase_kappa; - - for (int i = 0; i < n_mat; i++) { - phase_kappa.setZero(); - phase_stiffness[i] = Matrix::Zero(); - - phase_kappa.topLeftCorner(3, 3).setConstant(lambda[i]); - phase_kappa += 2 * mu[i] * Matrix::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 - { - double buf1 = lambda[mat_index] * (eps(i, 0) + eps(i + 1, 0) + eps(i + 2, 0)); - double buf2 = 2 * mu[mat_index]; - sigma(i + 0, 0) = buf1 + buf2 * eps(i + 0, 0); - sigma(i + 1, 0) = buf1 + buf2 * eps(i + 1, 0); - sigma(i + 2, 0) = buf1 + buf2 * eps(i + 2, 0); - sigma(i + 3, 0) = buf2 * eps(i + 3, 0); - sigma(i + 4, 0) = buf2 * eps(i + 4, 0); - sigma(i + 5, 0) = buf2 * eps(i + 5, 0); - } - - private: - vector bulk_modulus; - vector lambda; - vector mu; -}; - -#endif // LINEARELASTICISOTROPIC_H diff --git a/include/setup.h b/include/setup.h index 79d3581..beecd1c 100644 --- a/include/setup.h +++ b/include/setup.h @@ -5,7 +5,7 @@ #include "material_models/LinearThermalIsotropic.h" // Mechanical models -#include "material_models/LinearElasticIsotropic.h" +#include "material_models/LinearElastic.h" #include "material_models/PseudoPlastic.h" #include "material_models/J2Plasticity.h" @@ -28,6 +28,8 @@ Matmodel<3> *createMatmodel(const Reader &reader) // Linear Elastic models if (reader.matmodel == "LinearElasticIsotropic") { return new LinearElasticIsotropic(reader.l_e, reader.materialProperties); + } else if (reader.matmodel == "LinearElasticTriclinic") { + return new LinearElasticTriclinic(reader.l_e, reader.materialProperties); // Pseudo Plastic models } else if (reader.matmodel == "PseudoPlasticLinearHardening") {