KOKO\Mimi
1149 lines
47 KiB
C++
1149 lines
47 KiB
C++
#include "fesa/elements/mitc4_shell.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <utility>
|
|
|
|
#include "fesa/model/model_types.h"
|
|
#include "math/dense_blas_internal.h"
|
|
#include "math/finite_values.h"
|
|
#include "math/vector3_operations.h"
|
|
#include "elements/shell_kernel_values.h"
|
|
#include "elements/mitc4_drilling_internal.h"
|
|
|
|
namespace fesa {
|
|
|
|
Mitc4ShellDefinition::Mitc4ShellDefinition(
|
|
SourceEntityId source_id_value,
|
|
const ShellSourceElementType source_type_value,
|
|
std::array<EntityIndex, 4> node_indices_value,
|
|
const EntityIndex material_index_value,
|
|
const EntityIndex section_index_value, SourceLocation location_value)
|
|
: source_id{std::move(source_id_value)},
|
|
source_type{source_type_value},
|
|
node_indices{std::move(node_indices_value)},
|
|
material_index{material_index_value},
|
|
section_index{section_index_value},
|
|
location{std::move(location_value)},
|
|
property_index_{section_index_value},
|
|
node_indices_view_{node_indices.begin(), node_indices.end()} {}
|
|
|
|
ElementDefinitionKind Mitc4ShellDefinition::Kind() const noexcept {
|
|
return ElementDefinitionKind::kMitc4Shell;
|
|
}
|
|
|
|
const SourceEntityId& Mitc4ShellDefinition::SourceId() const noexcept {
|
|
return source_id;
|
|
}
|
|
|
|
std::string_view Mitc4ShellDefinition::SourceElementType() const noexcept {
|
|
return source_type == ShellSourceElementType::kS4 ? "S4" : "S4R";
|
|
}
|
|
|
|
const std::vector<EntityIndex>& Mitc4ShellDefinition::NodeIndices()
|
|
const noexcept {
|
|
return node_indices_view_;
|
|
}
|
|
|
|
EntityIndex Mitc4ShellDefinition::MaterialIndex() const noexcept {
|
|
return material_index;
|
|
}
|
|
|
|
EntityIndex Mitc4ShellDefinition::PropertyIndex() const noexcept {
|
|
return property_index_;
|
|
}
|
|
|
|
void Mitc4ShellDefinition::SetPropertyIndex(
|
|
const EntityIndex property_index) noexcept {
|
|
property_index_ = property_index;
|
|
}
|
|
|
|
void Mitc4ShellDefinition::SynchronizeNodeIndices() {
|
|
node_indices_view_.assign(node_indices.begin(), node_indices.end());
|
|
}
|
|
|
|
namespace {
|
|
|
|
using shell_kernel_internal::FrameComponent;
|
|
using shell_kernel_internal::LocalEngineeringComponents;
|
|
using shell_kernel_internal::ElementIdentity;
|
|
using math_internal::DerivativeSum;
|
|
using math_internal::WeightedSum;
|
|
using math_internal::IsFinite;
|
|
|
|
constexpr std::size_t kNodeCount = 4U;
|
|
constexpr std::size_t kPhysicalDofsPerNode = 5U;
|
|
constexpr std::size_t kGlobalDofsPerNode = 6U;
|
|
constexpr std::size_t kPhysicalDofCount = 20U;
|
|
constexpr std::size_t kGlobalDofCount = 24U;
|
|
constexpr double kShearCorrection = 5.0 / 6.0;
|
|
constexpr double kFrameTolerance = 1.0e-12;
|
|
|
|
bool SameCoordinates(const Vector3& left, const Vector3& right) {
|
|
return left == right;
|
|
}
|
|
|
|
/// @brief Selects each director's least-aligned global axis deterministically
|
|
/// and projects it into the tangent plane.
|
|
std::array<Vector3, kNodeCount> NodalTangentsA(
|
|
const std::array<Vector3, kNodeCount>& directors) {
|
|
const std::array<Vector3, 3> global_axes{
|
|
Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};
|
|
std::array<Vector3, kNodeCount> tangents{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
std::size_t selected = 0U;
|
|
double alignment = std::abs(global_axes[0].Dot(directors[node]));
|
|
for (std::size_t axis = 1U; axis < global_axes.size(); ++axis) {
|
|
const double candidate = std::abs(global_axes[axis].Dot(directors[node]));
|
|
if (candidate < alignment) {
|
|
selected = axis;
|
|
alignment = candidate;
|
|
}
|
|
}
|
|
const Vector3 tangent_candidate =
|
|
global_axes[selected] -
|
|
global_axes[selected].Dot(directors[node]) * directors[node];
|
|
tangents[node] = (1.0 / tangent_candidate.Norm()) * tangent_candidate;
|
|
}
|
|
return tangents;
|
|
}
|
|
|
|
/// @brief Completes each right-handed nodal director frame.
|
|
std::array<Vector3, kNodeCount> NodalTangentsB(
|
|
const std::array<Vector3, kNodeCount>& directors,
|
|
const std::array<Vector3, kNodeCount>& tangent_a) {
|
|
std::array<Vector3, kNodeCount> tangents{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
tangents[node] = directors[node].Cross(tangent_a[node]);
|
|
}
|
|
return tangents;
|
|
}
|
|
|
|
Result<Mitc4Shell> ModelFailure(std::string code,
|
|
const SourceLocation& location,
|
|
const std::string& identity,
|
|
std::string message) {
|
|
return Result<Mitc4Shell>::Failure(Status::Failure(
|
|
FailureCategory::kModel, {{Severity::kError, std::move(code), location,
|
|
"*ELEMENT", identity, std::move(message)}}));
|
|
}
|
|
|
|
/// @brief Forms one physical covariant strain column with thickness stretch
|
|
/// excluded.
|
|
std::array<std::array<double, 3>, 3> CovariantStrainColumn(
|
|
const std::array<Vector3, 3>& covariant,
|
|
const std::array<Vector3, 3>& derivatives) {
|
|
std::array<std::array<double, 3>, 3> strain{};
|
|
for (std::size_t first = 0U; first < 3U; ++first) {
|
|
for (std::size_t second = 0U; second < 3U; ++second) {
|
|
strain[first][second] = 0.5 * (covariant[first].Dot(derivatives[second]) +
|
|
covariant[second].Dot(derivatives[first]));
|
|
}
|
|
}
|
|
// Thickness stretch is excluded from the five-component shell law.
|
|
strain[2U][2U] = 0.0;
|
|
return strain;
|
|
}
|
|
|
|
/// @brief Solves C^T E C = E_cov for the actual strain tensor.
|
|
std::array<std::array<double, 3>, 3> ReconstructCartesianStrain(
|
|
const std::array<std::array<double, 3>, 3>& covariant_strain,
|
|
const dense_blas_internal::DenseFactor& factor) {
|
|
Matrix rhs{3U, 3U};
|
|
for (std::size_t i = 0; i < 3U; ++i)
|
|
for (std::size_t j = 0; j < 3U; ++j) rhs(i, j) = covariant_strain[i][j];
|
|
const auto left = factor.Solve(rhs, true);
|
|
if (!left.HasValue())
|
|
throw std::invalid_argument{"MITC4 Cartesian strain is nonfinite."};
|
|
const auto solved = factor.SolveRight(left.Value());
|
|
if (!solved.HasValue())
|
|
throw std::invalid_argument{"MITC4 Cartesian strain is nonfinite."};
|
|
std::array<std::array<double, 3>, 3> tensor{};
|
|
for (std::size_t i = 0; i < 3U; ++i)
|
|
for (std::size_t j = 0; j < 3U; ++j) tensor[i][j] = solved.Value()(i, j);
|
|
return tensor;
|
|
}
|
|
|
|
/// @brief Packs covariant columns without constructing reciprocal vectors.
|
|
Matrix CovariantMatrix(const std::array<Vector3, 3>& covariant) {
|
|
Matrix coefficients{3U, 3U};
|
|
for (std::size_t i = 0; i < 3U; ++i)
|
|
for (std::size_t j = 0; j < 3U; ++j) coefficients(i, j) = covariant[j][i];
|
|
return coefficients;
|
|
}
|
|
|
|
Matrix ScaledMatrix(const Matrix& source, double factor) {
|
|
Matrix result{source};
|
|
result.Scale(factor);
|
|
return result;
|
|
}
|
|
|
|
/// @brief Lifts a local matrix through MKL and mirrors the upper triangle.
|
|
Matrix Congruence(const Matrix& local, const Matrix& transformation) {
|
|
if (local.Rows() != local.Columns() ||
|
|
local.Rows() != transformation.Rows()) {
|
|
throw std::invalid_argument{
|
|
"MITC4 congruence dimensions are incompatible."};
|
|
}
|
|
Matrix result = dense_blas_internal::MultiplyTransposeLeft(
|
|
transformation, local.Multiply(transformation));
|
|
for (std::size_t row = 0U; row < result.Rows(); ++row) {
|
|
for (std::size_t column = row; column < result.Columns(); ++column) {
|
|
result(column, row) = result(row, column);
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
Result<Mitc4Stiffness> StiffnessFailure(const SourceLocation& location,
|
|
const std::string& identity,
|
|
std::string message) {
|
|
return Result<Mitc4Stiffness>::Failure(
|
|
Status::Failure(FailureCategory::kModel,
|
|
{{Severity::kError, "invalid-shell-stiffness", location,
|
|
"*ELEMENT", identity, std::move(message)}}));
|
|
}
|
|
|
|
Result<Mitc4Mass> MassFailure(const SourceLocation& location,
|
|
const std::string& identity,
|
|
std::string message) {
|
|
return Result<Mitc4Mass>::Failure(
|
|
Status::Failure(FailureCategory::kModel,
|
|
{{Severity::kError, "invalid-shell-mass", location,
|
|
"*ELEMENT", identity, std::move(message)}}));
|
|
}
|
|
|
|
/// @brief Returns S=-[director]_cross for theta cross director.
|
|
Matrix DirectorSkew(const Vector3& director) {
|
|
Matrix result{3U, 3U};
|
|
result(0U, 1U) = director.Z();
|
|
result(0U, 2U) = -director.Y();
|
|
result(1U, 0U) = -director.Z();
|
|
result(1U, 2U) = director.X();
|
|
result(2U, 0U) = director.Y();
|
|
result(2U, 1U) = -director.X();
|
|
return result;
|
|
}
|
|
|
|
Result<Mitc4PhysicalRecovery> RecoveryFailure(const SourceLocation& location,
|
|
const std::string& identity,
|
|
std::string message) {
|
|
return Result<Mitc4PhysicalRecovery>::Failure(
|
|
Status::Failure(FailureCategory::kModel,
|
|
{{Severity::kError, "invalid-shell-recovery", location,
|
|
"*ELEMENT", identity, std::move(message)}}));
|
|
}
|
|
|
|
/// @brief Creates a structured failure at the runtime element boundary.
|
|
template <class T>
|
|
Result<T> RuntimeFailure(const SourceLocation& location,
|
|
const SourceEntityId& source_id, std::string message) {
|
|
return Result<T>::Failure(Status::Failure(
|
|
FailureCategory::kModel,
|
|
{{Severity::kError, "invalid-runtime-element", location, "*ELEMENT",
|
|
source_id.source_label_text, std::move(message)}}));
|
|
}
|
|
|
|
} // namespace
|
|
|
|
Result<Mitc4Shell> Mitc4Shell::Create(
|
|
std::array<const Node*, 4> nodes,
|
|
std::array<std::array<double, 3>, 4> initial_directors,
|
|
const ShellSection& section, const LinearElasticMaterial& material) {
|
|
const std::string identity = ElementIdentity(nodes);
|
|
if (std::any_of(nodes.begin(), nodes.end(),
|
|
[](const Node* value) { return value == nullptr; })) {
|
|
return ModelFailure("invalid-shell-geometry", section.location, identity,
|
|
"MITC4 creation requires four valid node references.");
|
|
}
|
|
|
|
std::array<Vector3, kNodeCount> coordinates{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
coordinates[node] = Vector3{nodes[node]->coordinates};
|
|
if (!coordinates[node].IsFinite()) {
|
|
return ModelFailure("invalid-shell-geometry", nodes[node]->location,
|
|
identity, "MITC4 node coordinates must be finite.");
|
|
}
|
|
for (std::size_t previous = 0U; previous < node; ++previous) {
|
|
if (SameCoordinates(coordinates[node], coordinates[previous])) {
|
|
return ModelFailure("invalid-shell-geometry", nodes[node]->location,
|
|
identity,
|
|
"MITC4 nodes must have distinct coordinates.");
|
|
}
|
|
}
|
|
}
|
|
|
|
const std::array<Vector3, kNodeCount> directors{
|
|
Vector3{initial_directors[0U]}, Vector3{initial_directors[1U]},
|
|
Vector3{initial_directors[2U]}, Vector3{initial_directors[3U]}};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const double director_norm = directors[node].Norm();
|
|
if (!directors[node].IsFinite() || !std::isfinite(director_norm) ||
|
|
std::abs(director_norm - 1.0) > kFrameTolerance) {
|
|
return ModelFailure(
|
|
"invalid-shell-director", nodes[node]->location, identity,
|
|
"MITC4 initial directors must be finite unit vectors.");
|
|
}
|
|
}
|
|
|
|
if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
|
|
return ModelFailure("invalid-shell-section", section.location, identity,
|
|
"MITC4 shell thickness must be finite and positive.");
|
|
}
|
|
if (!std::isfinite(material.youngs_modulus) ||
|
|
!(material.youngs_modulus > 0.0) ||
|
|
!std::isfinite(material.poisson_ratio) ||
|
|
!(material.poisson_ratio > -1.0) || !(material.poisson_ratio < 0.5)) {
|
|
return ModelFailure(
|
|
"invalid-shell-material", material.location, identity,
|
|
"MITC4 isotropic material requires finite E>0 and -1<nu<0.5.");
|
|
}
|
|
|
|
const auto center_shape = ShapeFunctions(0.0, 0.0);
|
|
const Vector3 center_xi =
|
|
DerivativeSum(center_shape.xi_derivatives, coordinates);
|
|
const Vector3 center_eta =
|
|
DerivativeSum(center_shape.eta_derivatives, coordinates);
|
|
const Vector3 center_area = center_xi.Cross(center_eta);
|
|
const double center_measure = center_area.Norm();
|
|
if (!center_area.IsFinite() || !std::isfinite(center_measure) ||
|
|
!(center_measure > 0.0)) {
|
|
return ModelFailure(
|
|
"invalid-shell-geometry", nodes[0]->location, identity,
|
|
"MITC4 center surface basis must be finite and nonzero.");
|
|
}
|
|
const Vector3 normal_candidate = (1.0 / center_measure) * center_area;
|
|
if (std::any_of(directors.begin(), directors.end(),
|
|
[&normal_candidate](const Vector3& director) {
|
|
return !(normal_candidate.Dot(director) > 0.0);
|
|
})) {
|
|
return ModelFailure(
|
|
"invalid-shell-director", nodes[0]->location, identity,
|
|
"MITC4 directors must follow the source-order positive face.");
|
|
}
|
|
|
|
const auto tangent_a = NodalTangentsA(directors);
|
|
const auto tangent_b = NodalTangentsB(directors, tangent_a);
|
|
Mitc4Shell shell{coordinates,
|
|
directors,
|
|
tangent_a,
|
|
tangent_b,
|
|
normal_candidate,
|
|
section.thickness,
|
|
material.youngs_modulus,
|
|
material.poisson_ratio,
|
|
nodes[0]->location,
|
|
identity};
|
|
|
|
const double shear_modulus =
|
|
material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));
|
|
const double plane_stress_factor =
|
|
material.youngs_modulus /
|
|
(1.0 - material.poisson_ratio * material.poisson_ratio);
|
|
const double thickness_cubed =
|
|
section.thickness * section.thickness * section.thickness;
|
|
const std::array<double, 6> derived{
|
|
shear_modulus,
|
|
plane_stress_factor,
|
|
kShearCorrection * shear_modulus,
|
|
plane_stress_factor * section.thickness,
|
|
plane_stress_factor * thickness_cubed / 12.0,
|
|
kShearCorrection * shear_modulus * section.thickness};
|
|
if (std::any_of(derived.begin(), derived.end(), [](double value) {
|
|
return !std::isfinite(value) || !(value > 0.0);
|
|
})) {
|
|
return ModelFailure(
|
|
"invalid-shell-material", material.location, identity,
|
|
"Derived MITC4 constitutive coefficients must be finite and positive.");
|
|
}
|
|
|
|
GeometryData geometry{};
|
|
if (!shell.EvaluateGeometry(0.0, 0.0, 0.0, geometry)) {
|
|
return ModelFailure("invalid-shell-jacobian", nodes[0]->location, identity,
|
|
"MITC4 center frame or Jacobian is invalid.");
|
|
}
|
|
for (const auto& point : VolumeQuadrature()) {
|
|
if (!shell.EvaluateGeometry(point.natural_coordinates[0],
|
|
point.natural_coordinates[1],
|
|
point.natural_coordinates[2], geometry)) {
|
|
return ModelFailure("invalid-shell-jacobian", nodes[0]->location,
|
|
identity,
|
|
"MITC4 quadrature frame or Jacobian is invalid.");
|
|
}
|
|
}
|
|
constexpr std::array<Vector3, 4> tying_points{
|
|
Vector3{0.0, -1.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{-1.0, 0.0, 0.0},
|
|
Vector3{1.0, 0.0, 0.0}};
|
|
for (const auto& point : tying_points) {
|
|
if (!shell.EvaluateGeometry(point[0], point[1], point[2], geometry)) {
|
|
return ModelFailure("invalid-shell-jacobian", nodes[0]->location,
|
|
identity,
|
|
"MITC4 tying-point frame or Jacobian is invalid.");
|
|
}
|
|
}
|
|
|
|
return Result<Mitc4Shell>::Success(std::move(shell));
|
|
}
|
|
|
|
Mitc4ShapeFunctions Mitc4Shell::ShapeFunctions(double xi, double eta) noexcept {
|
|
constexpr std::array<double, kNodeCount> xi_signs{-1.0, 1.0, 1.0, -1.0};
|
|
constexpr std::array<double, kNodeCount> eta_signs{-1.0, -1.0, 1.0, 1.0};
|
|
Mitc4ShapeFunctions shape{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
shape.values[node] =
|
|
0.25 * (1.0 + xi_signs[node] * xi) * (1.0 + eta_signs[node] * eta);
|
|
shape.xi_derivatives[node] =
|
|
0.25 * xi_signs[node] * (1.0 + eta_signs[node] * eta);
|
|
shape.eta_derivatives[node] =
|
|
0.25 * eta_signs[node] * (1.0 + xi_signs[node] * xi);
|
|
}
|
|
return shape;
|
|
}
|
|
|
|
Mitc4TyingWeights Mitc4Shell::TyingWeights(double xi, double eta) noexcept {
|
|
return {{(1.0 - eta) * 0.5, (1.0 + eta) * 0.5},
|
|
{(1.0 - xi) * 0.5, (1.0 + xi) * 0.5}};
|
|
}
|
|
|
|
const std::array<Mitc4QuadraturePoint, 8>&
|
|
Mitc4Shell::VolumeQuadrature() noexcept {
|
|
static const std::array<Mitc4QuadraturePoint, 8> points = [] {
|
|
const double gauss = 1.0 / std::sqrt(3.0);
|
|
return std::array<Mitc4QuadraturePoint, 8>{
|
|
Mitc4QuadraturePoint{{-gauss, -gauss, -gauss}, 1.0},
|
|
{{-gauss, -gauss, gauss}, 1.0},
|
|
{{gauss, -gauss, -gauss}, 1.0},
|
|
{{gauss, -gauss, gauss}, 1.0},
|
|
{{gauss, gauss, -gauss}, 1.0},
|
|
{{gauss, gauss, gauss}, 1.0},
|
|
{{-gauss, gauss, -gauss}, 1.0},
|
|
{{-gauss, gauss, gauss}, 1.0}};
|
|
}();
|
|
return points;
|
|
}
|
|
|
|
Mitc4LocalFrame Mitc4Shell::LocalFrame(double xi, double eta) const {
|
|
GeometryData geometry{};
|
|
if (!EvaluateGeometry(xi, eta, 0.0, geometry)) {
|
|
throw std::invalid_argument{
|
|
"MITC4 local frame is invalid at the requested point."};
|
|
}
|
|
return geometry.frame;
|
|
}
|
|
|
|
Matrix Mitc4Shell::PhysicalTransformation20() const {
|
|
Matrix transformation{kPhysicalDofCount, kGlobalDofCount};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const std::size_t physical_offset = node * kPhysicalDofsPerNode;
|
|
const std::size_t global_offset = node * kGlobalDofsPerNode;
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
transformation(physical_offset + component, global_offset + component) =
|
|
1.0;
|
|
transformation(physical_offset + 3U, global_offset + 3U + component) =
|
|
tangent_a_[node][component];
|
|
transformation(physical_offset + 4U, global_offset + 3U + component) =
|
|
tangent_b_[node][component];
|
|
}
|
|
}
|
|
return transformation;
|
|
}
|
|
|
|
Matrix Mitc4Shell::DrillingTransformation4() const {
|
|
Matrix transformation{kNodeCount, kGlobalDofCount};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const std::size_t global_offset = node * kGlobalDofsPerNode;
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
transformation(node, global_offset + 3U + component) =
|
|
directors_[node][component];
|
|
}
|
|
}
|
|
return transformation;
|
|
}
|
|
|
|
Matrix Mitc4Shell::DirectStrainDisplacement20(double xi, double eta,
|
|
double zeta) const {
|
|
return StrainDisplacement(xi, eta, zeta, nullptr);
|
|
}
|
|
|
|
Matrix Mitc4Shell::CovariantTyingShearSamples20() const {
|
|
constexpr std::array<Vector3, 4> points{
|
|
Vector3{0.0, -1.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{-1.0, 0.0, 0.0},
|
|
Vector3{1.0, 0.0, 0.0}};
|
|
Matrix samples{4U, kPhysicalDofCount};
|
|
for (std::size_t point = 0U; point < points.size(); ++point) {
|
|
GeometryData geometry{};
|
|
if (!EvaluateGeometry(points[point][0], points[point][1], 0.0, geometry)) {
|
|
throw std::logic_error{"Validated MITC4 tying geometry became invalid."};
|
|
}
|
|
const auto derivatives =
|
|
BasisDerivatives(points[point][0], points[point][1], 0.0);
|
|
const std::size_t first = point < 2U ? 0U : 1U;
|
|
for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
|
|
const auto strain =
|
|
CovariantStrainColumn(geometry.covariant, derivatives[dof]);
|
|
samples(point, dof) = strain[first][2U];
|
|
}
|
|
}
|
|
return samples;
|
|
}
|
|
|
|
Matrix Mitc4Shell::StrainDisplacement20(double xi, double eta,
|
|
double zeta) const {
|
|
const Matrix samples = CovariantTyingShearSamples20();
|
|
return StrainDisplacement(xi, eta, zeta, &samples);
|
|
}
|
|
|
|
Matrix Mitc4Shell::PlaneStressConstitutive() const {
|
|
const double factor =
|
|
youngs_modulus_ / (1.0 - poisson_ratio_ * poisson_ratio_);
|
|
Matrix constitutive{3U, 3U};
|
|
constitutive(0U, 0U) = factor;
|
|
constitutive(0U, 1U) = factor * poisson_ratio_;
|
|
constitutive(1U, 0U) = factor * poisson_ratio_;
|
|
constitutive(1U, 1U) = factor;
|
|
constitutive(2U, 2U) = factor * (1.0 - poisson_ratio_) * 0.5;
|
|
return constitutive;
|
|
}
|
|
|
|
Matrix Mitc4Shell::MaterialConstitutive5() const {
|
|
Matrix constitutive{5U, 5U};
|
|
const Matrix plane_stress = PlaneStressConstitutive();
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
constitutive(row, column) = plane_stress(row, column);
|
|
}
|
|
}
|
|
const double shear_modulus = youngs_modulus_ / (2.0 * (1.0 + poisson_ratio_));
|
|
constitutive(3U, 3U) = kShearCorrection * shear_modulus;
|
|
constitutive(4U, 4U) = kShearCorrection * shear_modulus;
|
|
return constitutive;
|
|
}
|
|
|
|
Matrix Mitc4Shell::MembraneSectionMatrix() const {
|
|
return ScaledMatrix(PlaneStressConstitutive(), thickness_);
|
|
}
|
|
|
|
Matrix Mitc4Shell::BendingSectionMatrix() const {
|
|
return ScaledMatrix(PlaneStressConstitutive(),
|
|
thickness_ * thickness_ * thickness_ / 12.0);
|
|
}
|
|
|
|
Matrix Mitc4Shell::TransverseShearSectionMatrix() const {
|
|
const double shear_modulus = youngs_modulus_ / (2.0 * (1.0 + poisson_ratio_));
|
|
Matrix result{2U, 2U};
|
|
result(0U, 0U) = kShearCorrection * shear_modulus * thickness_;
|
|
result(1U, 1U) = result(0U, 0U);
|
|
return result;
|
|
}
|
|
|
|
Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const try {
|
|
const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
|
|
const Matrix constitutive = MaterialConstitutive5();
|
|
const Matrix tying_samples = CovariantTyingShearSamples20();
|
|
Matrix physical_local{kPhysicalDofCount, kPhysicalDofCount};
|
|
for (const auto& point : VolumeQuadrature()) {
|
|
GeometryData geometry{};
|
|
if (!EvaluateGeometry(point.natural_coordinates[0],
|
|
point.natural_coordinates[1],
|
|
point.natural_coordinates[2], geometry)) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"Validated MITC4 quadrature geometry became invalid.");
|
|
}
|
|
const Matrix strain = StrainDisplacement(
|
|
point.natural_coordinates[0], point.natural_coordinates[1],
|
|
point.natural_coordinates[2], &tying_samples);
|
|
const Matrix integrand = dense_blas_internal::MultiplyTransposeLeft(
|
|
strain, constitutive.Multiply(strain));
|
|
for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) {
|
|
for (std::size_t column = row; column < kPhysicalDofCount; ++column) {
|
|
const double contribution =
|
|
integrand(row, column) * geometry.jacobian * point.weight;
|
|
physical_local(row, column) += contribution;
|
|
if (row != column) {
|
|
physical_local(column, row) += contribution;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
double maximum_rotary_diagonal = 0.0;
|
|
std::array<double, 8> rotary_diagonals{};
|
|
std::size_t rotary_index = 0U;
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const std::size_t offset = node * kPhysicalDofsPerNode;
|
|
for (std::size_t rotation = 3U; rotation < 5U; ++rotation) {
|
|
const double diagonal =
|
|
physical_local(offset + rotation, offset + rotation);
|
|
if (!std::isfinite(diagonal) || diagonal < 0.0) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"MITC4 drilling stabilization requires finite nonnegative "
|
|
"physical tangent-rotation diagonals.");
|
|
}
|
|
rotary_diagonals[rotary_index++] = diagonal;
|
|
maximum_rotary_diagonal = (std::max)(maximum_rotary_diagonal, diagonal);
|
|
}
|
|
}
|
|
if (!(maximum_rotary_diagonal > 0.0)) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"MITC4 drilling stabilization requires a positive physical "
|
|
"tangent-rotation diagonal inventory.");
|
|
}
|
|
if (!IsFinite(physical_local)) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"MITC4 physical stiffness must contain only finite values.");
|
|
}
|
|
|
|
double normalized_diagonal_sum = 0.0;
|
|
for (const double diagonal : rotary_diagonals) {
|
|
normalized_diagonal_sum += diagonal / maximum_rotary_diagonal;
|
|
}
|
|
const double drilling_stiffness =
|
|
1.0e-3 * maximum_rotary_diagonal * (normalized_diagonal_sum / 8.0);
|
|
if (!std::isfinite(drilling_stiffness) || !(drilling_stiffness > 0.0)) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"MITC4 drilling stiffness must be finite and positive.");
|
|
}
|
|
|
|
Matrix physical_global =
|
|
Congruence(physical_local, PhysicalTransformation20());
|
|
Matrix drilling_local{kNodeCount, kNodeCount};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
drilling_local(node, node) = drilling_stiffness;
|
|
}
|
|
Matrix drilling_global =
|
|
Congruence(drilling_local, DrillingTransformation4());
|
|
Matrix stabilized_global{kGlobalDofCount, kGlobalDofCount};
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
stabilized_global(row, column) =
|
|
physical_global(row, column) + drilling_global(row, column);
|
|
}
|
|
}
|
|
if (!IsFinite(physical_global) || !IsFinite(drilling_global) ||
|
|
!IsFinite(stabilized_global)) {
|
|
return StiffnessFailure(
|
|
source_location_, identity_,
|
|
"MITC4 transformed stiffness must contain only finite values.");
|
|
}
|
|
|
|
return Result<Mitc4Stiffness>::Success(
|
|
Mitc4Stiffness{std::move(physical_local), std::move(physical_global),
|
|
std::move(drilling_global), std::move(stabilized_global),
|
|
drilling_stiffness});
|
|
} catch (const std::exception& error) {
|
|
return StiffnessFailure(source_location_, identity_, error.what());
|
|
}
|
|
|
|
Result<Mitc4Mass> Mitc4Shell::Mass(const double density) const {
|
|
if (!std::isfinite(density) || !(density > 0.0)) {
|
|
return MassFailure(source_location_, identity_,
|
|
"MITC4 mass requires finite density > 0.");
|
|
}
|
|
|
|
std::array<Matrix, 3> moments{Matrix{kNodeCount, kNodeCount},
|
|
Matrix{kNodeCount, kNodeCount},
|
|
Matrix{kNodeCount, kNodeCount}};
|
|
constexpr double kGaussCoordinate = 0.77459666924148337704;
|
|
constexpr std::array<double, 3> kGaussPoints{-kGaussCoordinate, 0.0,
|
|
kGaussCoordinate};
|
|
constexpr std::array<double, 3> kGaussWeights{5.0 / 9.0, 8.0 / 9.0,
|
|
5.0 / 9.0};
|
|
for (std::size_t xi_index = 0U; xi_index < kGaussPoints.size(); ++xi_index) {
|
|
for (std::size_t eta_index = 0U; eta_index < kGaussPoints.size();
|
|
++eta_index) {
|
|
const Mitc4ShapeFunctions shape =
|
|
ShapeFunctions(kGaussPoints[xi_index], kGaussPoints[eta_index]);
|
|
for (std::size_t zeta_index = 0U; zeta_index < kGaussPoints.size();
|
|
++zeta_index) {
|
|
GeometryData geometry{};
|
|
const double zeta = kGaussPoints[zeta_index];
|
|
if (!EvaluateGeometry(kGaussPoints[xi_index], kGaussPoints[eta_index],
|
|
zeta, geometry)) {
|
|
return MassFailure(source_location_, identity_,
|
|
"MITC4 mass quadrature geometry is invalid.");
|
|
}
|
|
const double z = 0.5 * thickness_ * zeta;
|
|
const double weight =
|
|
density * geometry.jacobian * kGaussWeights[xi_index] *
|
|
kGaussWeights[eta_index] * kGaussWeights[zeta_index];
|
|
if (!std::isfinite(z) || !std::isfinite(weight)) {
|
|
return MassFailure(
|
|
source_location_, identity_,
|
|
"MITC4 mass quadrature contribution is nonfinite.");
|
|
}
|
|
const std::array<double, 3> z_powers{1.0, z, z * z};
|
|
for (std::size_t row = 0U; row < kNodeCount; ++row) {
|
|
for (std::size_t column = row; column < kNodeCount; ++column) {
|
|
const double shape_product =
|
|
shape.values[row] * shape.values[column];
|
|
for (std::size_t moment = 0U; moment < moments.size(); ++moment) {
|
|
const double contribution =
|
|
weight * shape_product * z_powers[moment];
|
|
moments[moment](row, column) += contribution;
|
|
if (row != column) {
|
|
moments[moment](column, row) += contribution;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
for (const Matrix& moment : moments) {
|
|
if (!IsFinite(moment)) {
|
|
return MassFailure(source_location_, identity_,
|
|
"MITC4 consistent moments must be finite.");
|
|
}
|
|
}
|
|
|
|
Matrix consistent{kGlobalDofCount, kGlobalDofCount};
|
|
for (std::size_t row_node = 0U; row_node < kNodeCount; ++row_node) {
|
|
const Matrix row_skew = DirectorSkew(directors_[row_node]);
|
|
for (std::size_t column_node = 0U; column_node < kNodeCount;
|
|
++column_node) {
|
|
const Matrix column_skew = DirectorSkew(directors_[column_node]);
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
const std::size_t global_row = kGlobalDofsPerNode * row_node + row;
|
|
const std::size_t global_column =
|
|
kGlobalDofsPerNode * column_node + column;
|
|
consistent(global_row, global_column) =
|
|
moments[0U](row_node, column_node) * (row == column ? 1.0 : 0.0);
|
|
consistent(global_row, global_column + 3U) =
|
|
moments[1U](row_node, column_node) * column_skew(row, column);
|
|
consistent(global_row + 3U, global_column) =
|
|
moments[1U](row_node, column_node) * row_skew(column, row);
|
|
double rotary = 0.0;
|
|
for (std::size_t axis = 0U; axis < 3U; ++axis) {
|
|
rotary += row_skew(axis, row) * column_skew(axis, column);
|
|
}
|
|
consistent(global_row + 3U, global_column + 3U) =
|
|
moments[2U](row_node, column_node) * rotary;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Matrix lumped{kGlobalDofCount, kGlobalDofCount};
|
|
std::array<std::array<double, 3>, kNodeCount> nodal_moments{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
double mass = 0.0;
|
|
double first_moment = 0.0;
|
|
double second_moment = 0.0;
|
|
for (std::size_t column = 0U; column < kNodeCount; ++column) {
|
|
mass += moments[0U](node, column);
|
|
first_moment += moments[1U](node, column);
|
|
second_moment += moments[2U](node, column);
|
|
}
|
|
const double centered_inertia =
|
|
second_moment - first_moment * (first_moment / mass);
|
|
if (!std::isfinite(mass) || !std::isfinite(first_moment) ||
|
|
!std::isfinite(second_moment) || !std::isfinite(centered_inertia) ||
|
|
!(mass > 0.0) || !(second_moment > 0.0) || !(centered_inertia > 0.0)) {
|
|
return MassFailure(
|
|
source_location_, identity_,
|
|
"MITC4 nodal mass moments must be finite and positive.");
|
|
}
|
|
nodal_moments[node] = {mass, first_moment, second_moment};
|
|
const Matrix skew = DirectorSkew(directors_[node]);
|
|
const std::size_t offset = kGlobalDofsPerNode * node;
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
lumped(offset + row, offset + column) =
|
|
mass * (row == column ? 1.0 : 0.0);
|
|
lumped(offset + row, offset + column + 3U) =
|
|
first_moment * skew(row, column);
|
|
lumped(offset + row + 3U, offset + column) =
|
|
first_moment * skew(column, row);
|
|
double projector = row == column ? 1.0 : 0.0;
|
|
projector -= directors_[node][row] * directors_[node][column];
|
|
lumped(offset + row + 3U, offset + column + 3U) =
|
|
second_moment * projector;
|
|
}
|
|
}
|
|
}
|
|
if (!IsFinite(consistent) || !IsFinite(lumped)) {
|
|
return MassFailure(
|
|
source_location_, identity_,
|
|
"MITC4 transformed mass must contain only finite values.");
|
|
}
|
|
return Result<Mitc4Mass>::Success(
|
|
Mitc4Mass{std::move(moments), std::move(consistent), std::move(lumped),
|
|
std::move(nodal_moments)});
|
|
}
|
|
|
|
Result<Mitc4PhysicalRecovery> Mitc4Shell::RecoverPhysical(
|
|
const Vector& global_element_displacement24) const {
|
|
if (global_element_displacement24.Size() != kGlobalDofCount) {
|
|
return RecoveryFailure(
|
|
source_location_, identity_,
|
|
"MITC4 physical recovery requires exactly 24 global element DOFs.");
|
|
}
|
|
if (!IsFinite(global_element_displacement24)) {
|
|
return RecoveryFailure(
|
|
source_location_, identity_,
|
|
"MITC4 physical recovery displacement must be finite.");
|
|
}
|
|
|
|
const Vector physical_displacement =
|
|
PhysicalTransformation20().Multiply(global_element_displacement24);
|
|
const Matrix tying_samples = CovariantTyingShearSamples20();
|
|
const Matrix constitutive = MaterialConstitutive5();
|
|
const Matrix plane_stress = PlaneStressConstitutive();
|
|
const double gauss = 1.0 / std::sqrt(3.0);
|
|
const std::array<std::array<double, 2>, 4> surface_points{
|
|
std::array<double, 2>{-gauss, -gauss},
|
|
std::array<double, 2>{gauss, -gauss}, std::array<double, 2>{gauss, gauss},
|
|
std::array<double, 2>{-gauss, gauss}};
|
|
constexpr std::array<double, 2> thickness_points{-1.0, 1.0};
|
|
constexpr std::array<double, 3> section_positions{-1.0, 0.0, 1.0};
|
|
|
|
Mitc4PhysicalRecovery recovery{};
|
|
for (std::size_t surface = 0U; surface < surface_points.size(); ++surface) {
|
|
auto& point = recovery.points[surface];
|
|
point.natural_coordinates = surface_points[surface];
|
|
GeometryData midsurface_geometry{};
|
|
if (!EvaluateGeometry(point.natural_coordinates[0],
|
|
point.natural_coordinates[1], 0.0,
|
|
midsurface_geometry)) {
|
|
return RecoveryFailure(source_location_, identity_,
|
|
"MITC4 midsurface recovery geometry is invalid.");
|
|
}
|
|
point.local_frame = midsurface_geometry.frame;
|
|
|
|
for (double thickness_sign : thickness_points) {
|
|
const double zeta = thickness_sign * gauss;
|
|
const Matrix strain_matrix = StrainDisplacement(
|
|
point.natural_coordinates[0], point.natural_coordinates[1], zeta,
|
|
&tying_samples);
|
|
const Vector strain = strain_matrix.Multiply(physical_displacement);
|
|
const Vector stress = constitutive.Multiply(strain);
|
|
GeometryData geometry{};
|
|
if (!EvaluateGeometry(point.natural_coordinates[0],
|
|
point.natural_coordinates[1], zeta, geometry)) {
|
|
return RecoveryFailure(
|
|
source_location_, identity_,
|
|
"Validated MITC4 recovery geometry became invalid.");
|
|
}
|
|
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
point.generalized_strain[component] += 0.5 * strain[component];
|
|
point.generalized_strain[3U + component] +=
|
|
3.0 * zeta * strain[component] / thickness_;
|
|
point.section_resultant[component] +=
|
|
0.5 * thickness_ * stress[component];
|
|
point.section_resultant[3U + component] +=
|
|
0.25 * thickness_ * thickness_ * zeta * stress[component];
|
|
}
|
|
for (std::size_t component = 0U; component < 2U; ++component) {
|
|
point.generalized_strain[6U + component] +=
|
|
0.5 * strain[3U + component];
|
|
point.section_resultant[6U + component] +=
|
|
0.5 * thickness_ * stress[3U + component];
|
|
}
|
|
recovery.strain_energy += 0.5 * strain.Dot(stress) * geometry.jacobian;
|
|
}
|
|
|
|
for (std::size_t position = 0U; position < section_positions.size();
|
|
++position) {
|
|
GeometryData section_geometry{};
|
|
if (!EvaluateGeometry(point.natural_coordinates[0],
|
|
point.natural_coordinates[1],
|
|
section_positions[position], section_geometry)) {
|
|
return RecoveryFailure(
|
|
source_location_, identity_,
|
|
"MITC4 section-position recovery geometry is invalid.");
|
|
}
|
|
const Vector strain =
|
|
StrainDisplacement(point.natural_coordinates[0],
|
|
point.natural_coordinates[1],
|
|
section_positions[position], &tying_samples)
|
|
.Multiply(physical_displacement);
|
|
Vector in_plane_strain{3U};
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
in_plane_strain[component] = strain[component];
|
|
}
|
|
const Vector stress = plane_stress.Multiply(in_plane_strain);
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
point.in_plane_stress[position][component] = stress[component];
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!std::isfinite(recovery.strain_energy)) {
|
|
return RecoveryFailure(source_location_, identity_,
|
|
"MITC4 physical strain energy must be finite.");
|
|
}
|
|
for (const auto& point : recovery.points) {
|
|
const auto finite = [](const auto& values) {
|
|
return std::all_of(values.begin(), values.end(),
|
|
[](double value) { return std::isfinite(value); });
|
|
};
|
|
if (!finite(point.generalized_strain) || !finite(point.section_resultant) ||
|
|
std::any_of(
|
|
point.in_plane_stress.begin(), point.in_plane_stress.end(),
|
|
[&finite](const auto& stress) { return !finite(stress); })) {
|
|
return RecoveryFailure(source_location_, identity_,
|
|
"MITC4 physical recovery values must be finite.");
|
|
}
|
|
}
|
|
|
|
return Result<Mitc4PhysicalRecovery>::Success(std::move(recovery));
|
|
}
|
|
|
|
const ElementDofLayout& Mitc4Shell::DofLayout() const noexcept {
|
|
return dof_layout_;
|
|
}
|
|
|
|
Result<ElementStiffnessContribution> Mitc4Shell::ComputeStiffness() const {
|
|
if (dof_layout_.node_indices.size() != kNodeCount ||
|
|
dof_layout_.components_per_node.size() != kGlobalDofsPerNode) {
|
|
return RuntimeFailure<ElementStiffnessContribution>(
|
|
runtime_location_, dof_layout_.source_id,
|
|
"MITC4 runtime stiffness requires a bound four-node six-DOF layout.");
|
|
}
|
|
try {
|
|
auto stiffness = Stiffness();
|
|
if (!stiffness.HasValue()) {
|
|
return Result<ElementStiffnessContribution>::Failure(
|
|
stiffness.GetStatus());
|
|
}
|
|
return Result<ElementStiffnessContribution>::Success(
|
|
{dof_layout_, std::move(stiffness.Value().stabilized_global24)});
|
|
} catch (const std::exception& error) {
|
|
return RuntimeFailure<ElementStiffnessContribution>(
|
|
runtime_location_, dof_layout_.source_id, error.what());
|
|
}
|
|
}
|
|
|
|
Result<ElementResultBundle> Mitc4Shell::Recover(
|
|
const Vector& element_displacement) const {
|
|
if (dof_layout_.node_indices.size() != kNodeCount ||
|
|
dof_layout_.components_per_node.size() != kGlobalDofsPerNode ||
|
|
element_displacement.Size() != kGlobalDofCount) {
|
|
return RuntimeFailure<ElementResultBundle>(
|
|
runtime_location_, dof_layout_.source_id,
|
|
"MITC4 runtime recovery requires a bound four-node layout and 24 "
|
|
"DOFs.");
|
|
}
|
|
|
|
Result<Mitc4PhysicalRecovery> recovered = [&]() {
|
|
try {
|
|
return RecoverPhysical(element_displacement);
|
|
} catch (const std::exception& error) {
|
|
return RuntimeFailure<Mitc4PhysicalRecovery>(
|
|
runtime_location_, dof_layout_.source_id, error.what());
|
|
}
|
|
}();
|
|
if (!recovered.HasValue()) {
|
|
return Result<ElementResultBundle>::Failure(recovered.GetStatus());
|
|
}
|
|
|
|
constexpr std::array<ShellMidsurfaceLocation, kNodeCount> locations{
|
|
ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
|
|
ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
|
|
constexpr std::array<ShellSectionPosition, 3> positions{
|
|
ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
|
|
ShellSectionPosition::kTop};
|
|
constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
|
|
|
|
ShellElementResultRows rows{};
|
|
rows.rows.reserve(kNodeCount);
|
|
rows.physical_strain_energy = recovered.Value().strain_energy;
|
|
for (std::size_t point = 0U; point < recovered.Value().points.size();
|
|
++point) {
|
|
const auto& physical_point = recovered.Value().points[point];
|
|
ShellResultRow row{};
|
|
row.element = element_index_;
|
|
row.location = locations[point];
|
|
row.natural_coordinates = physical_point.natural_coordinates;
|
|
row.local_frame = {physical_point.local_frame.e1,
|
|
physical_point.local_frame.e2,
|
|
physical_point.local_frame.e3};
|
|
row.generalized_strain = physical_point.generalized_strain;
|
|
row.section_resultant = physical_point.section_resultant;
|
|
for (std::size_t position = 0U; position < positions.size(); ++position) {
|
|
row.stress[position] = {positions[position], zeta[position],
|
|
physical_point.in_plane_stress[position]};
|
|
}
|
|
rows.rows.push_back(std::move(row));
|
|
}
|
|
return Result<ElementResultBundle>::Success(
|
|
{dof_layout_.source_id, std::move(rows)});
|
|
}
|
|
|
|
void Mitc4Shell::BindRuntime(ElementDofLayout layout,
|
|
const EntityIndex element_index,
|
|
SourceLocation location) {
|
|
dof_layout_ = std::move(layout);
|
|
element_index_ = element_index;
|
|
runtime_location_ = std::move(location);
|
|
}
|
|
|
|
Mitc4Shell::Mitc4Shell(std::array<Vector3, 4> coordinates,
|
|
std::array<Vector3, 4> directors,
|
|
std::array<Vector3, 4> tangent_a,
|
|
std::array<Vector3, 4> tangent_b,
|
|
Vector3 normal_candidate, double thickness,
|
|
double youngs_modulus, double poisson_ratio,
|
|
SourceLocation source_location, std::string identity)
|
|
: coordinates_{std::move(coordinates)},
|
|
directors_{std::move(directors)},
|
|
tangent_a_{std::move(tangent_a)},
|
|
tangent_b_{std::move(tangent_b)},
|
|
normal_candidate_{std::move(normal_candidate)},
|
|
thickness_{thickness},
|
|
youngs_modulus_{youngs_modulus},
|
|
poisson_ratio_{poisson_ratio},
|
|
drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode},
|
|
source_location_{std::move(source_location)},
|
|
identity_{std::move(identity)} {}
|
|
|
|
bool Mitc4Shell::EvaluateGeometry(double xi, double eta, double zeta,
|
|
GeometryData& result) const noexcept {
|
|
const auto shape = ShapeFunctions(xi, eta);
|
|
const Vector3 midsurface_xi =
|
|
DerivativeSum(shape.xi_derivatives, coordinates_);
|
|
const Vector3 midsurface_eta =
|
|
DerivativeSum(shape.eta_derivatives, coordinates_);
|
|
const Vector3 director_xi = DerivativeSum(shape.xi_derivatives, directors_);
|
|
const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors_);
|
|
const Vector3 director_value = WeightedSum(shape.values, directors_);
|
|
const double half_thickness = 0.5 * thickness_;
|
|
|
|
result.covariant[0] = midsurface_xi + half_thickness * zeta * director_xi;
|
|
result.covariant[1] = midsurface_eta + half_thickness * zeta * director_eta;
|
|
result.covariant[2] = half_thickness * director_value;
|
|
result.jacobian =
|
|
result.covariant[0].Dot(result.covariant[1].Cross(result.covariant[2]));
|
|
if (!result.covariant[0].IsFinite() || !result.covariant[1].IsFinite() ||
|
|
!result.covariant[2].IsFinite() || !std::isfinite(result.jacobian) ||
|
|
!(result.jacobian > 0.0)) {
|
|
return false;
|
|
}
|
|
|
|
if (!std::isfinite(1.0 / result.jacobian) ||
|
|
!result.covariant[1].Cross(result.covariant[2]).IsFinite() ||
|
|
!result.covariant[2].Cross(result.covariant[0]).IsFinite() ||
|
|
!result.covariant[0].Cross(result.covariant[1]).IsFinite()) return false;
|
|
// Preserve the reciprocal finiteness gate, streaming each action directly
|
|
// into validation rather than retaining a reciprocal-basis matrix.
|
|
try {
|
|
const auto factor = dense_blas_internal::DenseFactor::Lu(
|
|
CovariantMatrix(result.covariant));
|
|
if (!factor.HasValue()) return false;
|
|
Matrix rhs{3U, 1U};
|
|
for (std::size_t axis = 0; axis < 3U; ++axis) {
|
|
if (axis != 0U) rhs(axis - 1U, 0U) = 0.0;
|
|
rhs(axis, 0U) = 1.0;
|
|
if (!factor.Value().Solve(rhs, true).HasValue()) return false;
|
|
}
|
|
} catch (const std::exception&) {
|
|
return false;
|
|
}
|
|
|
|
const Vector3 area = midsurface_xi.Cross(midsurface_eta);
|
|
const double director_norm = director_value.Norm();
|
|
if (!area.IsFinite() || !director_value.IsFinite() ||
|
|
!std::isfinite(director_norm) || !(director_norm > 0.0) ||
|
|
!(area.Dot(normal_candidate_) > 0.0)) {
|
|
return false;
|
|
}
|
|
const Vector3 e3 = (1.0 / director_norm) * director_value;
|
|
result.frame.e3 = e3.Components();
|
|
if (!(area.Dot(e3) > 0.0)) {
|
|
return false;
|
|
}
|
|
const Vector3 e1_candidate = midsurface_xi - midsurface_xi.Dot(e3) * e3;
|
|
const double e1_norm = e1_candidate.Norm();
|
|
if (!e1_candidate.IsFinite() || !std::isfinite(e1_norm) || !(e1_norm > 0.0)) {
|
|
return false;
|
|
}
|
|
const Vector3 e1 = (1.0 / e1_norm) * e1_candidate;
|
|
const Vector3 e2 = e3.Cross(e1);
|
|
result.frame.e1 = e1.Components();
|
|
result.frame.e2 = e2.Components();
|
|
return e2.IsFinite();
|
|
}
|
|
|
|
std::array<std::array<Vector3, 3>, 20> Mitc4Shell::BasisDerivatives(
|
|
double xi, double eta, double zeta) const noexcept {
|
|
const auto shape = ShapeFunctions(xi, eta);
|
|
const double half_thickness = 0.5 * thickness_;
|
|
std::array<std::array<Vector3, 3>, kPhysicalDofCount> derivatives{};
|
|
constexpr std::array<Vector3, 3> global_axes{
|
|
Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};
|
|
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const std::size_t offset = node * kPhysicalDofsPerNode;
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
derivatives[offset + component][0] =
|
|
shape.xi_derivatives[node] * global_axes[component];
|
|
derivatives[offset + component][1] =
|
|
shape.eta_derivatives[node] * global_axes[component];
|
|
}
|
|
|
|
const Vector3 alpha_direction = -half_thickness * tangent_b_[node];
|
|
derivatives[offset + 3U][0] =
|
|
zeta * shape.xi_derivatives[node] * alpha_direction;
|
|
derivatives[offset + 3U][1] =
|
|
zeta * shape.eta_derivatives[node] * alpha_direction;
|
|
derivatives[offset + 3U][2] = shape.values[node] * alpha_direction;
|
|
|
|
const Vector3 beta_direction = half_thickness * tangent_a_[node];
|
|
derivatives[offset + 4U][0] =
|
|
zeta * shape.xi_derivatives[node] * beta_direction;
|
|
derivatives[offset + 4U][1] =
|
|
zeta * shape.eta_derivatives[node] * beta_direction;
|
|
derivatives[offset + 4U][2] = shape.values[node] * beta_direction;
|
|
}
|
|
return derivatives;
|
|
}
|
|
|
|
Matrix Mitc4Shell::StrainDisplacement(double xi, double eta, double zeta,
|
|
const Matrix* tying_samples) const {
|
|
GeometryData geometry{};
|
|
if (!EvaluateGeometry(xi, eta, zeta, geometry)) {
|
|
throw std::invalid_argument{
|
|
"MITC4 strain geometry is invalid at the requested point."};
|
|
}
|
|
const auto factor = dense_blas_internal::DenseFactor::Lu(
|
|
CovariantMatrix(geometry.covariant));
|
|
if (!factor.HasValue())
|
|
throw std::invalid_argument{
|
|
"MITC4 strain geometry is invalid at the requested point."};
|
|
const auto derivatives = BasisDerivatives(xi, eta, zeta);
|
|
const Mitc4TyingWeights weights = TyingWeights(xi, eta);
|
|
Matrix result{5U, kPhysicalDofCount};
|
|
for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
|
|
auto covariant =
|
|
CovariantStrainColumn(geometry.covariant, derivatives[dof]);
|
|
if (tying_samples != nullptr) {
|
|
covariant[0U][2U] = weights.xi_zeta[0] * (*tying_samples)(0U, dof) +
|
|
weights.xi_zeta[1] * (*tying_samples)(1U, dof);
|
|
covariant[2U][0U] = covariant[0U][2U];
|
|
covariant[1U][2U] = weights.eta_zeta[0] * (*tying_samples)(2U, dof) +
|
|
weights.eta_zeta[1] * (*tying_samples)(3U, dof);
|
|
covariant[2U][1U] = covariant[1U][2U];
|
|
}
|
|
const auto engineering = LocalEngineeringComponents(
|
|
ReconstructCartesianStrain(covariant, factor.Value()),
|
|
geometry.frame);
|
|
for (std::size_t component = 0U; component < engineering.size();
|
|
++component) {
|
|
result(component, dof) = engineering[component];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
} // namespace fesa
|