feat(linear-static-mitc4-shell): step 6 - shell-dof-scatter

include/fesa/fem/dof_manager.hpp

@@ -38,6 +38,8 @@ public:
     std::optional<std::size_t> freeEquation(std::size_t fullDof) const;
     const std::array<std::size_t, 12>& elementScatter(
         EntityIndex element) const;
+    const std::array<std::size_t, 24>& shellElementScatter(
+        EntityIndex element) const;
     const std::vector<std::size_t>& freeDofs() const noexcept;
     const std::vector<std::size_t>& constrainedDofs() const noexcept;
     const Vector& prescribedValues() const noexcept;
@@ -48,6 +50,7 @@ private:
         std::size_t fullDofCount,
         std::vector<std::optional<std::size_t>> freeEquations,
         std::vector<std::array<std::size_t, 12>> elementScatters,
+        std::vector<std::array<std::size_t, 24>> shellElementScatters,
         std::vector<std::size_t> freeDofs,
         std::vector<std::size_t> constrainedDofs,
         Vector prescribedValues,
@@ -56,6 +59,7 @@ private:
     std::size_t fullDofCount_;
     std::vector<std::optional<std::size_t>> freeEquations_;
     std::vector<std::array<std::size_t, 12>> elementScatters_;
+    std::vector<std::array<std::size_t, 24>> shellElementScatters_;
     std::vector<std::size_t> freeDofs_;
     std::vector<std::size_t> constrainedDofs_;
     Vector prescribedValues_;

src/fesa/fem/dof_manager.cpp

@@ -54,17 +54,28 @@ std::vector<EntityIndex> expandBoundaryTarget(
     return {};
 }
 
+template <std::size_t scatterSize>
+void appendScatter(
+    std::vector<std::vector<std::size_t>>& columnsByRow,
+    const std::array<std::size_t, scatterSize>& scatter) {
+    for (const std::size_t row : scatter) {
+        auto& columns = columnsByRow[row];
+        columns.insert(columns.end(), scatter.begin(), scatter.end());
+    }
+}
+
 SparsePattern buildSparsePattern(
     std::size_t fullDofCount,
     const std::vector<EntityIndex>& activeElements,
-    const std::vector<std::array<std::size_t, 12>>& elementScatters) {
+    const std::vector<std::array<std::size_t, 12>>& elementScatters,
+    const std::vector<std::array<std::size_t, 24>>& shellElementScatters) {
     std::vector<std::vector<std::size_t>> columnsByRow(fullDofCount);
     for (const EntityIndex element : activeElements) {
-        const auto& scatter = elementScatters.at(element);
+        appendScatter(columnsByRow, elementScatters.at(element));
-        for (const std::size_t row : scatter) {
+    }
-            auto& columns = columnsByRow[row];
+    // Every shell in the approved single-step shell subset is active.
-            columns.insert(columns.end(), scatter.begin(), scatter.end());
+    for (const auto& scatter : shellElementScatters) {
-        }
+        appendScatter(columnsByRow, scatter);
     }
 
     SparsePattern pattern;
@@ -152,12 +163,36 @@ Result<DofManager> DofManager::create(const AnalysisModel& model) {
         }
     }
 
