diff --git a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md index 1a8202e..d674912 100644 --- a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md +++ b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md @@ -969,3 +969,55 @@ that the factorized CSR must be nonempty or that empty input is rejected. - concerns: none; the transient build contention was resolved before recorded RED/GREEN/VERIFY runs and no critical or upstream blocker remains. + +## Step 21 — load-assembly + +- task_id: `TASK-21` +- status: `completed` +- changed_files: `include/fesa/assembly/load_assembler.hpp`, + `src/fesa/assembly/load_assembler.cpp`, + `tests/unit/assembly/load_assembler_test.cpp`, `src/fesa/CMakeLists.txt`, + `tests/CMakeLists.txt`, + `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`, + `phases/linear-static-3d-euler-beam/index.json` +- requirement_ids: `FESA-REQ-LS3DEB-007`, `FESA-REQ-LS3DEB-011`, + `FESA-REQ-LS3DEB-027`, `FESA-REQ-LS3DEB-034` +- test_ids: `T21-LOAD-001`, `T21-LOAD-002`, `T21-LOAD-003`, + `T21-LOAD-004`, `T21-LOAD-005` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exactly five planned tests were registered before production and the public load API was absent | MSVC C1083 reported missing `fesa/assembly/load_assembler.hpp` from `load_assembler_test.cpp` | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | Minimum nodal load/effective RHS implementation and all five tests compiled and linked | `load_assembler.cpp`, its test, `fesa_solver.lib`, and `fesa_unit_tests.exe` built without a FESA warning under `/W4 /WX` | +| GREEN-test | `ctest --test-dir .harness/build -C Debug -R LoadAssembly --output-on-failure` | 0 | Node/set loads, source-order accumulation, nonzero prescribed RHS, rejection, and zero cases pass | 5/5 exact `LoadAssembly` tests passed | +| VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | Approved explicit-dependency MSVC x64 build tree generates | Windows SDK 10.0.26100.0, oneMKL 2026.1 ILP64/dynamic, oneTBB, and HDF5 resolved; configure/generate completed | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R LoadAssembly --output-on-failure` | 0 | Focused Step 21 suite remains green | 5/5 exact `LoadAssembly` tests passed | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated suite and all five exact names | 62 tests discovered, including 5 `LoadAssembly` tests | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 62/62 tests passed | +| VERIFY-contract-scans | Exact-test-count, public-backend, factorization/solve/distributed-load, whitespace, and reference diff/status scans | 0 | Load/RHS scope and backend boundary remain isolated | exact tests 5; public backend leaks 0; forbidden scope 0; whitespace clean; reference unchanged | + +- contract_checks: semantic node label/node-set targets resolve case-insensitively + to unique stable node indices. All six global components use DofManager full + order, constrained loads remain in the full vector, and duplicate CLOAD rows + accumulate strictly in active source order; a cancellation fixture proves + that the implementation does not reorder floating-point additions. +- contract_checks: model/full/free/constrained/prescribed dimensions, stable + DOF partition/equation order, active-load order, target indices, component + range, and finite individual inputs are validated before use. Nonfinite + source-order load sums, `Kfc*dc` products/row sums, and final subtraction fail + closed with one structured model diagnostic. +- contract_checks: effective RHS is exactly the stable gather of `Ff` minus + CSR row-order `Kfc*dc`; hand-computed nonzero `dc`, zero-load, zero-free, and + zero-constrained cases pass. No factorization/substitution call, distributed + load object, `*DLOAD`, or line-load behavior was added. +- contract_checks: the public header exposes only the two approved static + `LoadAssembler` functions and no MKL, oneTBB, HDF5, or PARDISO type. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`, + `.harness/build/tests/Debug/fesa_unit_tests.exe` +- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0 +- handoff: Step 22 can consume the full external force and effective free RHS; + Step 24 retains ownership of factorize-before-load orchestration, + substitution, reconstruction, recovery, and output sequencing. +- concerns: none; no critical implementation, environment, backend, or + upstream-contract conflict was found. diff --git a/include/fesa/assembly/load_assembler.hpp b/include/fesa/assembly/load_assembler.hpp new file mode 100644 index 0000000..3c14857 --- /dev/null +++ b/include/fesa/assembly/load_assembler.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/math/sparse_matrix.hpp" +#include "fesa/math/vector.hpp" + +namespace fesa { + +// Assembles only semantic nodal CLOAD records and forms the eliminated free +// right-hand side; stiffness factorization remains an analysis responsibility. +class LoadAssembler { +public: + static Result assembleFullNodalLoad( + const AnalysisModel& model, + const DofManager& dofs); + static Result effectiveFreeRhs( + const Vector& fullLoad, + const SparseMatrix& kfc, + const Vector& prescribedValues, + const DofManager& dofs); +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 679f27a..0157d1a 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -3,6 +3,7 @@ add_library( STATIC analysis/analysis_model.cpp analysis/analysis_state.cpp + assembly/load_assembler.cpp assembly/parallel_for.cpp assembly/sparse_assembler.cpp build_info.cpp diff --git a/src/fesa/assembly/load_assembler.cpp b/src/fesa/assembly/load_assembler.cpp new file mode 100644 index 0000000..695e911 --- /dev/null +++ b/src/fesa/assembly/load_assembler.cpp @@ -0,0 +1,418 @@ +#include "fesa/assembly/load_assembler.hpp" + +#include "fesa/constraints/essential_constraints.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr std::size_t dofsPerNode = 6U; + +Status loadFailure( + const std::string& code, + const SourceLocation& location, + const std::string& keyword, + const std::string& identity, + const std::string& message) { + return Status::failure( + FailureCategory::model, + {{Severity::error, code, location, keyword, identity, message}}); +} + +char asciiLower(const char value) { + if (value >= 'A' && value <= 'Z') { + return static_cast(value + ('a' - 'A')); + } + return value; +} + +bool equalName(const std::string& left, const std::string& right) { + return left.size() == right.size() && + std::equal( + left.begin(), + left.end(), + right.begin(), + [](const char leftValue, const char rightValue) { + return asciiLower(leftValue) == asciiLower(rightValue); + }); +} + +bool tryPositiveInteger(const std::string& text, std::int64_t& value) { + const char* const first = text.data(); + const char* const last = first + text.size(); + const auto parsed = std::from_chars(first, last, value); + return parsed.ec == std::errc{} && parsed.ptr == last && value > 0; +} + +bool isStrictlyIncreasing(const std::vector& values) { + return std::adjacent_find( + values.begin(), + values.end(), + [](const std::size_t left, const std::size_t right) { + return left >= right; + }) == values.end(); +} + +Status validateDofOrder( + const DofManager& dofs, + const std::size_t expectedFullCount, + const SourceLocation& location) { + const std::size_t fullCount = dofs.fullDofCount(); + const auto& freeDofs = dofs.freeDofs(); + const auto& constrainedDofs = dofs.constrainedDofs(); + if (fullCount != expectedFullCount || + freeDofs.size() != dofs.freeDofCount() || + constrainedDofs.size() != dofs.constrainedDofCount() || + dofs.prescribedValues().size() != constrainedDofs.size() || + constrainedDofs.size() > fullCount || + freeDofs.size() != fullCount - constrainedDofs.size()) { + return loadFailure( + "invalid-load-dimensions", + location, + "LOAD_ASSEMBLER", + std::to_string(fullCount), + "Full, free, constrained, prescribed, and model dimensions must agree."); + } + if (!isStrictlyIncreasing(freeDofs) || + !isStrictlyIncreasing(constrainedDofs)) { + return loadFailure( + "invalid-load-order", + location, + "LOAD_ASSEMBLER", + std::to_string(fullCount), + "Free and constrained DOFs must use stable increasing full-DOF order."); + } + + std::vector ownership(fullCount, 0U); + try { + for (std::size_t equation = 0U; + equation < freeDofs.size(); + ++equation) { + const std::size_t fullDof = freeDofs[equation]; + if (fullDof >= fullCount || ownership[fullDof] != 0U || + dofs.freeEquation(fullDof) != equation) { + return loadFailure( + "invalid-load-order", + location, + "LOAD_ASSEMBLER", + std::to_string(fullDof), + "Free equation numbering must match stable full-DOF order."); + } + ownership[fullDof] = 1U; + } + for (const std::size_t fullDof : constrainedDofs) { + if (fullDof >= fullCount || ownership[fullDof] != 0U || + dofs.freeEquation(fullDof).has_value()) { + return loadFailure( + "invalid-load-order", + location, + "LOAD_ASSEMBLER", + std::to_string(fullDof), + "Constrained DOFs must be unique and absent from free equations."); + } + ownership[fullDof] = 2U; + } + } catch (const std::out_of_range&) { + return loadFailure( + "invalid-load-dimensions", + location, + "LOAD_ASSEMBLER", + std::to_string(fullCount), + "DofManager equation storage must cover every full DOF."); + } + if (std::find(ownership.begin(), ownership.end(), 0U) != ownership.end()) { + return loadFailure( + "invalid-load-order", + location, + "LOAD_ASSEMBLER", + std::to_string(fullCount), + "Free and constrained DOFs must partition the full range."); + } + return Status::ok(); +} + +Result> resolveTarget( + const Domain& domain, + const NodalLoad& load) { + std::vector matchingSets; + for (const auto& set : domain.nodeSets()) { + if (equalName(set.name, load.target)) { + matchingSets.push_back(&set); + } + } + + std::vector matchingNodes; + std::int64_t label = 0; + if (tryPositiveInteger(load.target, label)) { + for (std::size_t index = 0U; index < domain.nodes().size(); ++index) { + if (domain.nodes()[index].sourceId.sourceLabel == label) { + matchingNodes.push_back(static_cast(index)); + } + } + } + + if (matchingSets.size() > 1U || matchingNodes.size() > 1U || + (!matchingSets.empty() && !matchingNodes.empty())) { + return Result>::failure(loadFailure( + "invalid-load-target", + load.location, + "CLOAD", + load.target, + "The load target must resolve unambiguously to one node or one expanded node set.")); + } + if (!matchingSets.empty()) { + const auto& nodes = matchingSets.front()->nodeIndices; + std::vector seen(domain.nodes().size(), 0U); + for (const EntityIndex node : nodes) { + if (node >= domain.nodes().size() || seen[node] != 0U) { + return Result>::failure(loadFailure( + "invalid-load-target", + load.location, + "CLOAD", + load.target, + "The expanded node set must contain unique in-range stable node identities.")); + } + seen[node] = 1U; + } + return Result>::success(nodes); + } + if (!matchingNodes.empty()) { + return Result>::success( + std::move(matchingNodes)); + } + return Result>::failure(loadFailure( + "invalid-load-target", + load.location, + "CLOAD", + load.target, + "The load target must resolve to one semantic node or node set.")); +} + +Status validateFiniteVector( + const Vector& values, + const SourceLocation& location, + const std::string& identity) { + for (std::size_t index = 0U; index < values.size(); ++index) { + if (!std::isfinite(values[index])) { + return loadFailure( + "nonfinite-load-value", + location, + "LOAD_ASSEMBLER", + identity + ":" + std::to_string(index), + "Load and prescribed displacement vectors must contain finite values."); + } + } + return Status::ok(); +} + +} // namespace + +Result LoadAssembler::assembleFullNodalLoad( + const AnalysisModel& model, + const DofManager& dofs) { + const Domain& domain = model.domain(); + if (domain.nodes().size() > + (std::numeric_limits::max)() / dofsPerNode) { + return Result::failure(loadFailure( + "invalid-load-dimensions", + {domain.sourcePath(), 0U}, + "LOAD_ASSEMBLER", + domain.sourceContentIdentity(), + "The semantic node count cannot be represented in full-DOF order.")); + } + const std::size_t expectedFullCount = + domain.nodes().size() * dofsPerNode; + const Status dofStatus = validateDofOrder( + dofs, expectedFullCount, {domain.sourcePath(), 0U}); + if (!dofStatus.isOk()) { + return Result::failure(dofStatus); + } + for (std::size_t node = 0U; node < domain.nodes().size(); ++node) { + for (std::size_t component = 0U; + component < dofsPerNode; + ++component) { + try { + if (dofs.fullDof( + static_cast(node), + static_cast(component)) != + node * dofsPerNode + component) { + return Result::failure(loadFailure( + "invalid-load-order", + domain.nodes()[node].location, + "LOAD_ASSEMBLER", + domain.nodes()[node].sourceId.sourceLabelText, + "DofManager node/component identity must match full-DOF order.")); + } + } catch (const std::out_of_range&) { + return Result::failure(loadFailure( + "invalid-load-dimensions", + domain.nodes()[node].location, + "LOAD_ASSEMBLER", + domain.nodes()[node].sourceId.sourceLabelText, + "DofManager must provide all six DOFs for every semantic node.")); + } + } + } + + const auto& activeLoads = model.activeLoads(); + const auto& loads = model.step().loads; + if (activeLoads.size() != loads.size()) { + return Result::failure(loadFailure( + "invalid-load-order", + model.step().location, + "CLOAD", + model.step().name, + "The active load view must include every sole-step load once.")); + } + + Vector fullLoad{expectedFullCount}; + // Active load indices are required to be the original source order; this + // loop is therefore also the fixed floating-point accumulation order. + for (std::size_t sourceOrder = 0U; + sourceOrder < activeLoads.size(); + ++sourceOrder) { + const EntityIndex loadIndex = activeLoads[sourceOrder]; + if (static_cast(loadIndex) != sourceOrder || + loadIndex >= loads.size()) { + return Result::failure(loadFailure( + "invalid-load-order", + model.step().location, + "CLOAD", + std::to_string(sourceOrder), + "Active loads must retain complete stable source order.")); + } + const auto& load = loads[loadIndex]; + if (load.dof < 1 || load.dof > static_cast(dofsPerNode)) { + return Result::failure(loadFailure( + "invalid-load-dof", + load.location, + "CLOAD", + load.target, + "A nodal load component must be in the range 1 through 6.")); + } + if (!std::isfinite(load.magnitude)) { + return Result::failure(loadFailure( + "nonfinite-load-value", + load.location, + "CLOAD", + load.target, + "A nodal load magnitude must be finite.")); + } + + auto target = resolveTarget(domain, load); + if (!target.hasValue()) { + return Result::failure(target.status()); + } + const auto component = static_cast(load.dof - 1); + for (const EntityIndex node : target.value()) { + const std::size_t fullDof = dofs.fullDof(node, component); + const double accumulated = fullLoad[fullDof] + load.magnitude; + if (!std::isfinite(accumulated)) { + return Result::failure(loadFailure( + "nonfinite-load-accumulation", + load.location, + "CLOAD", + load.target, + "Source-order load accumulation produced a nonfinite value.")); + } + fullLoad[fullDof] = accumulated; + } + } + return Result::success(std::move(fullLoad)); +} + +Result LoadAssembler::effectiveFreeRhs( + const Vector& fullLoad, + const SparseMatrix& kfc, + const Vector& prescribedValues, + const DofManager& dofs) { + const SourceLocation location{{}, 0U}; + const Status dofStatus = + validateDofOrder(dofs, fullLoad.size(), location); + if (!dofStatus.isOk()) { + return Result::failure(dofStatus); + } + if (kfc.rows() != dofs.freeDofCount() || + kfc.columns() != dofs.constrainedDofCount() || + prescribedValues.size() != dofs.constrainedDofCount()) { + return Result::failure(loadFailure( + "invalid-load-dimensions", + location, + "LOAD_ASSEMBLER", + std::to_string(kfc.rows()) + "x" + + std::to_string(kfc.columns()), + "Kfc rows/columns and prescribed values must match free/constrained order.")); + } + const Status matrixStatus = kfc.validate(); + if (!matrixStatus.isOk()) { + return Result::failure(matrixStatus); + } + const Status loadStatus = + validateFiniteVector(fullLoad, location, "full-load"); + if (!loadStatus.isOk()) { + return Result::failure(loadStatus); + } + const Status prescribedStatus = validateFiniteVector( + prescribedValues, location, "prescribed-values"); + if (!prescribedStatus.isOk()) { + return Result::failure(prescribedStatus); + } + + Vector correction{kfc.rows()}; + for (std::size_t row = 0U; row < kfc.rows(); ++row) { + double sum = 0.0; + for (std::size_t position = kfc.rowOffsets()[row]; + position < kfc.rowOffsets()[row + 1U]; + ++position) { + const double product = kfc.values()[position] * + prescribedValues[kfc.columnIndices()[position]]; + if (!std::isfinite(product)) { + return Result::failure(loadFailure( + "nonfinite-load-accumulation", + location, + "LOAD_ASSEMBLER", + std::to_string(row), + "Kfc times prescribed displacement produced a nonfinite product.")); + } + sum += product; + if (!std::isfinite(sum)) { + return Result::failure(loadFailure( + "nonfinite-load-accumulation", + location, + "LOAD_ASSEMBLER", + std::to_string(row), + "Kfc times prescribed displacement produced a nonfinite row sum.")); + } + } + correction[row] = sum; + } + + Vector rhs = EssentialConstraints::gatherFree(fullLoad, dofs); + // The constrained vector is already in DofManager order, so this is the + // approved elimination equation rhs = Ff - Kfc*dc without reordering dc. + for (std::size_t row = 0U; row < rhs.size(); ++row) { + const double value = rhs[row] - correction[row]; + if (!std::isfinite(value)) { + return Result::failure(loadFailure( + "nonfinite-load-accumulation", + location, + "LOAD_ASSEMBLER", + std::to_string(row), + "Effective RHS subtraction produced a nonfinite value.")); + } + rhs[row] = value; + } + return Result::success(std::move(rhs)); +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 36b02a9..7a93343 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,7 @@ add_executable( unit/analysis/analysis_model_test.cpp unit/analysis/analysis_state_test.cpp unit/assembly/parallel_for_test.cpp + unit/assembly/load_assembler_test.cpp unit/assembly/sparse_assembler_test.cpp unit/constraints/essential_constraints_test.cpp unit/core/diagnostic_test.cpp diff --git a/tests/unit/assembly/load_assembler_test.cpp b/tests/unit/assembly/load_assembler_test.cpp new file mode 100644 index 0000000..241d702 --- /dev/null +++ b/tests/unit/assembly/load_assembler_test.cpp @@ -0,0 +1,349 @@ +#include "fesa/assembly/load_assembler.hpp" + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/model/domain.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +struct LoadFixture { + std::unique_ptr domain; + std::unique_ptr model; + std::unique_ptr dofs; +}; + +LoadFixture makeFixture( + const std::size_t nodeCount, + std::vector nodeSets, + std::vector boundaries, + std::vector loads) { + const std::filesystem::path source{"models/load-assembly.inp"}; + fesa::ModelDefinition definition{}; + definition.sourcePath = source; + definition.sourceContentIdentity = "fnv1a64:abcdef0123456789"; + for (std::size_t index = 0U; index < nodeCount; ++index) { + const auto label = static_cast((index + 1U) * 10U); + definition.nodes.push_back({ + {"Beam-1", label, std::to_string(label)}, + {static_cast(index), 0.0, 0.0}, + {source, index + 2U}}); + } + definition.nodeSets = std::move(nodeSets); + definition.steps = {{ + "Step-1", + std::move(boundaries), + std::move(loads), + 0.1, + 1.0, + 0.01, + 1.0, + {source, 20U}}}; + + auto domainResult = fesa::Domain::create(std::move(definition)); + if (!domainResult.hasValue()) { + throw std::runtime_error{"Load fixture Domain construction failed."}; + } + auto domain = std::make_unique( + std::move(domainResult.value())); + + auto modelResult = fesa::AnalysisModel::create(*domain); + if (!modelResult.hasValue()) { + throw std::runtime_error{"Load fixture AnalysisModel construction failed."}; + } + auto model = std::make_unique( + std::move(modelResult.value())); + + auto dofResult = fesa::DofManager::create(*model); + if (!dofResult.hasValue()) { + throw std::runtime_error{"Load fixture DofManager construction failed."}; + } + auto dofs = std::make_unique( + std::move(dofResult.value())); + return {std::move(domain), std::move(model), std::move(dofs)}; +} + +fesa::SparseMatrix makeDenseSparse( + const std::size_t rows, + const std::size_t columns, + const std::vector& values) { + if (values.size() != rows * columns) { + throw std::invalid_argument{"Dense sparse fixture has the wrong value count."}; + } + + fesa::SparsePattern pattern; + std::vector contributions; + pattern.rowOffsets.reserve(rows + 1U); + pattern.rowOffsets.push_back(0U); + for (std::size_t row = 0U; row < rows; ++row) { + for (std::size_t column = 0U; column < columns; ++column) { + pattern.columnIndices.push_back(column); + contributions.push_back({ + row, + column, + values[row * columns + column], + row, + column}); + } + pattern.rowOffsets.push_back(pattern.columnIndices.size()); + } + + auto result = fesa::SparseMatrix::fromCoo( + rows, columns, std::move(contributions), pattern); + if (!result.hasValue()) { + throw std::runtime_error{"Sparse fixture construction failed."}; + } + return std::move(result.value()); +} + +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()[0U].code, code); +} + +} // namespace + +TEST(LoadAssembly, AssemblesNodeSetAndSixComponentLoads) { + const std::filesystem::path source{"models/load-assembly.inp"}; + auto fixture = makeFixture( + 2U, + {{"Pair", std::nullopt, {0U, 1U}, {source, 10U}}}, + {{"10", 1, 1, 0.0, {source, 21U}}}, + {{"pair", 1, 1.0, {source, 30U}}, + {"10", 2, 2.0, {source, 31U}}, + {"20", 3, 3.0, {source, 32U}}, + {"10", 4, -4.0, {source, 33U}}, + {"PAIR", 5, 5.0, {source, 34U}}, + {"20", 6, 6.0, {source, 35U}}}); + + auto result = fesa::LoadAssembler::assembleFullNodalLoad( + *fixture.model, *fixture.dofs); + ASSERT_TRUE(result.hasValue()); + ASSERT_EQ(result.value().size(), 12U); + EXPECT_EQ( + std::vector(result.value().data(), result.value().data() + 12U), + (std::vector{ + 1.0, 2.0, 0.0, -4.0, 5.0, 0.0, + 1.0, 0.0, 3.0, 0.0, 5.0, 6.0})); + EXPECT_EQ(fixture.dofs->constrainedDofs(), + (std::vector{0U})); + EXPECT_DOUBLE_EQ(result.value()[0U], 1.0); +} + +TEST(LoadAssembly, AccumulatesSignedLoadsInSourceOrder) { + const std::filesystem::path source{"models/load-assembly.inp"}; + auto firstOrder = makeFixture( + 1U, + {}, + {}, + {{"10", 1, 1.0e16, {source, 30U}}, + {"10", 1, -1.0e16, {source, 31U}}, + {"10", 1, 1.0, {source, 32U}}}); + auto secondOrder = makeFixture( + 1U, + {}, + {}, + {{"10", 1, 1.0e16, {source, 30U}}, + {"10", 1, 1.0, {source, 31U}}, + {"10", 1, -1.0e16, {source, 32U}}}); + + auto first = fesa::LoadAssembler::assembleFullNodalLoad( + *firstOrder.model, *firstOrder.dofs); + auto second = fesa::LoadAssembler::assembleFullNodalLoad( + *secondOrder.model, *secondOrder.dofs); + ASSERT_TRUE(first.hasValue()); + ASSERT_TRUE(second.hasValue()); + EXPECT_DOUBLE_EQ(first.value()[0U], 1.0); + EXPECT_DOUBLE_EQ(second.value()[0U], 0.0); +} + +TEST(LoadAssembly, FormsNonzeroPrescribedEffectiveRhs) { + const std::filesystem::path source{"models/load-assembly.inp"}; + auto fixture = makeFixture( + 1U, + {}, + {{"10", 2, 2, 2.0, {source, 21U}}, + {"10", 5, 5, -1.0, {source, 22U}}}, + {{"10", 1, 10.0, {source, 30U}}, + {"10", 2, 900.0, {source, 31U}}, + {"10", 3, 20.0, {source, 32U}}, + {"10", 4, 30.0, {source, 33U}}, + {"10", 5, 800.0, {source, 34U}}, + {"10", 6, 40.0, {source, 35U}}}); + auto full = fesa::LoadAssembler::assembleFullNodalLoad( + *fixture.model, *fixture.dofs); + ASSERT_TRUE(full.hasValue()); + const auto kfc = makeDenseSparse( + 4U, + 2U, + {1.0, 2.0, + 3.0, 4.0, + -2.0, 5.0, + 0.5, -1.0}); + + auto rhs = fesa::LoadAssembler::effectiveFreeRhs( + full.value(), kfc, fixture.dofs->prescribedValues(), *fixture.dofs); + ASSERT_TRUE(rhs.hasValue()); + ASSERT_EQ(rhs.value().size(), 4U); + EXPECT_EQ( + std::vector(rhs.value().data(), rhs.value().data() + 4U), + (std::vector{10.0, 18.0, 39.0, 38.0})); +} + +TEST(LoadAssembly, RejectsNonfiniteOrDimensionMismatch) { + const std::filesystem::path source{"models/load-assembly.inp"}; + const double maximum = (std::numeric_limits::max)(); + auto nonfinite = makeFixture( + 1U, + {}, + {}, + {{"10", 1, std::numeric_limits::quiet_NaN(), {source, 30U}}}); + expectFailureCode( + fesa::LoadAssembler::assembleFullNodalLoad( + *nonfinite.model, *nonfinite.dofs), + "nonfinite-load-value"); + + auto overflow = makeFixture( + 1U, + {}, + {}, + {{"10", 1, maximum, {source, 30U}}, + {"10", 1, maximum, {source, 31U}}}); + expectFailureCode( + fesa::LoadAssembler::assembleFullNodalLoad( + *overflow.model, *overflow.dofs), + "nonfinite-load-accumulation"); + + auto oneNode = makeFixture( + 1U, + {}, + {{"10", 2, 2, 2.0, {source, 21U}}, + {"10", 5, 5, -1.0, {source, 22U}}}, + {}); + auto twoNodes = makeFixture(2U, {}, {}, {}); + expectFailureCode( + fesa::LoadAssembler::assembleFullNodalLoad( + *twoNodes.model, *oneNode.dofs), + "invalid-load-dimensions"); + + const auto validKfc = makeDenseSparse(4U, 2U, std::vector(8U, 0.0)); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{5U}, + validKfc, + oneNode.dofs->prescribedValues(), + *oneNode.dofs), + "invalid-load-dimensions"); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{6U}, + makeDenseSparse(3U, 2U, std::vector(6U, 0.0)), + oneNode.dofs->prescribedValues(), + *oneNode.dofs), + "invalid-load-dimensions"); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{6U}, + makeDenseSparse(4U, 1U, std::vector(4U, 0.0)), + oneNode.dofs->prescribedValues(), + *oneNode.dofs), + "invalid-load-dimensions"); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{6U}, validKfc, fesa::Vector{1U}, *oneNode.dofs), + "invalid-load-dimensions"); + + fesa::Vector nonfiniteFull{6U}; + nonfiniteFull[0U] = std::numeric_limits::infinity(); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + nonfiniteFull, + validKfc, + oneNode.dofs->prescribedValues(), + *oneNode.dofs), + "nonfinite-load-value"); + + fesa::Vector nonfinitePrescribed{2U}; + nonfinitePrescribed[0U] = std::numeric_limits::quiet_NaN(); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{6U}, validKfc, nonfinitePrescribed, *oneNode.dofs), + "nonfinite-load-value"); + + const auto overflowingKfc = makeDenseSparse( + 4U, + 2U, + {maximum, 0.0, + 0.0, 0.0, + 0.0, 0.0, + 0.0, 0.0}); + expectFailureCode( + fesa::LoadAssembler::effectiveFreeRhs( + fesa::Vector{6U}, + overflowingKfc, + oneNode.dofs->prescribedValues(), + *oneNode.dofs), + "nonfinite-load-accumulation"); +} + +TEST(LoadAssembly, ZeroLoadsRemainZero) { + const std::filesystem::path source{"models/load-assembly.inp"}; + auto freeFixture = makeFixture( + 1U, + {}, + {}, + {{"10", 3, 0.0, {source, 30U}}}); + auto full = fesa::LoadAssembler::assembleFullNodalLoad( + *freeFixture.model, *freeFixture.dofs); + ASSERT_TRUE(full.hasValue()); + EXPECT_TRUE(std::all_of( + full.value().data(), + full.value().data() + full.value().size(), + [](const double value) { return value == 0.0; })); + const auto noConstrainedColumns = makeDenseSparse(6U, 0U, {}); + auto freeRhs = fesa::LoadAssembler::effectiveFreeRhs( + full.value(), + noConstrainedColumns, + freeFixture.dofs->prescribedValues(), + *freeFixture.dofs); + ASSERT_TRUE(freeRhs.hasValue()); + EXPECT_EQ(freeRhs.value().size(), 6U); + EXPECT_TRUE(std::all_of( + freeRhs.value().data(), + freeRhs.value().data() + freeRhs.value().size(), + [](const double value) { return value == 0.0; })); + + auto constrainedFixture = makeFixture( + 1U, + {}, + {{"10", 1, 6, 0.0, {source, 21U}}}, + {}); + auto constrainedFull = fesa::LoadAssembler::assembleFullNodalLoad( + *constrainedFixture.model, *constrainedFixture.dofs); + ASSERT_TRUE(constrainedFull.hasValue()); + const auto noFreeRows = makeDenseSparse(0U, 6U, {}); + auto constrainedRhs = fesa::LoadAssembler::effectiveFreeRhs( + constrainedFull.value(), + noFreeRows, + constrainedFixture.dofs->prescribedValues(), + *constrainedFixture.dofs); + ASSERT_TRUE(constrainedRhs.hasValue()); + EXPECT_EQ(constrainedRhs.value().size(), 0U); +}