feat(nonlinear-static-newton-raphson): step 4 - nonlinear-mitc4-residual-and-tangent
This commit is contained in:
KOKO\Mimi
@@ -4,6 +4,8 @@
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <cstddef>
|
||||
#include <limits>
|
||||
#include <memory>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
|
||||
@@ -343,19 +345,675 @@ bool IsFinite(const NonlinearMitc4Recovery& recovery) {
|
||||
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;
|
||||
};
|
||||
|
||||
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) {
|
||||
Jet result = Constant(value);
|
||||
result.first[index] = 1.0;
|
||||
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;
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
result.first[index] = left.first[index] + right.first[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;
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
result.first[index] = left.first[index] - right.first[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;
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
result.first[index] =
|
||||
left.first[index] * right.value + left.value * right.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] = 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;
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
result.first[index] = scalar * value.first[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) {
|
||||
Jet result{};
|
||||
result.value = value;
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
result.first[index] = first * input.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] = 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);
|
||||
}
|
||||
|
||||
/// @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));
|
||||
}
|
||||
|
||||
/// @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));
|
||||
}
|
||||
|
||||
/// @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));
|
||||
}
|
||||
|
||||
/// @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];
|
||||
}
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/// @brief Checks every coefficient in an automatic-differentiation quantity.
|
||||
bool IsFinite(const Jet& value) {
|
||||
if (!std::isfinite(value.value)) {
|
||||
return false;
|
||||
}
|
||||
return 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 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);
|
||||
return Add(Add(Identity(), Scale(sine_over_norm, skew)),
|
||||
Scale(one_minus_cosine_over_norm_squared, Multiply(skew, skew)));
|
||||
}
|
||||
|
||||
/// @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;
|
||||
};
|
||||
|
||||
/// @brief Lifts one validated nonlinear trial into independent Jet variables.
|
||||
JetTrial MakeJetTrial(const NonlinearMitc4TrialDofs& trial,
|
||||
const std::array<Matrix, kNodeCount>& initial_frames) {
|
||||
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);
|
||||
result.rotation[node][component] = Variable(
|
||||
trial.rotation_vector[node][component], 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[node] =
|
||||
Multiply(ExpSkew(result.rotation[node]), initial_director);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/// @brief Forms one current covariant basis from the exact trial variables.
|
||||
std::array<JetVector3, 3> CurrentCovariant(
|
||||
const std::array<Vector3, kNodeCount>& coordinates, const JetTrial& trial,
|
||||
const double thickness, const double xi, const double eta,
|
||||
const double zeta) {
|
||||
const Mitc4ShapeFunctions shape = Mitc4Shell::ShapeFunctions(xi, eta);
|
||||
JetVector3 displacement_xi{};
|
||||
JetVector3 displacement_eta{};
|
||||
JetVector3 director{};
|
||||
JetVector3 director_xi{};
|
||||
JetVector3 director_eta{};
|
||||
Vector3 reference_xi{};
|
||||
Vector3 reference_eta{};
|
||||
for (std::size_t node = 0U; node < kNodeCount; ++node) {
|
||||
reference_xi =
|
||||
reference_xi + shape.xi_derivatives[node] * coordinates[node];
|
||||
reference_eta =
|
||||
reference_eta + shape.eta_derivatives[node] * coordinates[node];
|
||||
displacement_xi = Add(displacement_xi, Scale(shape.xi_derivatives[node],
|
||||
trial.displacement[node]));
|
||||
displacement_eta = Add(displacement_eta, Scale(shape.eta_derivatives[node],
|
||||
trial.displacement[node]));
|
||||
director = Add(director, Scale(shape.values[node], trial.director[node]));
|
||||
director_xi = Add(director_xi,
|
||||
Scale(shape.xi_derivatives[node], trial.director[node]));
|
||||
director_eta = Add(
|
||||
director_eta, Scale(shape.eta_derivatives[node], trial.director[node]));
|
||||
}
|
||||
const double half_thickness = 0.5 * thickness;
|
||||
return {Add(Add(ConstantVector(reference_xi), displacement_xi),
|
||||
Scale(half_thickness * zeta, director_xi)),
|
||||
Add(Add(ConstantVector(reference_eta), displacement_eta),
|
||||
Scale(half_thickness * zeta, director_eta)),
|
||||
Scale(half_thickness, director)};
|
||||
}
|
||||
|
||||
/// @brief Forms the direct covariant Green--Lagrange strain tensor.
|
||||
JetTensor DirectStrain(const std::array<JetVector3, 3>& current,
|
||||
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(current[row], current[column]) -
|
||||
Constant(reference.reference_covariant[row].Dot(
|
||||
reference.reference_covariant[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(CurrentCovariant(coordinates, 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(CurrentCovariant(coordinates, 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 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) {
|
||||
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 =
|
||||
CurrentCovariant(coordinates, trial, thickness, 0.0, 0.0, 0.0);
|
||||
JetMatrix3 deformation{};
|
||||
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);
|
||||
energy = energy + 0.5 * drilling_stiffness * gamma * gamma;
|
||||
}
|
||||
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,
|
||||
SourceLocation source_location, std::string identity)
|
||||
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)} {}
|
||||
|
||||
@@ -448,12 +1106,31 @@ Result<NonlinearMitc4Shell> NonlinearMitc4Shell::Create(
|
||||
"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{});
|
||||
@@ -611,4 +1288,67 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
|
||||
return EvaluateKinematics(trial);
|
||||
}
|
||||
|
||||
Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
|
||||
const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial) {
|
||||
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_);
|
||||
const Jet physical =
|
||||
PhysicalEnergy(element.coordinates_, element.initial_directors_,
|
||||
jet_trial, element.thickness_, element.youngs_modulus_,
|
||||
element.poisson_ratio_, valid);
|
||||
const Jet drilling = DrillingEnergy(
|
||||
element.coordinates_, element.initial_directors_, jet_trial,
|
||||
element.thickness_, element.drilling_stiffness_, valid);
|
||||
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};
|
||||
for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
|
||||
equilibrium.internal_force[index] = equilibrium.internal_residual[index];
|
||||
}
|
||||
return Result<NonlinearMitc4InternalEquilibrium>::Success(
|
||||
std::move(equilibrium));
|
||||
}
|
||||
|
||||
} // namespace fesa
|
||||
|
||||
Reference in New Issue
Block a user