KOKO\Mimi
eaa5cbeb21
Integrate Hex8 mass verification and fixed-null-basis shell dynamics with numerical and physics handoff evidence. Preserve interrupted reference and full-validation status; full integration validation is deferred until after performance optimization.
1716 lines
72 KiB
C++
1716 lines
72 KiB
C++
#include "fesa/elements/nonlinear_mitc4_shell.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <limits>
|
|
#include <memory>
|
|
#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> covariant_increment;
|
|
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 Forms (R(phi)-I)D0 without adding and subtracting the baseline.
|
|
Vector3 DirectorIncrement(const Vector3& phi, const Vector3& initial) {
|
|
const double squared_norm = phi.Dot(phi);
|
|
double first = 1.0;
|
|
double second = 0.5;
|
|
if (squared_norm < 0.25) {
|
|
double power = 1.0;
|
|
double first_coefficient = 1.0;
|
|
double second_coefficient = 0.5;
|
|
for (unsigned int k = 1U; k <= 9U; ++k) {
|
|
power *= squared_norm;
|
|
first_coefficient /= -static_cast<double>((2U * k) * (2U * k + 1U));
|
|
second_coefficient /= -static_cast<double>((2U * k + 1U) * (2U * k + 2U));
|
|
first += first_coefficient * power;
|
|
second += second_coefficient * power;
|
|
}
|
|
} else {
|
|
const double norm = std::sqrt(squared_norm);
|
|
first = std::sin(norm) / norm;
|
|
second = (1.0 - std::cos(norm)) / squared_norm;
|
|
}
|
|
const Vector3 crossed = phi.Cross(initial);
|
|
return first * crossed + second * phi.Cross(crossed);
|
|
}
|
|
|
|
/// @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;
|
|
std::array<Vector3, kNodeCount> director_increment{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node)
|
|
director_increment[node] =
|
|
DirectorIncrement(trial.rotation_vector[node], initial_directors[node]);
|
|
geometry.covariant_increment[0U] =
|
|
displacement_xi + half_thickness * zeta *
|
|
DerivativeSum(shape.xi_derivatives, director_increment);
|
|
geometry.covariant_increment[1U] =
|
|
displacement_eta + half_thickness * zeta *
|
|
DerivativeSum(shape.eta_derivatives, director_increment);
|
|
geometry.covariant_increment[2U] =
|
|
half_thickness * WeightedSum(shape.values, director_increment);
|
|
|
|
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.reference_covariant[row].Dot(
|
|
geometry.covariant_increment[column]) +
|
|
geometry.covariant_increment[row].Dot(
|
|
geometry.reference_covariant[column]) +
|
|
geometry.covariant_increment[row].Dot(
|
|
geometry.covariant_increment[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 Distinguishes invalid geometry from arithmetic recovery failure.
|
|
enum class PointStatus { kValid, kInvalidGeometry, kNonfiniteRecovery };
|
|
|
|
/// @brief Evaluates direct strain at one mandatory edge-midpoint tying point.
|
|
PointStatus 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 PointStatus::kInvalidGeometry;
|
|
}
|
|
strain = DirectStrain(geometry);
|
|
return IsFinite(strain) ? PointStatus::kValid : PointStatus::kNonfiniteRecovery;
|
|
}
|
|
|
|
/// @brief Computes one MITC-projected local strain and material stress point.
|
|
PointStatus 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 PointStatus::kInvalidGeometry;
|
|
}
|
|
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)
|
|
? PointStatus::kValid
|
|
: PointStatus::kNonfiniteRecovery;
|
|
}
|
|
|
|
/// @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;
|
|
}
|
|
|
|
constexpr std::size_t kGlobalDofCount = 24U;
|
|
constexpr std::size_t kHessianSize = kGlobalDofCount * kGlobalDofCount;
|
|
|
|
/// @brief Carries exact first and second trial-coordinate variations.
|
|
class Hessian {
|
|
public:
|
|
Hessian() : values_{std::make_shared<std::array<double, kHessianSize>>()} {}
|
|
|
|
double& operator[](const std::size_t index) { return (*values_)[index]; }
|
|
|
|
const double& operator[](const std::size_t index) const {
|
|
return (*values_)[index];
|
|
}
|
|
|
|
auto begin() const { return values_->begin(); }
|
|
|
|
auto end() const { return values_->end(); }
|
|
|
|
private:
|
|
std::shared_ptr<std::array<double, kHessianSize>> values_;
|
|
};
|
|
|
|
/// @brief Carries exact first and second trial-coordinate variations.
|
|
struct Jet {
|
|
double value{};
|
|
std::array<double, kGlobalDofCount> first{};
|
|
Hessian second;
|
|
double directional_first{};
|
|
double directional_second{};
|
|
std::array<double, kGlobalDofCount> mixed_first{};
|
|
std::array<double, kGlobalDofCount> mixed_second{};
|
|
};
|
|
|
|
using JetVector3 = std::array<Jet, 3>;
|
|
using JetMatrix3 = std::array<std::array<Jet, 3>, 3>;
|
|
using JetTensor = std::array<std::array<Jet, 3>, 3>;
|
|
|
|
/// @brief Returns a constant second-order trial-coordinate quantity.
|
|
Jet Constant(const double value) {
|
|
Jet result{};
|
|
result.value = value;
|
|
return result;
|
|
}
|
|
|
|
/// @brief Returns one independent global trial coordinate.
|
|
Jet Variable(const double value, const std::size_t index,
|
|
const double direction = 0.0) {
|
|
Jet result = Constant(value);
|
|
result.first[index] = 1.0;
|
|
result.directional_first = direction;
|
|
return result;
|
|
}
|
|
|
|
/// @brief Adds two second-order trial-coordinate quantities.
|
|
Jet operator+(const Jet& left, const Jet& right) {
|
|
Jet result{};
|
|
result.value = left.value + right.value;
|
|
result.directional_first = left.directional_first + right.directional_first;
|
|
result.directional_second =
|
|
left.directional_second + right.directional_second;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] = left.first[index] + right.first[index];
|
|
result.mixed_first[index] =
|
|
left.mixed_first[index] + right.mixed_first[index];
|
|
result.mixed_second[index] =
|
|
left.mixed_second[index] + right.mixed_second[index];
|
|
}
|
|
for (std::size_t index = 0U; index < kHessianSize; ++index) {
|
|
result.second[index] = left.second[index] + right.second[index];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Subtracts two second-order trial-coordinate quantities.
|
|
Jet operator-(const Jet& left, const Jet& right) {
|
|
Jet result{};
|
|
result.value = left.value - right.value;
|
|
result.directional_first = left.directional_first - right.directional_first;
|
|
result.directional_second =
|
|
left.directional_second - right.directional_second;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] = left.first[index] - right.first[index];
|
|
result.mixed_first[index] =
|
|
left.mixed_first[index] - right.mixed_first[index];
|
|
result.mixed_second[index] =
|
|
left.mixed_second[index] - right.mixed_second[index];
|
|
}
|
|
for (std::size_t index = 0U; index < kHessianSize; ++index) {
|
|
result.second[index] = left.second[index] - right.second[index];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Multiplies two second-order trial-coordinate quantities.
|
|
Jet operator*(const Jet& left, const Jet& right) {
|
|
Jet result{};
|
|
result.value = left.value * right.value;
|
|
result.directional_first = left.directional_first * right.value +
|
|
left.value * right.directional_first;
|
|
result.directional_second =
|
|
left.directional_second * right.value +
|
|
2.0 * left.directional_first * right.directional_first +
|
|
left.value * right.directional_second;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] =
|
|
left.first[index] * right.value + left.value * right.first[index];
|
|
result.mixed_first[index] = left.mixed_first[index] * right.value +
|
|
left.first[index] * right.directional_first +
|
|
left.directional_first * right.first[index] +
|
|
left.value * right.mixed_first[index];
|
|
result.mixed_second[index] =
|
|
left.mixed_second[index] * right.value +
|
|
2.0 * left.mixed_first[index] * right.directional_first +
|
|
left.first[index] * right.directional_second +
|
|
left.directional_second * right.first[index] +
|
|
2.0 * left.directional_first * right.mixed_first[index] +
|
|
left.value * right.mixed_second[index];
|
|
}
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
const std::size_t index = row * kGlobalDofCount + column;
|
|
result.second[index] = left.second[index] * right.value +
|
|
left.first[row] * right.first[column] +
|
|
left.first[column] * right.first[row] +
|
|
left.value * right.second[index];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Scales a second-order trial-coordinate quantity.
|
|
Jet operator*(const double scalar, const Jet& value) {
|
|
Jet result{};
|
|
result.value = scalar * value.value;
|
|
result.directional_first = scalar * value.directional_first;
|
|
result.directional_second = scalar * value.directional_second;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] = scalar * value.first[index];
|
|
result.mixed_first[index] = scalar * value.mixed_first[index];
|
|
result.mixed_second[index] = scalar * value.mixed_second[index];
|
|
}
|
|
for (std::size_t index = 0U; index < kHessianSize; ++index) {
|
|
result.second[index] = scalar * value.second[index];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Scales a second-order trial-coordinate quantity.
|
|
Jet operator*(const Jet& value, const double scalar) { return scalar * value; }
|
|
|
|
/// @brief Applies a scalar unary function through its first two derivatives.
|
|
Jet Compose(const Jet& input, const double value, const double first,
|
|
const double second, const double third) {
|
|
Jet result{};
|
|
result.value = value;
|
|
result.directional_first = first * input.directional_first;
|
|
result.directional_second =
|
|
second * input.directional_first * input.directional_first +
|
|
first * input.directional_second;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] = first * input.first[index];
|
|
result.mixed_first[index] =
|
|
second * input.first[index] * input.directional_first +
|
|
first * input.mixed_first[index];
|
|
result.mixed_second[index] =
|
|
third * input.first[index] * input.directional_first *
|
|
input.directional_first +
|
|
second * (2.0 * input.mixed_first[index] * input.directional_first +
|
|
input.first[index] * input.directional_second) +
|
|
first * input.mixed_second[index];
|
|
}
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
const std::size_t index = row * kGlobalDofCount + column;
|
|
result.second[index] = first * input.second[index] +
|
|
second * input.first[row] * input.first[column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Forms an analytic reciprocal without a secant approximation.
|
|
Jet Reciprocal(const Jet& value) {
|
|
const double inverse = 1.0 / value.value;
|
|
return Compose(value, inverse, -inverse * inverse,
|
|
2.0 * inverse * inverse * inverse,
|
|
-6.0 * inverse * inverse * inverse * inverse);
|
|
}
|
|
|
|
/// @brief Divides two second-order trial-coordinate quantities.
|
|
Jet operator/(const Jet& left, const Jet& right) {
|
|
return left * Reciprocal(right);
|
|
}
|
|
|
|
/// @brief Computes the analytic square root variation.
|
|
Jet Sqrt(const Jet& value) {
|
|
const double root = std::sqrt(value.value);
|
|
return Compose(value, root, 0.5 / root, -0.25 / (root * root * root),
|
|
0.375 / (root * root * root * root * root));
|
|
}
|
|
|
|
/// @brief Computes the analytic sine variation.
|
|
Jet Sin(const Jet& value) {
|
|
return Compose(value, std::sin(value.value), std::cos(value.value),
|
|
-std::sin(value.value), -std::cos(value.value));
|
|
}
|
|
|
|
/// @brief Computes the analytic cosine variation.
|
|
Jet Cos(const Jet& value) {
|
|
return Compose(value, std::cos(value.value), -std::sin(value.value),
|
|
-std::cos(value.value), std::sin(value.value));
|
|
}
|
|
|
|
/// @brief Computes the two-argument angle and its analytic first/second map.
|
|
Jet Atan2(const Jet& y, const Jet& x) {
|
|
const double radius_squared = x.value * x.value + y.value * y.value;
|
|
Jet result{};
|
|
result.value = std::atan2(y.value, x.value);
|
|
const double derivative_y = x.value / radius_squared;
|
|
const double derivative_x = -y.value / radius_squared;
|
|
const double radius_fourth = radius_squared * radius_squared;
|
|
const double derivative_yy = -2.0 * x.value * y.value / radius_fourth;
|
|
const double derivative_yx =
|
|
(y.value * y.value - x.value * x.value) / radius_fourth;
|
|
const double derivative_xx = 2.0 * x.value * y.value / radius_fourth;
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
result.first[index] =
|
|
derivative_y * y.first[index] + derivative_x * x.first[index];
|
|
}
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
const std::size_t index = row * kGlobalDofCount + column;
|
|
result.second[index] = derivative_y * y.second[index] +
|
|
derivative_x * x.second[index] +
|
|
derivative_yy * y.first[row] * y.first[column] +
|
|
derivative_yx * (y.first[row] * x.first[column] +
|
|
x.first[row] * y.first[column]) +
|
|
derivative_xx * x.first[row] * x.first[column];
|
|
}
|
|
}
|
|
// Use a bounded ratio for the directional third derivative only. Retain
|
|
// the established value/gradient/Hessian reduction above for static runs.
|
|
const bool use_x = std::abs(x.value) >= std::abs(y.value);
|
|
const Jet ratio = use_x ? y / x : x / y;
|
|
const double denominator = 1.0 + ratio.value * ratio.value;
|
|
const Jet angle = Compose(ratio, std::atan(ratio.value), 1.0 / denominator,
|
|
-2.0 * ratio.value / (denominator * denominator),
|
|
(6.0 * ratio.value * ratio.value - 2.0) /
|
|
(denominator * denominator * denominator));
|
|
const double sign = use_x ? 1.0 : -1.0;
|
|
result.directional_first = sign * angle.directional_first;
|
|
result.directional_second = sign * angle.directional_second;
|
|
for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
|
|
result.mixed_first[i] = sign * angle.mixed_first[i];
|
|
result.mixed_second[i] = sign * angle.mixed_second[i];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Checks every coefficient in an automatic-differentiation quantity.
|
|
bool IsFinite(const Jet& value) {
|
|
if (!std::isfinite(value.value) || !std::isfinite(value.directional_first) ||
|
|
!std::isfinite(value.directional_second)) {
|
|
return false;
|
|
}
|
|
return std::all_of(value.mixed_first.begin(), value.mixed_first.end(),
|
|
[](const double coefficient) {
|
|
return std::isfinite(coefficient);
|
|
}) &&
|
|
std::all_of(value.mixed_second.begin(), value.mixed_second.end(),
|
|
[](const double coefficient) {
|
|
return std::isfinite(coefficient);
|
|
}) &&
|
|
std::all_of(value.first.begin(), value.first.end(),
|
|
[](const double coefficient) {
|
|
return std::isfinite(coefficient);
|
|
}) &&
|
|
std::all_of(value.second.begin(), value.second.end(),
|
|
[](const double coefficient) {
|
|
return std::isfinite(coefficient);
|
|
});
|
|
}
|
|
|
|
/// @brief Forms a Jet vector from a fixed reference vector.
|
|
JetVector3 ConstantVector(const Vector3& value) {
|
|
return {Constant(value.X()), Constant(value.Y()), Constant(value.Z())};
|
|
}
|
|
|
|
/// @brief Adds two Jet vectors.
|
|
JetVector3 Add(const JetVector3& left, const JetVector3& right) {
|
|
return {left[0U] + right[0U], left[1U] + right[1U], left[2U] + right[2U]};
|
|
}
|
|
|
|
/// @brief Subtracts two Jet vectors.
|
|
JetVector3 Subtract(const JetVector3& left, const JetVector3& right) {
|
|
return {left[0U] - right[0U], left[1U] - right[1U], left[2U] - right[2U]};
|
|
}
|
|
|
|
/// @brief Scales a Jet vector.
|
|
JetVector3 Scale(const Jet& scalar, const JetVector3& vector) {
|
|
return {scalar * vector[0U], scalar * vector[1U], scalar * vector[2U]};
|
|
}
|
|
|
|
/// @brief Scales a Jet vector by a fixed value.
|
|
JetVector3 Scale(const double scalar, const JetVector3& vector) {
|
|
return {scalar * vector[0U], scalar * vector[1U], scalar * vector[2U]};
|
|
}
|
|
|
|
/// @brief Computes the Jet dot product.
|
|
Jet Dot(const JetVector3& left, const JetVector3& right) {
|
|
return left[0U] * right[0U] + left[1U] * right[1U] + left[2U] * right[2U];
|
|
}
|
|
|
|
/// @brief Computes the Jet right-handed cross product.
|
|
JetVector3 Cross(const JetVector3& left, const JetVector3& right) {
|
|
return {left[1U] * right[2U] - left[2U] * right[1U],
|
|
left[2U] * right[0U] - left[0U] * right[2U],
|
|
left[0U] * right[1U] - left[1U] * right[0U]};
|
|
}
|
|
|
|
/// @brief Builds a Jet identity matrix.
|
|
JetMatrix3 Identity() {
|
|
JetMatrix3 result{};
|
|
for (std::size_t index = 0U; index < 3U; ++index) {
|
|
result[index][index] = Constant(1.0);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Forms a skew matrix from a Jet vector.
|
|
JetMatrix3 Skew(const JetVector3& value) {
|
|
JetMatrix3 result{};
|
|
result[0U][1U] = -1.0 * value[2U];
|
|
result[0U][2U] = value[1U];
|
|
result[1U][0U] = value[2U];
|
|
result[1U][2U] = -1.0 * value[0U];
|
|
result[2U][0U] = -1.0 * value[1U];
|
|
result[2U][1U] = value[0U];
|
|
return result;
|
|
}
|
|
|
|
/// @brief Multiplies two three-by-three Jet matrices.
|
|
JetMatrix3 Multiply(const JetMatrix3& left, const JetMatrix3& right) {
|
|
JetMatrix3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
for (std::size_t index = 0U; index < 3U; ++index) {
|
|
result[row][column] =
|
|
result[row][column] + left[row][index] * right[index][column];
|
|
}
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Multiplies a Jet matrix by a Jet vector.
|
|
JetVector3 Multiply(const JetMatrix3& matrix, const JetVector3& vector) {
|
|
JetVector3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row] = result[row] + matrix[row][column] * vector[column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Transposes a three-by-three Jet matrix.
|
|
JetMatrix3 Transpose(const JetMatrix3& matrix) {
|
|
JetMatrix3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] = matrix[column][row];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Returns a scalar multiple of a Jet matrix.
|
|
JetMatrix3 Scale(const double scalar, const JetMatrix3& matrix) {
|
|
JetMatrix3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] = scalar * matrix[row][column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Returns a trial-coordinate scalar multiple of a Jet matrix.
|
|
JetMatrix3 Scale(const Jet& scalar, const JetMatrix3& matrix) {
|
|
JetMatrix3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] = scalar * matrix[row][column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Adds two three-by-three Jet matrices.
|
|
JetMatrix3 Add(const JetMatrix3& left, const JetMatrix3& right) {
|
|
JetMatrix3 result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] = left[row][column] + right[row][column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Computes the analytic three-by-three inverse through its adjugate.
|
|
JetMatrix3 Inverse(const JetMatrix3& matrix) {
|
|
JetMatrix3 cofactor{};
|
|
cofactor[0U][0U] =
|
|
matrix[1U][1U] * matrix[2U][2U] - matrix[1U][2U] * matrix[2U][1U];
|
|
cofactor[0U][1U] =
|
|
matrix[1U][2U] * matrix[2U][0U] - matrix[1U][0U] * matrix[2U][2U];
|
|
cofactor[0U][2U] =
|
|
matrix[1U][0U] * matrix[2U][1U] - matrix[1U][1U] * matrix[2U][0U];
|
|
cofactor[1U][0U] =
|
|
matrix[0U][2U] * matrix[2U][1U] - matrix[0U][1U] * matrix[2U][2U];
|
|
cofactor[1U][1U] =
|
|
matrix[0U][0U] * matrix[2U][2U] - matrix[0U][2U] * matrix[2U][0U];
|
|
cofactor[1U][2U] =
|
|
matrix[0U][1U] * matrix[2U][0U] - matrix[0U][0U] * matrix[2U][1U];
|
|
cofactor[2U][0U] =
|
|
matrix[0U][1U] * matrix[1U][2U] - matrix[0U][2U] * matrix[1U][1U];
|
|
cofactor[2U][1U] =
|
|
matrix[0U][2U] * matrix[1U][0U] - matrix[0U][0U] * matrix[1U][2U];
|
|
cofactor[2U][2U] =
|
|
matrix[0U][0U] * matrix[1U][1U] - matrix[0U][1U] * matrix[1U][0U];
|
|
const Jet determinant = matrix[0U][0U] * cofactor[0U][0U] +
|
|
matrix[0U][1U] * cofactor[0U][1U] +
|
|
matrix[0U][2U] * cofactor[0U][2U];
|
|
JetMatrix3 result{};
|
|
const Jet inverse_determinant = Reciprocal(determinant);
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] = cofactor[column][row] * inverse_determinant;
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Evaluates Rodrigues' map with an analytic near-zero series branch.
|
|
JetMatrix3 ExpSkew(const JetVector3& rotation, const bool increment_only = false) {
|
|
const JetMatrix3 identity = increment_only ? JetMatrix3{} : Identity();
|
|
const Jet squared_norm = Dot(rotation, rotation);
|
|
Jet sine_over_norm{};
|
|
Jet one_minus_cosine_over_norm_squared{};
|
|
if (squared_norm.value < 1.0e-8) {
|
|
const Jet squared = squared_norm * squared_norm;
|
|
const Jet cubed = squared * squared_norm;
|
|
sine_over_norm = Constant(1.0) - squared_norm / Constant(6.0) +
|
|
squared / Constant(120.0) - cubed / Constant(5040.0);
|
|
one_minus_cosine_over_norm_squared =
|
|
Constant(0.5) - squared_norm / Constant(24.0) +
|
|
squared / Constant(720.0) - cubed / Constant(40320.0);
|
|
} else {
|
|
const Jet norm = Sqrt(squared_norm);
|
|
sine_over_norm = Sin(norm) / norm;
|
|
one_minus_cosine_over_norm_squared =
|
|
(Constant(1.0) - Cos(norm)) / squared_norm;
|
|
}
|
|
const JetMatrix3 skew = Skew(rotation);
|
|
JetMatrix3 result =
|
|
Add(Add(identity, Scale(sine_over_norm, skew)),
|
|
Scale(one_minus_cosine_over_norm_squared, Multiply(skew, skew)));
|
|
if (squared_norm.value < 0.25) {
|
|
// Near the original small-angle boundary, cancellation in 1-cos(theta)
|
|
// also corrupts the physical-force gradient and Hessian. Use the existing
|
|
// stable series consistently for the complete Jet and rotation tensor.
|
|
Jet sine_series = Constant(1.0);
|
|
Jet cosine_series = Constant(0.5);
|
|
Jet power = Constant(1.0);
|
|
double sine_coefficient = 1.0;
|
|
double cosine_coefficient = 0.5;
|
|
for (unsigned int k = 1U; k <= 9U; ++k) {
|
|
power = power * squared_norm;
|
|
sine_coefficient /= -static_cast<double>((2U * k) * (2U * k + 1U));
|
|
cosine_coefficient /= -static_cast<double>((2U * k + 1U) * (2U * k + 2U));
|
|
sine_series = sine_series + sine_coefficient * power;
|
|
cosine_series = cosine_series + cosine_coefficient * power;
|
|
}
|
|
const JetMatrix3 stable = Add(Add(identity, Scale(sine_series, skew)),
|
|
Scale(cosine_series, Multiply(skew, skew)));
|
|
result = stable;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Stores the complete Jet-valued global element trial coordinates.
|
|
struct JetTrial {
|
|
std::array<JetVector3, kNodeCount> displacement;
|
|
std::array<JetVector3, kNodeCount> rotation;
|
|
std::array<JetVector3, kNodeCount> director_increment;
|
|
};
|
|
|
|
/// @brief Lifts one validated nonlinear trial into independent Jet variables.
|
|
JetTrial MakeJetTrial(const NonlinearMitc4TrialDofs& trial,
|
|
const std::array<Matrix, kNodeCount>& initial_frames,
|
|
const Vector* direction) {
|
|
JetTrial result{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
result.displacement[node][component] = Variable(
|
|
trial.displacement[node][component], 6U * node + component,
|
|
direction == nullptr ? 0.0 : (*direction)[6U * node + component]);
|
|
result.rotation[node][component] = Variable(
|
|
trial.rotation_vector[node][component], 6U * node + 3U + component,
|
|
direction == nullptr ? 0.0
|
|
: (*direction)[6U * node + 3U + component]);
|
|
}
|
|
JetMatrix3 frame{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
frame[row][column] = Constant(initial_frames[node](row, column));
|
|
}
|
|
}
|
|
const JetVector3 initial_director{frame[0U][2U], frame[1U][2U],
|
|
frame[2U][2U]};
|
|
result.director_increment[node] =
|
|
Multiply(ExpSkew(result.rotation[node], true), initial_director);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Forms covariant increments directly, preserving exact zero motion.
|
|
std::array<JetVector3, 3> CovariantIncrement(
|
|
const std::array<Vector3, kNodeCount>& initial_directors,
|
|
const JetTrial& trial, const double thickness, const double xi,
|
|
const double eta, const double zeta) {
|
|
static_cast<void>(initial_directors);
|
|
const auto shape = Mitc4Shell::ShapeFunctions(xi, eta);
|
|
std::array<JetVector3, 3> increment{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const auto displacement =
|
|
Subtract(trial.displacement[node], trial.displacement[0U]);
|
|
const auto director = trial.director_increment[node];
|
|
const auto layer =
|
|
Add(displacement, Scale(0.5 * thickness * zeta, director));
|
|
increment[0U] =
|
|
Add(increment[0U], Scale(shape.xi_derivatives[node], layer));
|
|
increment[1U] =
|
|
Add(increment[1U], Scale(shape.eta_derivatives[node], layer));
|
|
increment[2U] = Add(increment[2U],
|
|
Scale(0.5 * thickness * shape.values[node], director));
|
|
}
|
|
return increment;
|
|
}
|
|
|
|
/// @brief Forms the direct covariant Green--Lagrange strain tensor.
|
|
JetTensor DirectStrain(const std::array<JetVector3, 3>& increment,
|
|
const Geometry& reference) {
|
|
JetTensor result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result[row][column] =
|
|
0.5 * (Dot(ConstantVector(reference.reference_covariant[row]),
|
|
increment[column]) +
|
|
Dot(increment[row],
|
|
ConstantVector(reference.reference_covariant[column])) +
|
|
Dot(increment[row], increment[column]));
|
|
}
|
|
}
|
|
result[2U][2U] = Constant(0.0);
|
|
return result;
|
|
}
|
|
|
|
/// @brief Reconstructs a fixed-reference Cartesian Jet strain tensor.
|
|
JetTensor CartesianStrain(const JetTensor& covariant,
|
|
const Geometry& reference) {
|
|
JetTensor result{};
|
|
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) {
|
|
result[row][column] =
|
|
result[row][column] + reference.reciprocal[first][row] *
|
|
reference.reciprocal[second][column] *
|
|
covariant[first][second];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Forms one fixed-frame engineering strain component.
|
|
Jet FrameComponent(const JetVector3& left, const JetTensor& tensor,
|
|
const JetVector3& right) {
|
|
Jet result{};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result = result + left[row] * tensor[row][column] * right[column];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Evaluates the physical strain energy through exact Jet variations.
|
|
Jet PhysicalEnergy(const std::array<Vector3, kNodeCount>& coordinates,
|
|
const std::array<Vector3, kNodeCount>& initial_directors,
|
|
const JetTrial& trial, const double thickness,
|
|
const double youngs, const double poisson, bool& valid) {
|
|
const NonlinearMitc4TrialDofs reference_trial{};
|
|
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<JetTensor, 4> tying{};
|
|
for (std::size_t point = 0U; point < tying_locations.size(); ++point) {
|
|
Geometry reference{};
|
|
if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
|
|
reference_trial, thickness,
|
|
tying_locations[point][0U],
|
|
tying_locations[point][1U], 0.0, reference)) {
|
|
valid = false;
|
|
return {};
|
|
}
|
|
tying[point] =
|
|
DirectStrain(CovariantIncrement(initial_directors, trial, thickness,
|
|
tying_locations[point][0U],
|
|
tying_locations[point][1U], 0.0),
|
|
reference);
|
|
}
|
|
|
|
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}};
|
|
const std::array<double, 2> thickness_points{-gauss, gauss};
|
|
const double plane_stress = youngs / (1.0 - poisson * poisson);
|
|
const double shear = kShearCorrection * youngs / (2.0 * (1.0 + poisson));
|
|
Jet energy{};
|
|
for (const auto& surface : surface_points) {
|
|
for (const double zeta : thickness_points) {
|
|
Geometry reference{};
|
|
if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
|
|
reference_trial, thickness, surface[0U],
|
|
surface[1U], zeta, reference)) {
|
|
valid = false;
|
|
return {};
|
|
}
|
|
JetTensor covariant =
|
|
DirectStrain(CovariantIncrement(initial_directors, trial, thickness,
|
|
surface[0U], surface[1U], zeta),
|
|
reference);
|
|
const Mitc4TyingWeights weights =
|
|
Mitc4Shell::TyingWeights(surface[0U], surface[1U]);
|
|
covariant[0U][2U] = weights.xi_zeta[0U] * tying[0U][0U][2U] +
|
|
weights.xi_zeta[1U] * tying[1U][0U][2U];
|
|
covariant[2U][0U] = covariant[0U][2U];
|
|
covariant[1U][2U] = weights.eta_zeta[0U] * tying[2U][1U][2U] +
|
|
weights.eta_zeta[1U] * tying[3U][1U][2U];
|
|
covariant[2U][1U] = covariant[1U][2U];
|
|
const JetTensor cartesian = CartesianStrain(covariant, reference);
|
|
const JetVector3 e1 = ConstantVector(Vector3{reference.local_frame.e1});
|
|
const JetVector3 e2 = ConstantVector(Vector3{reference.local_frame.e2});
|
|
const JetVector3 e3 = ConstantVector(Vector3{reference.local_frame.e3});
|
|
const std::array<Jet, 5> strain{FrameComponent(e1, cartesian, e1),
|
|
FrameComponent(e2, cartesian, e2),
|
|
2.0 * FrameComponent(e1, cartesian, e2),
|
|
2.0 * FrameComponent(e1, cartesian, e3),
|
|
2.0 * FrameComponent(e2, cartesian, e3)};
|
|
const std::array<Jet, 5> stress{
|
|
plane_stress * (strain[0U] + poisson * strain[1U]),
|
|
plane_stress * (poisson * strain[0U] + strain[1U]),
|
|
plane_stress * (1.0 - poisson) * 0.5 * strain[2U], shear * strain[3U],
|
|
shear * strain[4U]};
|
|
Jet density{};
|
|
for (std::size_t component = 0U; component < strain.size(); ++component) {
|
|
density = density + strain[component] * stress[component];
|
|
}
|
|
energy = energy + 0.5 * reference.reference_jacobian * density;
|
|
}
|
|
}
|
|
if (!IsFinite(energy)) {
|
|
valid = false;
|
|
}
|
|
return energy;
|
|
}
|
|
|
|
/// @brief Computes a fixed-iteration analytic polar factor for drill
|
|
/// corotation.
|
|
JetMatrix3 PolarRotation(const JetMatrix3& deformation) {
|
|
JetMatrix3 result = deformation;
|
|
for (std::size_t iteration = 0U; iteration < 12U; ++iteration) {
|
|
result = Scale(0.5, Add(result, Transpose(Inverse(result))));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// @brief Returns the unnormalized physical drill-null chart vector.
|
|
JetVector3 DrillNullPreimage(const JetVector3& phi, const Vector3& director) {
|
|
const Jet square = Dot(phi, phi);
|
|
Jet coefficient{};
|
|
if (square.value < 0.25) {
|
|
constexpr std::array<double, 9> series{1.0 / 12.0,
|
|
1.0 / 720.0,
|
|
1.0 / 30240.0,
|
|
1.0 / 1209600.0,
|
|
1.0 / 47900160.0,
|
|
691.0 / 1307674368000.0,
|
|
1.0 / 74724249600.0,
|
|
3617.0 / 10670622842880000.0,
|
|
43867.0 / 5109094217170944000.0};
|
|
coefficient = Constant(series.back());
|
|
for (std::size_t i = series.size() - 1U; i > 0U; --i)
|
|
coefficient = Constant(series[i - 1U]) + square * coefficient;
|
|
} else {
|
|
const Jet half = 0.5 * Sqrt(square);
|
|
coefficient = (Constant(1.0) - half * Cos(half) * Reciprocal(Sin(half))) *
|
|
Reciprocal(square);
|
|
}
|
|
const auto d0 = ConstantVector(director);
|
|
const auto cross = Cross(phi, d0);
|
|
const auto inverse_action =
|
|
Add(Add(d0, Scale(0.5, cross)), Scale(coefficient, Cross(phi, cross)));
|
|
return inverse_action;
|
|
}
|
|
|
|
/// @brief Returns the normalized physical drill-null chart vector.
|
|
JetVector3 DrillNull(const JetVector3& phi, const Vector3& director) {
|
|
const JetVector3 preimage = DrillNullPreimage(phi, director);
|
|
return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage);
|
|
}
|
|
|
|
/// @brief Returns the speed of physical spin for a unit chart null vector.
|
|
Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
|
|
const JetVector3 preimage = DrillNullPreimage(phi, director);
|
|
return Reciprocal(Sqrt(Dot(preimage, preimage)));
|
|
}
|
|
|
|
/// @brief Evaluates the objective numerical drilling potential and variations.
|
|
Jet DrillingEnergy(const std::array<Vector3, kNodeCount>& coordinates,
|
|
const std::array<Vector3, kNodeCount>& initial_directors,
|
|
const JetTrial& trial, const double thickness,
|
|
const double drilling_stiffness, bool& valid, Matrix* metric,
|
|
std::array<Jet, 4>* constraints) {
|
|
const NonlinearMitc4TrialDofs reference_trial{};
|
|
Geometry reference{};
|
|
if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
|
|
reference_trial, thickness, 0.0, 0.0, 0.0, reference)) {
|
|
valid = false;
|
|
return {};
|
|
}
|
|
const auto current =
|
|
CovariantIncrement(initial_directors, trial, thickness, 0.0, 0.0, 0.0);
|
|
JetMatrix3 deformation = Identity();
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
for (std::size_t alpha = 0U; alpha < 3U; ++alpha) {
|
|
deformation[row][column] =
|
|
deformation[row][column] +
|
|
current[alpha][row] * reference.reciprocal[alpha][column];
|
|
}
|
|
}
|
|
}
|
|
const JetMatrix3 element_rotation = PolarRotation(deformation);
|
|
Jet energy{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
const JetMatrix3 nodal_rotation = ExpSkew(trial.rotation[node]);
|
|
const JetMatrix3 corotated =
|
|
Multiply(Transpose(element_rotation), nodal_rotation);
|
|
const JetVector3 initial_director = ConstantVector(initial_directors[node]);
|
|
const JetVector3 swung_director = Multiply(corotated, initial_director);
|
|
const Jet alignment = Dot(initial_director, swung_director);
|
|
if (!std::isfinite(alignment.value) ||
|
|
!(alignment.value > -1.0 + kDirectorOppositionTolerance)) {
|
|
valid = false;
|
|
return {};
|
|
}
|
|
const JetVector3 swing_axis = Cross(initial_director, swung_director);
|
|
const JetMatrix3 swing_skew = Skew(swing_axis);
|
|
const JetMatrix3 swing = Add(Add(Identity(), swing_skew),
|
|
Scale(Reciprocal(Constant(1.0) + alignment),
|
|
Multiply(swing_skew, swing_skew)));
|
|
const JetMatrix3 twist = Multiply(Transpose(swing), corotated);
|
|
const JetVector3 twist_vex{twist[2U][1U] - twist[1U][2U],
|
|
twist[0U][2U] - twist[2U][0U],
|
|
twist[1U][0U] - twist[0U][1U]};
|
|
const Jet sine = 0.5 * Dot(initial_director, twist_vex);
|
|
const Jet cosine =
|
|
0.5 * (twist[0U][0U] + twist[1U][1U] + twist[2U][2U] - Constant(1.0));
|
|
const Jet gamma = Atan2(sine, cosine);
|
|
if (constraints != nullptr) {
|
|
// n^T grad(Uphys)=0 and Dgamma_j[n_i]=delta_ij/||Jl^-1 d||.
|
|
(*constraints)[node] =
|
|
drilling_stiffness *
|
|
DrillSpinRate(trial.rotation[node], initial_directors[node]) * gamma;
|
|
}
|
|
energy = energy + 0.5 * drilling_stiffness * gamma * gamma;
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
(*metric)(row, column) +=
|
|
drilling_stiffness * gamma.first[row] * gamma.first[column];
|
|
}
|
|
}
|
|
}
|
|
if (!IsFinite(energy)) {
|
|
valid = false;
|
|
}
|
|
return energy;
|
|
}
|
|
|
|
} // 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,
|
|
const double drilling_stiffness, 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},
|
|
drilling_stiffness_{drilling_stiffness},
|
|
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.");
|
|
}
|
|
|
|
std::array<std::array<double, 3>, kNodeCount> director_components{};
|
|
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
|
director_components[node] = directors[node].Components();
|
|
}
|
|
const auto linear_shell =
|
|
Mitc4Shell::Create(nodes, director_components, section, material);
|
|
if (!linear_shell.HasValue()) {
|
|
return Result<NonlinearMitc4Shell>::Failure(linear_shell.GetStatus());
|
|
}
|
|
const auto linear_stiffness = linear_shell.Value().Stiffness();
|
|
if (!linear_stiffness.HasValue() ||
|
|
!std::isfinite(linear_stiffness.Value().drilling_stiffness) ||
|
|
!(linear_stiffness.Value().drilling_stiffness > 0.0)) {
|
|
return Failure<NonlinearMitc4Shell>(
|
|
"invalid-nonlinear-shell-drilling", section.location, identity,
|
|
"Nonlinear MITC4 requires a finite positive fixed drilling stiffness.");
|
|
}
|
|
|
|
NonlinearMitc4Shell shell{coordinates,
|
|
std::move(frames),
|
|
directors,
|
|
section.thickness,
|
|
material.youngs_modulus,
|
|
material.poisson_ratio,
|
|
linear_stiffness.Value().drilling_stiffness,
|
|
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) {
|
|
const auto point_status = EvaluateTyingSample(
|
|
coordinates_, initial_directors_, directors, trial, thickness_,
|
|
tying_locations[point][0U], tying_locations[point][1U],
|
|
tying_samples[point]);
|
|
if (point_status != PointStatus::kValid) {
|
|
return Failure<NonlinearMitc4Recovery>(
|
|
point_status == PointStatus::kInvalidGeometry
|
|
? "invalid-nonlinear-shell-geometry"
|
|
: "invalid-nonlinear-shell-recovery",
|
|
source_location_, identity_,
|
|
"Nonlinear MITC4 tying-point geometry or strain 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{};
|
|
const auto point_status = 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);
|
|
if (point_status != PointStatus::kValid) {
|
|
return Failure<NonlinearMitc4Recovery>(
|
|
point_status == PointStatus::kInvalidGeometry
|
|
? "invalid-nonlinear-shell-geometry"
|
|
: "invalid-nonlinear-shell-recovery",
|
|
source_location_, identity_,
|
|
"Nonlinear MITC4 volume recovery geometry or values are 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_ * 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{};
|
|
const auto point_status = 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);
|
|
if (point_status != PointStatus::kValid) {
|
|
return Failure<NonlinearMitc4Recovery>(
|
|
point_status == PointStatus::kInvalidGeometry
|
|
? "invalid-nonlinear-shell-geometry"
|
|
: "invalid-nonlinear-shell-recovery",
|
|
source_location_, identity_,
|
|
"Nonlinear MITC4 section recovery geometry or values are 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);
|
|
}
|
|
|
|
Result<NonlinearMitc4ExternalDrillingConstraint>
|
|
NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
|
|
const std::size_t node, const Vector3& phi, const Vector3& phi_dot,
|
|
const Vector3& force, const Vector3& force_rate,
|
|
const Vector3& force_acceleration) const {
|
|
if (node >= kNodeCount || !phi_dot.IsFinite() || !force.IsFinite() ||
|
|
!force_rate.IsFinite() || !force_acceleration.IsFinite()) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-external-constraint", source_location_,
|
|
identity_, "External drilling constraint inputs must be finite and "
|
|
"identify one source-local node.");
|
|
}
|
|
const auto rotation = MakeFiniteRotationState(phi, initial_frames_[node]);
|
|
if (!rotation.HasValue()) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-rotation", source_location_, identity_,
|
|
"External drilling constraint rotations must remain in the principal "
|
|
"chart.");
|
|
}
|
|
|
|
JetVector3 jet_phi{};
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
jet_phi[component] = Variable(phi[component], component,
|
|
phi_dot[component]);
|
|
}
|
|
const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]);
|
|
const Jet scalar = -1.0 * Dot(null_direction, ConstantVector(force));
|
|
const Jet scalar_time =
|
|
-1.0 * Dot(null_direction, ConstantVector(force_rate));
|
|
const Jet scalar_acceleration =
|
|
-1.0 * Dot(null_direction, ConstantVector(force_acceleration));
|
|
if (!IsFinite(scalar) || !IsFinite(scalar_time) ||
|
|
!IsFinite(scalar_acceleration)) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-external-constraint", source_location_,
|
|
identity_, "External drilling constraint derivatives must remain "
|
|
"finite.");
|
|
}
|
|
|
|
NonlinearMitc4ExternalDrillingConstraint result{};
|
|
result.value = scalar.value;
|
|
result.gradient = Vector{3U};
|
|
result.hessian = Matrix{3U, 3U};
|
|
result.gradient_directional_first = Vector{3U};
|
|
result.gradient_directional_second = Vector{3U};
|
|
result.directional_first =
|
|
scalar.directional_first + scalar_time.value;
|
|
result.directional_second = scalar.directional_second +
|
|
2.0 * scalar_time.directional_first +
|
|
scalar_acceleration.value;
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
result.gradient[row] = scalar.first[row];
|
|
result.gradient_directional_first[row] =
|
|
scalar.mixed_first[row] + scalar_time.first[row];
|
|
result.gradient_directional_second[row] =
|
|
scalar.mixed_second[row] +
|
|
2.0 * scalar_time.mixed_first[row] + scalar_acceleration.first[row];
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
result.hessian(row, column) =
|
|
scalar.second[row * kGlobalDofCount + column];
|
|
}
|
|
}
|
|
if (!std::isfinite(result.value) || !std::isfinite(result.directional_first) ||
|
|
!std::isfinite(result.directional_second)) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-external-constraint", source_location_,
|
|
identity_, "External drilling constraint derivatives must remain "
|
|
"finite.");
|
|
}
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
if (!std::isfinite(result.gradient[row]) ||
|
|
!std::isfinite(result.gradient_directional_first[row]) ||
|
|
!std::isfinite(result.gradient_directional_second[row])) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-external-constraint", source_location_,
|
|
identity_, "External drilling constraint derivatives must remain "
|
|
"finite.");
|
|
}
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
if (!std::isfinite(result.hessian(row, column))) {
|
|
return Failure<NonlinearMitc4ExternalDrillingConstraint>(
|
|
"invalid-nonlinear-shell-external-constraint", source_location_,
|
|
identity_, "External drilling constraint derivatives must remain "
|
|
"finite.");
|
|
}
|
|
}
|
|
}
|
|
return Result<NonlinearMitc4ExternalDrillingConstraint>::Success(
|
|
std::move(result));
|
|
}
|
|
|
|
Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
|
|
const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial,
|
|
const Vector* coordinate_direction) {
|
|
if (coordinate_direction != nullptr) {
|
|
if (coordinate_direction->Size() != kGlobalDofCount) {
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_, "The directional derivative requires 24 DOFs.");
|
|
}
|
|
for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
|
|
if (!std::isfinite((*coordinate_direction)[i])) {
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_, "The directional derivative must be finite.");
|
|
}
|
|
}
|
|
}
|
|
const auto recovery = element.EvaluateKinematics(trial);
|
|
if (!recovery.HasValue()) {
|
|
return Result<NonlinearMitc4InternalEquilibrium>::Failure(
|
|
recovery.GetStatus());
|
|
}
|
|
|
|
bool valid = true;
|
|
const JetTrial jet_trial =
|
|
MakeJetTrial(trial, element.initial_frames_, coordinate_direction);
|
|
const Jet physical =
|
|
PhysicalEnergy(element.coordinates_, element.initial_directors_,
|
|
jet_trial, element.thickness_, element.youngs_modulus_,
|
|
element.poisson_ratio_, valid);
|
|
Matrix drilling_metric{kGlobalDofCount, kGlobalDofCount};
|
|
std::array<Jet, 4> constraints{};
|
|
const Jet drilling = DrillingEnergy(
|
|
element.coordinates_, element.initial_directors_, jet_trial,
|
|
element.thickness_, element.drilling_stiffness_, valid, &drilling_metric,
|
|
coordinate_direction != nullptr ? &constraints : nullptr);
|
|
const Jet internal_energy = physical + drilling;
|
|
if (!valid || !IsFinite(internal_energy)) {
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_,
|
|
"Nonlinear MITC4 internal residual and tangent must remain finite.");
|
|
}
|
|
|
|
Vector internal_residual{kGlobalDofCount};
|
|
Matrix internal_tangent{kGlobalDofCount, kGlobalDofCount};
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
internal_residual[row] = internal_energy.first[row];
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
internal_tangent(row, column) =
|
|
internal_energy.second[row * kGlobalDofCount + column];
|
|
}
|
|
}
|
|
for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
|
|
if (!std::isfinite(internal_residual[row])) {
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_,
|
|
"Nonlinear MITC4 internal residual and tangent must remain finite.");
|
|
}
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
if (!std::isfinite(internal_tangent(row, column))) {
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_,
|
|
"Nonlinear MITC4 internal residual and tangent must remain "
|
|
"finite.");
|
|
}
|
|
}
|
|
}
|
|
|
|
NonlinearMitc4InternalEquilibrium equilibrium{
|
|
std::move(internal_residual), std::move(internal_tangent),
|
|
Vector{kGlobalDofCount}, recovery.Value(),
|
|
recovery.Value().physical_strain_energy};
|
|
if (coordinate_direction != nullptr) {
|
|
equilibrium.drilling_constraint = Vector{4U};
|
|
equilibrium.drilling_constraint_gradient = Matrix{4U, kGlobalDofCount};
|
|
for (std::size_t node = 0U; node < 4U; ++node) {
|
|
if (!IsFinite(constraints[node]))
|
|
return Failure<NonlinearMitc4InternalEquilibrium>(
|
|
"invalid-nonlinear-shell-equilibrium", element.source_location_,
|
|
element.identity_,
|
|
"Drilling constraint derivatives must remain finite.");
|
|
equilibrium.drilling_constraint[node] = constraints[node].value;
|
|
equilibrium.drilling_constraint_hessian[node] =
|
|
Matrix{kGlobalDofCount, kGlobalDofCount};
|
|
equilibrium.drilling_constraint_third_directional_gradient[node] =
|
|
Vector{kGlobalDofCount};
|
|
for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
|
|
equilibrium.drilling_constraint_gradient(node, i) =
|
|
constraints[node].first[i];
|
|
equilibrium.drilling_constraint_third_directional_gradient[node][i] =
|
|
constraints[node].mixed_second[i];
|
|
for (std::size_t j = 0U; j < kGlobalDofCount; ++j)
|
|
equilibrium.drilling_constraint_hessian[node](i, j) =
|
|
constraints[node].second[i * kGlobalDofCount + j];
|
|
}
|
|
}
|
|
}
|
|
equilibrium.physical_force = Vector{kGlobalDofCount};
|
|
equilibrium.drilling_force = Vector{kGlobalDofCount};
|
|
equilibrium.drilling_tangent = Matrix{kGlobalDofCount, kGlobalDofCount};
|
|
equilibrium.drilling_metric = std::move(drilling_metric);
|
|
equilibrium.drilling_energy = drilling.value;
|
|
equilibrium.force_directional_second = Vector{kGlobalDofCount};
|
|
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
|
equilibrium.internal_force[index] = equilibrium.internal_residual[index];
|
|
equilibrium.physical_force[index] = physical.first[index];
|
|
equilibrium.drilling_force[index] = drilling.first[index];
|
|
equilibrium.force_directional_second[index] =
|
|
internal_energy.mixed_second[index];
|
|
for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
|
|
equilibrium.drilling_tangent(index, column) =
|
|
drilling.second[index * kGlobalDofCount + column];
|
|
}
|
|
}
|
|
return Result<NonlinearMitc4InternalEquilibrium>::Success(
|
|
std::move(equilibrium));
|
|
}
|
|
|
|
} // namespace fesa
|