|
| 1 | +#ifndef LARGESTRAINMECHMODEL_H |
| 2 | +#define LARGESTRAINMECHMODEL_H |
| 3 | + |
| 4 | +#include "matmodel.h" |
| 5 | + |
| 6 | +/** |
| 7 | + * @brief Base class for large strain mechanical material models |
| 8 | + * This class provides the kinematic framework for finite deformation mechanics. |
| 9 | + * Derived classes only need to implement the constitutive response (S from E). |
| 10 | + */ |
| 11 | +class LargeStrainMechModel : public Matmodel<3, 9> { |
| 12 | + public: |
| 13 | + LargeStrainMechModel(const Reader &reader) |
| 14 | + : Matmodel(reader) |
| 15 | + { |
| 16 | + Construct_B(); |
| 17 | + } |
| 18 | + |
| 19 | + protected: |
| 20 | + // ============================================================================ |
| 21 | + // Kinematic helper functions |
| 22 | + // ============================================================================ |
| 23 | + |
| 24 | + /** |
| 25 | + * @brief Extract deformation gradient from eps vector at Gauss point |
| 26 | + * @param i Index in eps array (should be n_str * gauss_point) |
| 27 | + * @return 3×3 deformation gradient matrix F |
| 28 | + */ |
| 29 | + inline Matrix3d extract_F(int i) const |
| 30 | + { |
| 31 | + Matrix3d F; |
| 32 | + F << eps(i), eps(i + 1), eps(i + 2), |
| 33 | + eps(i + 3), eps(i + 4), eps(i + 5), |
| 34 | + eps(i + 6), eps(i + 7), eps(i + 8); |
| 35 | + return F; |
| 36 | + } |
| 37 | + /** |
| 38 | + * @brief Store 1st Piola-Kirchhoff stress in sigma array |
| 39 | + * @param i Index in sigma array (should be n_str * gauss_point) |
| 40 | + * @param P 3×3 1st Piola-Kirchhoff stress tensor |
| 41 | + */ |
| 42 | + inline void store_P(int i, const Matrix3d &P) |
| 43 | + { |
| 44 | + sigma(i) = P(0, 0); |
| 45 | + sigma(i + 1) = P(0, 1); |
| 46 | + sigma(i + 2) = P(0, 2); |
| 47 | + sigma(i + 3) = P(1, 0); |
| 48 | + sigma(i + 4) = P(1, 1); |
| 49 | + sigma(i + 5) = P(1, 2); |
| 50 | + sigma(i + 6) = P(2, 0); |
| 51 | + sigma(i + 7) = P(2, 1); |
| 52 | + sigma(i + 8) = P(2, 2); |
| 53 | + } |
| 54 | + |
| 55 | + inline Matrix3d compute_C(const Matrix3d &F) const |
| 56 | + { |
| 57 | + return F.transpose() * F; // C = F^T F |
| 58 | + } |
| 59 | + inline Matrix3d compute_E(const Matrix3d &C) const |
| 60 | + { |
| 61 | + return 0.5 * (C - Matrix3d::Identity()); // E = 0.5 * (C - I) |
| 62 | + } |
| 63 | + inline Matrix3d compute_E_from_F(const Matrix3d &F) const |
| 64 | + { |
| 65 | + return compute_E(compute_C(F)); // E = 0.5 * (F^T F - I) |
| 66 | + } |
| 67 | + inline Matrix3d push_forward(const Matrix3d &F, const Matrix3d &S) const |
| 68 | + { |
| 69 | + return F * S; // P = F S |
| 70 | + } |
| 71 | + |
| 72 | + // ============================================================================ |
| 73 | + // Virtual functions for derived material models to implement |
| 74 | + // ============================================================================ |
| 75 | + |
| 76 | + /** |
| 77 | + * @brief Compute 2nd Piola-Kirchhoff stress and material tangent C = dS/dE from deformation gradient |
| 78 | + * |
| 79 | + * This is the core constitutive function that derived classes must implement. |
| 80 | + * It should compute S based on the material's constitutive law. |
| 81 | + * @param F Deformation gradient (3×3) |
| 82 | + * @param mat_index Material index |
| 83 | + * @param element_idx Element index |
| 84 | + * @return 2nd Piola-Kirchhoff stress S (symmetric 3×3) |
| 85 | + * @return the 6×6 material tangent C = dS/dE in Mandel notation: Ordering: (11, 22, 33, sqrt(2)·12, sqrt(2)·13, sqrt(2)·23) |
| 86 | + */ |
| 87 | + virtual Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx) = 0; |
| 88 | + virtual Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index) = 0; |
| 89 | + |
| 90 | + /** |
| 91 | + * @brief Convert material tangent C (dS/dE) to spatial tangent A (dP/dF) |
| 92 | + * |
| 93 | + * Computes the spatial tangent dP/dF from the material tangent dS/dE. |
| 94 | + * |
| 95 | + * The full expression is: |
| 96 | + * dP_iJ/dF_kL = δ_ik S_LJ + F_iM * dS_MJ/dE_PQ * dE_PQ/dF_kL |
| 97 | + * where: |
| 98 | + * dE_PQ/dF_kL = 0.5 * (F_kP δ_QL + F_kQ δ_PL) |
| 99 | + * |
| 100 | + * @param F Deformation gradient (3×3) |
| 101 | + * @param S 2nd Piola-Kirchhoff stress (3×3 symmetric) |
| 102 | + * @param C_mandel Material tangent dS/dE in Mandel notation (6×6) |
| 103 | + * @return Spatial tangent A = dP/dF in full notation (9×9) |
| 104 | + */ |
| 105 | + inline Matrix<double, 9, 9> compute_spatial_tangent(const Matrix3d &F, |
| 106 | + const Matrix3d &S, |
| 107 | + const Matrix<double, 6, 6> &C_mandel) const |
| 108 | + { |
| 109 | + Matrix<double, 9, 9> A = Matrix<double, 9, 9>::Zero(); |
| 110 | + |
| 111 | + // Mandel to tensor index mapping |
| 112 | + int mandel_to_ij[6][2] = {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}; |
| 113 | + |
| 114 | + // Compute A_iJkL = dP_iJ/dF_kL |
| 115 | + for (int i = 0; i < 3; ++i) { |
| 116 | + for (int J = 0; J < 3; ++J) { |
| 117 | + int row = 3 * i + J; // Index for P_iJ |
| 118 | + |
| 119 | + for (int k = 0; k < 3; ++k) { |
| 120 | + for (int L = 0; L < 3; ++L) { |
| 121 | + int col = 3 * k + L; // Index for F_kL |
| 122 | + |
| 123 | + // First term: δ_ik S_LJ |
| 124 | + if (i == k) { |
| 125 | + A(row, col) += S(L, J); |
| 126 | + } |
| 127 | + |
| 128 | + // Second term: F_iM * C_MJPQ * dE_PQ/dF_kL |
| 129 | + for (int M = 0; M < 3; ++M) { |
| 130 | + // Find MJ in Mandel notation |
| 131 | + int MJ_mandel = -1; |
| 132 | + for (int idx = 0; idx < 6; ++idx) { |
| 133 | + if ((mandel_to_ij[idx][0] == M && mandel_to_ij[idx][1] == J) || |
| 134 | + (mandel_to_ij[idx][0] == J && mandel_to_ij[idx][1] == M)) { |
| 135 | + MJ_mandel = idx; |
| 136 | + break; |
| 137 | + } |
| 138 | + } |
| 139 | + if (MJ_mandel < 0) |
| 140 | + continue; |
| 141 | + |
| 142 | + for (int P = 0; P < 3; ++P) { |
| 143 | + for (int Q = P; Q < 3; ++Q) { // Q >= P due to symmetry |
| 144 | + // Find PQ in Mandel notation |
| 145 | + int PQ_mandel = -1; |
| 146 | + for (int idx = 0; idx < 6; ++idx) { |
| 147 | + if ((mandel_to_ij[idx][0] == P && mandel_to_ij[idx][1] == Q) || |
| 148 | + (mandel_to_ij[idx][0] == Q && mandel_to_ij[idx][1] == P)) { |
| 149 | + PQ_mandel = idx; |
| 150 | + break; |
| 151 | + } |
| 152 | + } |
| 153 | + if (PQ_mandel < 0) |
| 154 | + continue; |
| 155 | + |
| 156 | + // Get C_MJPQ and account for Mandel factors |
| 157 | + double C_val = C_mandel(MJ_mandel, PQ_mandel); |
| 158 | + if (MJ_mandel >= 3) |
| 159 | + C_val /= sqrt(2.0); |
| 160 | + if (PQ_mandel >= 3) |
| 161 | + C_val /= sqrt(2.0); |
| 162 | + |
| 163 | + // Compute dE_PQ/dF_kL = 0.5 * (F_kP δ_QL + F_kQ δ_PL) |
| 164 | + double dE_dF = 0.0; |
| 165 | + if (Q == L) |
| 166 | + dE_dF += 0.5 * F(k, P); |
| 167 | + if (P == L) |
| 168 | + dE_dF += 0.5 * F(k, Q); |
| 169 | + |
| 170 | + A(row, col) += F(i, M) * C_val * dE_dF; |
| 171 | + } |
| 172 | + } |
| 173 | + } |
| 174 | + } |
| 175 | + } |
| 176 | + } |
| 177 | + } |
| 178 | + |
| 179 | + return A; |
| 180 | + } |
| 181 | + |
| 182 | + // ============================================================================ |
| 183 | + // Standard interface implementation |
| 184 | + // ============================================================================ |
| 185 | + |
| 186 | + /** |
| 187 | + * @brief Compute B matrix for deformation gradient (9×24) |
| 188 | + * Returns B_F matrix that relates nodal displacements to F components: |
| 189 | + */ |
| 190 | + Matrix<double, 9, 24> Compute_B(const double x, const double y, const double z) override; |
| 191 | + |
| 192 | + /** |
| 193 | + * @brief Compute stress at Gauss point (implements base class interface) |
| 194 | + * This function: |
| 195 | + * 1. Extracts F from eps |
| 196 | + * 2. Calls compute_S(F) - material model implements this |
| 197 | + * 3. Computes P = F S |
| 198 | + * 4. Stores P in sigma |
| 199 | + */ |
| 200 | + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override final |
| 201 | + { |
| 202 | + Matrix3d F = extract_F(i); // Extract deformation gradient |
| 203 | + Matrix3d S = compute_S(F, mat_index, element_idx); // Material model computes 2nd PK stress from F |
| 204 | + Matrix3d P = push_forward(F, S); // Push forward to 1st PK stress |
| 205 | + store_P(i, P); // Store in sigma array |
| 206 | + } |
| 207 | +}; |
| 208 | + |
| 209 | +inline Matrix<double, 9, 24> LargeStrainMechModel::Compute_B(const double x, const double y, const double z) |
| 210 | +{ |
| 211 | + Matrix<double, 9, 24> B_F = Matrix<double, 9, 24>::Zero(); |
| 212 | + Matrix<double, 3, 8> dN = Matmodel<3, 9>::Compute_basic_B(x, y, z); |
| 213 | + |
| 214 | + for (int i = 0; i < 3; ++i) { // Displacement component |
| 215 | + for (int J = 0; J < 3; ++J) { // Material coordinate derivative |
| 216 | + int row = 3 * i + J; // Row-major: F_iJ |
| 217 | + for (int node = 0; node < 8; ++node) { |
| 218 | + int col = 3 * node + i; // Column: DOF for u_i at node |
| 219 | + B_F(row, col) = dN(J, node); |
| 220 | + } |
| 221 | + } |
| 222 | + } |
| 223 | + |
| 224 | + return B_F; |
| 225 | +} |
| 226 | + |
| 227 | +#endif // LARGESTRAINMECHMODEL_H |
0 commit comments