feat(linear-static-mitc4-shell): step 7 - shell-sparse-assembly

src/fesa/assembly/sparse_assembler.cpp

@@ -3,6 +3,7 @@
 #include "fesa/analysis/analysis_model.hpp"
 #include "fesa/assembly/parallel_for.hpp"
 #include "fesa/elements/euler_beam_3d.hpp"
+#include "fesa/elements/mitc4_shell.hpp"
 #include "fesa/fem/dof_manager.hpp"
 
 #include <array>
@@ -17,11 +18,17 @@ namespace fesa {
 namespace {
 
 constexpr std::size_t kDofsPerNode = 6U;
-constexpr std::size_t kElementDofCount = 12U;
+constexpr std::size_t kBeamElementDofCount = 12U;
-constexpr std::size_t kContributionCount =
+constexpr std::size_t kBeamContributionCount =
-    kElementDofCount * kElementDofCount;
+    kBeamElementDofCount * kBeamElementDofCount;
+constexpr std::size_t kShellElementDofCount = 24U;
+constexpr std::size_t kShellContributionCount =
+    kShellElementDofCount * kShellElementDofCount;
 
-using ElementBuffer = std::array<CooContribution, kContributionCount>;
+using BeamElementBuffer =
+    std::array<CooContribution, kBeamContributionCount>;
+using ShellElementBuffer =
+    std::array<CooContribution, kShellContributionCount>;
 
 Result<SparseMatrix> assemblyFailure(
     const std::string& code,
@@ -54,8 +61,178 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
             std::to_string(dofs.fullDofCount()),
             "DofManager dimensions do not match the active model nodes.");
     }
+    if (!model.activeElements().empty() && !domain.shellElements().empty()) {
+        return assemblyFailure(
+            "unsupported-mixed-element-model",
+            {domain.sourcePath(), 0U},
+            "B33:FESA-MITC4",
+            "Sparse assembly does not support mixed beam and shell models.");
+    }
+
+    if (!domain.shellElements().empty()) {
+        if (domain.shellElements().size() >
+            (std::numeric_limits<std::size_t>::max)() /
+                kShellContributionCount) {
+            return assemblyFailure(
+                "invalid-assembly-dimensions",
+                {domain.sourcePath(), 0U},
+                std::to_string(domain.shellElements().size()),
+                "Shell contribution storage exceeds the addressable range.");
+        }
+
+        std::vector<std::optional<std::array<double, 3>>> directorsByNode(
+            domain.nodes().size());
+        for (const auto& frame : domain.shellNodeInitialFrames()) {
+            if (frame.nodeIndex >= directorsByNode.size() ||
+                directorsByNode[frame.nodeIndex]) {
+                return assemblyFailure(
+                    "invalid-assembly-element",
+                    {domain.sourcePath(), 0U},
+                    std::to_string(frame.nodeIndex),
+                    "Shell initial frames must map uniquely to model nodes.");
+            }
+            directorsByNode[frame.nodeIndex] = frame.director;
+        }
+
+        struct ShellInput {
+            std::array<const Node*, 4> nodes;
+            std::array<std::array<double, 3>, 4> directors;
+            const ShellSection* section;
+            const LinearElasticMaterial* material;
+            std::array<std::size_t, kShellElementDofCount> scatter;
+        };
+        std::vector<ShellInput> inputs;
+        inputs.reserve(domain.shellElements().size());
+        for (std::size_t elementOrder = 0U;
+             elementOrder < domain.shellElements().size();
+             ++elementOrder) {
+            const auto& element = domain.shellElements()[elementOrder];
+            if (element.materialIndex >= domain.materials().size() ||
+                element.sectionIndex >= domain.shellSections().size()) {
+                return assemblyFailure(
+                    "invalid-assembly-element",
+                    element.location,
+                    element.sourceId.sourceLabelText,
+                    "Shell element references an entity outside the Domain.");
+            }
+
+            ShellInput input{};
+            input.section = &domain.shellSections()[element.sectionIndex];
+            input.material = &domain.materials()[element.materialIndex];
+            try {
+                input.scatter = dofs.shellElementScatter(
+                    static_cast<EntityIndex>(elementOrder));
+            } catch (const std::out_of_range&) {
+                return assemblyFailure(
+                    "invalid-assembly-scatter",
+                    element.location,
+                    element.sourceId.sourceLabelText,
+                    "DofManager does not contain the active shell scatter.");
+            }
+            for (std::size_t nodePosition = 0U;
+                 nodePosition < element.nodeIndices.size();
+                 ++nodePosition) {
+                const EntityIndex nodeIndex = element.nodeIndices[nodePosition];
+                if (nodeIndex >= domain.nodes().size() ||
+                    !directorsByNode[nodeIndex]) {
+                    return assemblyFailure(
+                        "invalid-assembly-element",
+                        element.location,
+                        element.sourceId.sourceLabelText,
+                        "Shell element requires a valid node and initial director.");
+                }
+                input.nodes[nodePosition] = &domain.nodes()[nodeIndex];
+                input.directors[nodePosition] = *directorsByNode[nodeIndex];
+                for (std::size_t component = 0U;
+                     component < kDofsPerNode;
+                     ++component) {
+                    const std::size_t local =
+                        nodePosition * kDofsPerNode + component;
+                    const std::size_t expected =
+                        static_cast<std::size_t>(nodeIndex) * kDofsPerNode +
+                        component;
+                    if (input.scatter[local] != expected ||
+                        input.scatter[local] >= dofs.fullDofCount()) {
+                        return assemblyFailure(
+                            "invalid-assembly-scatter",
+                            element.location,
+                            element.sourceId.sourceLabelText,
+                            "Shell scatter does not match the active model topology.");
+                    }
+                }
+            }
+            inputs.push_back(input);
+        }
+
+        std::vector<ShellElementBuffer> localBuffers(inputs.size());
+        std::vector<std::optional<Status>> localFailures(inputs.size());
+        parallelFor.execute(
+            inputs.size(),
+            [&](const std::size_t elementOrder) {
+                const auto& input = inputs[elementOrder];
+                const auto shell = Mitc4Shell::create(
+                    input.nodes,
+                    input.directors,
+                    *input.section,
+                    *input.material);
+                if (!shell.hasValue()) {
+                    localFailures[elementOrder] = shell.status();
+                    return;
+                }
+                const auto stiffness = shell.value().stiffness();
+                if (!stiffness.hasValue()) {
+                    localFailures[elementOrder] = stiffness.status();
+                    return;
+                }
+
+                auto& buffer = localBuffers[elementOrder];
+                for (std::size_t localRow = 0U;
+                     localRow < kShellElementDofCount;
+                     ++localRow) {
+                    for (std::size_t localColumn = 0U;
+                         localColumn < kShellElementDofCount;
+                         ++localColumn) {
+                        const std::size_t localOrder =
+                            localRow * kShellElementDofCount + localColumn;
+                        buffer[localOrder] = {
+                            input.scatter[localRow],
+                            input.scatter[localColumn],
+                            stiffness.value().stabilizedGlobal24(
+                                localRow, localColumn),
+                            elementOrder,
+                            localOrder};
+                    }
+                }
+            });
+
+        for (std::size_t elementOrder = 0U;
+             elementOrder < localFailures.size();
+             ++elementOrder) {
+            if (localFailures[elementOrder]) {
+                return Result<SparseMatrix>::failure(
+                    *localFailures[elementOrder]);
+            }
+        }
+
+        std::vector<CooContribution> contributions;
+        contributions.reserve(
+            localBuffers.size() * kShellContributionCount);
+        // Flatten in source-element order after workers complete. The canonical
+        // COO reduction remains the sole writer of global CSR values.
+        for (const auto& buffer : localBuffers) {
+            contributions.insert(
+                contributions.end(), buffer.begin(), buffer.end());
+        }
+        return SparseMatrix::fromCoo(
+            dofs.fullDofCount(),
+            dofs.fullDofCount(),
+            std::move(contributions),
+            dofs.sparsePattern());
+    }
+
     if (model.activeElements().size() >
-        (std::numeric_limits<std::size_t>::max)() / kContributionCount) {
+        (std::numeric_limits<std::size_t>::max)() /
+            kBeamContributionCount) {
         return assemblyFailure(
             "invalid-assembly-dimensions",
             {domain.sourcePath(), 0U},
@@ -63,7 +240,7 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
             "Element contribution storage exceeds the addressable range.");
     }
 
-    std::vector<std::array<std::size_t, kElementDofCount>> scatters;
+    std::vector<std::array<std::size_t, kBeamElementDofCount>> scatters;
     scatters.reserve(model.activeElements().size());
     for (const EntityIndex elementIndex : model.activeElements()) {
         if (elementIndex >= domain.elements().size()) {
@@ -85,7 +262,7 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
                 "Element references an entity outside the Domain.");
         }
 
-        std::array<std::size_t, kElementDofCount> scatter{};
+        std::array<std::size_t, kBeamElementDofCount> scatter{};
         try {
             scatter = dofs.elementScatter(elementIndex);
         } catch (const std::out_of_range&) {
@@ -117,7 +294,7 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
         scatters.push_back(scatter);
     }
 
-    std::vector<ElementBuffer> localBuffers(model.activeElements().size());
+    std::vector<BeamElementBuffer> localBuffers(model.activeElements().size());
     std::vector<std::optional<Status>> localFailures(
         model.activeElements().size());
     parallelFor.execute(
@@ -139,13 +316,13 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
             auto& buffer = localBuffers[elementOrder];
             const auto& scatter = scatters[elementOrder];
             for (std::size_t localRow = 0U;
-                 localRow < kElementDofCount;
+                 localRow < kBeamElementDofCount;
                  ++localRow) {
                 for (std::size_t localColumn = 0U;
-                     localColumn < kElementDofCount;
+                     localColumn < kBeamElementDofCount;
                      ++localColumn) {
                     const std::size_t localOrder =
-                        localRow * kElementDofCount + localColumn;
+                        localRow * kBeamElementDofCount + localColumn;
                     buffer[localOrder] = {
                         scatter[localRow],
                         scatter[localColumn],
@@ -167,7 +344,7 @@ Result<SparseMatrix> SparseAssembler::assembleStiffness(
 
     std::vector<CooContribution> contributions;
     contributions.reserve(
-        localBuffers.size() * kContributionCount);
+        localBuffers.size() * kBeamContributionCount);
     // Flatten only after all workers complete; workers never share CSR state.
     for (const auto& buffer : localBuffers) {
         contributions.insert(

tests/unit/assembly/sparse_assembler_test.cpp

@@ -1,12 +1,14 @@
 #include "fesa/analysis/analysis_model.hpp"
 #include "fesa/assembly/parallel_for.hpp"
 #include "fesa/assembly/sparse_assembler.hpp"
+#include "fesa/elements/mitc4_shell.hpp"
 #include "fesa/fem/dof_manager.hpp"
 #include "fesa/model/domain.hpp"
 
 #include <gtest/gtest.h>
 
 #include <algorithm>
+#include <array>
 #include <cstring>
 #include <filesystem>
 #include <utility>
@@ -43,6 +45,103 @@ fesa::ModelDefinition makeDefinition() {
     return definition;
 }
 
+fesa::ModelDefinition makeShellDefinition(
+    const fesa::ShellSourceElementType sourceType,
+    const bool twoElements = false) {
+    const std::filesystem::path source{"models/shell-sparse-assembly.inp"};
+    fesa::ModelDefinition definition{};
+    definition.sourcePath = source;
+    definition.sourceContentIdentity = "fnv1a64:fedcba9876543210";
+    if (twoElements) {
+        definition.nodes = {
+            {{"Shell-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 10U}},
+            {{"Shell-1", 2, "2"}, {1.0, 0.0, 0.0}, {source, 11U}},
+            {{"Shell-1", 3, "3"}, {2.0, 0.0, 0.0}, {source, 12U}},
+            {{"Shell-1", 4, "4"}, {0.0, 1.0, 0.0}, {source, 13U}},
+            {{"Shell-1", 5, "5"}, {1.0, 1.0, 0.0}, {source, 14U}},
+            {{"Shell-1", 6, "6"}, {2.0, 1.0, 0.0}, {source, 15U}}};
+        for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+            definition.shellNodeInitialFrames.push_back({
+                static_cast<fesa::EntityIndex>(node),
+                {0.0, 0.0, 1.0},
+                {1.0, 0.0, 0.0},
+                {0.0, 1.0, 0.0}});
+        }
+        definition.shellElements = {
+            {{"Shell-1", 10, "10"}, sourceType, {0U, 1U, 4U, 3U},
+             0U, 0U, {source, 40U}},
+            {{"Shell-1", 20, "20"}, sourceType, {1U, 2U, 5U, 4U},
+             0U, 0U, {source, 41U}}};
+    } else {
+        // A YZ-plane fixture catches any accidental global-Z director assumption.
+        definition.nodes = {
+            {{"Shell-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 10U}},
+            {{"Shell-1", 2, "2"}, {0.0, 1.0, 0.0}, {source, 11U}},
+            {{"Shell-1", 3, "3"}, {0.0, 1.0, 1.0}, {source, 12U}},
+            {{"Shell-1", 4, "4"}, {0.0, 0.0, 1.0}, {source, 13U}}};
+        for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+            definition.shellNodeInitialFrames.push_back({
+                static_cast<fesa::EntityIndex>(node),
+                {1.0, 0.0, 0.0},
+                {0.0, 1.0, 0.0},
+                {0.0, 0.0, 1.0}});
+        }
+        definition.shellElements = {{
+            {"Shell-1", 10, "10"}, sourceType, {0U, 1U, 2U, 3U},
+            0U, 0U, {source, 40U}}};
+    }
+    definition.materials = {
+        {"Material", 120.0, 0.25, {source, 20U}}};
+    definition.shellSections = {
+        {"ShellSection", 0.2, 0U, {source, 30U}}};
+    definition.steps = {{
+        "Step-1", {}, {}, 0.1, 1.0, 0.01, 1.0, {source, 50U}}};
+    return definition;
+}
+
+fesa::Result<fesa::Mitc4Stiffness> directShellStiffness(
+    const fesa::Domain& domain,
+    const fesa::EntityIndex elementIndex) {
+    const auto& definition = domain.shellElements().at(elementIndex);
+    std::array<const fesa::Node*, 4> nodes{};
+    std::array<std::array<double, 3>, 4> directors{};
+    for (std::size_t node = 0U; node < definition.nodeIndices.size(); ++node) {
+        const fesa::EntityIndex nodeIndex = definition.nodeIndices[node];
+        nodes[node] = &domain.nodes().at(nodeIndex);
+        directors[node] = domain.shellNodeInitialFrames().at(nodeIndex).director;
+    }
+    auto shell = fesa::Mitc4Shell::create(
+        nodes,
+        directors,
+        domain.shellSections().at(definition.sectionIndex),
+        domain.materials().at(definition.materialIndex));
+    if (!shell.hasValue()) {
+        return fesa::Result<fesa::Mitc4Stiffness>::failure(shell.status());
+    }
+    return shell.value().stiffness();
+}
+
+fesa::Result<fesa::SparseMatrix> assembleShell(
+    const fesa::ShellSourceElementType sourceType,
+    const fesa::ParallelFor& parallelFor,
+    const bool twoElements = false) {
+    auto domain = fesa::Domain::create(
+        makeShellDefinition(sourceType, twoElements));
+    if (!domain.hasValue()) {
+        return fesa::Result<fesa::SparseMatrix>::failure(domain.status());
+    }
+    auto model = fesa::AnalysisModel::create(domain.value());
+    if (!model.hasValue()) {
+        return fesa::Result<fesa::SparseMatrix>::failure(model.status());
+    }
+    auto dofs = fesa::DofManager::create(model.value());
+    if (!dofs.hasValue()) {
+        return fesa::Result<fesa::SparseMatrix>::failure(dofs.status());
+    }
+    return fesa::SparseAssembler::assembleStiffness(
+        model.value(), dofs.value(), parallelFor);
+}
+
 template<class T>
 bool byteIdentical(const std::vector<T>& left, const std::vector<T>& right) {
     return left.size() == right.size() &&
@@ -155,4 +254,86 @@ TEST(SparseAssembly, SerialTbbAndRepeatedRunsAreByteIdentical) {
     EXPECT_NEAR(entry(serial.value(), 12U, 12U), 80.0, 1.0e-12);
 }
 
+TEST(
+    SparseAssembly,
+    AssemblesFourNodeTwentyFourDofKernelAndPreservesDiagonalSlots) {
+    auto domain = fesa::Domain::create(
+        makeShellDefinition(fesa::ShellSourceElementType::s4));
+    ASSERT_TRUE(domain.hasValue());
+    auto model = fesa::AnalysisModel::create(domain.value());
+    ASSERT_TRUE(model.hasValue());
+    auto dofs = fesa::DofManager::create(model.value());
+    ASSERT_TRUE(dofs.hasValue());
+    fesa::SerialParallelFor serialExecutor;
+
+    auto assembled = fesa::SparseAssembler::assembleStiffness(
+        model.value(), dofs.value(), serialExecutor);
+    auto expected = directShellStiffness(domain.value(), 0U);
+    ASSERT_TRUE(assembled.hasValue());
+    ASSERT_TRUE(expected.hasValue());
+
+    EXPECT_EQ(assembled.value().rows(), 24U);
+    EXPECT_EQ(assembled.value().columns(), 24U);
+    EXPECT_EQ(assembled.value().values().size(), 24U * 24U);
+    EXPECT_EQ(
+        assembled.value().rowOffsets(),
+        dofs.value().sparsePattern().rowOffsets);
+    EXPECT_EQ(
+        assembled.value().columnIndices(),
+        dofs.value().sparsePattern().columnIndices);
+    for (std::size_t row = 0U; row < 24U; ++row) {
+        const auto begin = assembled.value().columnIndices().begin() +
+                           assembled.value().rowOffsets()[row];
+        const auto end = assembled.value().columnIndices().begin() +
+                         assembled.value().rowOffsets()[row + 1U];
+        EXPECT_NE(std::lower_bound(begin, end, row), end);
+        for (std::size_t column = 0U; column < 24U; ++column) {
+            EXPECT_DOUBLE_EQ(
+                entry(assembled.value(), row, column),
+                expected.value().stabilizedGlobal24(row, column));
+        }
+    }
+}
+
+TEST(SparseAssembly, ShellSerialTbbReverseAndRepeatedRunsAreByteIdentical) {
+    fesa::SerialParallelFor serialExecutor;
+    fesa::TbbParallelFor tbbExecutor;
+    ReverseParallelFor reverseExecutor;
+    auto serial = assembleShell(
+        fesa::ShellSourceElementType::s4, serialExecutor, true);
+    auto tbb = assembleShell(
+        fesa::ShellSourceElementType::s4, tbbExecutor, true);
+    auto reversed = assembleShell(
+        fesa::ShellSourceElementType::s4, reverseExecutor, true);
+    ASSERT_TRUE(serial.hasValue());
+    ASSERT_TRUE(tbb.hasValue());
+    ASSERT_TRUE(reversed.hasValue());
+
+    EXPECT_EQ(reverseExecutor.calls(), 1U);
+    EXPECT_EQ(reverseExecutor.observedCount(), 2U);
+    expectByteIdentical(tbb.value(), serial.value());
+    expectByteIdentical(reversed.value(), serial.value());
+    for (std::size_t repetition = 0U; repetition < 8U; ++repetition) {
+        auto repeated = assembleShell(
+            fesa::ShellSourceElementType::s4, tbbExecutor, true);
+        ASSERT_TRUE(repeated.hasValue());
+        expectByteIdentical(repeated.value(), serial.value());
+    }
+}
+
+TEST(SparseAssembly, S4AndS4rSemanticFixturesAssembleIdenticalStiffness) {
+    fesa::SerialParallelFor serialExecutor;
+    auto s4 = assembleShell(
+        fesa::ShellSourceElementType::s4, serialExecutor);
+    auto s4r = assembleShell(
+        fesa::ShellSourceElementType::s4r, serialExecutor);
+    ASSERT_TRUE(s4.hasValue());
+    ASSERT_TRUE(s4r.hasValue());
+    EXPECT_TRUE(std::any_of(
+        s4.value().values().begin(),
+        s4.value().values().end(),
+        [](const double value) { return value != 0.0; }));
+    expectByteIdentical(s4r.value(), s4.value());
+}
+
 } // namespace