From c0a68ee964b9c5a2d757f5a143108610fe24705b Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 12 Aug 2026 19:28:54 +0900 Subject: [PATCH] feat(linear-static-mitc4-shell): step 2 - shell-director-geometry --- include/fesa/model/shell_geometry.hpp | 44 ++ src/fesa/CMakeLists.txt | 1 + src/fesa/io/abaqus/domain_mapper.cpp | 17 + src/fesa/model/shell_geometry.cpp | 491 ++++++++++++++++++++ tests/CMakeLists.txt | 1 + tests/unit/io/abaqus/domain_mapper_test.cpp | 12 + tests/unit/model/shell_geometry_test.cpp | 276 +++++++++++ 7 files changed, 842 insertions(+) create mode 100644 include/fesa/model/shell_geometry.hpp create mode 100644 src/fesa/model/shell_geometry.cpp create mode 100644 tests/unit/model/shell_geometry_test.cpp diff --git a/include/fesa/model/shell_geometry.hpp b/include/fesa/model/shell_geometry.hpp new file mode 100644 index 0000000..43729e5 --- /dev/null +++ b/include/fesa/model/shell_geometry.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include "fesa/core/status.hpp" +#include "fesa/model/model_types.hpp" + +#include +#include +#include + +namespace fesa { + +enum class ShellGeometryPointKind { + center, + stiffness, + tying, + recovery +}; + +struct ShellGeometryValidationPoint { + ShellGeometryPointKind kind; + std::size_t locationIndex; + std::array naturalCoordinates; +}; + +struct ShellElementGeometryData { + EntityIndex elementIndex; + std::array normalCandidate; + double surfaceAreaWeight; +}; + +struct ShellGeometry { + std::vector nodalFrames; + std::vector elementData; +}; + +const std::array& +shellGeometryValidationPoints() noexcept; + +Result preprocessShellGeometry( + const std::vector& nodes, + const std::vector& elements, + const std::vector& sections); + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index cbaa81b..f5bbcb7 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -21,6 +21,7 @@ add_library( math/sparse_matrix.cpp math/vector.cpp model/domain.cpp + model/shell_geometry.cpp results/result_recovery.cpp solvers/linear/mkl_pardiso_solver.cpp ) diff --git a/src/fesa/io/abaqus/domain_mapper.cpp b/src/fesa/io/abaqus/domain_mapper.cpp index 27dc04c..c7b0c63 100644 --- a/src/fesa/io/abaqus/domain_mapper.cpp +++ b/src/fesa/io/abaqus/domain_mapper.cpp @@ -1,5 +1,7 @@ #include "fesa/io/abaqus/domain_mapper.hpp" +#include "fesa/model/shell_geometry.hpp" + #include #include #include @@ -1768,6 +1770,21 @@ private: if (failure_) { return; } + if (!definition_.shellElements.empty()) { + auto geometry = preprocessShellGeometry( + definition_.nodes, + definition_.shellElements, + definition_.shellSections); + if (!geometry.hasValue()) { + const auto& status = geometry.status(); + failure_ = MappingFailure{ + status.failureCategory().value_or(FailureCategory::model), + status.diagnostics().front()}; + return; + } + definition_.shellNodeInitialFrames = + std::move(geometry.value().nodalFrames); + } expandAssemblySets(); if (failure_) { return; diff --git a/src/fesa/model/shell_geometry.cpp b/src/fesa/model/shell_geometry.cpp new file mode 100644 index 0000000..11f2163 --- /dev/null +++ b/src/fesa/model/shell_geometry.cpp @@ -0,0 +1,491 @@ +#include "fesa/model/shell_geometry.hpp" + +#include +#include +#include +#include + +namespace fesa { +namespace { + +using Vector3 = std::array; + +constexpr std::array kXiSigns{-1.0, 1.0, 1.0, -1.0}; +constexpr std::array kEtaSigns{-1.0, -1.0, 1.0, 1.0}; + +struct ShapeData { + std::array values; + std::array xiDerivatives; + std::array etaDerivatives; +}; + +struct ElementWork { + std::array coordinates; + Vector3 normal; + double areaWeight; +}; + +Vector3 add(const Vector3& left, const Vector3& right) { + return { + left[0] + right[0], left[1] + right[1], left[2] + right[2]}; +} + +Vector3 subtract(const Vector3& left, const Vector3& right) { + return { + left[0] - right[0], left[1] - right[1], left[2] - right[2]}; +} + +Vector3 scale(double factor, const Vector3& value) { + return {factor * value[0], factor * value[1], factor * value[2]}; +} + +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::hypot(value[0], value[1], value[2]); +} + +bool isFinite(const Vector3& value) { + return std::all_of(value.begin(), value.end(), [](double component) { + return std::isfinite(component); + }); +} + +ShapeData shapeData(double xi, double eta) { + ShapeData data{}; + for (std::size_t node = 0U; node < 4U; ++node) { + data.values[node] = + 0.25 * (1.0 + kXiSigns[node] * xi) * + (1.0 + kEtaSigns[node] * eta); + data.xiDerivatives[node] = + 0.25 * kXiSigns[node] * (1.0 + kEtaSigns[node] * eta); + data.etaDerivatives[node] = + 0.25 * kEtaSigns[node] * (1.0 + kXiSigns[node] * xi); + } + return data; +} + +Vector3 weightedSum( + const std::array& weights, + const std::array& values) { + Vector3 result{}; + for (std::size_t node = 0U; node < values.size(); ++node) { + result = add(result, scale(weights[node], values[node])); + } + return result; +} + +Vector3 derivativeSum( + const std::array& derivatives, + const std::array& coordinates) { + // Shape derivatives sum to zero, so translating by node 1 improves the + // numerical cancellation without changing the covariant tangent. + std::array relative{}; + for (std::size_t node = 0U; node < coordinates.size(); ++node) { + relative[node] = subtract(coordinates[node], coordinates[0]); + } + return weightedSum(derivatives, relative); +} + +Result geometryFailure( + std::string code, + const SourceLocation& location, + std::string keyword, + std::string identity, + std::string message) { + return Result::failure(Status::failure( + FailureCategory::model, + {{Severity::error, + std::move(code), + location, + std::move(keyword), + std::move(identity), + std::move(message)}})); +} + +bool sameCoordinates(const Vector3& left, const Vector3& right) { + return left == right; +} + +double orientation( + const Vector3& first, + const Vector3& second, + const Vector3& third, + const Vector3& normal) { + return dot(cross(subtract(second, first), subtract(third, first)), normal); +} + +bool hasOppositeSigns(double first, double second) { + return (first < 0.0 && second > 0.0) || + (first > 0.0 && second < 0.0); +} + +bool segmentsProperlyIntersect( + const Vector3& firstStart, + const Vector3& firstEnd, + const Vector3& secondStart, + const Vector3& secondEnd, + const Vector3& normal) { + const double firstSideStart = + orientation(firstStart, firstEnd, secondStart, normal); + const double firstSideEnd = + orientation(firstStart, firstEnd, secondEnd, normal); + const double secondSideStart = + orientation(secondStart, secondEnd, firstStart, normal); + const double secondSideEnd = + orientation(secondStart, secondEnd, firstEnd, normal); + return hasOppositeSigns(firstSideStart, firstSideEnd) && + hasOppositeSigns(secondSideStart, secondSideEnd); +} + +bool sourceIdentityLess( + const Mitc4ShellDefinition& left, + std::size_t leftIndex, + const Mitc4ShellDefinition& right, + std::size_t rightIndex) { + return std::tie( + left.sourceId.instanceName, + left.sourceId.sourceLabel, + left.sourceId.sourceLabelText, + leftIndex) < + std::tie( + right.sourceId.instanceName, + right.sourceId.sourceLabel, + right.sourceId.sourceLabelText, + rightIndex); +} + +} // namespace + +const std::array& +shellGeometryValidationPoints() noexcept { + static const std::array points = [] { + const double g = 1.0 / std::sqrt(3.0); + return std::array{ + ShellGeometryValidationPoint{ + ShellGeometryPointKind::center, 0U, {0.0, 0.0, 0.0}}, + {ShellGeometryPointKind::stiffness, 0U, {-g, -g, -g}}, + {ShellGeometryPointKind::stiffness, 1U, {-g, -g, g}}, + {ShellGeometryPointKind::stiffness, 2U, {g, -g, -g}}, + {ShellGeometryPointKind::stiffness, 3U, {g, -g, g}}, + {ShellGeometryPointKind::stiffness, 4U, {g, g, -g}}, + {ShellGeometryPointKind::stiffness, 5U, {g, g, g}}, + {ShellGeometryPointKind::stiffness, 6U, {-g, g, -g}}, + {ShellGeometryPointKind::stiffness, 7U, {-g, g, g}}, + {ShellGeometryPointKind::tying, 0U, {0.0, -1.0, 0.0}}, + {ShellGeometryPointKind::tying, 1U, {0.0, 1.0, 0.0}}, + {ShellGeometryPointKind::tying, 2U, {-1.0, 0.0, 0.0}}, + {ShellGeometryPointKind::tying, 3U, {1.0, 0.0, 0.0}}, + {ShellGeometryPointKind::recovery, 0U, {-g, -g, 0.0}}, + {ShellGeometryPointKind::recovery, 1U, {g, -g, 0.0}}, + {ShellGeometryPointKind::recovery, 2U, {g, g, 0.0}}, + {ShellGeometryPointKind::recovery, 3U, {-g, g, 0.0}}}; + }(); + return points; +} + +Result preprocessShellGeometry( + const std::vector& nodes, + const std::vector& elements, + const std::vector& sections) { + ShellGeometry geometry; + geometry.elementData.reserve(elements.size()); + std::vector work; + work.reserve(elements.size()); + + const double g = 1.0 / std::sqrt(3.0); + const std::array surfaceGaussPoints{ + Vector3{-g, -g, 0.0}, + Vector3{g, -g, 0.0}, + Vector3{g, g, 0.0}, + Vector3{-g, g, 0.0}}; + + for (std::size_t elementIndex = 0U; + elementIndex < elements.size(); + ++elementIndex) { + const auto& element = elements[elementIndex]; + ElementWork current{}; + for (std::size_t localNode = 0U; localNode < 4U; ++localNode) { + if (element.nodeIndices[localNode] >= nodes.size()) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell geometry references an unavailable internal node."); + } + current.coordinates[localNode] = + nodes[element.nodeIndices[localNode]].coordinates; + if (!isFinite(current.coordinates[localNode])) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell geometry contains a nonfinite source coordinate."); + } + } + for (std::size_t first = 0U; first < 4U; ++first) { + for (std::size_t second = first + 1U; second < 4U; ++second) { + if (element.nodeIndices[first] == element.nodeIndices[second] || + sameCoordinates( + current.coordinates[first], current.coordinates[second])) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell geometry contains duplicate nodes."); + } + } + } + + const ShapeData center = shapeData(0.0, 0.0); + const Vector3 centerXi = + derivativeSum(center.xiDerivatives, current.coordinates); + const Vector3 centerEta = + derivativeSum(center.etaDerivatives, current.coordinates); + const Vector3 centerCross = cross(centerXi, centerEta); + const double centerMeasure = norm(centerCross); + if (!isFinite(centerXi) || !isFinite(centerEta) || + !isFinite(centerCross) || !std::isfinite(centerMeasure) || + !(centerMeasure > 0.0)) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell center has no finite nonzero normal candidate."); + } + current.normal = scale(1.0 / centerMeasure, centerCross); + + if (segmentsProperlyIntersect( + current.coordinates[0], current.coordinates[1], + current.coordinates[2], current.coordinates[3], current.normal) || + segmentsProperlyIntersect( + current.coordinates[1], current.coordinates[2], + current.coordinates[3], current.coordinates[0], current.normal)) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell boundary is self-intersecting in the center-normal projection."); + } + + double areaWeight = 0.0; + for (const auto& point : surfaceGaussPoints) { + const ShapeData shape = shapeData(point[0], point[1]); + const Vector3 tangentXi = + derivativeSum(shape.xiDerivatives, current.coordinates); + const Vector3 tangentEta = + derivativeSum(shape.etaDerivatives, current.coordinates); + const Vector3 areaVector = cross(tangentXi, tangentEta); + const double measure = norm(areaVector); + if (!isFinite(tangentXi) || !isFinite(tangentEta) || + !isFinite(areaVector) || !std::isfinite(measure) || + !(measure > 0.0) || !(dot(areaVector, current.normal) > 0.0)) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell surface is zero-area or locally reversed at a required point."); + } + areaWeight += measure; + } + if (!std::isfinite(areaWeight) || !(areaWeight > 0.0)) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell surface-area weight is nonfinite or zero."); + } + current.areaWeight = areaWeight; + work.push_back(current); + geometry.elementData.push_back({ + static_cast(elementIndex), current.normal, areaWeight}); + } + + std::vector> incident(nodes.size()); + for (std::size_t elementIndex = 0U; + elementIndex < elements.size(); + ++elementIndex) { + for (const EntityIndex nodeIndex : elements[elementIndex].nodeIndices) { + incident[nodeIndex].push_back(elementIndex); + } + } + + geometry.nodalFrames.reserve(nodes.size()); + std::vector frameByNode(nodes.size(), nullptr); + for (std::size_t nodeIndex = 0U; nodeIndex < incident.size(); ++nodeIndex) { + auto& nodeIncident = incident[nodeIndex]; + if (nodeIncident.empty()) { + continue; + } + std::sort( + nodeIncident.begin(), nodeIncident.end(), + [&elements](std::size_t left, std::size_t right) { + return sourceIdentityLess( + elements[left], left, elements[right], right); + }); + for (std::size_t first = 0U; first < nodeIncident.size(); ++first) { + for (std::size_t second = first + 1U; + second < nodeIncident.size(); + ++second) { + const double pairDot = dot( + work[nodeIncident[first]].normal, + work[nodeIncident[second]].normal); + if (!std::isfinite(pairDot) || !(pairDot > 0.0)) { + return geometryFailure( + "opposed-incident-normal", nodes[nodeIndex].location, + "NODE", nodes[nodeIndex].sourceId.sourceLabelText, + "Incident shell normal candidates do not share a positive orientation hemisphere."); + } + } + } + + double maximumWeight = 0.0; + for (const std::size_t elementIndex : nodeIncident) { + maximumWeight = std::max(maximumWeight, work[elementIndex].areaWeight); + } + Vector3 directorSum{}; + for (const std::size_t elementIndex : nodeIncident) { + directorSum = add( + directorSum, + scale( + work[elementIndex].areaWeight / maximumWeight, + work[elementIndex].normal)); + } + const double directorNorm = norm(directorSum); + if (!isFinite(directorSum) || !std::isfinite(directorNorm) || + !(directorNorm > 0.0)) { + return geometryFailure( + "invalid-shell-director", nodes[nodeIndex].location, "NODE", + nodes[nodeIndex].sourceId.sourceLabelText, + "Area-weighted shell director is nonfinite or zero."); + } + const Vector3 director = scale(1.0 / directorNorm, directorSum); + + const std::array globalAxes{ + Vector3{1.0, 0.0, 0.0}, + Vector3{0.0, 1.0, 0.0}, + Vector3{0.0, 0.0, 1.0}}; + std::size_t selectedAxis = 0U; + double selectedAlignment = std::abs(dot(globalAxes[0], director)); + for (std::size_t axis = 1U; axis < globalAxes.size(); ++axis) { + const double alignment = std::abs(dot(globalAxes[axis], director)); + if (alignment < selectedAlignment) { + selectedAlignment = alignment; + selectedAxis = axis; + } + } + const Vector3 tangentCandidate = subtract( + globalAxes[selectedAxis], + scale(dot(globalAxes[selectedAxis], director), director)); + const double tangentNorm = norm(tangentCandidate); + if (!isFinite(tangentCandidate) || !std::isfinite(tangentNorm) || + !(tangentNorm > 0.0)) { + return geometryFailure( + "invalid-shell-director", nodes[nodeIndex].location, "NODE", + nodes[nodeIndex].sourceId.sourceLabelText, + "Least-aligned-axis tangent frame construction failed."); + } + const Vector3 tangentA = scale(1.0 / tangentNorm, tangentCandidate); + const Vector3 tangentB = cross(director, tangentA); + if (!isFinite(tangentB) || !(norm(tangentB) > 0.0)) { + return geometryFailure( + "invalid-shell-director", nodes[nodeIndex].location, "NODE", + nodes[nodeIndex].sourceId.sourceLabelText, + "Right-handed shell tangent frame construction failed."); + } + geometry.nodalFrames.push_back({ + static_cast(nodeIndex), director, tangentA, tangentB}); + } + + // Build lookup only after frame storage is complete so later code never + // observes a pointer invalidated by vector growth. + for (const auto& frame : geometry.nodalFrames) { + frameByNode[frame.nodeIndex] = &frame; + } + + for (std::size_t elementIndex = 0U; + elementIndex < elements.size(); + ++elementIndex) { + const auto& element = elements[elementIndex]; + if (element.sectionIndex >= sections.size()) { + return geometryFailure( + "invalid-shell-jacobian", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell geometry cannot resolve its thickness for Jacobian validation."); + } + const double thickness = sections[element.sectionIndex].thickness; + std::array directors{}; + for (std::size_t localNode = 0U; localNode < 4U; ++localNode) { + const auto* frame = frameByNode[element.nodeIndices[localNode]]; + if (frame == nullptr) { + return geometryFailure( + "invalid-shell-director", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell element is missing a nodal director."); + } + directors[localNode] = frame->director; + } + + for (const auto& point : shellGeometryValidationPoints()) { + const double xi = point.naturalCoordinates[0]; + const double eta = point.naturalCoordinates[1]; + const double zeta = point.naturalCoordinates[2]; + const ShapeData shape = shapeData(xi, eta); + const Vector3 midsurfaceXi = + derivativeSum(shape.xiDerivatives, work[elementIndex].coordinates); + const Vector3 midsurfaceEta = + derivativeSum(shape.etaDerivatives, work[elementIndex].coordinates); + const Vector3 areaVector = cross(midsurfaceXi, midsurfaceEta); + const double surfaceMeasure = norm(areaVector); + if (!isFinite(midsurfaceXi) || !isFinite(midsurfaceEta) || + !isFinite(areaVector) || !std::isfinite(surfaceMeasure) || + !(surfaceMeasure > 0.0) || + !(dot(areaVector, work[elementIndex].normal) > 0.0)) { + return geometryFailure( + "invalid-shell-geometry", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell surface basis is nonfinite, zero, or reversed at a required point."); + } + + const Vector3 directorXi = + weightedSum(shape.xiDerivatives, directors); + const Vector3 directorEta = + weightedSum(shape.etaDerivatives, directors); + const Vector3 directorValue = weightedSum(shape.values, directors); + const Vector3 covariantXi = add( + midsurfaceXi, scale(0.5 * thickness * zeta, directorXi)); + const Vector3 covariantEta = add( + midsurfaceEta, scale(0.5 * thickness * zeta, directorEta)); + const Vector3 covariantZeta = scale(0.5 * thickness, directorValue); + const double jacobian = + dot(covariantXi, cross(covariantEta, covariantZeta)); + if (!isFinite(covariantXi) || !isFinite(covariantEta) || + !isFinite(covariantZeta) || !std::isfinite(jacobian) || + !(jacobian > 0.0)) { + return geometryFailure( + "invalid-shell-jacobian", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell Jacobian is nonfinite or nonpositive at a required point."); + } + const Vector3 reciprocalXi = + scale(1.0 / jacobian, cross(covariantEta, covariantZeta)); + const Vector3 reciprocalEta = + scale(1.0 / jacobian, cross(covariantZeta, covariantXi)); + const Vector3 reciprocalZeta = + scale(1.0 / jacobian, cross(covariantXi, covariantEta)); + if (!isFinite(reciprocalXi) || !isFinite(reciprocalEta) || + !isFinite(reciprocalZeta)) { + return geometryFailure( + "invalid-shell-jacobian", element.location, "ELEMENT", + element.sourceId.sourceLabelText, + "Shell reciprocal basis is nonfinite at a required point."); + } + } + } + + return Result::success(std::move(geometry)); +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index dfced5c..ede415f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,6 +23,7 @@ add_executable( unit/io/hdf5/hdf5_results_writer_test.cpp unit/model/domain_test.cpp unit/model/model_types_test.cpp + unit/model/shell_geometry_test.cpp unit/results/result_records_test.cpp unit/results/result_recovery_test.cpp unit/results/results_writer_test.cpp diff --git a/tests/unit/io/abaqus/domain_mapper_test.cpp b/tests/unit/io/abaqus/domain_mapper_test.cpp index ca685f8..e1d90cd 100644 --- a/tests/unit/io/abaqus/domain_mapper_test.cpp +++ b/tests/unit/io/abaqus/domain_mapper_test.cpp @@ -517,6 +517,18 @@ TEST(InpDomainMapping, MapsS4AndS4rThroughOneMitc4Identity) { EXPECT_DOUBLE_EQ(domain.shellSections()[1].thickness, 0.2); EXPECT_EQ(domain.shellSections()[1].materialIndex, 1U); + ASSERT_EQ(domain.shellNodeInitialFrames().size(), domain.nodes().size()); + EXPECT_EQ(domain.shellNodeInitialFrames()[0].nodeIndex, 0U); + EXPECT_EQ( + domain.shellNodeInitialFrames()[0].director, + (std::array{0.0, 0.0, 1.0})); + EXPECT_EQ( + domain.shellNodeInitialFrames()[0].tangentA, + (std::array{1.0, 0.0, 0.0})); + EXPECT_EQ( + domain.shellNodeInitialFrames()[0].tangentB, + (std::array{0.0, 1.0, 0.0})); + ASSERT_EQ(domain.steps().size(), 1U); ASSERT_EQ(domain.steps()[0].boundaries.size(), 1U); EXPECT_EQ(domain.steps()[0].boundaries[0].lastDof, 6); diff --git a/tests/unit/model/shell_geometry_test.cpp b/tests/unit/model/shell_geometry_test.cpp new file mode 100644 index 0000000..a84860e --- /dev/null +++ b/tests/unit/model/shell_geometry_test.cpp @@ -0,0 +1,276 @@ +#include "fesa/model/shell_geometry.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using Vector3 = std::array; + +fesa::Node node(fesa::EntityIndex label, Vector3 coordinates) { + return { + {"Shell-Instance", static_cast(label), + std::to_string(label)}, + coordinates, + {"shell-geometry.inp", static_cast(label + 1U)}}; +} + +fesa::Mitc4ShellDefinition element( + fesa::EntityIndex label, + std::array nodeIndices) { + return { + {"Shell-Instance", static_cast(label), + std::to_string(label)}, + fesa::ShellSourceElementType::s4, + nodeIndices, + 0U, + 0U, + {"shell-geometry.inp", static_cast(100U + label)}}; +} + +std::vector sections(double thickness = 0.2) { + return {{"Section", thickness, 0U, {"shell-geometry.inp", 90U}}}; +} + +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 expectRightHandedFrame(const fesa::ShellNodeInitialFrame& frame) { + EXPECT_NEAR(norm(frame.director), 1.0, 1.0e-12); + EXPECT_NEAR(norm(frame.tangentA), 1.0, 1.0e-12); + EXPECT_NEAR(norm(frame.tangentB), 1.0, 1.0e-12); + EXPECT_NEAR(dot(frame.director, frame.tangentA), 0.0, 1.0e-12); + EXPECT_NEAR(dot(frame.director, frame.tangentB), 0.0, 1.0e-12); + EXPECT_NEAR(dot(frame.tangentA, frame.tangentB), 0.0, 1.0e-12); + expectVectorNear(cross(frame.tangentA, frame.tangentB), frame.director); +} + +const fesa::ShellNodeInitialFrame& frameFor( + const fesa::ShellGeometry& geometry, + fesa::EntityIndex nodeIndex) { + const auto found = std::find_if( + geometry.nodalFrames.begin(), + geometry.nodalFrames.end(), + [nodeIndex](const fesa::ShellNodeInitialFrame& frame) { + return frame.nodeIndex == nodeIndex; + }); + EXPECT_NE(found, geometry.nodalFrames.end()); + return *found; +} + +void expectFailureCode( + const fesa::Result& result, + const std::string& code) { + ASSERT_FALSE(result.hasValue()); + EXPECT_EQ(result.status().failureCategory(), fesa::FailureCategory::model); + ASSERT_EQ(result.status().diagnostics().size(), 1U); + EXPECT_EQ(result.status().diagnostics()[0].code, code); +} + +} // namespace + +// MITC4-GEO-001 +TEST(Mitc4Geometry, BuildsDeterministicFramesForPlanarRotatedAndWarpedElements) { + const std::vector planarNodes{ + node(0U, {0.0, 0.0, 0.0}), + node(1U, {1.0, 0.0, 0.0}), + node(2U, {1.0, 1.0, 0.0}), + node(3U, {0.0, 1.0, 0.0})}; + auto planar = fesa::preprocessShellGeometry( + planarNodes, {element(10U, {0U, 1U, 2U, 3U})}, sections()); + + ASSERT_TRUE(planar.hasValue()); + ASSERT_EQ(planar.value().elementData.size(), 1U); + expectVectorNear(planar.value().elementData[0].normalCandidate, {0.0, 0.0, 1.0}); + EXPECT_NEAR(planar.value().elementData[0].surfaceAreaWeight, 1.0, 1.0e-12); + ASSERT_EQ(planar.value().nodalFrames.size(), 4U); + for (const auto& frame : planar.value().nodalFrames) { + expectVectorNear(frame.director, {0.0, 0.0, 1.0}); + expectVectorNear(frame.tangentA, {1.0, 0.0, 0.0}); + expectVectorNear(frame.tangentB, {0.0, 1.0, 0.0}); + expectRightHandedFrame(frame); + } + + const std::vector rotatedNodes{ + node(0U, {0.0, 0.0, 0.0}), + node(1U, {0.0, 1.0, 0.0}), + node(2U, {0.0, 1.0, 1.0}), + node(3U, {0.0, 0.0, 1.0})}; + auto rotated = fesa::preprocessShellGeometry( + rotatedNodes, {element(11U, {0U, 1U, 2U, 3U})}, sections()); + + ASSERT_TRUE(rotated.hasValue()); + const auto& rotatedFrame = frameFor(rotated.value(), 0U); + expectVectorNear(rotatedFrame.director, {1.0, 0.0, 0.0}); + expectVectorNear(rotatedFrame.tangentA, {0.0, 1.0, 0.0}); + expectVectorNear(rotatedFrame.tangentB, {0.0, 0.0, 1.0}); + expectRightHandedFrame(rotatedFrame); + + const std::vector warpedNodes{ + node(0U, {0.0, 0.0, 0.0}), + node(1U, {2.0, 0.0, 0.0}), + node(2U, {2.0, 1.0, 0.2}), + node(3U, {0.0, 1.0, 0.0})}; + auto warped = fesa::preprocessShellGeometry( + warpedNodes, {element(12U, {0U, 1U, 2U, 3U})}, sections()); + + ASSERT_TRUE(warped.hasValue()); + EXPECT_GT(warped.value().elementData[0].surfaceAreaWeight, 2.0); + for (const auto& frame : warped.value().nodalFrames) { + expectRightHandedFrame(frame); + } +} + +// MITC4-GEO-002 +TEST(Mitc4Geometry, AreaWeightsSharedDirectorsInStableSourceIdentityOrder) { + const std::vector nodes{ + node(0U, {0.0, 0.0, 0.0}), + node(1U, {1.0, 0.0, 0.0}), + node(2U, {1.0, 1.0, 0.0}), + node(3U, {0.0, 1.0, 0.0}), + node(4U, {2.0, 0.0, 1.0}), + node(5U, {2.0, 1.0, 1.0})}; + const auto flat = element(10U, {0U, 1U, 2U, 3U}); + const auto tilted = element(20U, {1U, 4U, 5U, 2U}); + + auto first = fesa::preprocessShellGeometry( + nodes, {tilted, flat}, sections()); + auto second = fesa::preprocessShellGeometry( + nodes, {flat, tilted}, sections()); + + ASSERT_TRUE(first.hasValue()); + ASSERT_TRUE(second.hasValue()); + const Vector3 expectedSharedDirector{ + -1.0 / std::sqrt(5.0), 0.0, 2.0 / std::sqrt(5.0)}; + for (const auto sharedNode : {1U, 2U}) { + const auto& firstFrame = frameFor(first.value(), sharedNode); + const auto& secondFrame = frameFor(second.value(), sharedNode); + expectVectorNear(firstFrame.director, expectedSharedDirector); + expectVectorNear(firstFrame.director, secondFrame.director, 0.0); + expectVectorNear(firstFrame.tangentA, {0.0, 1.0, 0.0}); + expectRightHandedFrame(firstFrame); + } +} + +// MITC4-GEO-003 +TEST(Mitc4Geometry, RejectsInvalidSurfaceJacobianAndIncidentOrientationCases) { + const auto validElement = element(10U, {0U, 1U, 2U, 3U}); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), node(1U, {0.0, 0.0, 0.0}), + node(2U, {1.0, 1.0, 0.0}), node(3U, {0.0, 1.0, 0.0})}, + {validElement}, sections()), + "invalid-shell-geometry"); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), node(1U, {1.0, 1.0, 0.0}), + node(2U, {0.0, 1.0, 0.0}), node(3U, {1.0, 0.0, 0.0})}, + {validElement}, sections()), + "invalid-shell-geometry"); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), node(1U, {1.0, 0.0, 0.0}), + node(2U, {2.0, 0.0, 0.0}), node(3U, {3.0, 0.0, 0.0})}, + {validElement}, sections()), + "invalid-shell-geometry"); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), node(1U, {1.0, 0.0, 0.0}), + node(2U, {0.05, 0.05, 0.0}), node(3U, {0.0, 1.0, 0.0})}, + {validElement}, sections()), + "invalid-shell-geometry"); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), + node(1U, {1.0, 0.0, 0.0}), + node(2U, {1.0, std::numeric_limits::quiet_NaN(), 0.0}), + node(3U, {0.0, 1.0, 0.0})}, + {validElement}, sections()), + "invalid-shell-geometry"); + + expectFailureCode( + fesa::preprocessShellGeometry( + {node(0U, {0.0, 0.0, 0.0}), node(1U, {1.0, 0.0, 0.0}), + node(2U, {1.0, 1.0, 0.0}), node(3U, {0.0, 1.0, 0.0})}, + {validElement}, sections(0.0)), + "invalid-shell-jacobian"); + + const std::vector opposedNodes{ + node(0U, {0.0, 0.0, 0.0}), node(1U, {1.0, 0.0, 0.0}), + node(2U, {1.0, 1.0, 0.0}), node(3U, {0.0, 1.0, 0.0})}; + expectFailureCode( + fesa::preprocessShellGeometry( + opposedNodes, + {element(10U, {0U, 1U, 2U, 3U}), + element(20U, {0U, 3U, 2U, 1U})}, + sections()), + "opposed-incident-normal"); +} + +// MITC4-GEO-004 +TEST(Mitc4Geometry, ExposesTheCompleteRequiredValidationPointInventory) { + const auto& points = fesa::shellGeometryValidationPoints(); + ASSERT_EQ(points.size(), 17U); + EXPECT_EQ( + std::count_if(points.begin(), points.end(), [](const auto& point) { + return point.kind == fesa::ShellGeometryPointKind::center; + }), + 1); + EXPECT_EQ( + std::count_if(points.begin(), points.end(), [](const auto& point) { + return point.kind == fesa::ShellGeometryPointKind::stiffness; + }), + 8); + EXPECT_EQ( + std::count_if(points.begin(), points.end(), [](const auto& point) { + return point.kind == fesa::ShellGeometryPointKind::tying; + }), + 4); + EXPECT_EQ( + std::count_if(points.begin(), points.end(), [](const auto& point) { + return point.kind == fesa::ShellGeometryPointKind::recovery; + }), + 4); + + EXPECT_EQ(points.front().naturalCoordinates, (Vector3{0.0, 0.0, 0.0})); + const double gauss = 1.0 / std::sqrt(3.0); + EXPECT_EQ(points[1].naturalCoordinates, (Vector3{-gauss, -gauss, -gauss})); + EXPECT_EQ(points[8].naturalCoordinates, (Vector3{-gauss, gauss, gauss})); + EXPECT_EQ(points[9].naturalCoordinates, (Vector3{0.0, -1.0, 0.0})); + EXPECT_EQ(points[12].naturalCoordinates, (Vector3{1.0, 0.0, 0.0})); + EXPECT_EQ(points[13].naturalCoordinates, (Vector3{-gauss, -gauss, 0.0})); + EXPECT_EQ(points[16].naturalCoordinates, (Vector3{-gauss, gauss, 0.0})); +}