baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(cpp-object-oriented-modular-refactoring): step 9 - result-io-vector3

Browse Source
This commit is contained in:
KOKO\Mimi
2026-08-16 07:40:16 +09:00
parent 1cb1f26cdc
commit 0be8d1bd89
6 changed files with 140 additions and 148 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/fesa/io/abaqus/domain_mapper.cpp
+20 -39
View File
@@ -14,6 +14,7 @@
#include <utility>
#include <vector>
#include "fesa/math/vector3.h"
#include "fesa/model/shell_geometry.h"
namespace fesa {
@@ -1800,61 +1801,41 @@ class MappingContext {
bool ValidateGeometry(const RawElement& raw, const Node& first,
const Node& second, const GeneralBeamSection& section) {
const auto norm = [](const std::array<double, 3>& vector) {
return std::hypot(vector[0], vector[1], vector[2]);
};
const auto maximum_absolute = [](const std::array<double, 3>& vector) {
const auto maximum_absolute = [](const Vector3& vector) {
return std::max(
{std::abs(vector[0]), std::abs(vector[1]), std::abs(vector[2])});
};
const Vector3 first_position{first.coordinates};
const Vector3 second_position{second.coordinates};
// Compare both approved inequalities after a common scaling. This
// preserves the exact ratios while avoiding overflow in x*x and in
// subtraction between large finite coordinates.
const double global_coordinate_scale =
std::max({1.0, maximum_absolute(first.coordinates),
maximum_absolute(second.coordinates)});
std::array<double, 3> first_scaled{};
std::array<double, 3> second_scaled{};
std::array<double, 3> delta_scaled{};
for (std::size_t coordinate = 0U; coordinate < 3U; ++coordinate) {
first_scaled[coordinate] =
first.coordinates[coordinate] / global_coordinate_scale;
second_scaled[coordinate] =
second.coordinates[coordinate] / global_coordinate_scale;
delta_scaled[coordinate] =
second_scaled[coordinate] - first_scaled[coordinate];
}
const double length_ratio = norm(delta_scaled);
std::max({1.0, maximum_absolute(first_position),
maximum_absolute(second_position)});
const Vector3 first_scaled = first_position / global_coordinate_scale;
const Vector3 second_scaled = second_position / global_coordinate_scale;
const Vector3 delta_scaled = second_scaled - first_scaled;
const double length_ratio = delta_scaled.Norm();
const double coordinate_norm_ratio =
std::max({1.0 / global_coordinate_scale, norm(first_scaled),
norm(second_scaled)});
std::max({1.0 / global_coordinate_scale, first_scaled.Norm(),
second_scaled.Norm()});
if (!(length_ratio > 1.0e-12 * coordinate_norm_ratio)) {
return ModelFailure(
"invalid-beam-length", raw.location, "ELEMENT", raw.label_text,
"Beam length fails the approved scale-aware threshold.");
}
std::array<double, 3> tangent{delta_scaled[0] / length_ratio,
delta_scaled[1] / length_ratio,
delta_scaled[2] / length_ratio};
const Vector3 tangent = delta_scaled / length_ratio;
const double global_guide_scale =
std::max(1.0, maximum_absolute(section.first_axis));
std::array<double, 3> guide_scaled{};
for (std::size_t coordinate = 0U; coordinate < 3U; ++coordinate) {
guide_scaled[coordinate] =
section.first_axis[coordinate] / global_guide_scale;
}
const double projection = guide_scaled[0] * tangent[0] +
guide_scaled[1] * tangent[1] +
guide_scaled[2] * tangent[2];
std::array<double, 3> perpendicular{
guide_scaled[0] - projection * tangent[0],
guide_scaled[1] - projection * tangent[1],
guide_scaled[2] - projection * tangent[2]};
const Vector3 guide{section.first_axis};
const double global_guide_scale = std::max(1.0, maximum_absolute(guide));
const Vector3 guide_scaled = guide / global_guide_scale;
const double projection = guide_scaled.Dot(tangent);
const Vector3 perpendicular = guide_scaled - projection * tangent;
const double guide_norm_ratio =
std::max(1.0 / global_guide_scale, norm(guide_scaled));
if (!(norm(perpendicular) > 1.0e-12 * guide_norm_ratio)) {
std::max(1.0 / global_guide_scale, guide_scaled.Norm());
if (!(perpendicular.Norm() > 1.0e-12 * guide_norm_ratio)) {
return ModelFailure(
"invalid-beam-guide-vector", section.location, "BEAM GENERAL SECTION",
raw.label_text,
src/fesa/io/hdf5/hdf5_results_writer.cpp
+29 -43
View File
@@ -20,6 +20,7 @@
#include "fesa/analysis/analysis_model.h"
#include "fesa/build_info.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"
namespace fesa {
namespace {
@@ -163,11 +164,6 @@ Status OutputFailure(const std::string& code, const std::string& message) {
{{Severity::kError, code, {}, "", "", message}});
}
bool IsFinite(const std::array<double, 3>& values) {
return std::all_of(values.begin(), values.end(),
[](const double value) { return std::isfinite(value); });
}
template <std::size_t Size>
bool IsFinite(const std::array<double, Size>& values) {
return std::all_of(values.begin(), values.end(),
@@ -244,26 +240,21 @@ const char* ShellSourceTypeName(const ShellSourceElementType type) {
bool IsOrthonormalRightHanded(
const std::array<std::array<double, 3>, 3>& frame) {
constexpr double kTolerance = 1.0e-12;
const auto dot = [](const std::array<double, 3>& left,
const std::array<double, 3>& right) {
return left[0U] * right[0U] + left[1U] * right[1U] + left[2U] * right[2U];
};
for (const auto& axis : frame) {
if (!IsFinite(axis) || std::abs(dot(axis, axis) - 1.0) > kTolerance) {
const std::array<Vector3, 3> axes = {Vector3{frame[0U]}, Vector3{frame[1U]},
Vector3{frame[2U]}};
for (const Vector3& axis : axes) {
if (!axis.IsFinite() || std::abs(axis.Dot(axis) - 1.0) > kTolerance) {
return false;
}
}
if (std::abs(dot(frame[0U], frame[1U])) > kTolerance ||
std::abs(dot(frame[0U], frame[2U])) > kTolerance ||
std::abs(dot(frame[1U], frame[2U])) > kTolerance) {
if (std::abs(axes[0U].Dot(axes[1U])) > kTolerance ||
std::abs(axes[0U].Dot(axes[2U])) > kTolerance ||
std::abs(axes[1U].Dot(axes[2U])) > kTolerance) {
return false;
}
const std::array<double, 3> cross = {
frame[0U][1U] * frame[1U][2U] - frame[0U][2U] * frame[1U][1U],
frame[0U][2U] * frame[1U][0U] - frame[0U][0U] * frame[1U][2U],
frame[0U][0U] * frame[1U][1U] - frame[0U][1U] * frame[1U][0U]};
return dot(cross, frame[2U]) > 0.0 &&
std::abs(dot(cross, frame[2U]) - 1.0) <= kTolerance;
const Vector3 cross = axes[0U].Cross(axes[1U]);
const double handedness = cross.Dot(axes[2U]);
return handedness > 0.0 && std::abs(handedness - 1.0) <= kTolerance;
}
bool ComputeLocalAxes(const Domain& domain,
@@ -275,30 +266,25 @@ bool ComputeLocalAxes(const Domain& domain,
}
const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates;
const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates;
const auto& guide = domain.Sections()[element.section_index].first_axis;
const std::array<double, 3> delta = {
second[0U] - first[0U], second[1U] - first[1U], second[2U] - first[2U]};
const double length = std::hypot(delta[0U], delta[1U], delta[2U]);
if (!IsFinite(first) || !IsFinite(second) || !IsFinite(guide) ||
!IsFinite(delta) || !std::isfinite(length) || !(length > 0.0)) {
const Vector3 first_position{first};
const Vector3 second_position{second};
const Vector3 guide{domain.Sections()[element.section_index].first_axis};
const Vector3 delta = second_position - first_position;
const double length = delta.Norm();
if (!first_position.IsFinite() || !second_position.IsFinite() ||
!guide.IsFinite() || !delta.IsFinite() || !std::isfinite(length) ||
!(length > 0.0)) {
return false;
}
const std::array<double, 3> x = {delta[0U] / length, delta[1U] / length,
delta[2U] / length};
const double projection =
guide[0U] * x[0U] + guide[1U] * x[1U] + guide[2U] * x[2U];
const std::array<double, 3> y_trial = {guide[0U] - projection * x[0U],
guide[1U] - projection * x[1U],
guide[2U] - projection * x[2U]};
const double y_norm = std::hypot(y_trial[0U], y_trial[1U], y_trial[2U]);
if (!IsFinite(y_trial) || !std::isfinite(y_norm) || !(y_norm > 0.0)) {
const Vector3 x = delta / length;
const double projection = guide.Dot(x);
const Vector3 y_trial = guide - projection * x;
const double y_norm = y_trial.Norm();
if (!y_trial.IsFinite() || !std::isfinite(y_norm) || !(y_norm > 0.0)) {
return false;
}
const std::array<double, 3> y = {y_trial[0U] / y_norm, y_trial[1U] / y_norm,
y_trial[2U] / y_norm};
const std::array<double, 3> z = {x[1U] * y[2U] - x[2U] * y[1U],
x[2U] * y[0U] - x[0U] * y[2U],
x[0U] * y[1U] - x[1U] * y[0U]};
const Vector3 y = y_trial / y_norm;
const Vector3 z = x.Cross(y);
axes = {x[0U], x[1U], x[2U], y[0U], y[1U], y[2U], z[0U], z[1U], z[2U]};
return IsFinite(axes);
}
@@ -417,7 +403,7 @@ Status ValidateShellWriterInput(const Domain& domain,
++position) {
if (row.stress[position].position != positions[position] ||
row.stress[position].zeta != kZeta[position] ||
!IsFinite(row.stress[position].components)) {
!Vector3{row.stress[position].components}.IsFinite()) {
return OutputFailure("invalid-result-rows",
"Shell stress rows must preserve BOTTOM, "
"MIDDLE, TOP identity.");
@@ -426,7 +412,7 @@ Status ValidateShellWriterInput(const Domain& domain,
}
if (!std::isfinite(state.PhysicalStrainEnergy()) ||
!IsFinite(state.Equilibrium()) ||
!IsFinite(state.VerificationMetrics())) {
!Vector3{state.VerificationMetrics()}.IsFinite()) {
return OutputFailure(
"invalid-result-rows",
"Shell energy, equilibrium, and verification metrics must be finite.");
@@ -493,7 +479,7 @@ Status ValidateWriterInput(const std::filesystem::path& output_path,
node.source_id.source_label_text.empty() ||
!IsValidUtf8(node.source_id.instance_name) ||
!IsValidUtf8(node.source_id.source_label_text) ||
!IsFinite(node.coordinates)) {
!Vector3{node.coordinates}.IsFinite()) {
return OutputFailure(
"invalid-result-identity",
"Every node requires finite coordinates and UTF-8 source identity.");
src/fesa/results/result_recovery.cpp
+31 -51
View File
@@ -16,6 +16,7 @@
#include "fesa/elements/euler_beam_3d.h"
#include "fesa/elements/mitc4_shell.h"
#include "fesa/math/vector3.h"
namespace fesa {
namespace {
@@ -389,27 +390,12 @@ Result<std::vector<EntityIndex>> ResolveLoadTarget(const Domain& domain,
"A station-eligibility load target must resolve to a node or node set.");
}
std::array<double, 3> Cross(const std::array<double, 3>& left,
const std::array<double, 3>& right) {
return {left[1U] * right[2U] - left[2U] * right[1U],
left[2U] * right[0U] - left[0U] * right[2U],
left[0U] * right[1U] - left[1U] * right[0U]};
}
double Dot(const std::array<double, 3>& left,
const std::array<double, 3>& right) {
return left[0U] * right[0U] + left[1U] * right[1U] + left[2U] * right[2U];
}
double Norm(const std::array<double, 3>& value) {
return std::hypot(value[0U], value[1U], value[2U]);
}
bool AccumulateVectorAndScale(std::array<double, 3>& total, double& scale,
const std::array<double, 3>& contribution) {
const double magnitude = Norm(contribution);
const Vector3 contribution_vector{contribution};
const double magnitude = contribution_vector.Norm();
const double accumulated_scale = scale + magnitude;
if (!IsFinite(contribution) || !std::isfinite(magnitude) ||
if (!contribution_vector.IsFinite() || !std::isfinite(magnitude) ||
!std::isfinite(accumulated_scale)) {
return false;
}
@@ -427,7 +413,7 @@ bool AccumulateVectorAndScale(std::array<double, 3>& total, double& scale,
double NormalizedBalance(const std::array<double, 3>& balance,
const double scale) {
const double balance_norm = Norm(balance);
const double balance_norm = Vector3{balance}.Norm();
if (!std::isfinite(balance_norm) || !std::isfinite(scale)) {
return (std::numeric_limits<double>::infinity)();
}
@@ -474,10 +460,11 @@ Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs,
}
}
const auto applied_force_moment =
Cross(domain.Nodes()[node].coordinates, nodal_applied_force);
const auto reaction_force_moment =
Cross(domain.Nodes()[node].coordinates, nodal_reaction_force);
const Vector3 position{domain.Nodes()[node].coordinates};
const Vector3 applied_force_moment =
position.Cross(Vector3{nodal_applied_force});
const Vector3 reaction_force_moment =
position.Cross(Vector3{nodal_reaction_force});
for (std::size_t component = 0U; component < 3U; ++component) {
nodal_applied_moment[component] += applied_force_moment[component];
nodal_reaction_moment[component] += reaction_force_moment[component];
@@ -498,24 +485,24 @@ Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs,
}
}
std::array<double, 3> force_balance{};
std::array<double, 3> moment_balance{};
const Vector3 force_balance =
Vector3{applied_force} + Vector3{reaction_force};
const Vector3 moment_balance =
Vector3{applied_moment} + Vector3{reaction_moment};
for (std::size_t component = 0U; component < 3U; ++component) {
force_balance[component] =
applied_force[component] + reaction_force[component];
moment_balance[component] =
applied_moment[component] + reaction_moment[component];
candidate.equilibrium[component] = force_balance[component];
candidate.equilibrium[3U + component] = moment_balance[component];
}
const double force_metric = NormalizedBalance(
force_balance, (std::max)(applied_force_scale, reaction_force_scale));
const double force_metric =
NormalizedBalance(force_balance.Components(),
(std::max)(applied_force_scale, reaction_force_scale));
const double moment_metric = NormalizedBalance(
moment_balance, (std::max)(applied_moment_scale, reaction_moment_scale));
moment_balance.Components(),
(std::max)(applied_moment_scale, reaction_moment_scale));
candidate.verification_metrics = {normalized_residual, force_metric,
moment_metric};
if (!IsFinite(candidate.equilibrium) ||
!IsFinite(candidate.verification_metrics)) {
!Vector3{candidate.verification_metrics}.IsFinite()) {
return RecoveryFailure("nonfinite-recovery-value",
{domain.SourcePath(), 0U}, "global-equilibrium",
"Global equilibrium values and their physical "
@@ -535,34 +522,27 @@ std::optional<AxisSet> LocalAxes(const Domain& domain,
const EulerBeam3DDefinition& element) {
const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates;
const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates;
const std::array<double, 3> delta = {
second[0U] - first[0U], second[1U] - first[1U], second[2U] - first[2U]};
const double length = Norm(delta);
const Vector3 delta = Vector3{second} - Vector3{first};
const double length = delta.Norm();
if (!std::isfinite(length) || !(length > 0.0)) {
return std::nullopt;
}
const std::array<double, 3> ex = {delta[0U] / length, delta[1U] / length,
delta[2U] / length};
const auto& guide = domain.Sections()[element.section_index].first_axis;
const double projection = Dot(guide, ex);
const std::array<double, 3> ey_trial = {guide[0U] - projection * ex[0U],
guide[1U] - projection * ex[1U],
guide[2U] - projection * ex[2U]};
const double ey_norm = Norm(ey_trial);
const Vector3 ex = delta / length;
const Vector3 guide{domain.Sections()[element.section_index].first_axis};
const double projection = guide.Dot(ex);
const Vector3 ey_trial = guide - projection * ex;
const double ey_norm = ey_trial.Norm();
if (!std::isfinite(ey_norm) || !(ey_norm > 0.0)) {
return std::nullopt;
}
const std::array<double, 3> ey = {
ey_trial[0U] / ey_norm, ey_trial[1U] / ey_norm, ey_trial[2U] / ey_norm};
const std::array<double, 3> ez = Cross(ex, ey);
const AxisSet axes = {ex, ey, ez};
const Vector3 ey = ey_trial / ey_norm;
const Vector3 ez = ex.Cross(ey);
const AxisSet axes = {ex.Components(), ey.Components(), ez.Components()};
for (const auto& axis : axes) {
for (const double component : axis) {
if (!std::isfinite(component)) {
if (!Vector3{axis}.IsFinite()) {
return std::nullopt;
}
}
}
return axes;
}
tests/unit/io/abaqus/domain_mapper_test.cpp
+14
View File
@@ -3,6 +3,7 @@
#include <gtest/gtest.h>
#include <algorithm>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iterator>
@@ -13,6 +14,7 @@
#include <vector>
#include "fesa/io/abaqus/input_reader.h"
#include "fesa/math/vector3.h"
namespace {
@@ -239,6 +241,7 @@ RootAssembly, 6, -12.5
} // namespace
// C-DUP-002
TEST(InpDomainMapping, MapsEverySupportedKeywordAndLegacyDeck) {
auto result = MapText("supported-inventory", SupportedInventoryDeck(true));
ASSERT_TRUE(result.HasValue());
@@ -425,6 +428,17 @@ OnlySecond, 2, 5.
"2, 1., 0., 0.", "2, 1e308, 1e297, 0."),
"0., 1., 0.", "1e308, 0., 0."));
ASSERT_TRUE(large_finite.HasValue());
const auto& large_domain = large_finite.Value();
ASSERT_EQ(large_domain.Nodes().size(), 2U);
const fesa::Vector3 large_delta =
fesa::Vector3{large_domain.Nodes()[1U].coordinates} -
fesa::Vector3{large_domain.Nodes()[0U].coordinates};
EXPECT_TRUE(large_delta.IsFinite());
EXPECT_TRUE(std::isfinite(large_delta.Norm()));
EXPECT_DOUBLE_EQ(large_delta.Y(), 1.0e297);
const fesa::Vector3 large_guide{large_domain.Sections()[0U].first_axis};
EXPECT_TRUE(large_guide.IsFinite());
EXPECT_TRUE(std::isfinite(large_guide.Norm()));
}
TEST(InpDomainMapping, KeepsNodeAndElementSetNamesInSeparateNamespaces) {
tests/unit/io/hdf5/hdf5_results_writer_test.cpp
+38 -3
View File
@@ -24,6 +24,7 @@
#include "fesa/analysis/analysis_state.h"
#include "fesa/build_info.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"
#include "fesa/model/domain.h"
namespace {
@@ -1114,7 +1115,7 @@ TEST(Hdf5ResultsWriter, SuccessfullyReplacesExistingFinal) {
EXPECT_EQ(ReadUint64Attribute(metadata.Get(), "schema_version"), 0U);
}
// MITC4-H5-001
// MITC4-H5-001, C-DUP-002
TEST(Hdf5ResultsWriter, WritesExactShellMetadataAndModelIdentity) {
TempDirectory directory{"shell-model"};
const auto source = directory.Path() / "shell.inp";
@@ -1150,9 +1151,15 @@ TEST(Hdf5ResultsWriter, WritesExactShellMetadataAndModelIdentity) {
EXPECT_EQ(DatasetDimensions(file.Get(), "/model/shell/nodal_director"),
std::vector<hsize_t>({4U, 3U}));
constexpr fesa::Vector3 director{0.0, 0.0, 1.0};
std::vector<double> expected_directors;
for (std::size_t node = 0U; node < 4U; ++node) {
expected_directors.insert(expected_directors.end(),
director.Components().begin(),
director.Components().end());
}
EXPECT_EQ(ReadDoubleDataset(file.Get(), "/model/shell/nodal_director"),
std::vector<double>(
{0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0}));
expected_directors);
EXPECT_EQ(DatasetDimensions(file.Get(), "/model/shell/nodal_frame"),
std::vector<hsize_t>({4U, 3U, 3U}));
ExpectCompoundMemberNames(file.Get(), "/model/shell/materials",
@@ -1305,3 +1312,31 @@ TEST(Hdf5ResultsWriter, InvalidShellInventoryPreservesExistingFinal) {
EXPECT_EQ(ReadBytes(final), sentinel);
EXPECT_EQ(EntryCount(directory.Path()), 1U);
}
// C-DUP-002
TEST(Hdf5ResultsWriter,
Vector3ShellSerializationRejectsNonfiniteWithoutReplacingFinal) {
TempDirectory directory{"shell-vector3-atomic"};
auto fixture = MakeShellFixture(directory.Path() / "shell.inp");
const auto final = directory.Path() / "results.h5";
fesa::Hdf5ResultsWriter writer;
ASSERT_TRUE(writer.Write(final, *fixture.domain, *fixture.state, {}).IsOk());
const std::vector<char> valid_bytes = ReadBytes(final);
const auto file = OpenFile(final);
const auto directors =
ReadDoubleDataset(file.Get(), "/model/shell/nodal_director");
ASSERT_EQ(directors.size(), 12U);
for (std::size_t offset = 0U; offset < directors.size(); offset += 3U) {
const fesa::Vector3 director{directors[offset], directors[offset + 1U],
directors[offset + 2U]};
EXPECT_TRUE(director.IsFinite());
EXPECT_EQ(director, (fesa::Vector3{0.0, 0.0, 1.0}));
}
fixture.state->Displacement()[0U] = std::numeric_limits<double>::quiet_NaN();
ExpectOutputFailure(writer.Write(final, *fixture.domain, *fixture.state, {}),
"invalid-result-state");
EXPECT_EQ(ReadBytes(final), valid_bytes);
EXPECT_EQ(EntryCount(directory.Path()), 1U);
}
tests/unit/results/result_recovery_test.cpp
+7 -11
View File
@@ -20,6 +20,7 @@
#include "fesa/assembly/parallel_for.h"
#include "fesa/assembly/sparse_assembler.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"
#include "fesa/model/domain.h"
namespace {
@@ -719,14 +720,14 @@ TEST(ResultRecovery, KeepsFullResidualAndComputesGlobalShellEquilibrium) {
}
}
// MITC4-REC-004
// MITC4-REC-004, C-DUP-002
TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) {
auto centered_definition = MakeShellDefinition();
auto translated_definition = centered_definition;
constexpr std::array<double, 3> translation{7.0, 11.0, 0.0};
constexpr fesa::Vector3 translation{7.0, 11.0, 0.0};
for (auto& node : translated_definition.nodes) {
for (std::size_t component = 0U; component < translation.size();
++component) {
for (std::size_t component = 0U;
component < translation.Components().size(); ++component) {
node.coordinates[component] += translation[component];
}
}
@@ -753,13 +754,8 @@ TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) {
EXPECT_NEAR(translated.Equilibrium()[component], centered_force[component],
1.0e-12);
}
const std::array<double, 3> translated_moment_delta{
translation[1U] * centered_force[2U] -
translation[2U] * centered_force[1U],
translation[2U] * centered_force[0U] -
translation[0U] * centered_force[2U],
translation[0U] * centered_force[1U] -
translation[1U] * centered_force[0U]};
const fesa::Vector3 translated_moment_delta =
translation.Cross(fesa::Vector3{centered_force});
for (std::size_t component = 0U; component < 3U; ++component) {
EXPECT_NEAR(translated.Equilibrium()[3U + component] -
centered.Equilibrium()[3U + component],