feat: add mitc4 material integration scaffolding

2026-05-04 22:34:53 +09:00
parent d550e13e77
commit 2e1bd99d05
5 changed files with 433 additions and 6 deletions
@@ -1304,6 +1304,7 @@ struct MITC4IntegrationBasis {
 using MITC4ElementDofVector = std::array<Real, 24>;
 using MITC4StrainVector = std::array<Real, 6>;
 using MITC4StrainRow = std::array<Real, 24>;
+using MITC4MaterialMatrix = std::array<std::array<Real, 6>, 6>;

 enum class MITC4StrainComponent {
  Eps11 = 0,
@@ -1359,6 +1360,54 @@ struct MITC4StrainRows {
  }
 };

+struct MITC4MaterialMatrixEvaluation {
+  MITC4MaterialMatrix matrix{};
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
+struct MITC4StrainTransform {
+  MITC4MaterialMatrix matrix{};
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
+struct MITC4IntegrationPoint {
+  Real xi = 0.0;
+  Real eta = 0.0;
+  Real zeta = 0.0;
+  Real weight = 0.0;
+};
+
+struct MITC4MaterialIntegrationSample {
+  MITC4IntegrationPoint point;
+  MITC4IntegrationBasis basis;
+  MITC4StrainRows strain_rows;
+  MITC4MaterialMatrix local_material{};
+  MITC4MaterialMatrix strain_transform{};
+  MITC4MaterialMatrix convected_material{};
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
+struct MITC4MaterialIntegrationData {
+  std::vector<MITC4MaterialIntegrationSample> samples;
+  std::vector<Diagnostic> diagnostics;
+
+  bool ok() const {
+    return !hasError(diagnostics);
+  }
+};
+
 inline Vec3 globalEX() {
  return {1.0, 0.0, 0.0};
 }
@@ -1624,6 +1673,209 @@ inline MITC4StrainRows mitc4TiedCovariantStrainRows(const MITC4Geometry& geometr
  return result;
 }

+inline std::array<MITC4IntegrationPoint, 8> mitc4GaussQuadrature2x2x2() {
+  const Real gauss = 1.0 / std::sqrt(3.0);
+  const std::array<Real, 2> points = {-gauss, gauss};
+  std::array<MITC4IntegrationPoint, 8> integration_points{};
+  std::size_t index = 0;
+  for (Real xi : points) {
+    for (Real eta : points) {
+      for (Real zeta : points) {
+        integration_points[index++] = {xi, eta, zeta, 1.0};
+      }
+    }
+  }
+  return integration_points;
+}
+
+inline MITC4MaterialMatrixEvaluation mitc4PlaneStressMaterialMatrix(Real elastic_modulus, Real poisson_ratio,
+                                                                    Real shear_correction = 5.0 / 6.0,
+                                                                    Real tolerance = 1.0e-12) {
+  MITC4MaterialMatrixEvaluation result;
+  if (!isFinite(elastic_modulus) || elastic_modulus <= tolerance) {
+    result.diagnostics.push_back(
+        mitc4Diagnostic("FESA-MITC4-MATERIAL", "MITC4 elastic modulus must be positive and finite"));
+  }
+  if (!isFinite(poisson_ratio) || poisson_ratio <= -1.0 || poisson_ratio >= 0.5) {
+    result.diagnostics.push_back(
+        mitc4Diagnostic("FESA-MITC4-POISSON", "MITC4 isotropic Poisson ratio must satisfy -1 < nu < 0.5"));
+  }
+  if (!isFinite(shear_correction) || shear_correction <= tolerance) {
+    result.diagnostics.push_back(
+        mitc4Diagnostic("FESA-MITC4-SHEAR-CORRECTION", "MITC4 shear correction factor must be positive and finite"));
+  }
+  if (hasError(result.diagnostics)) {
+    return result;
+  }
+
+  const Real scale = elastic_modulus / (1.0 - poisson_ratio * poisson_ratio);
+  const Real shear_modulus = elastic_modulus / (2.0 * (1.0 + poisson_ratio));
+  const std::size_t eps11 = strainComponentIndex(MITC4StrainComponent::Eps11);
+  const std::size_t eps22 = strainComponentIndex(MITC4StrainComponent::Eps22);
+  const std::size_t gamma23 = strainComponentIndex(MITC4StrainComponent::Gamma23);
+  const std::size_t gamma13 = strainComponentIndex(MITC4StrainComponent::Gamma13);
+  const std::size_t gamma12 = strainComponentIndex(MITC4StrainComponent::Gamma12);
+
+  result.matrix[eps11][eps11] = scale;
+  result.matrix[eps11][eps22] = poisson_ratio * scale;
+  result.matrix[eps22][eps11] = poisson_ratio * scale;
+  result.matrix[eps22][eps22] = scale;
+  result.matrix[gamma23][gamma23] = shear_correction * shear_modulus;
+  result.matrix[gamma13][gamma13] = shear_correction * shear_modulus;
+  result.matrix[gamma12][gamma12] = shear_modulus;
+  return result;
+}
+
+inline MITC4StrainVector multiplyMITC4MaterialMatrix(const MITC4MaterialMatrix& matrix, const MITC4StrainVector& vector) {
+  MITC4StrainVector result{};
+  for (std::size_t row = 0; row < 6; ++row) {
+    for (std::size_t col = 0; col < 6; ++col) {
+      result[row] += matrix[row][col] * vector[col];
+    }
+  }
+  return result;
+}
+
+inline Real dotMITC4Vector(const MITC4StrainVector& a, const MITC4StrainVector& b) {
+  Real result = 0.0;
+  for (std::size_t i = 0; i < 6; ++i) {
+    result += a[i] * b[i];
+  }
+  return result;
+}
+
+inline MITC4MaterialMatrix mitc4TransformMaterialMatrix(const MITC4MaterialMatrix& local_material,
+                                                        const MITC4MaterialMatrix& covariant_to_local) {
+  MITC4MaterialMatrix result{};
+  for (std::size_t i = 0; i < 6; ++i) {
+    for (std::size_t j = 0; j < 6; ++j) {
+      Real value = 0.0;
+      for (std::size_t a = 0; a < 6; ++a) {
+        for (std::size_t b = 0; b < 6; ++b) {
+          value += covariant_to_local[a][i] * local_material[a][b] * covariant_to_local[b][j];
+        }
+      }
+      result[i][j] = value;
+    }
+  }
+  return result;
+}
+
+inline std::array<std::array<Real, 3>, 3> mitc4TensorFromEngineeringComponent(std::size_t component) {
+  std::array<std::array<Real, 3>, 3> tensor{};
+  switch (static_cast<MITC4StrainComponent>(component)) {
+    case MITC4StrainComponent::Eps11:
+      tensor[0][0] = 1.0;
+      break;
+    case MITC4StrainComponent::Eps22:
+      tensor[1][1] = 1.0;
+      break;
+    case MITC4StrainComponent::Eps33:
+      tensor[2][2] = 1.0;
+      break;
+    case MITC4StrainComponent::Gamma23:
+      tensor[1][2] = 0.5;
+      tensor[2][1] = 0.5;
+      break;
+    case MITC4StrainComponent::Gamma13:
+      tensor[0][2] = 0.5;
+      tensor[2][0] = 0.5;
+      break;
+    case MITC4StrainComponent::Gamma12:
+      tensor[0][1] = 0.5;
+      tensor[1][0] = 0.5;
+      break;
+  }
+  return tensor;
+}
+
+inline MITC4StrainVector mitc4EngineeringVectorFromTensor(const std::array<std::array<Real, 3>, 3>& tensor) {
+  MITC4StrainVector vector{};
+  vector[strainComponentIndex(MITC4StrainComponent::Eps11)] = tensor[0][0];
+  vector[strainComponentIndex(MITC4StrainComponent::Eps22)] = tensor[1][1];
+  vector[strainComponentIndex(MITC4StrainComponent::Eps33)] = tensor[2][2];
+  vector[strainComponentIndex(MITC4StrainComponent::Gamma23)] = 2.0 * tensor[1][2];
+  vector[strainComponentIndex(MITC4StrainComponent::Gamma13)] = 2.0 * tensor[0][2];
+  vector[strainComponentIndex(MITC4StrainComponent::Gamma12)] = 2.0 * tensor[0][1];
+  return vector;
+}
+
+inline MITC4StrainTransform mitc4CovariantToLocalStrainTransform(const MITC4IntegrationBasis& basis,
+                                                                 Real tolerance = 1.0e-12) {
+  MITC4StrainTransform result;
+  result.diagnostics = basis.diagnostics;
+  const Real jacobian = dot(cross(basis.g1, basis.g2), basis.g3);
+  if (!isFinite(jacobian) || std::fabs(jacobian) <= tolerance) {
+    result.diagnostics.push_back(
+        mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 material transform Jacobian is near zero"));
+    return result;
+  }
+  if (hasError(result.diagnostics)) {
+    return result;
+  }
+
+  const std::array<Vec3, 3> contravariant = {
+      (1.0 / jacobian) * cross(basis.g2, basis.g3),
+      (1.0 / jacobian) * cross(basis.g3, basis.g1),
+      (1.0 / jacobian) * cross(basis.g1, basis.g2)};
+  const std::array<Vec3, 3> local = {basis.local.e1, basis.local.e2, basis.local.e3};
+  std::array<std::array<Real, 3>, 3> direction_cosines{};
+  for (std::size_t local_axis = 0; local_axis < 3; ++local_axis) {
+    for (std::size_t convected_axis = 0; convected_axis < 3; ++convected_axis) {
+      direction_cosines[local_axis][convected_axis] = dot(local[local_axis], contravariant[convected_axis]);
+    }
+  }
+
+  for (std::size_t column = 0; column < 6; ++column) {
+    const auto covariant_tensor = mitc4TensorFromEngineeringComponent(column);
+    std::array<std::array<Real, 3>, 3> local_tensor{};
+    for (std::size_t a = 0; a < 3; ++a) {
+      for (std::size_t b = 0; b < 3; ++b) {
+        for (std::size_t i = 0; i < 3; ++i) {
+          for (std::size_t j = 0; j < 3; ++j) {
+            local_tensor[a][b] += direction_cosines[a][i] * direction_cosines[b][j] * covariant_tensor[i][j];
+          }
+        }
+      }
+    }
+    const auto local_vector = mitc4EngineeringVectorFromTensor(local_tensor);
+    for (std::size_t row = 0; row < 6; ++row) {
+      result.matrix[row][column] = local_vector[row];
+    }
+  }
+  return result;
+}
+
+inline MITC4MaterialIntegrationData mitc4BuildMaterialIntegrationData(const MITC4Geometry& geometry, Real elastic_modulus,
+                                                                      Real poisson_ratio,
+                                                                      Real shear_correction = 5.0 / 6.0) {
+  MITC4MaterialIntegrationData data;
+  const auto material = mitc4PlaneStressMaterialMatrix(elastic_modulus, poisson_ratio, shear_correction);
+  appendDiagnostics(data.diagnostics, material.diagnostics);
+  if (hasError(data.diagnostics)) {
+    return data;
+  }
+
+  for (const MITC4IntegrationPoint& point : mitc4GaussQuadrature2x2x2()) {
+    MITC4MaterialIntegrationSample sample;
+    sample.point = point;
+    sample.local_material = material.matrix;
+    sample.basis = computeMITC4IntegrationBasis(geometry, point.xi, point.eta, point.zeta);
+    appendDiagnostics(sample.diagnostics, sample.basis.diagnostics);
+    const auto transform = mitc4CovariantToLocalStrainTransform(sample.basis);
+    sample.strain_transform = transform.matrix;
+    appendDiagnostics(sample.diagnostics, transform.diagnostics);
+    sample.strain_rows = mitc4TiedCovariantStrainRows(geometry, point.xi, point.eta, point.zeta);
+    appendDiagnostics(sample.diagnostics, sample.strain_rows.diagnostics);
+    if (!hasError(sample.diagnostics)) {
+      sample.convected_material = mitc4TransformMaterialMatrix(sample.local_material, sample.strain_transform);
+    }
+    appendDiagnostics(data.diagnostics, sample.diagnostics);
+    data.samples.push_back(sample);
+  }
+  return data;
+}
+
 inline LocalBasis computeLocalBasis(const std::array<Vec3, 4>& coordinates) {
  const MITC4Geometry geometry = buildMITC4Geometry(coordinates, 1.0);
  if (!geometry.ok()) {