KOKO\Mimi
977 lines
41 KiB
C++
977 lines
41 KiB
C++
#include "fesa/elements/mitc4_shell.h"
|
|
|
|
#include <gtest/gtest.h>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <cstdint>
|
|
#include <limits>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <vector>
|
|
|
|
#include "fesa/math/vector3.h"
|
|
#include "fesa/model/model_types.h"
|
|
|
|
namespace {
|
|
|
|
using Vector3 = std::array<double, 3>;
|
|
|
|
fesa::Node Node(std::int64_t label, Vector3 coordinates) {
|
|
return {{"Shell-Instance", label, std::to_string(label)},
|
|
coordinates,
|
|
{"mitc4-shell.inp", static_cast<std::size_t>(label + 1)}};
|
|
}
|
|
|
|
std::array<const fesa::Node*, 4> NodePointers(
|
|
const std::array<fesa::Node, 4>& nodes) {
|
|
return {&nodes[0], &nodes[1], &nodes[2], &nodes[3]};
|
|
}
|
|
|
|
fesa::ShellSection Section(double thickness = 2.0) {
|
|
return {"Section", thickness, 0U, {"mitc4-shell.inp", 20U}};
|
|
}
|
|
|
|
fesa::LinearElasticMaterial Material(double youngs_modulus = 120.0,
|
|
double poisson_ratio = 0.25) {
|
|
return {"Material", youngs_modulus, poisson_ratio, {"mitc4-shell.inp", 21U}};
|
|
}
|
|
|
|
std::array<Vector3, 4> Directors(Vector3 director = {0.0, 0.0, 1.0}) {
|
|
return {director, director, director, director};
|
|
}
|
|
|
|
double Dot(const Vector3& left, const Vector3& right) {
|
|
return left[0] * right[0] + left[1] * right[1] + left[2] * right[2];
|
|
}
|
|
|
|
Vector3 Cross(const Vector3& left, const Vector3& right) {
|
|
return {left[1] * right[2] - left[2] * right[1],
|
|
left[2] * right[0] - left[0] * right[2],
|
|
left[0] * right[1] - left[1] * right[0]};
|
|
}
|
|
|
|
double Norm(const Vector3& value) { return std::sqrt(Dot(value, value)); }
|
|
|
|
void ExpectVectorNear(const Vector3& actual, const Vector3& expected,
|
|
double tolerance = 1.0e-12) {
|
|
for (std::size_t component = 0U; component < actual.size(); ++component) {
|
|
EXPECT_NEAR(actual[component], expected[component], tolerance);
|
|
}
|
|
}
|
|
|
|
void ExpectOrthonormalRightHanded(const fesa::Mitc4LocalFrame& frame) {
|
|
EXPECT_NEAR(Norm(frame.e1), 1.0, 1.0e-12);
|
|
EXPECT_NEAR(Norm(frame.e2), 1.0, 1.0e-12);
|
|
EXPECT_NEAR(Norm(frame.e3), 1.0, 1.0e-12);
|
|
EXPECT_NEAR(Dot(frame.e1, frame.e2), 0.0, 1.0e-12);
|
|
EXPECT_NEAR(Dot(frame.e1, frame.e3), 0.0, 1.0e-12);
|
|
EXPECT_NEAR(Dot(frame.e2, frame.e3), 0.0, 1.0e-12);
|
|
ExpectVectorNear(Cross(frame.e1, frame.e2), frame.e3);
|
|
}
|
|
|
|
void ExpectMatrixNear(const fesa::Matrix& actual, const fesa::Matrix& expected,
|
|
double tolerance = 1.0e-12) {
|
|
ASSERT_EQ(actual.Rows(), expected.Rows());
|
|
ASSERT_EQ(actual.Columns(), expected.Columns());
|
|
for (std::size_t row = 0U; row < actual.Rows(); ++row) {
|
|
for (std::size_t column = 0U; column < actual.Columns(); ++column) {
|
|
EXPECT_NEAR(actual(row, column), expected(row, column), tolerance)
|
|
<< "at (" << row << ", " << column << ")";
|
|
}
|
|
}
|
|
}
|
|
|
|
void ExpectSymmetric(const fesa::Matrix& matrix) {
|
|
ASSERT_EQ(matrix.Rows(), matrix.Columns());
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
|
|
EXPECT_NEAR(matrix(row, column), matrix(column, row), 1.0e-12);
|
|
}
|
|
}
|
|
}
|
|
|
|
bool HasPositiveCholeskyPivots(const fesa::Matrix& matrix) {
|
|
if (matrix.Rows() != matrix.Columns()) {
|
|
return false;
|
|
}
|
|
fesa::Matrix lower{matrix.Rows(), matrix.Columns()};
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
for (std::size_t column = 0U; column <= row; ++column) {
|
|
double value = matrix(row, column);
|
|
for (std::size_t inner = 0U; inner < column; ++inner) {
|
|
value -= lower(row, inner) * lower(column, inner);
|
|
}
|
|
if (row == column) {
|
|
if (!std::isfinite(value) || !(value > 0.0)) {
|
|
return false;
|
|
}
|
|
lower(row, column) = std::sqrt(value);
|
|
} else {
|
|
lower(row, column) = value / lower(column, column);
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
double FrobeniusNorm(const fesa::Matrix& matrix) {
|
|
double squared_norm = 0.0;
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
|
|
squared_norm += matrix(row, column) * matrix(row, column);
|
|
}
|
|
}
|
|
return std::sqrt(squared_norm);
|
|
}
|
|
|
|
double ScaledSymmetryError(const fesa::Matrix& matrix,
|
|
std::size_t dofs_per_node, double element_length) {
|
|
fesa::Matrix difference{matrix.Rows(), matrix.Columns()};
|
|
fesa::Matrix scaled{matrix.Rows(), matrix.Columns()};
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
const double row_scale = row % dofs_per_node < 3U ? element_length : 1.0;
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
|
|
const double column_scale =
|
|
column % dofs_per_node < 3U ? element_length : 1.0;
|
|
scaled(row, column) = row_scale * matrix(row, column) * column_scale;
|
|
difference(row, column) = row_scale *
|
|
(matrix(row, column) - matrix(column, row)) *
|
|
column_scale;
|
|
}
|
|
}
|
|
return FrobeniusNorm(difference) / FrobeniusNorm(scaled);
|
|
}
|
|
|
|
fesa::Matrix ScaledStiffness(const fesa::Matrix& matrix,
|
|
std::size_t dofs_per_node, double element_length) {
|
|
fesa::Matrix scaled{matrix.Rows(), matrix.Columns()};
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
const double row_scale = row % dofs_per_node < 3U ? element_length : 1.0;
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
|
|
const double column_scale =
|
|
column % dofs_per_node < 3U ? element_length : 1.0;
|
|
scaled(row, column) = row_scale * matrix(row, column) * column_scale;
|
|
}
|
|
}
|
|
return scaled;
|
|
}
|
|
|
|
std::vector<double> SymmetricEigenvalues(fesa::Matrix matrix) {
|
|
if (matrix.Rows() != matrix.Columns()) {
|
|
throw std::invalid_argument{
|
|
"Symmetric eigensolve requires a square matrix."};
|
|
}
|
|
const std::size_t size = matrix.Rows();
|
|
double matrix_scale = 0.0;
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
for (std::size_t column = 0U; column < size; ++column) {
|
|
matrix_scale = (std::max)(matrix_scale, std::abs(matrix(row, column)));
|
|
}
|
|
}
|
|
if (matrix_scale != 0.0) {
|
|
const double convergence_tolerance = 1.0e-14 * matrix_scale;
|
|
const std::size_t iteration_limit = 100U * size * size;
|
|
for (std::size_t iteration = 0U; iteration < iteration_limit; ++iteration) {
|
|
std::size_t pivot_row = 0U;
|
|
std::size_t pivot_column = 0U;
|
|
double largest_off_diagonal = 0.0;
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
for (std::size_t column = row + 1U; column < size; ++column) {
|
|
const double candidate = std::abs(matrix(row, column));
|
|
if (candidate > largest_off_diagonal) {
|
|
largest_off_diagonal = candidate;
|
|
pivot_row = row;
|
|
pivot_column = column;
|
|
}
|
|
}
|
|
}
|
|
if (largest_off_diagonal <= convergence_tolerance) {
|
|
break;
|
|
}
|
|
|
|
const double pivot = matrix(pivot_row, pivot_column);
|
|
const double tau =
|
|
(matrix(pivot_column, pivot_column) - matrix(pivot_row, pivot_row)) /
|
|
(2.0 * pivot);
|
|
const double tangent = tau >= 0.0
|
|
? 1.0 / (tau + std::sqrt(1.0 + tau * tau))
|
|
: -1.0 / (-tau + std::sqrt(1.0 + tau * tau));
|
|
const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent);
|
|
const double sine = tangent * cosine;
|
|
const double row_diagonal = matrix(pivot_row, pivot_row);
|
|
const double column_diagonal = matrix(pivot_column, pivot_column);
|
|
matrix(pivot_row, pivot_row) = row_diagonal - tangent * pivot;
|
|
matrix(pivot_column, pivot_column) = column_diagonal + tangent * pivot;
|
|
matrix(pivot_row, pivot_column) = 0.0;
|
|
matrix(pivot_column, pivot_row) = 0.0;
|
|
for (std::size_t index = 0U; index < size; ++index) {
|
|
if (index == pivot_row || index == pivot_column) {
|
|
continue;
|
|
}
|
|
const double row_value = matrix(index, pivot_row);
|
|
const double column_value = matrix(index, pivot_column);
|
|
const double rotated_row = cosine * row_value - sine * column_value;
|
|
const double rotated_column = sine * row_value + cosine * column_value;
|
|
matrix(index, pivot_row) = rotated_row;
|
|
matrix(pivot_row, index) = rotated_row;
|
|
matrix(index, pivot_column) = rotated_column;
|
|
matrix(pivot_column, index) = rotated_column;
|
|
}
|
|
}
|
|
}
|
|
|
|
std::vector<double> eigenvalues(size);
|
|
for (std::size_t index = 0U; index < size; ++index) {
|
|
eigenvalues[index] = matrix(index, index);
|
|
}
|
|
return eigenvalues;
|
|
}
|
|
|
|
std::size_t NumericalRank(const fesa::Matrix& scaled_matrix) {
|
|
const auto eigenvalues = SymmetricEigenvalues(scaled_matrix);
|
|
double spectral_scale = 0.0;
|
|
for (double eigenvalue : eigenvalues) {
|
|
spectral_scale = (std::max)(spectral_scale, std::abs(eigenvalue));
|
|
}
|
|
return static_cast<std::size_t>(
|
|
std::count_if(eigenvalues.begin(), eigenvalues.end(),
|
|
[spectral_scale](double eigenvalue) {
|
|
return std::abs(eigenvalue) > 1.0e-9 * spectral_scale;
|
|
}));
|
|
}
|
|
|
|
double SymmetricOperatorNorm(const fesa::Matrix& matrix) {
|
|
const auto eigenvalues = SymmetricEigenvalues(matrix);
|
|
double result = 0.0;
|
|
for (double eigenvalue : eigenvalues) {
|
|
result = (std::max)(result, std::abs(eigenvalue));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
double QuadraticEnergy(const fesa::Matrix& stiffness,
|
|
const fesa::Vector& vector) {
|
|
return 0.5 * vector.Dot(stiffness.Multiply(vector));
|
|
}
|
|
|
|
fesa::Vector PhysicalField(const std::array<std::array<double, 5>, 4>& values) {
|
|
fesa::Vector result{20U};
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
for (std::size_t component = 0U; component < 5U; ++component) {
|
|
result[5U * node_index + component] = values[node_index][component];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
void ExpectStrain(const fesa::Mitc4Shell& shell, const fesa::Vector& field,
|
|
double xi, double eta, double zeta,
|
|
const std::array<double, 5>& expected) {
|
|
const auto actual = shell.StrainDisplacement20(xi, eta, zeta).Multiply(field);
|
|
for (std::size_t component = 0U; component < expected.size(); ++component) {
|
|
EXPECT_NEAR(actual[component], expected[component], 1.0e-12)
|
|
<< "component " << component;
|
|
}
|
|
}
|
|
|
|
fesa::Vector PhysicalRigidMode(const std::array<fesa::Node, 4>& nodes,
|
|
const Vector3& translation,
|
|
const Vector3& rotation) {
|
|
fesa::Vector mode{24U};
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
const Vector3 rotational_translation =
|
|
Cross(rotation, nodes[node_index].coordinates);
|
|
const std::size_t offset = 6U * node_index;
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
mode[offset + component] =
|
|
translation[component] + rotational_translation[component];
|
|
}
|
|
// Remove the director-parallel component: it is numerical drilling,
|
|
// not part of the five-DOF physical rigid motion.
|
|
mode[offset + 3U] = rotation[0];
|
|
mode[offset + 4U] = rotation[1];
|
|
mode[offset + 5U] = 0.0;
|
|
}
|
|
return mode;
|
|
}
|
|
|
|
std::array<fesa::Node, 4> PlanarNodes() {
|
|
return {Node(1, {-1.0, -1.0, 0.0}), Node(2, {1.0, -1.0, 0.0}),
|
|
Node(3, {1.0, 1.0, 0.0}), Node(4, {-1.0, 1.0, 0.0})};
|
|
}
|
|
|
|
} // namespace
|
|
|
|
// MITC4-KIN-001
|
|
TEST(Mitc4ShellKinematics, ShapeFunctionsSatisfyNodalAndDerivativeIdentities) {
|
|
constexpr std::array<Vector3, 4> kNaturalNodes{
|
|
Vector3{-1.0, -1.0, 0.0}, Vector3{1.0, -1.0, 0.0}, Vector3{1.0, 1.0, 0.0},
|
|
Vector3{-1.0, 1.0, 0.0}};
|
|
|
|
for (std::size_t point = 0U; point < kNaturalNodes.size(); ++point) {
|
|
const auto shape = fesa::Mitc4Shell::ShapeFunctions(
|
|
kNaturalNodes[point][0], kNaturalNodes[point][1]);
|
|
for (std::size_t node_index = 0U; node_index < kNaturalNodes.size();
|
|
++node_index) {
|
|
EXPECT_DOUBLE_EQ(shape.values[node_index],
|
|
point == node_index ? 1.0 : 0.0);
|
|
}
|
|
}
|
|
|
|
const auto shape = fesa::Mitc4Shell::ShapeFunctions(0.25, -0.5);
|
|
double value_sum = 0.0;
|
|
double xi_derivative_sum = 0.0;
|
|
double eta_derivative_sum = 0.0;
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
value_sum += shape.values[node_index];
|
|
xi_derivative_sum += shape.xi_derivatives[node_index];
|
|
eta_derivative_sum += shape.eta_derivatives[node_index];
|
|
}
|
|
EXPECT_DOUBLE_EQ(value_sum, 1.0);
|
|
EXPECT_DOUBLE_EQ(xi_derivative_sum, 0.0);
|
|
EXPECT_DOUBLE_EQ(eta_derivative_sum, 0.0);
|
|
EXPECT_EQ(shape.values,
|
|
(std::array<double, 4>{0.28125, 0.46875, 0.15625, 0.09375}));
|
|
}
|
|
|
|
// MITC4-KIN-002
|
|
TEST(Mitc4ShellKinematics,
|
|
BuildsRightHandedFramesAndSeparatePhysicalDrillingMaps) {
|
|
const std::array<fesa::Node, 4> nodes{
|
|
Node(1, {0.0, -1.0, -1.0}), Node(2, {0.0, 1.0, -1.0}),
|
|
Node(3, {0.0, 1.0, 1.0}), Node(4, {0.0, -1.0, 1.0})};
|
|
const auto candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors({1.0, 0.0, 0.0}), Section(), Material());
|
|
ASSERT_TRUE(candidate.HasValue());
|
|
const auto& shell = candidate.Value();
|
|
|
|
const auto frame = shell.LocalFrame(0.0, 0.0);
|
|
ExpectVectorNear(frame.e1, {0.0, 1.0, 0.0});
|
|
ExpectVectorNear(frame.e2, {0.0, 0.0, 1.0});
|
|
ExpectVectorNear(frame.e3, {1.0, 0.0, 0.0});
|
|
ExpectOrthonormalRightHanded(frame);
|
|
|
|
const auto physical = shell.PhysicalTransformation20();
|
|
const auto drilling = shell.DrillingTransformation4();
|
|
ASSERT_EQ(physical.Rows(), 20U);
|
|
ASSERT_EQ(physical.Columns(), 24U);
|
|
ASSERT_EQ(drilling.Rows(), 4U);
|
|
ASSERT_EQ(drilling.Columns(), 24U);
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
const std::size_t physical_offset = 5U * node_index;
|
|
const std::size_t global_offset = 6U * node_index;
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
EXPECT_DOUBLE_EQ(
|
|
physical(physical_offset + component, global_offset + component),
|
|
1.0);
|
|
}
|
|
EXPECT_DOUBLE_EQ(physical(physical_offset + 3U, global_offset + 4U), 1.0);
|
|
EXPECT_DOUBLE_EQ(physical(physical_offset + 4U, global_offset + 5U), 1.0);
|
|
EXPECT_DOUBLE_EQ(drilling(node_index, global_offset + 3U), 1.0);
|
|
|
|
for (std::size_t global_dof = 0U; global_dof < 24U; ++global_dof) {
|
|
if (global_dof != global_offset + 4U) {
|
|
EXPECT_DOUBLE_EQ(physical(physical_offset + 3U, global_dof), 0.0);
|
|
}
|
|
if (global_dof != global_offset + 5U) {
|
|
EXPECT_DOUBLE_EQ(physical(physical_offset + 4U, global_dof), 0.0);
|
|
}
|
|
if (global_dof != global_offset + 3U) {
|
|
EXPECT_DOUBLE_EQ(drilling(node_index, global_dof), 0.0);
|
|
}
|
|
}
|
|
}
|
|
|
|
auto invalid_directors = Directors({1.0, 0.0, 0.0});
|
|
invalid_directors[2] = {0.0, 0.0, 0.0};
|
|
EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), invalid_directors,
|
|
Section(), Material())
|
|
.HasValue());
|
|
}
|
|
|
|
// MITC4-KIN-003
|
|
TEST(Mitc4ShellKinematics, FormsDirectColumnsAndAllCovariantTyingSamples) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(candidate.HasValue());
|
|
const auto& shell = candidate.Value();
|
|
|
|
const auto direct = shell.DirectStrainDisplacement20(0.0, 0.0, 0.5);
|
|
ASSERT_EQ(direct.Rows(), 5U);
|
|
ASSERT_EQ(direct.Columns(), 20U);
|
|
EXPECT_DOUBLE_EQ(direct(0U, 0U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(0U, 4U), -0.125);
|
|
EXPECT_DOUBLE_EQ(direct(1U, 1U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(1U, 3U), 0.125);
|
|
EXPECT_DOUBLE_EQ(direct(2U, 0U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(2U, 1U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(2U, 3U), 0.125);
|
|
EXPECT_DOUBLE_EQ(direct(2U, 4U), -0.125);
|
|
EXPECT_DOUBLE_EQ(direct(3U, 2U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(3U, 4U), 0.25);
|
|
EXPECT_DOUBLE_EQ(direct(4U, 2U), -0.25);
|
|
EXPECT_DOUBLE_EQ(direct(4U, 3U), -0.25);
|
|
|
|
const auto samples = shell.CovariantTyingShearSamples20();
|
|
ASSERT_EQ(samples.Rows(), 4U);
|
|
ASSERT_EQ(samples.Columns(), 20U);
|
|
EXPECT_DOUBLE_EQ(samples(0U, 2U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(0U, 4U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(0U, 7U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(0U, 9U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(1U, 12U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(1U, 14U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(1U, 17U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(1U, 19U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(2U, 2U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(2U, 3U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(2U, 17U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(2U, 18U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(3U, 7U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(3U, 8U), -0.25);
|
|
EXPECT_DOUBLE_EQ(samples(3U, 12U), 0.25);
|
|
EXPECT_DOUBLE_EQ(samples(3U, 13U), -0.25);
|
|
|
|
const auto weights = fesa::Mitc4Shell::TyingWeights(0.25, -0.5);
|
|
EXPECT_EQ(weights.xi_zeta, (std::array<double, 2>{0.75, 0.25}));
|
|
EXPECT_EQ(weights.eta_zeta, (std::array<double, 2>{0.375, 0.625}));
|
|
|
|
const auto tied = shell.StrainDisplacement20(0.0, 0.0, 0.0);
|
|
EXPECT_DOUBLE_EQ(tied(3U, 4U),
|
|
2.0 * (0.5 * samples(0U, 4U) + 0.5 * samples(1U, 4U)));
|
|
EXPECT_DOUBLE_EQ(tied(4U, 3U),
|
|
2.0 * (0.5 * samples(2U, 3U) + 0.5 * samples(3U, 3U)));
|
|
}
|
|
|
|
// MITC4-KIN-004
|
|
TEST(Mitc4ShellConstitutive,
|
|
BuildsExactPositiveDefiniteSectionMatricesAndRescalesUnits) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(candidate.HasValue());
|
|
const auto& shell = candidate.Value();
|
|
|
|
const auto cps = shell.PlaneStressConstitutive();
|
|
const auto c5 = shell.MaterialConstitutive5();
|
|
const auto a = shell.MembraneSectionMatrix();
|
|
const auto d = shell.BendingSectionMatrix();
|
|
const auto as = shell.TransverseShearSectionMatrix();
|
|
EXPECT_EQ(cps.Rows(), 3U);
|
|
EXPECT_EQ(cps.Columns(), 3U);
|
|
EXPECT_EQ(c5.Rows(), 5U);
|
|
EXPECT_EQ(c5.Columns(), 5U);
|
|
EXPECT_EQ(a.Rows(), 3U);
|
|
EXPECT_EQ(d.Rows(), 3U);
|
|
EXPECT_EQ(as.Rows(), 2U);
|
|
EXPECT_DOUBLE_EQ(cps(0U, 0U), 128.0);
|
|
EXPECT_DOUBLE_EQ(cps(0U, 1U), 32.0);
|
|
EXPECT_DOUBLE_EQ(cps(2U, 2U), 48.0);
|
|
EXPECT_DOUBLE_EQ(c5(3U, 3U), 40.0);
|
|
EXPECT_DOUBLE_EQ(c5(4U, 4U), 40.0);
|
|
EXPECT_DOUBLE_EQ(a(0U, 0U), 256.0);
|
|
EXPECT_NEAR(d(0U, 0U), 256.0 / 3.0, 1.0e-12);
|
|
EXPECT_DOUBLE_EQ(as(0U, 0U), 80.0);
|
|
ExpectSymmetric(cps);
|
|
ExpectSymmetric(c5);
|
|
EXPECT_TRUE(HasPositiveCholeskyPivots(cps));
|
|
EXPECT_TRUE(HasPositiveCholeskyPivots(c5));
|
|
EXPECT_TRUE(HasPositiveCholeskyPivots(a));
|
|
EXPECT_TRUE(HasPositiveCholeskyPivots(d));
|
|
EXPECT_TRUE(HasPositiveCholeskyPivots(as));
|
|
|
|
constexpr double kForceScale = 7.0;
|
|
constexpr double kLengthScale = 3.0;
|
|
const auto scaled_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(2.0 * kLengthScale),
|
|
Material(120.0 * kForceScale / (kLengthScale * kLengthScale), 0.25));
|
|
ASSERT_TRUE(scaled_candidate.HasValue());
|
|
const auto& scaled = scaled_candidate.Value();
|
|
fesa::Matrix expected_cps{3U, 3U};
|
|
fesa::Matrix expected_c5{5U, 5U};
|
|
fesa::Matrix expected_a{3U, 3U};
|
|
fesa::Matrix expected_d{3U, 3U};
|
|
fesa::Matrix expected_as{2U, 2U};
|
|
for (std::size_t row = 0U; row < 3U; ++row) {
|
|
for (std::size_t column = 0U; column < 3U; ++column) {
|
|
expected_cps(row, column) =
|
|
cps(row, column) * kForceScale / (kLengthScale * kLengthScale);
|
|
expected_a(row, column) = a(row, column) * kForceScale / kLengthScale;
|
|
expected_d(row, column) = d(row, column) * kForceScale * kLengthScale;
|
|
}
|
|
}
|
|
for (std::size_t row = 0U; row < 2U; ++row) {
|
|
for (std::size_t column = 0U; column < 2U; ++column) {
|
|
expected_as(row, column) = as(row, column) * kForceScale / kLengthScale;
|
|
}
|
|
}
|
|
for (std::size_t row = 0U; row < 5U; ++row) {
|
|
for (std::size_t column = 0U; column < 5U; ++column) {
|
|
expected_c5(row, column) =
|
|
c5(row, column) * kForceScale / (kLengthScale * kLengthScale);
|
|
}
|
|
}
|
|
ExpectMatrixNear(scaled.PlaneStressConstitutive(), expected_cps);
|
|
ExpectMatrixNear(scaled.MaterialConstitutive5(), expected_c5);
|
|
ExpectMatrixNear(scaled.MembraneSectionMatrix(), expected_a);
|
|
ExpectMatrixNear(scaled.BendingSectionMatrix(), expected_d);
|
|
ExpectMatrixNear(scaled.TransverseShearSectionMatrix(), expected_as);
|
|
|
|
EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
|
|
Section(0.0), Material())
|
|
.HasValue());
|
|
EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
|
|
Section(), Material(0.0, 0.25))
|
|
.HasValue());
|
|
EXPECT_FALSE(fesa::Mitc4Shell::Create(NodePointers(nodes), Directors(),
|
|
Section(), Material(120.0, 0.5))
|
|
.HasValue());
|
|
}
|
|
|
|
// MITC4-KIN-005
|
|
TEST(Mitc4ShellKinematics, UsesOneFixedTwoByTwoByTwoQuadratureOrder) {
|
|
const auto& points = fesa::Mitc4Shell::VolumeQuadrature();
|
|
ASSERT_EQ(points.size(), 8U);
|
|
const double g = 1.0 / std::sqrt(3.0);
|
|
const std::array<Vector3, 8> expected{Vector3{-g, -g, -g}, Vector3{-g, -g, g},
|
|
Vector3{g, -g, -g}, Vector3{g, -g, g},
|
|
Vector3{g, g, -g}, Vector3{g, g, g},
|
|
Vector3{-g, g, -g}, Vector3{-g, g, g}};
|
|
for (std::size_t point = 0U; point < points.size(); ++point) {
|
|
EXPECT_EQ(points[point].natural_coordinates, expected[point]);
|
|
EXPECT_DOUBLE_EQ(points[point].weight, 1.0);
|
|
}
|
|
}
|
|
|
|
// MITC4-KERNEL-001
|
|
TEST(Mitc4ShellKernel,
|
|
FormsFiniteScaledSymmetricPhysicalAndStabilizedStiffness) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
|
|
const auto stiffness_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
EXPECT_EQ(stiffness.physical_local20.Rows(), 20U);
|
|
EXPECT_EQ(stiffness.physical_local20.Columns(), 20U);
|
|
EXPECT_EQ(stiffness.physical_global24.Rows(), 24U);
|
|
EXPECT_EQ(stiffness.drilling_global24.Rows(), 24U);
|
|
EXPECT_EQ(stiffness.stabilized_global24.Rows(), 24U);
|
|
for (const fesa::Matrix* matrix :
|
|
{&stiffness.physical_local20, &stiffness.physical_global24,
|
|
&stiffness.drilling_global24, &stiffness.stabilized_global24}) {
|
|
for (std::size_t row = 0U; row < matrix->Rows(); ++row) {
|
|
for (std::size_t column = 0U; column < matrix->Columns(); ++column) {
|
|
EXPECT_TRUE(std::isfinite((*matrix)(row, column)));
|
|
}
|
|
}
|
|
}
|
|
EXPECT_LE(ScaledSymmetryError(stiffness.physical_local20, 5U, 2.0), 1.0e-12);
|
|
EXPECT_LE(ScaledSymmetryError(stiffness.physical_global24, 6U, 2.0), 1.0e-12);
|
|
EXPECT_LE(ScaledSymmetryError(stiffness.drilling_global24, 6U, 2.0), 1.0e-12);
|
|
EXPECT_LE(ScaledSymmetryError(stiffness.stabilized_global24, 6U, 2.0),
|
|
1.0e-12);
|
|
|
|
const auto repeated_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(repeated_candidate.HasValue());
|
|
const auto& repeated = repeated_candidate.Value();
|
|
ExpectMatrixNear(repeated.physical_local20, stiffness.physical_local20, 0.0);
|
|
ExpectMatrixNear(repeated.physical_global24, stiffness.physical_global24,
|
|
0.0);
|
|
ExpectMatrixNear(repeated.drilling_global24, stiffness.drilling_global24,
|
|
0.0);
|
|
ExpectMatrixNear(repeated.stabilized_global24, stiffness.stabilized_global24,
|
|
0.0);
|
|
EXPECT_DOUBLE_EQ(repeated.drilling_stiffness, stiffness.drilling_stiffness);
|
|
}
|
|
|
|
TEST(Mitc4ShellKernel,
|
|
PreservesExactStiffnessRecoveryAndPatchAcrossVector3Migration) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto& shell = shell_candidate.Value();
|
|
const auto stiffness_candidate = shell.Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
const fesa::Vector3 frame_e3{shell.LocalFrame(0.0, 0.0).e3};
|
|
EXPECT_DOUBLE_EQ(frame_e3.X(), 0.0);
|
|
EXPECT_DOUBLE_EQ(frame_e3.Y(), 0.0);
|
|
EXPECT_DOUBLE_EQ(frame_e3.Z(), 1.0);
|
|
|
|
constexpr std::array<double, 8> kGeneralized{0.1, -0.05, 0.2, 0.3,
|
|
-0.15, 0.25, 0.4, -0.3};
|
|
fesa::Vector global_field{24U};
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
const double x = nodes[node_index].coordinates[0];
|
|
const double y = nodes[node_index].coordinates[1];
|
|
const std::size_t offset = 6U * node_index;
|
|
global_field[offset] = kGeneralized[0] * x + 0.5 * kGeneralized[2] * y;
|
|
global_field[offset + 1U] = kGeneralized[1] * y + 0.5 * kGeneralized[2] * x;
|
|
global_field[offset + 2U] = kGeneralized[6] * x + kGeneralized[7] * y -
|
|
0.5 * kGeneralized[5] * x * y;
|
|
global_field[offset + 3U] =
|
|
-kGeneralized[4] * y - 0.5 * kGeneralized[5] * x;
|
|
global_field[offset + 4U] = kGeneralized[3] * x + 0.5 * kGeneralized[5] * y;
|
|
}
|
|
const auto recovery_candidate = shell.RecoverPhysical(global_field);
|
|
ASSERT_TRUE(recovery_candidate.HasValue());
|
|
const auto& recovery = recovery_candidate.Value();
|
|
|
|
std::array<std::array<double, 5>, 4> e11_values{};
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
e11_values[node_index][0U] = 0.2 * nodes[node_index].coordinates[0U];
|
|
}
|
|
const double gauss = 1.0 / std::sqrt(3.0);
|
|
const fesa::Vector patch = shell.StrainDisplacement20(gauss, -gauss, gauss)
|
|
.Multiply(PhysicalField(e11_values));
|
|
|
|
EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 0U), 0x1.d555555555554p+6);
|
|
EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 4U), 0.0);
|
|
EXPECT_DOUBLE_EQ(stiffness.physical_global24(2U, 2U), 0x1.aaaaaaaaaaaadp+5);
|
|
EXPECT_DOUBLE_EQ(stiffness.drilling_stiffness, 0x1.0d6cffc5beeb5p-4);
|
|
EXPECT_DOUBLE_EQ(recovery.strain_energy, 0x1.20a3d70a3d70bp+6);
|
|
constexpr std::array<double, 8> kExpectedStrain{
|
|
0x1.999999999999ap-4, -0x1.999999999999bp-5, 0x1.9999999999998p-3,
|
|
0x1.3333333333335p-2, -0x1.3333333333335p-3, 0x1.0000000000000p-2,
|
|
0x1.999999999999cp-2, -0x1.3333333333334p-2};
|
|
EXPECT_EQ(recovery.points[0U].generalized_strain, kExpectedStrain);
|
|
EXPECT_DOUBLE_EQ(patch[0U], 0x1.9999999999999p-3);
|
|
for (std::size_t component = 1U; component < patch.Size(); ++component) {
|
|
EXPECT_DOUBLE_EQ(patch[component], 0.0);
|
|
}
|
|
}
|
|
|
|
// MITC4-KERNEL-002
|
|
TEST(Mitc4ShellKernel,
|
|
PreservesPhysicalEnergyUnderTwentyToTwentyFourCongruence) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto stiffness_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
|
|
fesa::Vector global_field{24U};
|
|
for (std::size_t index = 0U; index < global_field.Size(); ++index) {
|
|
global_field[index] =
|
|
0.125 * static_cast<double>(static_cast<int>(index % 7U) - 3);
|
|
}
|
|
const auto physical_field20 =
|
|
shell_candidate.Value().PhysicalTransformation20().Multiply(global_field);
|
|
const double local_energy =
|
|
QuadraticEnergy(stiffness.physical_local20, physical_field20);
|
|
const double global_energy =
|
|
QuadraticEnergy(stiffness.physical_global24, global_field);
|
|
ASSERT_NE(local_energy, 0.0);
|
|
ASSERT_NE(global_energy, 0.0);
|
|
EXPECT_LE(std::abs(global_energy - local_energy) /
|
|
(std::abs(global_energy) + std::abs(local_energy)),
|
|
1.0e-12);
|
|
}
|
|
|
|
// MITC4-KERNEL-003
|
|
TEST(Mitc4ShellKernel,
|
|
RetainsSixRigidModesAndHasExpectedPhysicalAndStabilizedRank) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto stiffness_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
|
|
const auto scaled_physical20 =
|
|
ScaledStiffness(stiffness.physical_local20, 5U, 2.0);
|
|
const auto scaled_physical24 =
|
|
ScaledStiffness(stiffness.physical_global24, 6U, 2.0);
|
|
const auto scaled_stabilized24 =
|
|
ScaledStiffness(stiffness.stabilized_global24, 6U, 2.0);
|
|
EXPECT_EQ(NumericalRank(scaled_physical20), 14U);
|
|
EXPECT_EQ(NumericalRank(scaled_stabilized24), 18U);
|
|
|
|
const std::array<fesa::Vector, 6> rigid_modes{
|
|
PhysicalRigidMode(nodes, {1.0, 0.0, 0.0}, {}),
|
|
PhysicalRigidMode(nodes, {0.0, 1.0, 0.0}, {}),
|
|
PhysicalRigidMode(nodes, {0.0, 0.0, 1.0}, {}),
|
|
PhysicalRigidMode(nodes, {}, {1.0, 0.0, 0.0}),
|
|
PhysicalRigidMode(nodes, {}, {0.0, 1.0, 0.0}),
|
|
PhysicalRigidMode(nodes, {}, {0.0, 0.0, 1.0})};
|
|
const double physical_norm = SymmetricOperatorNorm(scaled_physical24);
|
|
const double stabilized_norm = SymmetricOperatorNorm(scaled_stabilized24);
|
|
ASSERT_GT(physical_norm, 0.0);
|
|
ASSERT_GT(stabilized_norm, 0.0);
|
|
for (const auto& rigidMode : rigid_modes) {
|
|
fesa::Vector scaledMode = rigidMode;
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
scaledMode[6U * node_index + component] /= 2.0;
|
|
}
|
|
}
|
|
const double modeNorm = scaledMode.Norm();
|
|
ASSERT_GT(modeNorm, 0.0);
|
|
EXPECT_LE(scaled_physical24.Multiply(scaledMode).Norm() /
|
|
(physical_norm * modeNorm),
|
|
1.0e-10);
|
|
EXPECT_LE(scaled_stabilized24.Multiply(scaledMode).Norm() /
|
|
(stabilized_norm * modeNorm),
|
|
1.0e-10);
|
|
}
|
|
}
|
|
|
|
// MITC4-KERNEL-004
|
|
TEST(Mitc4ShellPatch, ReproducesIndependentMembraneBendingShearAndTwistFields) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto& shell = shell_candidate.Value();
|
|
const auto stiffness_candidate = shell.Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value().physical_local20;
|
|
constexpr double kMagnitude = 0.2;
|
|
const double gauss = 1.0 / std::sqrt(3.0);
|
|
|
|
std::array<std::array<double, 5>, 4> e11_values{};
|
|
std::array<std::array<double, 5>, 4> e22_values{};
|
|
std::array<std::array<double, 5>, 4> g12_values{};
|
|
std::array<std::array<double, 5>, 4> k11_values{};
|
|
std::array<std::array<double, 5>, 4> k22_values{};
|
|
std::array<std::array<double, 5>, 4> g13_values{};
|
|
std::array<std::array<double, 5>, 4> g23_values{};
|
|
std::array<std::array<double, 5>, 4> k12_values{};
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
const double x = nodes[node_index].coordinates[0];
|
|
const double y = nodes[node_index].coordinates[1];
|
|
e11_values[node_index][0] = kMagnitude * x;
|
|
e22_values[node_index][1] = kMagnitude * y;
|
|
g12_values[node_index][0] = 0.5 * kMagnitude * y;
|
|
g12_values[node_index][1] = 0.5 * kMagnitude * x;
|
|
k11_values[node_index][4] = kMagnitude * x;
|
|
k22_values[node_index][3] = -kMagnitude * y;
|
|
g13_values[node_index][2] = kMagnitude * x;
|
|
g23_values[node_index][2] = kMagnitude * y;
|
|
k12_values[node_index][2] = -0.5 * kMagnitude * x * y;
|
|
k12_values[node_index][3] = -0.5 * kMagnitude * x;
|
|
k12_values[node_index][4] = 0.5 * kMagnitude * y;
|
|
}
|
|
|
|
const std::array<fesa::Vector, 8> fields{
|
|
PhysicalField(e11_values), PhysicalField(e22_values),
|
|
PhysicalField(g12_values), PhysicalField(k11_values),
|
|
PhysicalField(k22_values), PhysicalField(g13_values),
|
|
PhysicalField(g23_values), PhysicalField(k12_values)};
|
|
ExpectStrain(shell, fields[0], gauss, -gauss, gauss,
|
|
{kMagnitude, 0.0, 0.0, 0.0, 0.0});
|
|
ExpectStrain(shell, fields[1], gauss, -gauss, gauss,
|
|
{0.0, kMagnitude, 0.0, 0.0, 0.0});
|
|
ExpectStrain(shell, fields[2], gauss, -gauss, gauss,
|
|
{0.0, 0.0, kMagnitude, 0.0, 0.0});
|
|
ExpectStrain(shell, fields[3], gauss, -gauss, gauss,
|
|
{gauss * kMagnitude, 0.0, 0.0, 0.0, 0.0});
|
|
ExpectStrain(shell, fields[4], gauss, -gauss, gauss,
|
|
{0.0, gauss * kMagnitude, 0.0, 0.0, 0.0});
|
|
ExpectStrain(shell, fields[5], gauss, -gauss, gauss,
|
|
{0.0, 0.0, 0.0, kMagnitude, 0.0});
|
|
ExpectStrain(shell, fields[6], gauss, -gauss, gauss,
|
|
{0.0, 0.0, 0.0, 0.0, kMagnitude});
|
|
ExpectStrain(shell, fields[7], gauss, -gauss, gauss,
|
|
{0.0, 0.0, gauss * kMagnitude, 0.0, 0.0});
|
|
for (const auto& field : fields) {
|
|
EXPECT_GT(QuadraticEnergy(stiffness, field), 0.0);
|
|
}
|
|
}
|
|
|
|
// MITC4-KERNEL-005
|
|
TEST(Mitc4ShellDrilling,
|
|
UsesOnlyEightPositivePhysicalRotationDiagonalsAndFixedFactor) {
|
|
const std::array<fesa::Node, 4> nodes{
|
|
Node(1, {-50.0, -50.0, 0.0}), Node(2, {50.0, -50.0, 0.0}),
|
|
Node(3, {50.0, 50.0, 0.0}), Node(4, {-50.0, 50.0, 0.0})};
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(0.1), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto stiffness_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
|
|
double expected_reference = (std::numeric_limits<double>::max)();
|
|
double all_diagonal_minimum = (std::numeric_limits<double>::max)();
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
for (std::size_t rotation = 3U; rotation < 5U; ++rotation) {
|
|
const std::size_t index = 5U * node_index + rotation;
|
|
const double diagonal = stiffness.physical_local20(index, index);
|
|
ASSERT_TRUE(std::isfinite(diagonal));
|
|
ASSERT_GT(diagonal, 0.0);
|
|
expected_reference = (std::min)(expected_reference, diagonal);
|
|
}
|
|
for (std::size_t component = 0U; component < 5U; ++component) {
|
|
const std::size_t index = 5U * node_index + component;
|
|
const double diagonal = stiffness.physical_local20(index, index);
|
|
if (std::isfinite(diagonal) && diagonal > 0.0) {
|
|
all_diagonal_minimum = (std::min)(all_diagonal_minimum, diagonal);
|
|
}
|
|
}
|
|
}
|
|
EXPECT_DOUBLE_EQ(stiffness.drilling_stiffness, 1.0e-3 * expected_reference);
|
|
EXPECT_LT(all_diagonal_minimum, expected_reference);
|
|
EXPECT_NE(stiffness.drilling_stiffness, 1.0e-3 * all_diagonal_minimum);
|
|
}
|
|
|
|
// MITC4-KERNEL-006
|
|
TEST(Mitc4ShellDrilling,
|
|
FailsNonfiniteReferenceAndStabilizesEachPureDrillCoordinate) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto stiffness_candidate = shell_candidate.Value().Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
const auto& stiffness = stiffness_candidate.Value();
|
|
|
|
for (std::size_t node_index = 0U; node_index < 4U; ++node_index) {
|
|
fesa::Vector pure_drill{24U};
|
|
pure_drill[6U * node_index + 5U] = 1.0;
|
|
EXPECT_DOUBLE_EQ(stiffness.physical_global24.Multiply(pure_drill).Norm(),
|
|
0.0);
|
|
const auto drill_action = stiffness.drilling_global24.Multiply(pure_drill);
|
|
EXPECT_DOUBLE_EQ(drill_action[6U * node_index + 5U],
|
|
stiffness.drilling_stiffness);
|
|
EXPECT_GT(QuadraticEnergy(stiffness.drilling_global24, pure_drill), 0.0);
|
|
}
|
|
|
|
const std::array<fesa::Node, 4> extreme_nodes{
|
|
Node(1, {-5.0e9, -5.0e9, 0.0}), Node(2, {5.0e9, -5.0e9, 0.0}),
|
|
Node(3, {5.0e9, 5.0e9, 0.0}), Node(4, {-5.0e9, 5.0e9, 0.0})};
|
|
const auto extreme_shell =
|
|
fesa::Mitc4Shell::Create(NodePointers(extreme_nodes), Directors(),
|
|
Section(1.0), Material(1.0e300));
|
|
ASSERT_TRUE(extreme_shell.HasValue());
|
|
const auto failure = extreme_shell.Value().Stiffness();
|
|
ASSERT_FALSE(failure.HasValue());
|
|
ASSERT_EQ(failure.GetStatus().Diagnostics().size(), 1U);
|
|
EXPECT_EQ(failure.GetStatus().Diagnostics()[0].code,
|
|
"invalid-shell-stiffness");
|
|
const auto repeated_failure = extreme_shell.Value().Stiffness();
|
|
ASSERT_FALSE(repeated_failure.HasValue());
|
|
ASSERT_EQ(repeated_failure.GetStatus().Diagnostics().size(), 1U);
|
|
EXPECT_EQ(repeated_failure.GetStatus().Diagnostics()[0].code,
|
|
failure.GetStatus().Diagnostics()[0].code);
|
|
EXPECT_EQ(repeated_failure.GetStatus().Diagnostics()[0].message,
|
|
failure.GetStatus().Diagnostics()[0].message);
|
|
}
|
|
|
|
// MITC4-KERNEL-007
|
|
TEST(Mitc4ShellDrilling, ExcludesPureDrillFromPhysicalRecoveryAndEnergy) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto& shell = shell_candidate.Value();
|
|
const auto stiffness_candidate = shell.Stiffness();
|
|
ASSERT_TRUE(stiffness_candidate.HasValue());
|
|
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
fesa::Vector pure_drill{24U};
|
|
pure_drill[6U * node_index + 5U] = 1.0;
|
|
EXPECT_GT(stiffness_candidate.Value()
|
|
.stabilized_global24.Multiply(pure_drill)
|
|
.Norm(),
|
|
0.0);
|
|
|
|
const auto recovery_candidate = shell.RecoverPhysical(pure_drill);
|
|
ASSERT_TRUE(recovery_candidate.HasValue());
|
|
const auto& recovery = recovery_candidate.Value();
|
|
EXPECT_DOUBLE_EQ(recovery.strain_energy, 0.0);
|
|
for (const auto& point : recovery.points) {
|
|
for (double value : point.generalized_strain) {
|
|
EXPECT_DOUBLE_EQ(value, 0.0);
|
|
}
|
|
for (double value : point.section_resultant) {
|
|
EXPECT_DOUBLE_EQ(value, 0.0);
|
|
}
|
|
for (const auto& stress : point.in_plane_stress) {
|
|
for (double value : stress) {
|
|
EXPECT_DOUBLE_EQ(value, 0.0);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// MITC4-PHYSREC-001
|
|
TEST(Mitc4ShellPhysicalRecovery,
|
|
RecoversHandFieldAtFixedLocationsAndSectionPositions) {
|
|
const auto nodes = PlanarNodes();
|
|
const auto shell_candidate = fesa::Mitc4Shell::Create(
|
|
NodePointers(nodes), Directors(), Section(), Material());
|
|
ASSERT_TRUE(shell_candidate.HasValue());
|
|
const auto& shell = shell_candidate.Value();
|
|
|
|
constexpr std::array<double, 8> kGeneralized{0.1, -0.05, 0.2, 0.3,
|
|
-0.15, 0.25, 0.4, -0.3};
|
|
fesa::Vector global_field{24U};
|
|
for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
|
|
const double x = nodes[node_index].coordinates[0];
|
|
const double y = nodes[node_index].coordinates[1];
|
|
const std::size_t offset = 6U * node_index;
|
|
global_field[offset] = kGeneralized[0] * x + 0.5 * kGeneralized[2] * y;
|
|
global_field[offset + 1U] = kGeneralized[1] * y + 0.5 * kGeneralized[2] * x;
|
|
global_field[offset + 2U] = kGeneralized[6] * x + kGeneralized[7] * y -
|
|
0.5 * kGeneralized[5] * x * y;
|
|
global_field[offset + 3U] =
|
|
-kGeneralized[4] * y - 0.5 * kGeneralized[5] * x;
|
|
global_field[offset + 4U] = kGeneralized[3] * x + 0.5 * kGeneralized[5] * y;
|
|
}
|
|
|
|
const auto recovery_candidate = shell.RecoverPhysical(global_field);
|
|
ASSERT_TRUE(recovery_candidate.HasValue());
|
|
const auto& recovery = recovery_candidate.Value();
|
|
const double gauss = 1.0 / std::sqrt(3.0);
|
|
const std::array<std::array<double, 2>, 4> expected_coordinates{
|
|
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, 8> kExpectedResultant{22.4, -6.4, 19.2, 22.4,
|
|
-6.4, 8.0, 32.0, -24.0};
|
|
constexpr std::array<std::array<double, 3>, 3> kExpectedStress{
|
|
std::array<double, 3>{-22.4, 6.4, -2.4},
|
|
std::array<double, 3>{11.2, -3.2, 9.6},
|
|
std::array<double, 3>{44.8, -12.8, 21.6}};
|
|
|
|
ASSERT_EQ(recovery.points.size(), expected_coordinates.size());
|
|
for (std::size_t point_index = 0U; point_index < recovery.points.size();
|
|
++point_index) {
|
|
const auto& point = recovery.points[point_index];
|
|
EXPECT_EQ(point.natural_coordinates, expected_coordinates[point_index]);
|
|
ExpectOrthonormalRightHanded(point.local_frame);
|
|
ExpectVectorNear(point.local_frame.e1, {1.0, 0.0, 0.0});
|
|
ExpectVectorNear(point.local_frame.e2, {0.0, 1.0, 0.0});
|
|
ExpectVectorNear(point.local_frame.e3, {0.0, 0.0, 1.0});
|
|
for (std::size_t component = 0U; component < kGeneralized.size();
|
|
++component) {
|
|
EXPECT_NEAR(point.generalized_strain[component], kGeneralized[component],
|
|
1.0e-12);
|
|
EXPECT_NEAR(point.section_resultant[component],
|
|
kExpectedResultant[component], 1.0e-12);
|
|
}
|
|
for (std::size_t position = 0U; position < kExpectedStress.size();
|
|
++position) {
|
|
for (std::size_t component = 0U;
|
|
component < kExpectedStress[position].size(); ++component) {
|
|
EXPECT_NEAR(point.in_plane_stress[position][component],
|
|
kExpectedStress[position][component], 1.0e-12);
|
|
}
|
|
}
|
|
}
|
|
EXPECT_NEAR(recovery.strain_energy, 72.16, 1.0e-12);
|
|
}
|