feat(linear-static-mitc4-shell): step 2 - shell-director-geometry

include/fesa/model/shell_geometry.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "fesa/core/status.hpp"
+#include "fesa/model/model_types.hpp"
+
+#include <array>
+#include <cstddef>
+#include <vector>
+
+namespace fesa {
+
+enum class ShellGeometryPointKind {
+    center,
+    stiffness,
+    tying,
+    recovery
+};
+
+struct ShellGeometryValidationPoint {
+    ShellGeometryPointKind kind;
+    std::size_t locationIndex;
+    std::array<double, 3> naturalCoordinates;
+};
+
+struct ShellElementGeometryData {
+    EntityIndex elementIndex;
+    std::array<double, 3> normalCandidate;
+    double surfaceAreaWeight;
+};
+
+struct ShellGeometry {
+    std::vector<ShellNodeInitialFrame> nodalFrames;
+    std::vector<ShellElementGeometryData> elementData;
+};
+
+const std::array<ShellGeometryValidationPoint, 17>&
+shellGeometryValidationPoints() noexcept;
+
+Result<ShellGeometry> preprocessShellGeometry(
+    const std::vector<Node>& nodes,
+    const std::vector<Mitc4ShellDefinition>& elements,
+    const std::vector<ShellSection>& sections);
+
+} // namespace fesa

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
 )

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 <algorithm>
 #include <array>
 #include <cerrno>
@@ -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;

src/fesa/model/shell_geometry.cpp

@@ -0,0 +1,491 @@
+#include "fesa/model/shell_geometry.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <tuple>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+using Vector3 = std::array<double, 3>;
+
+constexpr std::array<double, 4> kXiSigns{-1.0, 1.0, 1.0, -1.0};
+constexpr std::array<double, 4> kEtaSigns{-1.0, -1.0, 1.0, 1.0};
+
+struct ShapeData {
+    std::array<double, 4> values;
+    std::array<double, 4> xiDerivatives;
+    std::array<double, 4> etaDerivatives;
+};
+
+struct ElementWork {
+    std::array<Vector3, 4> 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<double, 4>& weights,
+    const std::array<Vector3, 4>& 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<double, 4>& derivatives,
+    const std::array<Vector3, 4>& coordinates) {
+    // Shape derivatives sum to zero, so translating by node 1 improves the
+    // numerical cancellation without changing the covariant tangent.
+    std::array<Vector3, 4> relative{};
+    for (std::size_t node = 0U; node < coordinates.size(); ++node) {
+        relative[node] = subtract(coordinates[node], coordinates[0]);
+    }
+    return weightedSum(derivatives, relative);
+}
+
+Result<ShellGeometry> geometryFailure(
+    std::string code,
+    const SourceLocation& location,
+    std::string keyword,
+    std::string identity,
+    std::string message) {
+    return Result<ShellGeometry>::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<ShellGeometryValidationPoint, 17>&
+shellGeometryValidationPoints() noexcept {
+    static const std::array<ShellGeometryValidationPoint, 17> points = [] {
+        const double g = 1.0 / std::sqrt(3.0);
+        return std::array<ShellGeometryValidationPoint, 17>{
+            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<ShellGeometry> preprocessShellGeometry(
+    const std::vector<Node>& nodes,
+    const std::vector<Mitc4ShellDefinition>& elements,
+    const std::vector<ShellSection>& sections) {
+    ShellGeometry geometry;
+    geometry.elementData.reserve(elements.size());
+    std::vector<ElementWork> work;
+    work.reserve(elements.size());
+
+    const double g = 1.0 / std::sqrt(3.0);
+    const std::array<Vector3, 4> 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<EntityIndex>(elementIndex), current.normal, areaWeight});
+    }
+
+    std::vector<std::vector<std::size_t>> 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<const ShellNodeInitialFrame*> 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<Vector3, 3> 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<EntityIndex>(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<Vector3, 4> 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<ShellGeometry>::success(std::move(geometry));
+}
+
+} // namespace fesa

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

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<double, 3>{0.0, 0.0, 1.0}));
+    EXPECT_EQ(
+        domain.shellNodeInitialFrames()[0].tangentA,
+        (std::array<double, 3>{1.0, 0.0, 0.0}));
+    EXPECT_EQ(
+        domain.shellNodeInitialFrames()[0].tangentB,
+        (std::array<double, 3>{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);

tests/unit/model/shell_geometry_test.cpp

@@ -0,0 +1,276 @@
+#include "fesa/model/shell_geometry.hpp"
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <string>
+#include <vector>
+
+namespace {
+
+using Vector3 = std::array<double, 3>;
+
+fesa::Node node(fesa::EntityIndex label, Vector3 coordinates) {
+    return {
+        {"Shell-Instance", static_cast<std::int64_t>(label),
+         std::to_string(label)},
+        coordinates,
+        {"shell-geometry.inp", static_cast<std::size_t>(label + 1U)}};
+}
+
+fesa::Mitc4ShellDefinition element(
+    fesa::EntityIndex label,
+    std::array<fesa::EntityIndex, 4> nodeIndices) {
+    return {
+        {"Shell-Instance", static_cast<std::int64_t>(label),
+         std::to_string(label)},
+        fesa::ShellSourceElementType::s4,
+        nodeIndices,
+        0U,
+        0U,
+        {"shell-geometry.inp", static_cast<std::size_t>(100U + label)}};
+}
+
+std::vector<fesa::ShellSection> 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<fesa::ShellGeometry>& 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<fesa::Node> 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<fesa::Node> 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<fesa::Node> 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<fesa::Node> 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<double>::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<fesa::Node> 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}));
+}