+    std::vector<std::array<std::size_t, 24>> shellElementScatters(
+        domain.shellElements().size());
+    for (std::size_t elementIndex = 0U;
+         elementIndex < domain.shellElements().size();
+         ++elementIndex) {
+        const auto& element = domain.shellElements()[elementIndex];
+        auto& scatter = shellElementScatters[elementIndex];
+        for (std::size_t nodePosition = 0U;
+             nodePosition < element.nodeIndices.size();
+             ++nodePosition) {
+            const std::size_t node = element.nodeIndices[nodePosition];
+            for (std::size_t component = 0U;
+                 component < dofsPerNode;
+                 ++component) {
+                scatter[nodePosition * dofsPerNode + component] =
+                    node * dofsPerNode + component;
+            }
+        }
+    }
+
     auto pattern = buildSparsePattern(
-        fullCount, model.activeElements(), elementScatters);
+        fullCount,
+        model.activeElements(),
+        elementScatters,
+        shellElementScatters);
     return Result<DofManager>::success(DofManager{
         fullCount,
         std::move(freeEquations),
         std::move(elementScatters),
+        std::move(shellElementScatters),
         std::move(freeDofs),
         std::move(constrainedDofs),
         std::move(prescribedValues),
@@ -195,6 +230,11 @@ const std::array<std::size_t, 12>& DofManager::elementScatter(
     return elementScatters_.at(element);
 }
 
+const std::array<std::size_t, 24>& DofManager::shellElementScatter(
+    EntityIndex element) const {
+    return shellElementScatters_.at(element);
+}
+
 const std::vector<std::size_t>& DofManager::freeDofs() const noexcept {
     return freeDofs_;
 }
@@ -215,6 +255,7 @@ DofManager::DofManager(
     std::size_t fullDofCount,
     std::vector<std::optional<std::size_t>> freeEquations,
     std::vector<std::array<std::size_t, 12>> elementScatters,
+    std::vector<std::array<std::size_t, 24>> shellElementScatters,
     std::vector<std::size_t> freeDofs,
     std::vector<std::size_t> constrainedDofs,
     Vector prescribedValues,
@@ -222,6 +263,7 @@ DofManager::DofManager(
     : fullDofCount_{fullDofCount},
       freeEquations_{std::move(freeEquations)},
       elementScatters_{std::move(elementScatters)},
+      shellElementScatters_{std::move(shellElementScatters)},
       freeDofs_{std::move(freeDofs)},
       constrainedDofs_{std::move(constrainedDofs)},
       prescribedValues_{std::move(prescribedValues)},

tests/unit/constraints/essential_constraints_test.cpp

@@ -42,6 +42,47 @@ fesa::DofManager makeDofs(
     return std::move(dofs.value());
 }
 
+fesa::DofManager makeShellSizedDofs(
+    std::vector<fesa::BoundaryCondition> boundaries) {
+    const std::filesystem::path source{"models/shell-essential-constraints.inp"};
+    fesa::ModelDefinition definition{};
+    definition.sourcePath = source;
+    definition.sourceContentIdentity = "fnv1a64:8877665544332211";
+    definition.nodes = {
+        {{"Shell-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 2U}},
+        {{"Shell-1", 2, "2"}, {1.0, 0.0, 0.0}, {source, 3U}},
+        {{"Shell-1", 3, "3"}, {1.0, 1.0, 0.0}, {source, 4U}},
+        {{"Shell-1", 4, "4"}, {0.0, 1.0, 0.0}, {source, 5U}}};
+    definition.materials = {
+        {"Material", 1000.0, 0.25, {source, 6U}}};
+    definition.shellSections = {
+        {"ShellSection", 0.1, 0U, {source, 7U}}};
+    definition.shellElements = {{
+        {"Shell-1", 1, "1"},
+        fesa::ShellSourceElementType::s4,
+        {0U, 1U, 2U, 3U},
+        0U,
+        0U,
+        {source, 8U}}};
+    definition.steps = {{
+        "Step-1",
+        std::move(boundaries),
+        {},
+        0.1,
+        1.0,
+        0.01,
+        1.0,
+        {source, 10U}}};
+
+    auto domain = fesa::Domain::create(std::move(definition));
+    EXPECT_TRUE(domain.hasValue());
+    auto model = fesa::AnalysisModel::create(domain.value());
+    EXPECT_TRUE(model.hasValue());
+    auto dofs = fesa::DofManager::create(model.value());
+    EXPECT_TRUE(dofs.hasValue());
+    return std::move(dofs.value());
+}
+
 fesa::SparseMatrix makeMatrix(
     const std::size_t rows,
     const std::size_t columns,
@@ -172,6 +213,73 @@ TEST(EssentialConstraints, HandlesNoAllAndMixedConstraints) {
     expectShape(mixed.value().kcc, 2U, 2U);
 }
 
+// MITC4-DOF-003
+TEST(EssentialConstraints, PreservesShellSizedNoMixedAndAllConstraintRoundTrips) {
+    const auto full = makeMatrix(24U, 24U, sequentialDense(24U));
+
+    const auto noConstraints = makeShellSizedDofs({});
+    auto none = fesa::EssentialConstraints::partition(full, noConstraints);
+    ASSERT_TRUE(none.hasValue());
+    expectShape(none.value().kff, 24U, 24U);
+    expectShape(none.value().kfc, 24U, 0U);
+    expectShape(none.value().kcf, 0U, 24U);
+    expectShape(none.value().kcc, 0U, 0U);
+
+    const auto mixedConstraints = makeShellSizedDofs({
+        {"1", 1, 6, 0.0, {{}, 12U}},
+        {"4", 2, 2, 2.5, {{}, 13U}}});
+    auto mixed = fesa::EssentialConstraints::partition(full, mixedConstraints);
+    ASSERT_TRUE(mixed.hasValue());
+    expectShape(mixed.value().kff, 17U, 17U);
+    expectShape(mixed.value().kfc, 17U, 7U);
+    expectShape(mixed.value().kcf, 7U, 17U);
+    expectShape(mixed.value().kcc, 7U, 7U);
+
+    fesa::Vector mixedFull{24U};
+    for (std::size_t index = 0U; index < mixedFull.size(); ++index) {
+        mixedFull[index] = static_cast<double>(index) + 0.5;
+    }
+    for (std::size_t index = 0U;
+         index < mixedConstraints.constrainedDofCount();
+         ++index) {
+        mixedFull[mixedConstraints.constrainedDofs()[index]] =
+            mixedConstraints.prescribedValues()[index];
+    }
+    const auto mixedFree =
+        fesa::EssentialConstraints::gatherFree(mixedFull, mixedConstraints);
+    const auto mixedReconstructed =
+        fesa::EssentialConstraints::reconstructFull(
+            mixedFree, mixedConstraints.prescribedValues(), mixedConstraints);
+    ASSERT_EQ(mixedReconstructed.size(), mixedFull.size());
+    for (std::size_t index = 0U; index < mixedFull.size(); ++index) {
+        EXPECT_DOUBLE_EQ(mixedReconstructed[index], mixedFull[index]);
+    }
+
+    const auto allConstraints = makeShellSizedDofs({
+        {"1", 1, 6, 1.0, {{}, 14U}},
+        {"2", 1, 6, 2.0, {{}, 15U}},
+        {"3", 1, 6, 3.0, {{}, 16U}},
+        {"4", 1, 6, 4.0, {{}, 17U}}});
+    auto all = fesa::EssentialConstraints::partition(full, allConstraints);
+    ASSERT_TRUE(all.hasValue());
+    expectShape(all.value().kff, 0U, 0U);
+    expectShape(all.value().kfc, 0U, 24U);
+    expectShape(all.value().kcf, 24U, 0U);
+    expectShape(all.value().kcc, 24U, 24U);
+
+    const auto allReconstructed =
+        fesa::EssentialConstraints::reconstructFull(
+            fesa::Vector{0U},
+            allConstraints.prescribedValues(),
+            allConstraints);
+    ASSERT_EQ(allReconstructed.size(), 24U);
+    for (std::size_t node = 0U; node < 4U; ++node) {
+        for (std::size_t component = 0U; component < 6U; ++component) {
+            EXPECT_DOUBLE_EQ(allReconstructed[node * 6U + component], node + 1.0);
+        }
+    }
+}
+
 TEST(EssentialConstraints, ReconstructsNonzeroPrescribedValues) {
     const auto dofs = makeDofs({
         {"1", 2, 2, 2.5, {{}, 12U}},

tests/unit/fem/dof_manager_test.cpp

@@ -46,6 +46,32 @@ fesa::ModelDefinition makeDefinition() {
     return definition;
 }
 
+fesa::ModelDefinition makeShellDefinition() {
+    const std::filesystem::path source{"models/shell-dof-manager.inp"};
+    fesa::ModelDefinition definition{};
+    definition.sourcePath = source;
+    definition.sourceContentIdentity = "fnv1a64:1122334455667788";
+    definition.nodes = {
+        {{"Shell-1", 10, "10"}, {0.0, 0.0, 0.0}, {source, 10U}},
+        {{"Shell-1", 20, "20"}, {1.0, 0.0, 0.0}, {source, 11U}},
+        {{"Shell-1", 30, "30"}, {1.0, 1.0, 0.0}, {source, 12U}},
+        {{"Shell-1", 40, "40"}, {0.0, 1.0, 0.0}, {source, 13U}}};
+    definition.materials = {
+        {"Material", 1000.0, 0.25, {source, 20U}}};
+    definition.shellSections = {
+        {"ShellSection", 0.1, 0U, {source, 30U}}};
+    definition.shellElements = {{
+        {"Shell-1", 100, "100"},
+        fesa::ShellSourceElementType::s4,
+        {2U, 0U, 3U, 1U},
+        0U,
+        0U,
+        {source, 40U}}};
+    definition.steps = {{
+        "Step-1", {}, {}, 0.1, 1.0, 0.01, 1.0, {source, 50U}}};
+    return definition;
+}
+
 struct DofFixture {
     fesa::DofManager dofs;
 };
@@ -161,6 +187,42 @@ TEST(DofManager, BuildsTwelveDofScatterAndSortedUniquePattern) {
     }
 }
 
+// MITC4-DOF-001
+TEST(DofManager, BuildsShellScatterInSourceNodeAndComponentOrder) {
+    const auto fixture = makeDofFixture(makeShellDefinition());
+    const auto& dofs = fixture.dofs;
+
+    EXPECT_EQ(dofs.fullDofCount(), 24U);
+    EXPECT_EQ(
+        dofs.shellElementScatter(0U),
+        (std::array<std::size_t, 24>{
+            12U, 13U, 14U, 15U, 16U, 17U,
+            0U, 1U, 2U, 3U, 4U, 5U,
+            18U, 19U, 20U, 21U, 22U, 23U,
+            6U, 7U, 8U, 9U, 10U, 11U}));
+}
+
+// MITC4-DOF-002
+TEST(DofManager, IncludesShellConnectivityInSortedUniqueSparsePattern) {
+    const auto fixture = makeDofFixture(makeShellDefinition());
+    const auto& dofs = fixture.dofs;
+    const auto& pattern = dofs.sparsePattern();
+
+    ASSERT_EQ(pattern.rowOffsets.size(), 25U);
+    ASSERT_EQ(pattern.columnIndices.size(), 24U * 24U);
+    for (std::size_t row = 0U; row < dofs.fullDofCount(); ++row) {
+        EXPECT_EQ(pattern.rowOffsets[row], row * 24U);
+        EXPECT_EQ(pattern.rowOffsets[row + 1U], (row + 1U) * 24U);
+        const auto columns = rowColumns(pattern, row);
+        ASSERT_EQ(columns.size(), 24U);
+        for (std::size_t column = 0U; column < columns.size(); ++column) {
+            EXPECT_EQ(columns[column], column);
+        }
+        EXPECT_TRUE(std::binary_search(columns.begin(), columns.end(), row));
+        EXPECT_EQ(std::adjacent_find(columns.begin(), columns.end()), columns.end());
+    }
+}
+
 TEST(DofManager, ReconstructsFullReducedRoundTrip) {
     const auto fixture = makeDofFixture();
     const auto& dofs = fixture.dofs;