KOKO\Mimi
eaa5cbeb21
Integrate Hex8 mass verification and fixed-null-basis shell dynamics with numerical and physics handoff evidence. Preserve interrupted reference and full-validation status; full integration validation is deferred until after performance optimization.
615 lines
28 KiB
C++
615 lines
28 KiB
C++
#include <algorithm>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <optional>
|
|
#include <string>
|
|
#include <utility>
|
|
|
|
#include "analysis/hht_dynamic_internal.h"
|
|
#include "fesa/assembly/parallel_for.h"
|
|
#include "fesa/elements/finite_rotation.h"
|
|
#include "fesa/loads/load.h"
|
|
#include "fesa/model/source_target_resolver.h"
|
|
|
|
namespace fesa::hht_internal {
|
|
namespace {
|
|
|
|
/// @brief Reports a dynamic assembly or equilibrium failure.
|
|
Status Failure(const std::string& message) {
|
|
return Status::Failure(FailureCategory::kSolver, {{Severity::kError,
|
|
"dynamic-invalid-state",
|
|
{},
|
|
"DYNAMIC",
|
|
{},
|
|
message}});
|
|
}
|
|
|
|
/// @brief Uses exact zero semantics for a dimensionally consistent ratio.
|
|
double Ratio(const double numerator, const double denominator) {
|
|
return denominator > 0.0 ? numerator / denominator
|
|
: numerator == 0.0 ? 0.0
|
|
: (std::numeric_limits<double>::infinity)();
|
|
}
|
|
|
|
/// @brief Gathers one source-node Cartesian component triple.
|
|
Vector3 Triple(const Vector& vector, const std::size_t first) {
|
|
return {vector[first], vector[first + 1U], vector[first + 2U]};
|
|
}
|
|
|
|
/// @brief Computes one matrix-vector spatial transform.
|
|
Vector3 Transform(const Matrix& matrix, const Vector3& vector) {
|
|
return {matrix(0U, 0U) * vector[0U] + matrix(0U, 1U) * vector[1U] +
|
|
matrix(0U, 2U) * vector[2U],
|
|
matrix(1U, 0U) * vector[0U] + matrix(1U, 1U) * vector[1U] +
|
|
matrix(1U, 2U) * vector[2U],
|
|
matrix(2U, 0U) * vector[0U] + matrix(2U, 1U) * vector[1U] +
|
|
matrix(2U, 2U) * vector[2U]};
|
|
}
|
|
|
|
/// @brief Gathers a local 24-vector in the DOF-owner scatter order.
|
|
Vector Gather(const Vector& full, const RuntimeShell& shell) {
|
|
Vector local(24U);
|
|
for (std::size_t i = 0U; i < 24U; ++i) local[i] = full[shell.scatter[i]];
|
|
return local;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
Result<std::vector<RuntimeShell>> BuildRuntimeShells(const Domain& domain,
|
|
const DofManager& dofs) {
|
|
std::vector<RuntimeShell> shells;
|
|
for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) {
|
|
const auto* definition =
|
|
dynamic_cast<const Mitc4ShellDefinition*>(&domain.Elements()[index]);
|
|
if (definition == nullptr ||
|
|
definition->section_index >= domain.Properties().Size() ||
|
|
definition->material_index >= domain.Materials().Size()) {
|
|
return Result<std::vector<RuntimeShell>>::Failure(
|
|
Failure("HHT requires valid MITC4 shell definitions."));
|
|
}
|
|
const auto* section = dynamic_cast<const ShellSection*>(
|
|
&domain.Properties()[definition->section_index]);
|
|
const auto* material = dynamic_cast<const LinearElasticMaterial*>(
|
|
&domain.Materials()[definition->material_index]);
|
|
if (section == nullptr || material == nullptr) {
|
|
return Result<std::vector<RuntimeShell>>::Failure(
|
|
Failure("HHT requires homogeneous elastic shell properties."));
|
|
}
|
|
double density = 0.0;
|
|
for (const auto& record : domain.Densities()) {
|
|
if (record.material_index == definition->material_index)
|
|
density = record.rho;
|
|
}
|
|
std::array<const Node*, 4> nodes{};
|
|
std::array<ShellNodeInitialFrame, 4> frames{};
|
|
for (std::size_t local = 0U; local < 4U; ++local) {
|
|
const auto node = definition->node_indices[local];
|
|
if (node >= domain.Nodes().size() ||
|
|
node >= domain.ShellNodeInitialFrames().size()) {
|
|
return Result<std::vector<RuntimeShell>>::Failure(
|
|
Failure("A dynamic shell node has no initial frame."));
|
|
}
|
|
nodes[local] = &domain.Nodes()[node];
|
|
frames[local] = domain.ShellNodeInitialFrames()[node];
|
|
}
|
|
auto internal =
|
|
NonlinearMitc4Shell::Create(nodes, frames, *section, *material);
|
|
if (!internal.HasValue())
|
|
return Result<std::vector<RuntimeShell>>::Failure(internal.GetStatus());
|
|
auto inertia = NonlinearMitc4Inertia::Create(nodes, frames, *section,
|
|
*material, density);
|
|
if (!inertia.HasValue())
|
|
return Result<std::vector<RuntimeShell>>::Failure(inertia.GetStatus());
|
|
shells.push_back({static_cast<EntityIndex>(index),
|
|
dofs.ShellElementScatter(static_cast<EntityIndex>(index)),
|
|
std::move(internal.Value()), std::move(inertia.Value())});
|
|
}
|
|
return Result<std::vector<RuntimeShell>>::Success(std::move(shells));
|
|
}
|
|
|
|
Result<double> ReferenceLength(const Domain& domain) {
|
|
double length = 0.0;
|
|
for (std::size_t i = 0U; i < domain.ShellElements().Size(); ++i) {
|
|
const auto& element = domain.ShellElements()[i];
|
|
for (std::size_t node = 0U; node < 4U; ++node) {
|
|
const Vector3 a(domain.Nodes()[element.node_indices[node]].coordinates);
|
|
const Vector3 b(
|
|
domain.Nodes()[element.node_indices[(node + 1U) % 4U]].coordinates);
|
|
length = (std::max)(length, (a - b).Norm());
|
|
}
|
|
}
|
|
if (!std::isfinite(length) || length <= 0.0)
|
|
return Result<double>::Failure(
|
|
Failure("The dynamic reference length must be positive."));
|
|
return Result<double>::Success(length);
|
|
}
|
|
|
|
Result<std::array<Vector, 2>> AssembleLoads(const Domain& domain,
|
|
const DofManager& dofs) {
|
|
const auto index = SourceTargetIndex::FromDomain(domain);
|
|
const SourceTargetResolver resolver(index);
|
|
const LoadContext context{domain, dofs, resolver};
|
|
std::array<Vector, 2> result{Vector(dofs.FullDofCount()),
|
|
Vector(dofs.FullDofCount())};
|
|
for (const auto& load : domain.Steps()[0U].Loads()) {
|
|
const auto contributions = load.get().ComputeContributions(context);
|
|
if (!contributions.HasValue())
|
|
return Result<std::array<Vector, 2>>::Failure(contributions.GetStatus());
|
|
for (const auto& entry : contributions.Value()) {
|
|
result[0U][entry.full_dof_index] += entry.value;
|
|
result[1U][entry.full_dof_index] += std::abs(entry.value);
|
|
if (!std::isfinite(result[0U][entry.full_dof_index]) ||
|
|
!std::isfinite(result[1U][entry.full_dof_index]))
|
|
return Result<std::array<Vector, 2>>::Failure(
|
|
Failure("Dynamic load reduction overflowed."));
|
|
}
|
|
}
|
|
return Result<std::array<Vector, 2>>::Success(std::move(result));
|
|
}
|
|
|
|
Result<SparseMatrix> MakeSparse(const std::size_t rows,
|
|
const std::size_t columns,
|
|
std::vector<CooContribution> contributions) {
|
|
std::vector<std::vector<std::size_t>> row_columns(rows);
|
|
for (const auto& entry : contributions) {
|
|
if (entry.row >= rows || entry.column >= columns)
|
|
return Result<SparseMatrix>::Failure(
|
|
Failure("Projected sparse index is invalid."));
|
|
row_columns[entry.row].push_back(entry.column);
|
|
}
|
|
SparsePattern pattern;
|
|
pattern.row_offsets.push_back(0U);
|
|
for (auto& row : row_columns) {
|
|
std::sort(row.begin(), row.end());
|
|
row.erase(std::unique(row.begin(), row.end()), row.end());
|
|
pattern.column_indices.insert(pattern.column_indices.end(), row.begin(),
|
|
row.end());
|
|
pattern.row_offsets.push_back(pattern.column_indices.size());
|
|
}
|
|
return SparseMatrix::FromCoo(rows, columns, std::move(contributions),
|
|
pattern);
|
|
}
|
|
|
|
Result<SparseMatrix> ScaledFree(const SparseMatrix& matrix,
|
|
const DofManager& dofs, const double length) {
|
|
std::vector<CooContribution> entries;
|
|
for (std::size_t i = 0U; i < dofs.FreeDofCount(); ++i) {
|
|
const auto row = dofs.FreeDofs()[i];
|
|
for (auto offset = matrix.RowOffsets()[row];
|
|
offset < matrix.RowOffsets()[row + 1U]; ++offset) {
|
|
const auto column = matrix.ColumnIndices()[offset];
|
|
const auto free_column = dofs.FreeEquation(column);
|
|
if (free_column.has_value()) {
|
|
const double factor =
|
|
(row % 6U < 3U ? length : 1.0) * (column % 6U < 3U ? length : 1.0);
|
|
entries.push_back(
|
|
{i, *free_column, factor * matrix.Values()[offset], row, offset});
|
|
}
|
|
}
|
|
}
|
|
return MakeSparse(dofs.FreeDofCount(), dofs.FreeDofCount(),
|
|
std::move(entries));
|
|
}
|
|
|
|
Status CheckLinearSolve(const SparseMatrix& matrix, const Vector& rhs,
|
|
const Vector& solution) {
|
|
if (matrix.Rows() != rhs.Size() || matrix.Columns() != solution.Size())
|
|
return Failure("Dynamic linear solve dimensions differ.");
|
|
double numerator = 0.0;
|
|
double denominator = 0.0;
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
double residual = -rhs[row];
|
|
double scale = std::abs(rhs[row]);
|
|
for (auto offset = matrix.RowOffsets()[row];
|
|
offset < matrix.RowOffsets()[row + 1U]; ++offset) {
|
|
const double term =
|
|
matrix.Values()[offset] * solution[matrix.ColumnIndices()[offset]];
|
|
residual += term;
|
|
scale += std::abs(term);
|
|
}
|
|
if (!std::isfinite(residual) || !std::isfinite(scale))
|
|
return Failure("Dynamic linear solve returned nonfinite data.");
|
|
numerator = (std::max)(numerator, std::abs(residual));
|
|
denominator = (std::max)(denominator, scale);
|
|
}
|
|
if (Ratio(numerator, denominator) > 1e-12) {
|
|
return Status::Failure(FailureCategory::kSolver,
|
|
{{Severity::kError,
|
|
"dynamic-linear-solve",
|
|
{},
|
|
"DYNAMIC",
|
|
{},
|
|
"Dynamic linear backward error exceeds 1e-12."}});
|
|
}
|
|
return Status::Ok();
|
|
}
|
|
|
|
Result<DynamicAssembly> EvaluateNewtonStart(
|
|
const Domain& domain, const DofManager& dofs,
|
|
const std::vector<RuntimeShell>& shells, DynamicAcceptedFrameState frame,
|
|
const HhtPredictor& predictor, const Vector& old_acceleration,
|
|
const Vector& load_scale, const double time_increment,
|
|
const HhtParameters& parameters, const ParallelFor& parallel_for) {
|
|
const Vector accepted_coordinates = frame.generalized_coordinates;
|
|
auto evaluate = [&](const Vector& coordinates) -> Result<DynamicAssembly> {
|
|
frame.generalized_coordinates = coordinates;
|
|
auto rates = EvaluateHhtRates(coordinates, predictor, old_acceleration,
|
|
time_increment, parameters);
|
|
if (!rates.HasValue())
|
|
return Result<DynamicAssembly>::Failure(rates.GetStatus());
|
|
frame.velocity = std::move(rates.Value().velocity);
|
|
frame.acceleration = std::move(rates.Value().acceleration);
|
|
frame.algorithmic_acceleration =
|
|
std::move(rates.Value().algorithmic_acceleration);
|
|
frame.iteration_count = 0U;
|
|
return EvaluateAssembly(domain, dofs, shells, frame, load_scale,
|
|
time_increment, parameters, parallel_for, false);
|
|
};
|
|
auto start = evaluate(predictor.coordinates);
|
|
if (start.HasValue()) {
|
|
bool same_coordinates =
|
|
accepted_coordinates.Size() == predictor.coordinates.Size();
|
|
for (std::size_t i = 0U; same_coordinates && i < accepted_coordinates.Size();
|
|
++i)
|
|
same_coordinates = accepted_coordinates[i] == predictor.coordinates[i];
|
|
if (same_coordinates) return start;
|
|
auto accepted_start = evaluate(accepted_coordinates);
|
|
if (!accepted_start.HasValue()) return accepted_start;
|
|
const auto length = ReferenceLength(domain);
|
|
if (!length.HasValue())
|
|
return Result<DynamicAssembly>::Failure(length.GetStatus());
|
|
double predictor_residual = 0.0;
|
|
double accepted_residual = 0.0;
|
|
for (const auto full : dofs.FreeDofs()) {
|
|
const double scale = full % 6U < 3U ? length.Value() : 1.0;
|
|
predictor_residual =
|
|
std::max(predictor_residual,
|
|
scale * std::abs(start.Value().frame.residual[full]));
|
|
accepted_residual =
|
|
std::max(accepted_residual,
|
|
scale * std::abs(accepted_start.Value().frame.residual[full]));
|
|
}
|
|
if (!std::isfinite(predictor_residual) || !std::isfinite(accepted_residual))
|
|
return Result<DynamicAssembly>::Failure(
|
|
Failure("The scaled initial Newton residual must remain finite."));
|
|
// A geometrically valid large DIRECT predictor may still be far from the
|
|
// endpoint. Compare only these two initial trials, in common energy units.
|
|
return predictor_residual <= accepted_residual ? std::move(start)
|
|
: std::move(accepted_start);
|
|
}
|
|
if (start.GetStatus().Diagnostics().empty()) return start;
|
|
const auto& code = start.GetStatus().Diagnostics().front().code;
|
|
bool invalid_geometry = code == "invalid-nonlinear-shell-geometry" ||
|
|
code == "invalid-nonlinear-shell-director" ||
|
|
code == "invalid-nonlinear-inertia-geometry" ||
|
|
code == "invalid-nonlinear-inertia-director";
|
|
if (code == "invalid-dynamic-rotation") {
|
|
// This diagnostic also covers rate overflow. Only a finite coordinate
|
|
// outside the existing principal margin permits another initial geometry.
|
|
const double pi = std::acos(-1.0);
|
|
const double margin = std::sqrt(std::numeric_limits<double>::epsilon());
|
|
for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
|
|
const std::size_t offset = 6U * node + 3U;
|
|
const auto phi = Triple(predictor.coordinates, offset);
|
|
const auto rotation = EvaluateDynamicRotation(
|
|
phi, Vector3(domain.ShellNodeInitialFrames()[node].director),
|
|
Triple(frame.velocity, offset), Triple(frame.acceleration, offset));
|
|
if (rotation.HasValue()) continue;
|
|
// Match assembly's first failing node; an unrelated later chart error
|
|
// must not turn an earlier valid-chart rate failure into a fallback.
|
|
const double angle = phi.Norm();
|
|
invalid_geometry =
|
|
phi.IsFinite() &&
|
|
(!std::isfinite(angle) || !(pi - angle > margin));
|
|
break;
|
|
}
|
|
}
|
|
if (!invalid_geometry) return start;
|
|
// This is initial-iterate selection only: neither a Newton correction nor
|
|
// the frozen HHT predictor is changed, and accepted history is never mutated.
|
|
return evaluate(accepted_coordinates);
|
|
}
|
|
|
|
Result<DynamicAssembly> EvaluateAssembly(
|
|
const Domain& domain, const DofManager& dofs,
|
|
const std::vector<RuntimeShell>& shells, DynamicAcceptedFrameState frame,
|
|
const Vector& load_scale, const double time_increment,
|
|
const HhtParameters& parameters, const ParallelFor& parallel_for,
|
|
const bool directional_third) {
|
|
const std::size_t count = dofs.FullDofCount();
|
|
for (const auto* values : std::array<const Vector*, 6>{
|
|
&frame.generalized_coordinates, &frame.velocity, &frame.acceleration,
|
|
&frame.algorithmic_acceleration, &frame.external_force,
|
|
&load_scale}) {
|
|
if (values->Size() != count)
|
|
return Result<DynamicAssembly>::Failure(
|
|
Failure("Dynamic trial vector size is invalid."));
|
|
for (std::size_t i = 0U; i < count; ++i)
|
|
if (!std::isfinite((*values)[i]))
|
|
return Result<DynamicAssembly>::Failure(
|
|
Failure("Dynamic trial contains nonfinite values."));
|
|
}
|
|
frame.nodal_internal_order.clear();
|
|
frame.nodal_source_ids.clear();
|
|
frame.nodal_triad.clear();
|
|
frame.nodal_director.clear();
|
|
frame.angular_velocity = Vector(3U * domain.Nodes().size());
|
|
frame.angular_acceleration = Vector(3U * domain.Nodes().size());
|
|
std::vector<DynamicRotationKinematics> rotations;
|
|
for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
|
|
const auto phi = Triple(frame.generalized_coordinates, 6U * node + 3U);
|
|
const auto& baseline = domain.ShellNodeInitialFrames()[node];
|
|
auto rotation = EvaluateDynamicRotation(
|
|
phi, Vector3(baseline.director), Triple(frame.velocity, 6U * node + 3U),
|
|
Triple(frame.acceleration, 6U * node + 3U));
|
|
if (!rotation.HasValue())
|
|
return Result<DynamicAssembly>::Failure(rotation.GetStatus());
|
|
const Matrix r = ExpSkew(phi);
|
|
frame.nodal_internal_order.push_back(static_cast<EntityIndex>(node));
|
|
frame.nodal_source_ids.push_back(domain.Nodes()[node].source_id);
|
|
frame.nodal_director.push_back(rotation.Value().director.Components());
|
|
frame.nodal_triad.push_back(
|
|
{Transform(r, Vector3(baseline.tangent_a)).Components(),
|
|
Transform(r, Vector3(baseline.tangent_b)).Components(),
|
|
rotation.Value().director.Components()});
|
|
for (std::size_t c = 0U; c < 3U; ++c) {
|
|
frame.angular_velocity[3U * node + c] =
|
|
rotation.Value().angular_velocity[c];
|
|
frame.angular_acceleration[3U * node + c] =
|
|
rotation.Value().angular_acceleration[c];
|
|
}
|
|
rotations.push_back(std::move(rotation.Value()));
|
|
}
|
|
std::vector<std::optional<NonlinearMitc4InternalEquilibrium>> internal(
|
|
shells.size());
|
|
std::vector<std::optional<NonlinearMitc4InertiaEvaluation>> inertia(
|
|
shells.size());
|
|
std::vector<std::optional<Status>> failures(shells.size());
|
|
parallel_for.Execute(shells.size(), [&](const std::size_t e) {
|
|
const auto q = Gather(frame.generalized_coordinates, shells[e]);
|
|
const auto v = Gather(frame.velocity, shells[e]);
|
|
const auto a = Gather(frame.acceleration, shells[e]);
|
|
NonlinearMitc4TrialDofs trial{};
|
|
for (std::size_t node = 0U; node < 4U; ++node) {
|
|
trial.displacement[node] = Triple(q, 6U * node);
|
|
trial.rotation_vector[node] = Triple(q, 6U * node + 3U);
|
|
}
|
|
auto f = EvaluateInternalResidualAndTangent(
|
|
shells[e].internal, trial, directional_third ? &v : nullptr);
|
|
if (!f.HasValue()) {
|
|
failures[e] = f.GetStatus();
|
|
return;
|
|
}
|
|
auto m = shells[e].inertia.Evaluate(trial, v, a);
|
|
if (!m.HasValue()) {
|
|
failures[e] = m.GetStatus();
|
|
return;
|
|
}
|
|
internal[e] = std::move(f.Value());
|
|
inertia[e] = std::move(m.Value());
|
|
});
|
|
frame.inertia_force = Vector(count);
|
|
frame.internal_force = Vector(count);
|
|
frame.residual = Vector(count);
|
|
frame.reaction = Vector(count);
|
|
frame.residual_scale = load_scale;
|
|
frame.shell_results.clear();
|
|
frame.kinetic_energy = frame.physical_strain_energy = frame.drilling_energy =
|
|
0.0;
|
|
frame.linear_momentum = {};
|
|
frame.angular_momentum = {};
|
|
frame.linear_momentum_rate = {};
|
|
frame.angular_momentum_rate = {};
|
|
std::vector<std::map<std::pair<std::size_t, std::size_t>, double>>
|
|
constraint_hessians(domain.Nodes().size());
|
|
std::vector<std::map<std::size_t, double>>
|
|
constraint_gradients(domain.Nodes().size());
|
|
std::vector<bool> external_done(domain.Nodes().size(), false);
|
|
Vector third(count);
|
|
Vector third_scale(count);
|
|
std::array<std::vector<CooContribution>, 6> entries;
|
|
for (auto& values : entries) values.reserve(shells.size() * 24U * 24U);
|
|
double momentum_scale = 0.0;
|
|
double angular_scale = 0.0;
|
|
const double a_coefficient =
|
|
1.0 / ((1.0 + parameters.alpha) * parameters.beta * time_increment *
|
|
time_increment);
|
|
const double v_coefficient =
|
|
parameters.gamma / (parameters.beta * time_increment);
|
|
for (std::size_t e = 0U; e < shells.size(); ++e) {
|
|
if (failures[e].has_value())
|
|
return Result<DynamicAssembly>::Failure(*failures[e]);
|
|
const auto& f = *internal[e];
|
|
const auto& m = *inertia[e];
|
|
if (directional_third) {
|
|
for (std::size_t local = 0U; local < 4U; ++local) {
|
|
const auto node = shells[e].scatter[6U * local] / 6U;
|
|
auto& gradient = constraint_gradients[node];
|
|
auto& hessian = constraint_hessians[node];
|
|
for (std::size_t i = 0U; i < 24U; ++i)
|
|
gradient[shells[e].scatter[i]] +=
|
|
f.drilling_constraint_gradient(local, i);
|
|
for (std::size_t i = 0U; i < 24U; ++i)
|
|
for (std::size_t j = 0U; j < 24U; ++j) {
|
|
hessian[{shells[e].scatter[i], shells[e].scatter[j]}] +=
|
|
f.drilling_constraint_hessian[local](i, j);
|
|
}
|
|
if (!external_done[node]) {
|
|
external_done[node] = true;
|
|
const auto first = 6U * node + 3U;
|
|
const auto external =
|
|
shells[e].internal.EvaluateExternalDrillingConstraint(
|
|
local, Triple(frame.generalized_coordinates, first),
|
|
Triple(frame.velocity, first),
|
|
Triple(frame.external_force, first), {}, {});
|
|
if (!external.HasValue())
|
|
return Result<DynamicAssembly>::Failure(external.GetStatus());
|
|
for (std::size_t i = 0U; i < 3U; ++i) {
|
|
gradient[first + i] += external.Value().gradient[i];
|
|
for (std::size_t j = 0U; j < 3U; ++j)
|
|
hessian[{first + i, first + j}] +=
|
|
external.Value().hessian(i, j);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
const auto a = Gather(frame.acceleration, shells[e]);
|
|
for (std::size_t row = 0U; row < 24U; ++row) {
|
|
const auto global = shells[e].scatter[row];
|
|
frame.internal_force[global] += f.internal_force[row];
|
|
frame.inertia_force[global] += m.force[row];
|
|
third[global] += f.force_directional_second[row];
|
|
third_scale[global] += std::abs(f.force_directional_second[row]);
|
|
frame.residual_scale[global] += std::abs(m.convective_force[row]) +
|
|
std::abs(f.physical_force[row]) +
|
|
std::abs(f.drilling_force[row]);
|
|
for (std::size_t column = 0U; column < 24U; ++column) {
|
|
frame.residual_scale[global] +=
|
|
std::abs(m.mass(row, column) * a[column]);
|
|
const std::array<double, 6> values{
|
|
m.mass(row, column), m.position_tangent(row, column),
|
|
m.velocity_tangent(row, column), f.internal_tangent(row, column),
|
|
f.internal_tangent(row, column) + m.position_tangent(row, column) +
|
|
v_coefficient * m.velocity_tangent(row, column) +
|
|
a_coefficient * m.mass(row, column),
|
|
f.drilling_metric(row, column)};
|
|
for (std::size_t k = 0U; k < 6U; ++k)
|
|
entries[k].push_back({global, shells[e].scatter[column], values[k], e,
|
|
row * 24U + column});
|
|
}
|
|
}
|
|
for (const auto& point : f.recovery.gauss_points) {
|
|
frame.shell_results.push_back(
|
|
{shells[e].element_index,
|
|
point.location,
|
|
point.natural_coordinates,
|
|
{point.local_frame.e1, point.local_frame.e2, point.local_frame.e3},
|
|
point.generalized_strain,
|
|
point.section_resultant,
|
|
point.stress});
|
|
}
|
|
frame.kinetic_energy += m.kinetic_energy;
|
|
frame.physical_strain_energy += f.physical_strain_energy;
|
|
frame.drilling_energy += f.drilling_energy;
|
|
momentum_scale += m.momentum_rate_scale;
|
|
angular_scale += m.angular_momentum_rate_scale;
|
|
for (std::size_t c = 0U; c < 3U; ++c) {
|
|
frame.linear_momentum[c] += m.momentum[c];
|
|
frame.angular_momentum[c] += m.angular_momentum[c];
|
|
frame.linear_momentum_rate[c] += m.momentum_rate[c];
|
|
frame.angular_momentum_rate[c] += m.angular_momentum_rate[c];
|
|
}
|
|
}
|
|
for (std::size_t i = 0U; i < count; ++i) {
|
|
frame.residual[i] = frame.inertia_force[i] + frame.internal_force[i] -
|
|
frame.external_force[i];
|
|
if (!std::isfinite(frame.residual[i]) ||
|
|
!std::isfinite(frame.residual_scale[i]) || !std::isfinite(third[i]) ||
|
|
!std::isfinite(third_scale[i])) {
|
|
return Result<DynamicAssembly>::Failure(
|
|
Failure("Dynamic force reduction overflowed."));
|
|
}
|
|
}
|
|
for (const auto i : dofs.ConstrainedDofs())
|
|
frame.reaction[i] = frame.residual[i];
|
|
std::array<double, 2> residual_norm{};
|
|
std::array<double, 2> scale_norm{};
|
|
for (const auto i : dofs.FreeDofs()) {
|
|
const auto channel = i % 6U < 3U ? 0U : 1U;
|
|
residual_norm[channel] =
|
|
(std::max)(residual_norm[channel], std::abs(frame.residual[i]));
|
|
scale_norm[channel] =
|
|
(std::max)(scale_norm[channel], frame.residual_scale[i]);
|
|
}
|
|
for (std::size_t c = 0U; c < 2U; ++c)
|
|
frame.free_residual_metrics[c] = Ratio(residual_norm[c], scale_norm[c]);
|
|
auto mass = SparseMatrix::FromCoo(count, count, std::move(entries[0U]),
|
|
dofs.GetSparsePattern());
|
|
auto position = SparseMatrix::FromCoo(count, count, std::move(entries[1U]),
|
|
dofs.GetSparsePattern());
|
|
auto velocity = SparseMatrix::FromCoo(count, count, std::move(entries[2U]),
|
|
dofs.GetSparsePattern());
|
|
auto stiffness = SparseMatrix::FromCoo(count, count, std::move(entries[3U]),
|
|
dofs.GetSparsePattern());
|
|
auto effective = SparseMatrix::FromCoo(count, count, std::move(entries[4U]),
|
|
dofs.GetSparsePattern());
|
|
auto metric = SparseMatrix::FromCoo(count, count, std::move(entries[5U]),
|
|
dofs.GetSparsePattern());
|
|
for (const auto* result : {&mass, &position, &velocity, &stiffness,
|
|
&effective, &metric})
|
|
if (!result->HasValue())
|
|
return Result<DynamicAssembly>::Failure(result->GetStatus());
|
|
Vector hidden_second(domain.Nodes().size());
|
|
Vector hidden_second_scale(domain.Nodes().size());
|
|
for (std::size_t node = 0U; node < constraint_hessians.size(); ++node) {
|
|
for (const auto& entry : constraint_hessians[node]) {
|
|
const double term = entry.second * frame.velocity[entry.first.first] *
|
|
frame.velocity[entry.first.second];
|
|
hidden_second[node] += term;
|
|
hidden_second_scale[node] += std::abs(term);
|
|
}
|
|
}
|
|
std::optional<SparseMatrix> hidden_gradient;
|
|
if (directional_third) {
|
|
std::vector<CooContribution> gradient_entries;
|
|
for (std::size_t node = 0U; node < constraint_gradients.size(); ++node)
|
|
for (const auto& entry : constraint_gradients[node])
|
|
gradient_entries.push_back(
|
|
{node, entry.first, entry.second, node, gradient_entries.size()});
|
|
auto gradient = MakeSparse(domain.Nodes().size(), count,
|
|
std::move(gradient_entries));
|
|
if (!gradient.HasValue())
|
|
return Result<DynamicAssembly>::Failure(gradient.GetStatus());
|
|
hidden_gradient = std::move(gradient.Value());
|
|
}
|
|
return Result<DynamicAssembly>::Success(
|
|
{std::move(frame), std::move(mass.Value()), std::move(position.Value()),
|
|
std::move(velocity.Value()), std::move(stiffness.Value()),
|
|
std::move(effective.Value()), std::move(metric.Value()),
|
|
std::move(third), std::move(third_scale), std::move(rotations),
|
|
momentum_scale, angular_scale, std::move(hidden_second),
|
|
std::move(hidden_second_scale), std::move(hidden_gradient)});
|
|
}
|
|
|
|
Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs,
|
|
const DynamicAssembly& assembly,
|
|
DynamicAcceptedFrameState* frame) {
|
|
static_cast<void>(dofs);
|
|
Vector3 total_force;
|
|
Vector3 total_moment;
|
|
double force_scale = assembly.momentum_rate_scale;
|
|
double moment_scale = assembly.angular_momentum_rate_scale;
|
|
for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
|
|
const auto first = 6U * node;
|
|
const Vector3 arm = Vector3(domain.Nodes()[node].coordinates) +
|
|
Triple(frame->generalized_coordinates, first);
|
|
const Matrix transform = InverseTransposeLeftJacobianSO3(
|
|
Triple(frame->generalized_coordinates, first + 3U));
|
|
for (const Vector* values : {&frame->external_force, &frame->reaction}) {
|
|
const Vector3 force = Triple(*values, first);
|
|
const Vector3 moment = Transform(transform, Triple(*values, first + 3U));
|
|
total_force = total_force + force;
|
|
total_moment = total_moment + arm.Cross(force) + moment;
|
|
force_scale += force.Norm();
|
|
moment_scale += arm.Cross(force).Norm() + moment.Norm();
|
|
}
|
|
}
|
|
const Vector3 force_error =
|
|
total_force - Vector3(frame->linear_momentum_rate);
|
|
const Vector3 moment_error =
|
|
total_moment - Vector3(frame->angular_momentum_rate);
|
|
for (std::size_t c = 0U; c < 3U; ++c) {
|
|
frame->global_equilibrium[c] = force_error[c];
|
|
frame->global_equilibrium[3U + c] = moment_error[c];
|
|
}
|
|
if (!force_error.IsFinite() || !moment_error.IsFinite() ||
|
|
!std::isfinite(force_scale) || !std::isfinite(moment_scale) ||
|
|
Ratio(force_error.Norm(), force_scale) > 1e-8 ||
|
|
Ratio(moment_error.Norm(), moment_scale) > 1e-8)
|
|
return Failure("Dynamic physical force or moment balance exceeds 1e-8.");
|
|
return Status::Ok();
|
|
}
|
|
|
|
} // namespace fesa::hht_internal
|