feat(nonlinear-static-newton-raphson): step 3 - nonlinear-mitc4-kinematics-and-recovery
This commit is contained in:
KOKO\Mimi
@@ -0,0 +1,614 @@
|
||||
#include "fesa/elements/nonlinear_mitc4_shell.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <cstddef>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
|
||||
#include "fesa/elements/finite_rotation.h"
|
||||
|
||||
namespace fesa {
|
||||
namespace {
|
||||
|
||||
constexpr std::size_t kNodeCount = 4U;
|
||||
constexpr double kFrameTolerance = 1.0e-12;
|
||||
constexpr double kDirectorOppositionTolerance = 1.0e-8;
|
||||
constexpr double kShearCorrection = 5.0 / 6.0;
|
||||
|
||||
using Tensor = std::array<std::array<double, 3>, 3>;
|
||||
using EngineeringStrain = std::array<double, 5>;
|
||||
using EngineeringStress = std::array<double, 5>;
|
||||
|
||||
struct Geometry {
|
||||
std::array<Vector3, 3> reference_covariant;
|
||||
std::array<Vector3, 3> current_covariant;
|
||||
std::array<Vector3, 3> reciprocal;
|
||||
Mitc4LocalFrame local_frame;
|
||||
double reference_jacobian;
|
||||
};
|
||||
|
||||
struct PointEvaluation {
|
||||
Geometry geometry;
|
||||
EngineeringStrain strain;
|
||||
EngineeringStress stress;
|
||||
};
|
||||
|
||||
/// @brief Builds one deterministic structured nonlinear-shell failure.
|
||||
template <class T>
|
||||
Result<T> Failure(const std::string& code, const SourceLocation& location,
|
||||
const std::string& identity, std::string message) {
|
||||
return Result<T>::Failure(Status::Failure(
|
||||
FailureCategory::kModel,
|
||||
{{Severity::kError, code, location, "*ELEMENT", identity,
|
||||
std::move(message)}}));
|
||||
}
|
||||
|
||||
/// @brief Builds the stable source label concatenation used by diagnostics.
|
||||
std::string ElementIdentity(const std::array<const Node*, kNodeCount>& nodes) {
|
||||
std::string identity;
|
||||
for (const Node* node : nodes) {
|
||||
if (node == nullptr) {
|
||||
continue;
|
||||
}
|
||||
if (!identity.empty()) {
|
||||
identity += "-";
|
||||
}
|
||||
identity += node->source_id.source_label_text;
|
||||
}
|
||||
return identity;
|
||||
}
|
||||
|
||||
/// @brief Forms a weighted vector sum in source-node order.
|
||||
Vector3 WeightedSum(const std::array<double, kNodeCount>& weights,
|
||||
const std::array<Vector3, kNodeCount>& values) {
|
||||
Vector3 result{};
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
result = result + weights[node] * values[node];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/// @brief Forms a stable derivative sum after cancelling the constant term.
|
||||
Vector3 DerivativeSum(const std::array<double, kNodeCount>& derivatives,
|
||||
const std::array<Vector3, kNodeCount>& values) {
|
||||
std::array<Vector3, kNodeCount> relative{};
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
relative[node] = values[node] - values[0U];
|
||||
}
|
||||
return WeightedSum(derivatives, relative);
|
||||
}
|
||||
|
||||
/// @brief Returns a frame matrix whose columns are a, b and the director.
|
||||
Matrix FrameMatrix(const ShellNodeInitialFrame& frame) {
|
||||
Matrix result{3U, 3U};
|
||||
for (std::size_t row = 0U; row < 3U; ++row) {
|
||||
result(row, 0U) = frame.tangent_a[row];
|
||||
result(row, 1U) = frame.tangent_b[row];
|
||||
result(row, 2U) = frame.director[row];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/// @brief Checks every scalar in a tensor before recovery commits it.
|
||||
bool IsFinite(const Tensor& tensor) {
|
||||
for (const auto& row : tensor) {
|
||||
for (const double value : row) {
|
||||
if (!std::isfinite(value)) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/// @brief Checks every scalar in one fixed-size engineering quantity.
|
||||
template <std::size_t N>
|
||||
bool IsFinite(const std::array<double, N>& values) {
|
||||
return std::all_of(values.begin(), values.end(),
|
||||
[](const double value) { return std::isfinite(value); });
|
||||
}
|
||||
|
||||
/// @brief Evaluates a determinant from covariant basis columns.
|
||||
double Determinant(const std::array<Vector3, 3>& covariant) {
|
||||
return covariant[0U].Dot(covariant[1U].Cross(covariant[2U]));
|
||||
}
|
||||
|
||||
/// @brief Builds one proper fixed reference local frame at a midsurface point.
|
||||
bool MakeLocalFrame(const Vector3& midsurface_xi, const Vector3& midsurface_eta,
|
||||
const Vector3& director, Mitc4LocalFrame& frame) {
|
||||
const Vector3 area = midsurface_xi.Cross(midsurface_eta);
|
||||
const double director_norm = director.Norm();
|
||||
if (!area.IsFinite() || !director.IsFinite() ||
|
||||
!(director_norm > kDirectorOppositionTolerance) ||
|
||||
!std::isfinite(director_norm)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
const Vector3 e3 = director / director_norm;
|
||||
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 = e1_candidate / e1_norm;
|
||||
const Vector3 e2 = e3.Cross(e1);
|
||||
if (!e2.IsFinite()) {
|
||||
return false;
|
||||
}
|
||||
frame = {e1.Components(), e2.Components(), e3.Components()};
|
||||
return true;
|
||||
}
|
||||
|
||||
/// @brief Evaluates fixed-reference and trial covariant geometry together.
|
||||
bool EvaluateGeometry(const std::array<Vector3, kNodeCount>& coordinates,
|
||||
const std::array<Vector3, kNodeCount>& initial_directors,
|
||||
const std::array<Vector3, kNodeCount>& directors,
|
||||
const NonlinearMitc4TrialDofs& trial,
|
||||
const double thickness, const double xi, const double eta,
|
||||
const double zeta, Geometry& geometry) {
|
||||
const Mitc4ShapeFunctions shape = Mitc4Shell::ShapeFunctions(xi, eta);
|
||||
const Vector3 reference_xi =
|
||||
DerivativeSum(shape.xi_derivatives, coordinates);
|
||||
const Vector3 reference_eta =
|
||||
DerivativeSum(shape.eta_derivatives, coordinates);
|
||||
const Vector3 initial_director_xi =
|
||||
DerivativeSum(shape.xi_derivatives, initial_directors);
|
||||
const Vector3 initial_director_eta =
|
||||
DerivativeSum(shape.eta_derivatives, initial_directors);
|
||||
const Vector3 initial_director = WeightedSum(shape.values, initial_directors);
|
||||
const Vector3 displacement_xi =
|
||||
DerivativeSum(shape.xi_derivatives, trial.displacement);
|
||||
const Vector3 displacement_eta =
|
||||
DerivativeSum(shape.eta_derivatives, trial.displacement);
|
||||
const Vector3 director_xi = DerivativeSum(shape.xi_derivatives, directors);
|
||||
const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors);
|
||||
const Vector3 director = WeightedSum(shape.values, directors);
|
||||
const double half_thickness = 0.5 * thickness;
|
||||
|
||||
geometry.reference_covariant[0U] =
|
||||
reference_xi + half_thickness * zeta * initial_director_xi;
|
||||
geometry.reference_covariant[1U] =
|
||||
reference_eta + half_thickness * zeta * initial_director_eta;
|
||||
geometry.reference_covariant[2U] = half_thickness * initial_director;
|
||||
geometry.current_covariant[0U] =
|
||||
reference_xi + displacement_xi + half_thickness * zeta * director_xi;
|
||||
geometry.current_covariant[1U] =
|
||||
reference_eta + displacement_eta + half_thickness * zeta * director_eta;
|
||||
geometry.current_covariant[2U] = half_thickness * director;
|
||||
geometry.reference_jacobian = Determinant(geometry.reference_covariant);
|
||||
const double current_jacobian = Determinant(geometry.current_covariant);
|
||||
|
||||
if (!reference_xi.IsFinite() || !reference_eta.IsFinite() ||
|
||||
!initial_director.IsFinite() || !director.IsFinite() ||
|
||||
!(director.Norm() > kDirectorOppositionTolerance) ||
|
||||
!std::isfinite(geometry.reference_jacobian) ||
|
||||
!(geometry.reference_jacobian > 0.0) ||
|
||||
!std::isfinite(current_jacobian) || !(current_jacobian > 0.0)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
geometry.reciprocal[0U] =
|
||||
geometry.reference_covariant[1U].Cross(geometry.reference_covariant[2U]) /
|
||||
geometry.reference_jacobian;
|
||||
geometry.reciprocal[1U] =
|
||||
geometry.reference_covariant[2U].Cross(geometry.reference_covariant[0U]) /
|
||||
geometry.reference_jacobian;
|
||||
geometry.reciprocal[2U] =
|
||||
geometry.reference_covariant[0U].Cross(geometry.reference_covariant[1U]) /
|
||||
geometry.reference_jacobian;
|
||||
if (!geometry.reciprocal[0U].IsFinite() ||
|
||||
!geometry.reciprocal[1U].IsFinite() ||
|
||||
!geometry.reciprocal[2U].IsFinite()) {
|
||||
return false;
|
||||
}
|
||||
return MakeLocalFrame(reference_xi, reference_eta, initial_director,
|
||||
geometry.local_frame);
|
||||
}
|
||||
|
||||
/// @brief Forms the direct Green-Lagrange covariant strain with E33 excluded.
|
||||
Tensor DirectStrain(const Geometry& geometry) {
|
||||
Tensor strain{};
|
||||
for (std::size_t row = 0U; row < 3U; ++row) {
|
||||
for (std::size_t column = 0U; column < 3U; ++column) {
|
||||
strain[row][column] =
|
||||
0.5 * (geometry.current_covariant[row].Dot(
|
||||
geometry.current_covariant[column]) -
|
||||
geometry.reference_covariant[row].Dot(
|
||||
geometry.reference_covariant[column]));
|
||||
}
|
||||
}
|
||||
strain[2U][2U] = 0.0;
|
||||
return strain;
|
||||
}
|
||||
|
||||
/// @brief Reconstructs a fixed-reference Cartesian strain from covariant data.
|
||||
Tensor CartesianStrain(const Tensor& covariant,
|
||||
const std::array<Vector3, 3>& reciprocal) {
|
||||
Tensor cartesian{};
|
||||
for (std::size_t row = 0U; row < 3U; ++row) {
|
||||
for (std::size_t column = 0U; column < 3U; ++column) {
|
||||
for (std::size_t first = 0U; first < 3U; ++first) {
|
||||
for (std::size_t second = 0U; second < 3U; ++second) {
|
||||
cartesian[row][column] += covariant[first][second] *
|
||||
reciprocal[first][row] *
|
||||
reciprocal[second][column];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return cartesian;
|
||||
}
|
||||
|
||||
/// @brief Evaluates a tensor component in a fixed local reference frame.
|
||||
double FrameComponent(const Vector3& left, const Tensor& tensor,
|
||||
const Vector3& right) {
|
||||
double value = 0.0;
|
||||
for (std::size_t row = 0U; row < 3U; ++row) {
|
||||
for (std::size_t column = 0U; column < 3U; ++column) {
|
||||
value += left[row] * tensor[row][column] * right[column];
|
||||
}
|
||||
}
|
||||
return value;
|
||||
}
|
||||
|
||||
/// @brief Converts fixed-reference Cartesian strain to engineering components.
|
||||
EngineeringStrain LocalEngineering(const Tensor& tensor,
|
||||
const Mitc4LocalFrame& frame) {
|
||||
const Vector3 e1{frame.e1};
|
||||
const Vector3 e2{frame.e2};
|
||||
const Vector3 e3{frame.e3};
|
||||
return {FrameComponent(e1, tensor, e1), FrameComponent(e2, tensor, e2),
|
||||
2.0 * FrameComponent(e1, tensor, e2),
|
||||
2.0 * FrameComponent(e1, tensor, e3),
|
||||
2.0 * FrameComponent(e2, tensor, e3)};
|
||||
}
|
||||
|
||||
/// @brief Applies the approved homogeneous plane-stress plus shear law.
|
||||
EngineeringStress Stress(const EngineeringStrain& strain, const double youngs,
|
||||
const double poisson) {
|
||||
const double factor = youngs / (1.0 - poisson * poisson);
|
||||
const double shear = kShearCorrection * youngs / (2.0 * (1.0 + poisson));
|
||||
return {factor * (strain[0U] + poisson * strain[1U]),
|
||||
factor * (poisson * strain[0U] + strain[1U]),
|
||||
factor * (1.0 - poisson) * 0.5 * strain[2U],
|
||||
shear * strain[3U], shear * strain[4U]};
|
||||
}
|
||||
|
||||
/// @brief Evaluates direct strain at one mandatory edge-midpoint tying point.
|
||||
bool EvaluateTyingSample(const std::array<Vector3, kNodeCount>& coordinates,
|
||||
const std::array<Vector3, kNodeCount>& initial_directors,
|
||||
const std::array<Vector3, kNodeCount>& directors,
|
||||
const NonlinearMitc4TrialDofs& trial,
|
||||
const double thickness, const double xi,
|
||||
const double eta, Tensor& strain) {
|
||||
Geometry geometry{};
|
||||
if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
|
||||
thickness, xi, eta, 0.0, geometry)) {
|
||||
return false;
|
||||
}
|
||||
strain = DirectStrain(geometry);
|
||||
return IsFinite(strain);
|
||||
}
|
||||
|
||||
/// @brief Computes one MITC-projected local strain and material stress point.
|
||||
bool EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
|
||||
const std::array<Vector3, kNodeCount>& initial_directors,
|
||||
const std::array<Vector3, kNodeCount>& directors,
|
||||
const NonlinearMitc4TrialDofs& trial,
|
||||
const std::array<Tensor, 4>& tying_samples,
|
||||
const double thickness, const double youngs,
|
||||
const double poisson, const double xi, const double eta,
|
||||
const double zeta, PointEvaluation& point) {
|
||||
if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
|
||||
thickness, xi, eta, zeta, point.geometry)) {
|
||||
return false;
|
||||
}
|
||||
Tensor covariant = DirectStrain(point.geometry);
|
||||
const Mitc4TyingWeights weights = Mitc4Shell::TyingWeights(xi, eta);
|
||||
covariant[0U][2U] = weights.xi_zeta[0U] * tying_samples[0U][0U][2U] +
|
||||
weights.xi_zeta[1U] * tying_samples[1U][0U][2U];
|
||||
covariant[2U][0U] = covariant[0U][2U];
|
||||
covariant[1U][2U] = weights.eta_zeta[0U] * tying_samples[2U][1U][2U] +
|
||||
weights.eta_zeta[1U] * tying_samples[3U][1U][2U];
|
||||
covariant[2U][1U] = covariant[1U][2U];
|
||||
const Tensor cartesian = CartesianStrain(covariant, point.geometry.reciprocal);
|
||||
point.strain = LocalEngineering(cartesian, point.geometry.local_frame);
|
||||
point.stress = Stress(point.strain, youngs, poisson);
|
||||
return IsFinite(covariant) && IsFinite(cartesian) && IsFinite(point.strain) &&
|
||||
IsFinite(point.stress);
|
||||
}
|
||||
|
||||
/// @brief Tests the fixed-order physical recovery inventory for finite values.
|
||||
bool IsFinite(const NonlinearMitc4Recovery& recovery) {
|
||||
if (!std::isfinite(recovery.physical_strain_energy)) {
|
||||
return false;
|
||||
}
|
||||
for (const auto& point : recovery.gauss_points) {
|
||||
if (!IsFinite(point.generalized_strain) ||
|
||||
!IsFinite(point.section_resultant)) {
|
||||
return false;
|
||||
}
|
||||
for (const auto& stress : point.stress) {
|
||||
if (!IsFinite(stress.components)) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
NonlinearMitc4Shell::NonlinearMitc4Shell(
|
||||
std::array<Vector3, 4> coordinates, std::array<Matrix, 4> initial_frames,
|
||||
std::array<Vector3, 4> initial_directors, const double thickness,
|
||||
const double youngs_modulus, const double poisson_ratio,
|
||||
SourceLocation source_location, std::string identity)
|
||||
: coordinates_{std::move(coordinates)},
|
||||
initial_frames_{std::move(initial_frames)},
|
||||
initial_directors_{std::move(initial_directors)},
|
||||
thickness_{thickness},
|
||||
youngs_modulus_{youngs_modulus},
|
||||
poisson_ratio_{poisson_ratio},
|
||||
source_location_{std::move(source_location)},
|
||||
identity_{std::move(identity)} {}
|
||||
|
||||
Result<NonlinearMitc4Shell> NonlinearMitc4Shell::Create(
|
||||
std::array<const Node*, 4> nodes,
|
||||
std::array<ShellNodeInitialFrame, 4> initial_frames,
|
||||
const ShellSection& section, const LinearElasticMaterial& material) {
|
||||
const std::string identity = ElementIdentity(nodes);
|
||||
if (std::any_of(nodes.begin(), nodes.end(),
|
||||
[](const Node* node) { return node == nullptr; })) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-geometry", section.location, identity,
|
||||
"Nonlinear 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 Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-geometry", nodes[node]->location, identity,
|
||||
"Nonlinear MITC4 node coordinates must be finite.");
|
||||
}
|
||||
for (std::size_t previous = 0U; previous < node; ++previous) {
|
||||
if (coordinates[node] == coordinates[previous]) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-geometry", nodes[node]->location,
|
||||
identity, "Nonlinear MITC4 nodes must have distinct coordinates.");
|
||||
}
|
||||
if (initial_frames[node].node_index == initial_frames[previous].node_index) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-frame", nodes[node]->location, identity,
|
||||
"Nonlinear MITC4 initial frames must identify four distinct nodes.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::array<Matrix, kNodeCount> frames{Matrix{3U, 3U}, Matrix{3U, 3U},
|
||||
Matrix{3U, 3U}, Matrix{3U, 3U}};
|
||||
std::array<Vector3, kNodeCount> directors{};
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
frames[node] = FrameMatrix(initial_frames[node]);
|
||||
const auto state = MakeFiniteRotationState(Vector3{}, frames[node]);
|
||||
if (!state.HasValue()) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-frame", nodes[node]->location, identity,
|
||||
"Nonlinear MITC4 initial frames must be finite proper triads.");
|
||||
}
|
||||
directors[node] = state.Value().director;
|
||||
if ((directors[node] - Vector3{initial_frames[node].director}).Norm() >
|
||||
kFrameTolerance) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-frame", nodes[node]->location, identity,
|
||||
"Nonlinear MITC4 initial frame director must be its third axis.");
|
||||
}
|
||||
}
|
||||
|
||||
const Mitc4ShapeFunctions center = Mitc4Shell::ShapeFunctions(0.0, 0.0);
|
||||
const Vector3 center_xi = DerivativeSum(center.xi_derivatives, coordinates);
|
||||
const Vector3 center_eta = DerivativeSum(center.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 Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-geometry", nodes[0U]->location, identity,
|
||||
"Nonlinear MITC4 center surface basis must be finite and nonzero.");
|
||||
}
|
||||
const Vector3 positive_face = center_area / center_measure;
|
||||
if (std::any_of(directors.begin(), directors.end(),
|
||||
[&positive_face](const Vector3& director) {
|
||||
return !(positive_face.Dot(director) > 0.0);
|
||||
})) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-frame", nodes[0U]->location, identity,
|
||||
"Nonlinear MITC4 initial directors must follow the source positive face.");
|
||||
}
|
||||
|
||||
if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
|
||||
return Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-section", section.location, identity,
|
||||
"Nonlinear 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 Failure<NonlinearMitc4Shell>(
|
||||
"invalid-nonlinear-shell-material", material.location, identity,
|
||||
"Nonlinear MITC4 requires finite E>0 and -1<nu<0.5.");
|
||||
}
|
||||
|
||||
NonlinearMitc4Shell shell{coordinates,
|
||||
std::move(frames),
|
||||
directors,
|
||||
section.thickness,
|
||||
material.youngs_modulus,
|
||||
material.poisson_ratio,
|
||||
nodes[0U]->location,
|
||||
identity};
|
||||
const auto reference = shell.EvaluateKinematics(NonlinearMitc4TrialDofs{});
|
||||
if (!reference.HasValue()) {
|
||||
return Result<NonlinearMitc4Shell>::Failure(reference.GetStatus());
|
||||
}
|
||||
return Result<NonlinearMitc4Shell>::Success(std::move(shell));
|
||||
}
|
||||
|
||||
Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(
|
||||
const NonlinearMitc4TrialDofs& trial) const {
|
||||
std::array<Vector3, kNodeCount> directors{};
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
if (!trial.displacement[node].IsFinite() ||
|
||||
!trial.rotation_vector[node].IsFinite()) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-trial", source_location_, identity_,
|
||||
"Nonlinear MITC4 trial displacements and rotations must be finite.");
|
||||
}
|
||||
const auto state =
|
||||
MakeFiniteRotationState(trial.rotation_vector[node], initial_frames_[node]);
|
||||
if (!state.HasValue()) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-rotation", source_location_, identity_,
|
||||
"Nonlinear MITC4 trial rotations must remain in the principal chart.");
|
||||
}
|
||||
directors[node] = state.Value().director;
|
||||
}
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
for (std::size_t other = 0U; other < node; ++other) {
|
||||
if (!(directors[node].Dot(directors[other]) >
|
||||
-1.0 + kDirectorOppositionTolerance)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-director", source_location_, identity_,
|
||||
"Nonlinear MITC4 trial directors must not be opposed.");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Geometry center{};
|
||||
if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
|
||||
thickness_, 0.0, 0.0, 0.0, center)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-geometry", source_location_, identity_,
|
||||
"Nonlinear MITC4 center reference or trial geometry is invalid.");
|
||||
}
|
||||
|
||||
constexpr std::array<std::array<double, 2>, 4> tying_locations{
|
||||
std::array<double, 2>{0.0, -1.0}, std::array<double, 2>{0.0, 1.0},
|
||||
std::array<double, 2>{-1.0, 0.0}, std::array<double, 2>{1.0, 0.0}};
|
||||
std::array<Tensor, 4> tying_samples{};
|
||||
for (std::size_t point = 0U; point < tying_locations.size(); ++point) {
|
||||
if (!EvaluateTyingSample(coordinates_, initial_directors_, directors, trial,
|
||||
thickness_, tying_locations[point][0U],
|
||||
tying_locations[point][1U], tying_samples[point])) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-geometry", source_location_, identity_,
|
||||
"Nonlinear MITC4 tying-point reference or trial geometry is invalid.");
|
||||
}
|
||||
}
|
||||
|
||||
const double gauss = 1.0 / std::sqrt(3.0);
|
||||
constexpr std::array<ShellMidsurfaceLocation, 4> locations{
|
||||
ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
|
||||
ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
|
||||
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_gauss{-1.0, 1.0};
|
||||
constexpr std::array<double, 3> stress_zeta{-1.0, 0.0, 1.0};
|
||||
constexpr std::array<ShellSectionPosition, 3> stress_position{
|
||||
ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
|
||||
ShellSectionPosition::kTop};
|
||||
|
||||
NonlinearMitc4Recovery recovery{};
|
||||
for (std::size_t surface = 0U; surface < surface_points.size(); ++surface) {
|
||||
auto& row = recovery.gauss_points[surface];
|
||||
row.location = locations[surface];
|
||||
row.natural_coordinates = surface_points[surface];
|
||||
Geometry midpoint{};
|
||||
if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
|
||||
thickness_, row.natural_coordinates[0U],
|
||||
row.natural_coordinates[1U], 0.0, midpoint)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-geometry", source_location_, identity_,
|
||||
"Nonlinear MITC4 midsurface recovery geometry is invalid.");
|
||||
}
|
||||
row.local_frame = midpoint.local_frame;
|
||||
|
||||
for (const double thickness_sign : thickness_gauss) {
|
||||
PointEvaluation point{};
|
||||
if (!EvaluatePoint(coordinates_, initial_directors_, directors, trial,
|
||||
tying_samples, thickness_, youngs_modulus_,
|
||||
poisson_ratio_, row.natural_coordinates[0U],
|
||||
row.natural_coordinates[1U], thickness_sign * gauss,
|
||||
point)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-geometry", source_location_, identity_,
|
||||
"Nonlinear MITC4 volume recovery geometry is invalid.");
|
||||
}
|
||||
for (std::size_t component = 0U; component < 3U; ++component) {
|
||||
row.generalized_strain[component] += 0.5 * point.strain[component];
|
||||
row.generalized_strain[3U + component] +=
|
||||
3.0 * thickness_sign * gauss * point.strain[component] /
|
||||
thickness_;
|
||||
row.section_resultant[component] +=
|
||||
0.5 * thickness_ * point.stress[component];
|
||||
row.section_resultant[3U + component] +=
|
||||
0.25 * thickness_ * thickness_sign * gauss * point.stress[component];
|
||||
}
|
||||
for (std::size_t component = 0U; component < 2U; ++component) {
|
||||
row.generalized_strain[6U + component] +=
|
||||
0.5 * point.strain[3U + component];
|
||||
row.section_resultant[6U + component] +=
|
||||
0.5 * thickness_ * point.stress[3U + component];
|
||||
}
|
||||
double density = 0.0;
|
||||
for (std::size_t component = 0U; component < point.strain.size();
|
||||
++component) {
|
||||
density += point.strain[component] * point.stress[component];
|
||||
}
|
||||
recovery.physical_strain_energy +=
|
||||
0.5 * density * point.geometry.reference_jacobian;
|
||||
}
|
||||
|
||||
for (std::size_t position = 0U; position < stress_zeta.size(); ++position) {
|
||||
PointEvaluation point{};
|
||||
if (!EvaluatePoint(coordinates_, initial_directors_, directors, trial,
|
||||
tying_samples, thickness_, youngs_modulus_,
|
||||
poisson_ratio_, row.natural_coordinates[0U],
|
||||
row.natural_coordinates[1U], stress_zeta[position],
|
||||
point)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-geometry", source_location_, identity_,
|
||||
"Nonlinear MITC4 section-position recovery geometry is invalid.");
|
||||
}
|
||||
row.stress[position] = {stress_position[position], stress_zeta[position],
|
||||
{point.stress[0U], point.stress[1U],
|
||||
point.stress[2U]}};
|
||||
}
|
||||
}
|
||||
|
||||
if (!IsFinite(recovery)) {
|
||||
return Failure<NonlinearMitc4Recovery>(
|
||||
"invalid-nonlinear-shell-recovery", source_location_, identity_,
|
||||
"Nonlinear MITC4 physical recovery must remain finite.");
|
||||
}
|
||||
return Result<NonlinearMitc4Recovery>::Success(std::move(recovery));
|
||||
}
|
||||
|
||||
Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
|
||||
const NonlinearMitc4TrialDofs& trial) const {
|
||||
return EvaluateKinematics(trial);
|
||||
}
|
||||
|
||||
} // namespace fesa
|
||||
Reference in New Issue
Block a user