diff --git a/include/fesa/assembly/load_assembler.h b/include/fesa/assembly/load_assembler.h index 9cca81a..de4802f 100644 --- a/include/fesa/assembly/load_assembler.h +++ b/include/fesa/assembly/load_assembler.h @@ -13,17 +13,21 @@ namespace fesa { class LoadAssembler { public: /// @brief Accumulates all active CLOAD rows in full-DOF space. + /// @param step_time Actual finite nonnegative dynamic endpoint for amplitude. /// @return A finite full load vector or a structured model failure. static Result AssembleFullNodalLoad(const AnalysisModel& model, - const DofManager& dofs); + const DofManager& dofs, + double step_time = 0.0); /// @brief Accumulates explicitly supplied loads in their view order. /// @param loads Non-owning loads whose contribution source orders must match /// their view positions. + /// @param step_time Actual finite nonnegative dynamic endpoint for amplitude. /// @return A candidate committed only after all contributions validate. static Result AssembleFullNodalLoad(const AnalysisModel& model, const DofManager& dofs, - const LoadView& loads); + const LoadView& loads, + double step_time = 0.0); /// @brief Forms Ff-Kfc*dc in stable free/constrained order. /// @note This operation neither factorizes nor invokes a solver. diff --git a/include/fesa/elements/element_definition.h b/include/fesa/elements/element_definition.h index ed832a2..5de8bb3 100644 --- a/include/fesa/elements/element_definition.h +++ b/include/fesa/elements/element_definition.h @@ -9,7 +9,12 @@ namespace fesa { /// @brief Identifies the supported semantic element-definition kinds. -enum class ElementDefinitionKind { kEulerBeam3D, kMitc4Shell, kDs4Thermal }; +enum class ElementDefinitionKind { + kEulerBeam3D, + kMitc4Shell, + kDs4Thermal, + kHex8 +}; /// @brief Provides immutable source identity and topology for one element. /// @note Numerical stiffness, recovery, and equation ids are intentionally diff --git a/include/fesa/elements/hex8_definition.h b/include/fesa/elements/hex8_definition.h new file mode 100644 index 0000000..80f27c5 --- /dev/null +++ b/include/fesa/elements/hex8_definition.h @@ -0,0 +1,58 @@ +#ifndef FESA_ELEMENTS_HEX8_DEFINITION_H_ +#define FESA_ELEMENTS_HEX8_DEFINITION_H_ + +#include +#include +#include + +#include "fesa/elements/element_definition.h" + +namespace fesa { +/// @brief Owns source-safe eight-node Solid topology without numerical state. +class Hex8Definition final : public ElementDefinition { + public: + /// @brief Constructs a semantic candidate in source connectivity order. + Hex8Definition(SourceEntityId source, std::array nodes, + EntityIndex material, EntityIndex section, + SourceLocation source_location) + : source_id{std::move(source)}, + node_indices{nodes}, + material_index{material}, + section_index{section}, + location{std::move(source_location)}, + property_index_{section}, + node_indices_view_{nodes.begin(), nodes.end()} {} + ElementDefinitionKind Kind() const noexcept override { + return ElementDefinitionKind::kHex8; + } + const SourceEntityId& SourceId() const noexcept override { return source_id; } + std::string_view SourceElementType() const noexcept override { + return "C3D8"; + } + const std::vector& NodeIndices() const noexcept override { + return node_indices_view_; + } + EntityIndex MaterialIndex() const noexcept override { return material_index; } + EntityIndex PropertyIndex() const noexcept override { + return property_index_; + } + + SourceEntityId source_id; + std::array node_indices; + EntityIndex material_index; + EntityIndex section_index; + SourceLocation location; + + private: + friend class Domain; + /// @brief Synchronizes candidate topology with the immutable base view. + void SynchronizeNodeIndices() { + node_indices_view_.assign(node_indices.begin(), node_indices.end()); + } + /// @brief Binds the local Solid section to the unified property collection. + void SetPropertyIndex(EntityIndex index) noexcept { property_index_ = index; } + EntityIndex property_index_; + std::vector node_indices_view_; +}; +} // namespace fesa +#endif // FESA_ELEMENTS_HEX8_DEFINITION_H_ diff --git a/include/fesa/fem/dof_manager.h b/include/fesa/fem/dof_manager.h index 24a92f5..afb3703 100644 --- a/include/fesa/fem/dof_manager.h +++ b/include/fesa/fem/dof_manager.h @@ -16,8 +16,8 @@ namespace fesa { class DofManagerTestAccess; -/// @brief Distinguishes structural six-component and thermal five-station maps. -enum class DofProfile { kStructural, kThermal }; +/// @brief Distinguishes six-component, Solid translation, and thermal maps. +enum class DofProfile { kStructural, kThermal, kSolidStructural }; /// @brief Stores the stable full-space CSR pattern. struct SparsePattern { @@ -59,6 +59,8 @@ class DofManager { /// @brief Reports which physical component profile owns this mapping. DofProfile Profile() const noexcept; + /// @brief Returns the physical component count for the active profile. + std::size_t ComponentsPerNode() const noexcept; /// @brief Returns the full node-by-component DOF count. std::size_t FullDofCount() const noexcept; diff --git a/include/fesa/loads/concentrated_nodal_load.h b/include/fesa/loads/concentrated_nodal_load.h index ec1f960..ca5c564 100644 --- a/include/fesa/loads/concentrated_nodal_load.h +++ b/include/fesa/loads/concentrated_nodal_load.h @@ -3,6 +3,8 @@ #include #include +#include +#include #include "fesa/core/diagnostic.h" #include "fesa/loads/load.h" @@ -21,7 +23,10 @@ class ConcentratedNodalLoad final : public Load { /// @brief Creates one parsed CLOAD component while preserving diagnostics. ConcentratedNodalLoad(SourceTargetQuery target, int source_dof, double magnitude, std::size_t source_order, - SourceLocation location); + SourceLocation location, + std::optional amplitude_index = {}, + std::string amplitude_name_raw = {}, + SourceLocation amplitude_location = {}); /// @brief Computes target-major, component-minor full-DOF contributions. Result> ComputeContributions( @@ -39,12 +44,28 @@ class ConcentratedNodalLoad final : public Load { /// @brief Returns the source location used by structured diagnostics. const SourceLocation& Location() const noexcept; + /// @brief Returns a model-level amplitude index when one was specified. + std::optional AmplitudeIndex() const noexcept { + return amplitude_index_; + } + /// @brief Returns preserved CLOAD amplitude option spelling. + const std::string& AmplitudeName() const noexcept { + return amplitude_name_raw_; + } + /// @brief Returns the source location of the CLOAD amplitude option. + const SourceLocation& AmplitudeLocation() const noexcept { + return amplitude_location_; + } + private: SourceTargetQuery target_; std::array global_components_{}; std::size_t source_order_; SourceLocation location_{}; int source_dof_{0}; + std::optional amplitude_index_; + std::string amplitude_name_raw_; + SourceLocation amplitude_location_; }; } // namespace fesa diff --git a/include/fesa/loads/load.h b/include/fesa/loads/load.h index f9deba9..e212dff 100644 --- a/include/fesa/loads/load.h +++ b/include/fesa/loads/load.h @@ -26,6 +26,8 @@ struct LoadContext { const Domain& domain; const DofManager& dof_manager; const SourceTargetResolver& target_resolver; + /// @brief Supplies the actual dynamic endpoint time for amplitude sampling. + double step_time{0.0}; }; /// @brief Produces local ordered load contributions without global mutation. diff --git a/include/fesa/model/domain.h b/include/fesa/model/domain.h index b7afd9b..43f8a85 100644 --- a/include/fesa/model/domain.h +++ b/include/fesa/model/domain.h @@ -112,6 +112,13 @@ class StepDefinition { /// @brief Returns the static maximum-increment provenance value. double MaximumIncrement() const noexcept; + /// @brief Reports the source STATIC DIRECT fixed-schedule selection. + bool IsDirectStatic() const noexcept; + /// @brief Returns the fixed Solid static increment count, or zero otherwise. + std::uint32_t IncrementCount() const noexcept; + /// @brief Returns the source static INC limit. + std::uint32_t MaximumIncrements() const noexcept; + /// @brief Returns the source location of the step declaration. const SourceLocation& Location() const noexcept; @@ -145,6 +152,9 @@ class StepDefinition { SourceLocation location_; std::optional frequency_; std::optional dynamic_; + bool direct_static_{false}; + std::uint32_t increment_count_{0U}; + std::uint32_t maximum_increments_{100U}; }; /// @brief Owns the complete immutable semantic model definition. @@ -176,6 +186,15 @@ class Domain { const DomainCollectionView& ShellElements() const noexcept; + /// @brief Returns C3D8 definitions in stable source order. + const DomainCollectionView& SolidElements() const noexcept; + /// @brief Returns homogeneous Solid section assignments. + const DomainCollectionView& SolidSections() const noexcept; + /// @brief Reports whether this model has exclusively Solid elements. + bool IsSolid() const noexcept; + /// @brief Returns immutable model-level amplitude curves in source order. + const std::vector& Amplitudes() const noexcept; + /// @brief Returns all materials in stable Domain index order. const DomainCollectionView& Materials() const noexcept; @@ -262,6 +281,8 @@ class Domain { DomainCollectionView elements_view_; DomainCollectionView beam_elements_view_; DomainCollectionView shell_elements_view_; + DomainCollectionView solid_elements_view_; + DomainCollectionView solid_sections_view_; DomainCollectionView properties_view_; DomainCollectionView sections_view_; DomainCollectionView shell_sections_view_; diff --git a/include/fesa/model/model_types.h b/include/fesa/model/model_types.h index fd5845e..84e4bb9 100644 --- a/include/fesa/model/model_types.h +++ b/include/fesa/model/model_types.h @@ -1,6 +1,7 @@ #ifndef FESA_MODEL_MODEL_TYPES_H_ #define FESA_MODEL_MODEL_TYPES_H_ +#include #include #include #include @@ -14,11 +15,13 @@ #include "fesa/core/source_identity.h" #include "fesa/elements/ds4_thermal.h" #include "fesa/elements/euler_beam_3d.h" +#include "fesa/elements/hex8_definition.h" #include "fesa/elements/mitc4_shell.h" #include "fesa/materials/isotropic_linear_elastic_material.h" #include "fesa/materials/isotropic_thermal_material.h" #include "fesa/properties/general_beam_section.h" #include "fesa/properties/shell_section.h" +#include "fesa/properties/solid_section.h" namespace fesa { @@ -47,12 +50,37 @@ struct PrescribedDisplacementDefinition { SourceLocation location; }; +/// @brief Owns a validated TABULAR STEP TIME RELATIVE load curve. +struct TabularAmplitudeDefinition { + std::string name_raw; + std::string name_canonical; + SourceLocation location; + std::vector> points; + + /// @brief Evaluates exact knots, linear interiors, and endpoint holds. + /// @pre Time is finite; points are finite, nonempty, strictly increasing. + double Evaluate(double time) const { + if (time <= points.front()[0]) return points.front()[1]; + if (time >= points.back()[0]) return points.back()[1]; + const auto upper = std::lower_bound( + points.begin(), points.end(), time, + [](const auto& point, double t) { return point[0] < t; }); + if ((*upper)[0] == time) return (*upper)[1]; + const auto& lower = *(upper - 1); + const double weight = (time - lower[0]) / ((*upper)[0] - lower[0]); + return (1.0 - weight) * lower[1] + weight * (*upper)[1]; + } +}; + /// @brief Stores one concentrated nodal load component. struct NodalLoad { std::string target; int dof; double magnitude; SourceLocation location; + std::optional amplitude_index; + std::string amplitude_name_raw; + SourceLocation amplitude_location; }; /// @brief Identifies the selected exterior shell surface side. @@ -169,6 +197,9 @@ struct StaticStepDefinition { std::vector boundaries; std::vector loads; std::vector surface_pressure_loads; + bool direct{false}; + std::uint32_t increment_count{0U}; + std::uint32_t maximum_increments{100U}; StaticProcedureKind procedure_kind; bool nlgeom; double initial_increment; @@ -226,6 +257,7 @@ struct DynamicStepDefinition { double gamma{0.55}; std::vector boundaries; std::vector loads; + bool nlgeom{true}; }; /// @brief Stores homogeneous isotropic thermal material source data. @@ -374,10 +406,13 @@ struct ModelDefinition { std::vector nodes; std::vector elements; std::vector shell_elements; + std::vector solid_elements; std::vector materials; std::vector thermal_materials; std::vector sections; std::vector shell_sections; + std::vector solid_sections; + std::vector amplitudes; std::vector thermal_sections; std::vector thermal_elements; std::vector shell_node_initial_frames; diff --git a/include/fesa/properties/element_property.h b/include/fesa/properties/element_property.h index ce889aa..0ad40ba 100644 --- a/include/fesa/properties/element_property.h +++ b/include/fesa/properties/element_property.h @@ -6,7 +6,11 @@ namespace fesa { /// @brief Identifies the supported concrete element-property semantics. -enum class ElementPropertyKind { kGeneralBeamSection, kShellSection }; +enum class ElementPropertyKind { + kGeneralBeamSection, + kShellSection, + kSolidSection +}; /// @brief Provides stable identity for a Domain-owned element property. class ElementProperty { diff --git a/include/fesa/properties/solid_section.h b/include/fesa/properties/solid_section.h new file mode 100644 index 0000000..fafe643 --- /dev/null +++ b/include/fesa/properties/solid_section.h @@ -0,0 +1,39 @@ +#ifndef FESA_PROPERTIES_SOLID_SECTION_H_ +#define FESA_PROPERTIES_SOLID_SECTION_H_ + +#include +#include + +#include "fesa/properties/element_property.h" + +namespace fesa { +/// @brief Owns a homogeneous Solid material assignment without thickness. +class SolidSection final : public ElementProperty { + public: + /// @brief Constructs a parser-validated ELSET and material assignment. + SolidSection(SourceEntityId source, std::string name, + EntityIndex material_index, SourceLocation location) + : source_id_{std::move(source)}, + name_{std::move(name)}, + material_index_{material_index}, + location_{std::move(location)} {} + ElementPropertyKind Kind() const noexcept override { + return ElementPropertyKind::kSolidSection; + } + const SourceEntityId& SourceId() const noexcept override { + return source_id_; + } + const SourceLocation& Location() const noexcept override { return location_; } + /// @brief Returns the source ELSET spelling. + const std::string& Name() const noexcept { return name_; } + /// @brief Returns the stable Domain material index. + EntityIndex MaterialIndex() const noexcept { return material_index_; } + + private: + SourceEntityId source_id_; + std::string name_; + EntityIndex material_index_; + SourceLocation location_; +}; +} // namespace fesa +#endif // FESA_PROPERTIES_SOLID_SECTION_H_ diff --git a/src/fesa/assembly/load_assembler.cpp b/src/fesa/assembly/load_assembler.cpp index 6b945d2..7ad0ed6 100644 --- a/src/fesa/assembly/load_assembler.cpp +++ b/src/fesa/assembly/load_assembler.cpp @@ -98,7 +98,8 @@ Status ValidateShellMoments(const Domain& domain, const Vector& full_load) { } // namespace Result LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model, - const DofManager& dofs) { + const DofManager& dofs, + double step_time) { const auto& active_loads = model.ActiveLoads(); const auto& owned_loads = model.Step().Loads(); if (active_loads.size() != owned_loads.size()) { @@ -122,21 +123,31 @@ Result LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model, } loads.push_back(owned_loads[load_index]); } - return AssembleFullNodalLoad(model, dofs, loads); + return AssembleFullNodalLoad(model, dofs, loads, step_time); } Result LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model, const DofManager& dofs, - const LoadView& loads) { + const LoadView& loads, + double step_time) { const Domain& domain = model.GetDomain(); + const std::size_t components = dofs.ComponentsPerNode(); + if (dofs.Profile() == DofProfile::kThermal || + (domain.IsSolid() != (dofs.Profile() == DofProfile::kSolidStructural)) || + !std::isfinite(step_time) || step_time < 0.0) { + return Result::Failure(LoadFailure( + "invalid-load-context", {domain.SourcePath(), 0U}, "LOAD_ASSEMBLER", "", + "Nodal loads require the matching structural profile and finite " + "nonnegative endpoint time.")); + } if (domain.Nodes().size() > - (std::numeric_limits::max)() / kDofsPerNode) { + (std::numeric_limits::max)() / components) { 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 expected_full_count = domain.Nodes().size() * kDofsPerNode; + const std::size_t expected_full_count = domain.Nodes().size() * components; if (dofs.FullDofCount() != expected_full_count) { return Result::Failure(LoadFailure( "invalid-load-dimensions", {domain.SourcePath(), 0U}, "LOAD_ASSEMBLER", @@ -148,11 +159,11 @@ Result LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model, return Result::Failure(dof_status); } for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { - for (std::size_t component = 0U; component < kDofsPerNode; ++component) { + for (std::size_t component = 0U; component < components; ++component) { try { if (dofs.FullDof(static_cast(node), static_cast(component)) != - node * kDofsPerNode + component) { + node * components + component) { return Result::Failure(LoadFailure( "invalid-load-order", domain.Nodes()[node].location, "LOAD_ASSEMBLER", @@ -170,7 +181,7 @@ Result LoadAssembler::AssembleFullNodalLoad(const AnalysisModel& model, const SourceTargetIndex target_index = SourceTargetIndex::FromDomain(domain); const SourceTargetResolver target_resolver{target_index}; - const LoadContext context{domain, dofs, target_resolver}; + const LoadContext context{domain, dofs, target_resolver, step_time}; std::vector> contributions_by_load; contributions_by_load.reserve(loads.size()); diff --git a/src/fesa/fem/dof_manager.cpp b/src/fesa/fem/dof_manager.cpp index f532553..17ac70a 100644 --- a/src/fesa/fem/dof_manager.cpp +++ b/src/fesa/fem/dof_manager.cpp @@ -116,8 +116,12 @@ Result DofManager::Create(const AnalysisModel& model) { "An active element index is outside the Domain.")); } const auto& definition = domain.Elements()[element_index]; - layouts.push_back({definition.SourceId(), definition.NodeIndices(), - FullNodeComponents()}); + layouts.push_back( + {definition.SourceId(), definition.NodeIndices(), + domain.IsSolid() + ? std::vector{DofComponent::kUx, DofComponent::kUy, + DofComponent::kUz} + : FullNodeComponents()}); } DofManager dofs; @@ -305,6 +309,19 @@ Status DofManager::BuildLayouts(const AnalysisModel& analysis_model, ++source_order) { const EntityIndex definition_index = analysis_model.ActiveElements()[source_order]; + if (domain.IsSolid()) { + const std::vector translations{ + DofComponent::kUx, DofComponent::kUy, DofComponent::kUz}; + if (layouts[source_order].components_per_node != translations || + definition_index >= domain.Elements().Size() || + layouts[source_order].node_indices != + domain.Elements()[definition_index].NodeIndices()) { + return DofFailure( + "invalid-solid-element-layout", {domain.SourcePath(), 0U}, + std::to_string(source_order), + "Solid layouts require source-ordered eight nodes and UX/UY/UZ."); + } + } if (definition_index >= domain.Elements().Size() || !SameSourceIdentity(layouts[source_order].source_id, domain.Elements()[definition_index].SourceId())) { @@ -318,10 +335,13 @@ Status DofManager::BuildLayouts(const AnalysisModel& analysis_model, const SourceTargetIndex target_index = SourceTargetIndex::FromDomain(domain); const SourceTargetResolver target_resolver{target_index}; - const std::size_t full_count = domain.Nodes().size() * kDofsPerNode; + const std::size_t component_count = domain.IsSolid() ? 3U : kDofsPerNode; + const std::size_t full_count = domain.Nodes().size() * component_count; DofManager full_dof_map; full_dof_map.full_dof_count_ = full_count; + full_dof_map.profile_ = + domain.IsSolid() ? DofProfile::kSolidStructural : DofProfile::kStructural; const BoundaryConditionContext context{domain, full_dof_map, target_resolver}; std::vector> definitions_by_boundary; definitions_by_boundary.reserve(boundaries.size()); @@ -417,6 +437,7 @@ Status DofManager::BuildLayouts(const AnalysisModel& analysis_model, std::move(constrained_dofs), std::move(prescribed_values), {}}; + candidate.profile_ = full_dof_map.profile_; candidate.element_scatters_.reserve(layouts.size()); for (const auto& layout : layouts) { auto scatter = candidate.ElementScatter(layout); @@ -441,6 +462,10 @@ std::size_t DofManager::FullDofCount() const noexcept { } DofProfile DofManager::Profile() const noexcept { return profile_; } +std::size_t DofManager::ComponentsPerNode() const noexcept { + if (profile_ == DofProfile::kThermal) return kThermalDofsPerNode; + return profile_ == DofProfile::kSolidStructural ? 3U : kDofsPerNode; +} std::size_t DofManager::FreeDofCount() const noexcept { return free_dofs_.size(); @@ -452,15 +477,16 @@ std::size_t DofManager::ConstrainedDofCount() const noexcept { std::size_t DofManager::FullDof(const EntityIndex node, const DofComponent component) const { - if (profile_ != DofProfile::kStructural) { - throw std::logic_error{"Structural components are invalid for thermal DOFs."}; + if (profile_ == DofProfile::kThermal) { + throw std::logic_error{ + "Structural components are invalid for thermal DOFs."}; } const std::size_t component_index = static_cast(component); - if (node >= full_dof_count_ / kDofsPerNode || - component_index >= kDofsPerNode) { + if (node >= full_dof_count_ / ComponentsPerNode() || + component_index >= ComponentsPerNode()) { throw std::out_of_range{"Node or DOF component is out of range."}; } - return static_cast(node) * kDofsPerNode + component_index; + return static_cast(node) * ComponentsPerNode() + component_index; } std::size_t DofManager::FullThermalDof( @@ -481,7 +507,7 @@ std::optional DofManager::FreeEquation( Result> DofManager::ElementScatter( const ElementDofLayout& layout) const { - if (profile_ != DofProfile::kStructural) { + if (profile_ == DofProfile::kThermal) { return Result>::Failure(DofFailure( "invalid-dof-profile", {}, layout.source_id.source_label_text, "Structural layouts are invalid for a thermal mapping.")); @@ -501,8 +527,8 @@ Result> DofManager::ElementScatter( for (const EntityIndex node : layout.node_indices) { for (const DofComponent component : layout.components_per_node) { const std::size_t component_index = static_cast(component); - if (node >= full_dof_count_ / kDofsPerNode || - component_index >= kDofsPerNode) { + if (node >= full_dof_count_ / ComponentsPerNode() || + component_index >= ComponentsPerNode()) { return Result>::Failure(DofFailure( "invalid-element-dof-layout", {}, layout.source_id.source_label_text, @@ -510,7 +536,8 @@ Result> DofManager::ElementScatter( "DOF range.")); } const std::size_t full_dof = - static_cast(node) * kDofsPerNode + component_index; + static_cast(node) * ComponentsPerNode() + + component_index; if (std::find(scatter.begin(), scatter.end(), full_dof) != scatter.end()) { return Result>::Failure(DofFailure( diff --git a/src/fesa/io/abaqus/domain_mapping_context.cpp b/src/fesa/io/abaqus/domain_mapping_context.cpp index 9a6c235..6784497 100644 --- a/src/fesa/io/abaqus/domain_mapping_context.cpp +++ b/src/fesa/io/abaqus/domain_mapping_context.cpp @@ -50,7 +50,7 @@ struct RawNode { }; struct RawElement { - enum class Type { kB33, kS4, kS4r, kDs4 }; + enum class Type { kB33, kS4, kS4r, kDs4, kC3d8 }; std::int64_t label; std::string label_text; @@ -59,7 +59,7 @@ struct RawElement { SourceLocation location; }; -enum class ElementFamily { kBeam, kShell, kThermalShell }; +enum class ElementFamily { kBeam, kShell, kThermalShell, kSolid }; struct RawSet { std::string name; @@ -83,6 +83,12 @@ struct RawShellSection { SourceLocation location; }; +struct RawSolidSection { + std::string element_set_name; + std::string material_name; + SourceLocation location; +}; + struct RawPart { std::string name; SourceLocation location; @@ -92,6 +98,7 @@ struct RawPart { std::vector element_sets; std::vector sections; std::vector shell_sections; + std::vector solid_sections; }; struct RawInstance { @@ -170,6 +177,8 @@ struct RawStep { std::string name_raw; std::string name_canonical; bool has_static{false}; + bool static_direct{false}; + std::uint32_t static_increment_count{0U}; std::optional dynamic; bool perturbation_explicit{false}; bool nlgeom_explicit{false}; @@ -746,6 +755,10 @@ class DomainMappingContext::Impl { ParseAssemblyStart(block); return; } + if (block.canonical_name == "AMPLITUDE") { + ParseAmplitude(block); + return; + } if (block.canonical_name == "MATERIAL") { ParseMaterial(block); return; @@ -841,7 +854,8 @@ class DomainMappingContext::Impl { model_element_family_ == ElementFamily::kShell || model_element_family_ == ElementFamily::kThermalShell; if (part.nodes.empty() || part.elements.empty() || - (!shell_part && part.sections.empty())) { + (!shell_part && model_element_family_ != ElementFamily::kSolid && + part.sections.empty())) { InvalidKeywordLocation(block, "A part closes only after node, element, " "and matching section blocks."); @@ -917,7 +931,8 @@ class DomainMappingContext::Impl { return; } if (model_element_family_ == ElementFamily::kShell || - model_element_family_ == ElementFamily::kThermalShell) { + model_element_family_ == ElementFamily::kThermalShell || + model_element_family_ == ElementFamily::kSolid) { InputFailure( "unsupported-mixed-element-model", block.location, block.canonical_name, "", @@ -934,7 +949,8 @@ class DomainMappingContext::Impl { block, "SHELL SECTION follows the part mesh and optional sets."); return; } - if (model_element_family_ == ElementFamily::kBeam) { + if (model_element_family_ == ElementFamily::kBeam || + model_element_family_ == ElementFamily::kSolid) { InputFailure( "unsupported-mixed-element-model", block.location, block.canonical_name, "", @@ -943,6 +959,17 @@ class DomainMappingContext::Impl { } ParseShellSection(block, part); part_sections_seen_ = !failure_; + } else if (block.canonical_name == "SOLID SECTION") { + beam_section_context_active_ = false; + if (!part_elements_seen_ || + model_element_family_ != ElementFamily::kSolid) { + InputFailure("invalid-solid-section", block.location, + block.canonical_name, "", + "SOLID SECTION requires an exclusive C3D8 mesh."); + return; + } + ParseSolidSection(block, part); + part_sections_seen_ = !failure_; } else if (block.canonical_name == "SECTION POINTS") { ParseSectionPoints(block, part); } else if (block.canonical_name == "TRANSVERSE SHEAR STIFFNESS") { @@ -1011,10 +1038,12 @@ class DomainMappingContext::Impl { const bool thermal_block = type_parameter && type_parameter->value && AsciiCaseInsensitiveEquals(*type_parameter->value, "DS4"); - if (!ValidateParameters(block, - thermal_block - ? std::vector{"TYPE", "ELSET"} - : std::vector{"TYPE"})) { + if (!ValidateParameters( + block, (thermal_block || (type_parameter && type_parameter->value && + AsciiCaseInsensitiveEquals( + *type_parameter->value, "C3D8"))) + ? std::vector{"TYPE", "ELSET"} + : std::vector{"TYPE"})) { return; } const auto* type = RequiredParameterValue(block, "TYPE"); @@ -1036,6 +1065,10 @@ class DomainMappingContext::Impl { element_type = RawElement::Type::kS4r; element_family = ElementFamily::kShell; expected_field_count = 5U; + } else if (AsciiCaseInsensitiveEquals(*type, "C3D8")) { + element_type = RawElement::Type::kC3d8; + element_family = ElementFamily::kSolid; + expected_field_count = 9U; } else if (AsciiCaseInsensitiveEquals(*type, "DS4")) { element_type = RawElement::Type::kDs4; element_family = ElementFamily::kThermalShell; @@ -1043,15 +1076,16 @@ class DomainMappingContext::Impl { } else { InputFailure("unsupported-element-formulation", block.location, block.canonical_name, *type, - "Only TYPE=B33, TYPE=S4, TYPE=S4R, and TYPE=DS4 belong to " + "Only B33, S4, S4R, DS4, and C3D8 belong to " "the approved " "element subsets."); return; } if (model_element_family_ && *model_element_family_ != element_family) { - InputFailure("unsupported-mixed-element-model", block.location, - block.canonical_name, *type, - "Beam and shell elements cannot be mixed in one model."); + InputFailure( + "unsupported-mixed-element-model", block.location, + block.canonical_name, *type, + "Beam, Shell, thermal, and Solid families cannot mix in one model."); return; } model_element_family_ = element_family; @@ -1067,9 +1101,11 @@ class DomainMappingContext::Impl { inline_set = RawSet{*name, {}, block.location}; } if (block.data.empty()) { - InputFailure(element_family == ElementFamily::kShell - ? "invalid-shell-connectivity" - : "invalid-data-arity", + InputFailure(element_family == ElementFamily::kSolid + ? "invalid-solid-connectivity" + : (element_family == ElementFamily::kShell + ? "invalid-shell-connectivity" + : "invalid-data-arity"), block.location, block.canonical_name, "", "ELEMENT requires at least one row with the approved " "connectivity arity."); @@ -1078,9 +1114,11 @@ class DomainMappingContext::Impl { for (const auto& row : block.data) { if (row.fields.size() != expected_field_count) { InputFailure( - element_family == ElementFamily::kShell - ? "invalid-shell-connectivity" - : "invalid-data-arity", + element_family == ElementFamily::kSolid + ? "invalid-solid-connectivity" + : (element_family == ElementFamily::kShell + ? "invalid-shell-connectivity" + : "invalid-data-arity"), row.location, block.canonical_name, "", element_family == ElementFamily::kShell ? "S4 and S4R rows require a label and exactly four nodes." @@ -1108,10 +1146,12 @@ class DomainMappingContext::Impl { const std::set distinct_nodes{element.node_labels.begin(), element.node_labels.end()}; if (distinct_nodes.size() != element.node_labels.size()) { - InputFailure( - "invalid-shell-connectivity", row.location, block.canonical_name, - element.label_text, - "Shell connectivity requires four distinct source nodes."); + InputFailure(element_family == ElementFamily::kSolid + ? "invalid-solid-connectivity" + : "invalid-shell-connectivity", + row.location, block.canonical_name, element.label_text, + "Connectivity requires the specified number of distinct " + "source nodes."); return; } } @@ -1132,6 +1172,101 @@ class DomainMappingContext::Impl { if (inline_set) part.element_sets.push_back(std::move(*inline_set)); } + /// @brief Parses material-only Solid assignment without thickness defaults. + void ParseSolidSection(const KeywordBlock& block, RawPart& part) { + if (!ValidateParameters(block, {"ELSET", "MATERIAL"})) return; + const auto* elset = RequiredParameterValue(block, "ELSET"); + const auto* material = RequiredParameterValue(block, "MATERIAL"); + if (!elset || !material) return; + if (block.data.size() > 1U || + (!block.data.empty() && + std::any_of(block.data[0].fields.begin(), block.data[0].fields.end(), + [](const auto& field) { return !field.empty(); }))) { + InputFailure("invalid-solid-section", block.location, + block.canonical_name, *elset, + "SOLID SECTION permits omitted or blank data only."); + return; + } + part.solid_sections.push_back({*elset, *material, block.location}); + } + + /// @brief Validates and owns one model-level supported load amplitude. + void ParseAmplitude(const KeywordBlock& block) { + const bool dynamic_source = std::any_of( + input_.blocks.begin(), input_.blocks.end(), [](const auto& candidate) { + return candidate.canonical_name == "DYNAMIC"; + }); + if (model_element_family_ != ElementFamily::kSolid || !dynamic_source) { + InputFailure("unsupported-load-amplitude", block.location, + block.canonical_name, "", + "Amplitude is supported only in Solid dynamic models."); + return; + } + if (!ValidateParameters(block, {"NAME", "DEFINITION", "TIME", "VALUE"})) + return; + const auto* name = RequiredParameterValue(block, "NAME"); + if (!name) return; + for (const auto& option : + std::array, 3>{ + {{"DEFINITION", "TABULAR"}, + {"TIME", "STEP TIME"}, + {"VALUE", "RELATIVE"}}}) { + const auto* parameter = Parameter(block, option.first); + if (parameter && + (!parameter->value || + !AsciiCaseInsensitiveEquals(*parameter->value, option.second))) { + InputFailure("unsupported-amplitude-option", block.location, + block.canonical_name, option.first, + "Only TABULAR STEP TIME RELATIVE amplitude is supported."); + return; + } + } + const std::string canonical = CanonicalizeAscii(*name); + if (std::any_of(definition_.amplitudes.begin(), + definition_.amplitudes.end(), [&canonical](const auto& a) { + return a.name_canonical == canonical; + })) { + InputFailure("duplicate-amplitude", block.location, block.canonical_name, + *name, + "Amplitude names must be unique ignoring ASCII case."); + return; + } + TabularAmplitudeDefinition amplitude{*name, canonical, block.location, {}}; + for (const auto& row : block.data) { + if (row.fields.empty() || row.fields.size() % 2U != 0U) { + InputFailure("invalid-amplitude-data", row.location, + block.canonical_name, *name, + "Amplitude rows require complete time/value pairs."); + return; + } + for (std::size_t field = 0U; field < row.fields.size(); field += 2U) { + double time = 0.0; + double value = 0.0; + if (!ParseDouble(row.fields[field], time, row.location, + block.canonical_name) || + !ParseDouble(row.fields[field + 1U], value, row.location, + block.canonical_name)) + return; + if (time < 0.0 || (!amplitude.points.empty() && + !(time > amplitude.points.back()[0]))) { + InputFailure( + "invalid-amplitude-data", row.location, block.canonical_name, + *name, + "Amplitude times must be nonnegative and strictly increasing."); + return; + } + amplitude.points.push_back({time, value}); + } + } + if (amplitude.points.empty()) { + InputFailure("invalid-amplitude-data", block.location, + block.canonical_name, *name, + "Amplitude requires at least one finite time/value pair."); + return; + } + definition_.amplitudes.push_back(std::move(amplitude)); + } + void ParseShellSection(const KeywordBlock& block, RawPart& part) { for (const auto& candidate : block.parameters) { if (candidate.name != "ELSET" && candidate.name != "MATERIAL") { @@ -1735,13 +1870,17 @@ class DomainMappingContext::Impl { } if (!ParseModelDouble(block.data[0].fields[0], material.youngs_modulus, block.data[0].location, block.canonical_name, - model_element_family_ == ElementFamily::kShell + model_element_family_ == ElementFamily::kSolid + ? "invalid-solid-material" + : model_element_family_ == ElementFamily::kShell ? "invalid-shell-material" : "invalid-beam-property", "Elastic material values must be finite.") || !ParseModelDouble(block.data[0].fields[1], material.poisson_ratio, block.data[0].location, block.canonical_name, - model_element_family_ == ElementFamily::kShell + model_element_family_ == ElementFamily::kSolid + ? "invalid-solid-material" + : model_element_family_ == ElementFamily::kShell ? "invalid-shell-material" : "invalid-beam-property", "Elastic material values must be finite.")) { @@ -2006,6 +2145,8 @@ class DomainMappingContext::Impl { } double value = 0.0; if (row.fields.size() == 4U && + !(model_element_family_ == ElementFamily::kSolid && + row.fields[3].empty()) && !ParseDouble(row.fields[3], value, row.location, block.canonical_name)) { return; @@ -2016,19 +2157,46 @@ class DomainMappingContext::Impl { } void ParseCload(const KeywordBlock& block, std::vector& loads) { - if (step_.dynamic && !block.parameters.empty()) { + if (step_.dynamic && model_element_family_ != ElementFamily::kSolid && + !block.parameters.empty()) { InputFailure("unsupported-dynamic-load", block.location, block.canonical_name, "", "Dynamic CLOAD does not support keyword parameters."); return; } - if (!ValidateParameters(block, {}) || block.data.empty()) { + const bool solid_dynamic = + step_.dynamic && model_element_family_ == ElementFamily::kSolid; + if (!ValidateParameters( + block, solid_dynamic ? std::vector{"AMPLITUDE"} + : std::vector{}) || + block.data.empty()) { if (!failure_) { InputFailure("invalid-data-arity", block.location, block.canonical_name, "", "CLOAD requires one or more rows."); } return; } + std::optional amplitude_index; + std::string amplitude_name; + if (Parameter(block, "AMPLITUDE")) { + const auto* name = RequiredParameterValue(block, "AMPLITUDE"); + if (!name) return; + amplitude_name = *name; + for (std::size_t index = 0U; index < definition_.amplitudes.size(); + ++index) { + if (AsciiCaseInsensitiveEquals(definition_.amplitudes[index].name_raw, + *name)) { + amplitude_index = static_cast(index); + break; + } + } + if (!amplitude_index) { + InputFailure("unresolved-amplitude", block.location, + block.canonical_name, *name, + "CLOAD amplitude must resolve to a model-level curve."); + return; + } + } for (const auto& row : block.data) { if (row.fields.size() != 3U || row.fields[0].empty()) { InputFailure("invalid-data-arity", row.location, block.canonical_name, @@ -2048,15 +2216,17 @@ class DomainMappingContext::Impl { row.fields[0], "CLOAD DOF must be in 1..6."); return; } - if (step_.dynamic && dof > 3) { + if ((step_.dynamic || model_element_family_ == ElementFamily::kSolid) && + dof > 3) { InputFailure("unsupported-dynamic-load", row.location, block.canonical_name, row.fields[1], "Dynamic CLOAD supports global translational DOF 1 " "through 3 only."); return; } - loads.push_back( - {row.fields[0], static_cast(dof), magnitude, row.location}); + loads.push_back({row.fields[0], static_cast(dof), magnitude, + row.location, amplitude_index, amplitude_name, + block.location}); } } @@ -2148,10 +2318,11 @@ class DomainMappingContext::Impl { return; } if (AsciiCaseInsensitiveEquals(*nlgeom_parameter->value, "YES")) { - if (model_element_family_ != ElementFamily::kShell) { + if (model_element_family_ != ElementFamily::kShell && + model_element_family_ != ElementFamily::kSolid) { InputFailure("unsupported-analysis-procedure", block.location, block.canonical_name, *nlgeom_parameter->value, - "NLGEOM=YES requires an S4 or S4R shell model."); + "NLGEOM=YES requires a supported Shell or Solid model."); return; } procedure_kind = StaticProcedureKind::kNonlinearFullNewton; @@ -2353,7 +2524,8 @@ class DomainMappingContext::Impl { return; } step_no_op_seen_ = true; - if (step_.frequency || step_.has_heat_transfer) { + if (step_.frequency || step_.has_heat_transfer || + (step_.dynamic && model_element_family_ == ElementFamily::kSolid)) { const auto* write = Parameter(block, "WRITE"); const auto* frequency = Parameter(block, "FREQUENCY"); if (!ValidateParameters(block, {"WRITE", "FREQUENCY"}) || @@ -2553,18 +2725,26 @@ class DomainMappingContext::Impl { } void ParseDynamic(const KeywordBlock& block) { - if (step_.perturbation_explicit || !step_.nlgeom || - model_element_family_ != ElementFamily::kShell || + if (step_.perturbation_explicit || + (model_element_family_ != ElementFamily::kSolid && + (!step_.nlgeom || model_element_family_ != ElementFamily::kShell)) || !ValidateParameters(block, {"DIRECT"}, "unsupported-dynamic-parameter")) { if (!failure_) { - InputFailure( - "unsupported-analysis-procedure", block.location, - block.canonical_name, "", - "DYNAMIC requires the approved nonlinear shell configuration."); + InputFailure("unsupported-analysis-procedure", block.location, + block.canonical_name, "", + "DYNAMIC requires nonlinear Shell or linear/nonlinear " + "Solid geometry."); } return; } + if (model_element_family_ == ElementFamily::kSolid && + block.data.size() == 1U && block.data[0].fields.size() > 3U) { + InputFailure("invalid-dynamic-data", block.data[0].location, + block.canonical_name, "", + "Solid DYNAMIC requires h,T with at most one blank field."); + return; + } const auto* direct = Parameter(block, "DIRECT"); if (direct == nullptr || direct->value || block.data.size() != 1U || block.data[0].fields.size() < 2U || @@ -2622,6 +2802,7 @@ class DomainMappingContext::Impl { "DYNAMIC,DIRECT increment count cannot exceed STEP INC."); return; } + dynamic.nlgeom = step_.nlgeom; step_.dynamic = std::move(dynamic); } @@ -2716,6 +2897,46 @@ class DomainMappingContext::Impl { } void ParseStatic(const KeywordBlock& block) { + if (model_element_family_ == ElementFamily::kSolid && + Parameter(block, "DIRECT")) { + if (step_.perturbation_explicit || step_.has_static || step_.frequency || + !ValidateParameters(block, {"DIRECT"})) + return; + const auto* direct = Parameter(block, "DIRECT"); + if (direct->value || block.data.size() != 1U || + (block.data[0].fields.size() < 2U || + block.data[0].fields.size() > 3U || + std::any_of(block.data[0].fields.begin() + 2, + block.data[0].fields.end(), + [](const auto& field) { return !field.empty(); }))) { + InputFailure("invalid-static-data", block.location, + block.canonical_name, "DIRECT", + "STATIC DIRECT requires exactly h,T."); + return; + } + double h = 0.0; + double period = 0.0; + if (!ParseDouble(block.data[0].fields[0], h, block.data[0].location, + "STATIC") || + !ParseDouble(block.data[0].fields[1], period, block.data[0].location, + "STATIC")) + return; + const double count = std::round(period / h); + if (!(h > 0.0) || !(period > 0.0) || !std::isfinite(count) || + count < 1.0 || + count > static_cast(step_.maximum_increments) || + h * count != period) { + InputFailure("invalid-static-schedule", block.data[0].location, + block.canonical_name, "DIRECT", + "STATIC DIRECT requires positive exact N*h=T and N<=INC."); + return; + } + step_.static_values = {h, period, h, h}; + step_.static_direct = true; + step_.static_increment_count = static_cast(count); + step_.has_static = true; + return; + } if (step_.perturbation_explicit || step_.frequency) { InputFailure("unsupported-analysis-procedure", block.location, block.canonical_name, "PERTURBATION", @@ -2811,6 +3032,17 @@ class DomainMappingContext::Impl { "OUTPUT must select exactly FIELD or HISTORY."); return; } + if (model_element_family_ == ElementFamily::kSolid && + ((variable && + (!variable->value || + !AsciiCaseInsensitiveEquals(*variable->value, "PRESELECT"))) || + !RequireNoData(block))) { + if (!failure_) + InputFailure("unsupported-keyword", block.location, + block.canonical_name, "VARIABLE", + "Solid OUTPUT accepts only PRESELECT requests."); + return; + } if (step_.nlgeom && (variable == nullptr || !variable->value || !AsciiCaseInsensitiveEquals(*variable->value, "PRESELECT"))) { @@ -2865,6 +3097,50 @@ class DomainMappingContext::Impl { "Output variable keywords require an active OUTPUT request."); return; } + if (model_element_family_ == ElementFamily::kSolid && step_.dynamic) { + InputFailure( + "unsupported-keyword", block.location, block.canonical_name, "", + "Solid dynamic output children are outside the approved subset."); + return; + } + if (model_element_family_ == ElementFamily::kSolid && step_.has_static) { + const auto* variable = Parameter(block, "VARIABLE"); + const auto* directions = Parameter(block, "DIRECTIONS"); + const bool node = block.canonical_name == "NODE OUTPUT"; + const bool element = block.canonical_name == "ELEMENT OUTPUT"; + const bool contact = block.canonical_name == "CONTACT OUTPUT"; + if (!active_output_is_field_ || (!node && !element && !contact) || + !ValidateParameters( + block, + element ? std::vector{"VARIABLE", "DIRECTIONS"} + : std::vector{"VARIABLE"}) || + (directions && + (!directions->value || + !AsciiCaseInsensitiveEquals(*directions->value, "YES"))) || + (variable && + (!variable->value || + !AsciiCaseInsensitiveEquals(*variable->value, "PRESELECT"))) || + (contact && (!variable || !block.data.empty()))) { + if (!failure_) + InputFailure("unsupported-keyword", block.location, + block.canonical_name, "", + "Solid static output child is unsupported."); + return; + } + for (const auto& row : block.data) { + for (const auto& field : row.fields) { + const std::string token = CanonicalizeAscii(field); + if (!(node ? (token == "U" || token == "RF") + : (token == "S" || token == "E"))) { + InputFailure("unsupported-keyword", row.location, + block.canonical_name, field, + "Solid output quantities must be physical supported " + "components."); + return; + } + } + } + } if (step_.frequency) { const std::string_view required_variable = block.canonical_name == "NODE OUTPUT" ? "PRESELECT" : "ALL"; @@ -3001,16 +3277,32 @@ class DomainMappingContext::Impl { "A material must own exactly one ELASTIC row."); return; } - if (model_element_family_ == ElementFamily::kShell) { + if (model_element_family_ == ElementFamily::kShell || + model_element_family_ == ElementFamily::kSolid) { if (!(material.youngs_modulus > 0.0) || !(material.poisson_ratio > -1.0) || !(material.poisson_ratio < 0.5)) { - ModelFailure( - "invalid-shell-material", material.elastic_location, "ELASTIC", - material.name, - "Shell isotropic elasticity requires E>0 and -10 and -1(index), + *materials_[index].density, + materials_[index].density_location}); + } + return; + } if (!step_.frequency && !step_.dynamic) { if (std::any_of(materials_.begin(), materials_.end(), [](const RawMaterial& material) { @@ -3090,10 +3392,11 @@ class DomainMappingContext::Impl { } } } - if (model_element_family_ != ElementFamily::kShell) { + if (model_element_family_ != ElementFamily::kShell && + model_element_family_ != ElementFamily::kSolid) { InputFailure("unsupported-analysis-procedure", step_.location, step_.frequency ? "FREQUENCY" : "DYNAMIC", step_.name_raw, - "The procedure requires the approved shell model family."); + "The procedure requires an approved Shell or Solid model."); return; } if (!element_surfaces_.empty()) { @@ -3112,6 +3415,10 @@ class DomainMappingContext::Impl { break; } } + for (const auto& section : definition_.solid_sections) { + if (section.MaterialIndex() == static_cast(index)) + active = true; + } if (active) { ModelFailure( "missing-material-density", material.location, "MATERIAL", @@ -3196,7 +3503,8 @@ class DomainMappingContext::Impl { ++part_index) { const RawPart& part = parts_[part_index]; if (part.nodes.empty() || part.elements.empty() || - (!shell_model && part.sections.empty())) { + (!shell_model && model_element_family_ != ElementFamily::kSolid && + part.sections.empty())) { InputFailure("invalid-model-cardinality", part.location, "PART", part.name, "A part requires nodes, elements, and its approved " @@ -3245,6 +3553,45 @@ class DomainMappingContext::Impl { } definition_.parts.push_back(std::move(part_definition)); + if (model_element_family_ == ElementFamily::kSolid) { + auto& assignments = part_shell_section_assignments_[part_index]; + for (const auto& section : part.solid_sections) { + const auto* element_set = + FindSet(part.element_sets, section.element_set_name); + const auto material = FindMaterialIndex(section.material_name); + if (!element_set || !material) { + InputFailure("unresolved-solid-section", section.location, + "SOLID SECTION", section.element_set_name, + "Solid ELSET and MATERIAL must resolve."); + return; + } + const auto index = + static_cast(definition_.solid_sections.size()); + definition_.solid_sections.emplace_back( + SourceEntityId{part.name, static_cast(index), + section.element_set_name}, + section.element_set_name, *material, section.location); + for (const auto label : element_set->members) { + if (!assignments.emplace(label, std::make_pair(index, *material)) + .second) { + InputFailure( + "invalid-solid-section-assignment", section.location, + "SOLID SECTION", std::to_string(label), + "Solid elements require exactly one section assignment."); + return; + } + } + } + for (const auto& element : part.elements) { + if (assignments.find(element.label) == assignments.end()) { + InputFailure("invalid-solid-section-assignment", element.location, + "ELEMENT", element.label_text, + "A Solid section assignment is missing."); + return; + } + } + continue; + } if (shell_model) { auto& assignments = part_shell_section_assignments_[part_index]; for (const auto& section : part.shell_sections) { @@ -3475,6 +3822,36 @@ class DomainMappingContext::Impl { assignment->second.second, assignment->second.first, raw_element.location}); + } else if (raw_element.type == RawElement::Type::kC3d8) { + std::array connected_nodes{}; + for (std::size_t node = 0U; node < connected_nodes.size(); ++node) { + const auto found = node_indices.find(raw_element.node_labels[node]); + if (found == node_indices.end()) { + InputFailure("invalid-solid-connectivity", raw_element.location, + "ELEMENT", raw_element.label_text, + "Solid nodes must resolve."); + return; + } + connected_nodes[node] = found->second; + } + const auto& assignments = + part_shell_section_assignments_[PartIndex(*part)]; + const auto assignment = assignments.find(raw_element.label); + if (assignment == assignments.end() || + definition_.solid_elements.size() > + std::numeric_limits::max()) { + InputFailure("invalid-solid-section-assignment", + raw_element.location, "ELEMENT", + raw_element.label_text, + "Solid assignment must resolve."); + return; + } + index = static_cast(definition_.solid_elements.size()); + definition_.solid_elements.emplace_back( + SourceEntityId{raw_instance.name, raw_element.label, + raw_element.label_text}, + connected_nodes, assignment->second.second, + assignment->second.first, raw_element.location); } else { std::array connected_nodes{}; for (std::size_t node = 0U; node < connected_nodes.size(); ++node) { @@ -4048,6 +4425,47 @@ class DomainMappingContext::Impl { const SourceTargetIndex target_index = BuildSourceTargetIndex(); const SourceTargetResolver target_resolver{target_index}; + if (model_element_family_ == ElementFamily::kSolid) { + std::vector physical; + for (auto boundary : boundaries) { + if (!ResolveNodeTarget(target_resolver, boundary.target, + boundary.location, "BOUNDARY")) + return; + if (boundary.last_dof > 3) { + if (boundary.value != 0.0) { + InputFailure("unsupported-solid-rotation-boundary", + boundary.location, "BOUNDARY", boundary.target, + "Solid rotational boundaries must be exact zero to be " + "inactive."); + return; + } + definition_.warnings.push_back( + {Severity::kWarning, "inactive-solid-rotation-boundary", + boundary.location, "BOUNDARY", boundary.target, + "Exact-zero rotational boundary DOFs " + + std::to_string(std::max(4, boundary.first_dof)) + ".." + + std::to_string(boundary.last_dof) + + " are inactive; Solid equations contain translations " + "only."}); + boundary.last_dof = 3; + } + if (boundary.first_dof <= boundary.last_dof) + physical.push_back(std::move(boundary)); + } + boundaries = std::move(physical); + std::vector active(definition_.nodes.size(), false); + for (const auto& element : definition_.solid_elements) + for (const auto node : element.node_indices) active[node] = true; + for (std::size_t node = 0U; node < active.size(); ++node) { + if (!active[node]) { + ModelFailure("unused-solid-node", definition_.nodes[node].location, + "NODE", + definition_.nodes[node].source_id.source_label_text, + "Every Solid node must connect to a C3D8 element."); + return; + } + } + } std::map, double> prescribed_values; for (const auto& boundary : boundaries) { auto target = ResolveNodeTarget(target_resolver, boundary.target, @@ -4115,10 +4533,15 @@ class DomainMappingContext::Impl { // Static time-control values are provenance only: V0 creates exactly // the canonical final frame 0 and performs no increment loop here. definition_.steps.push_back( - {"Step-1", std::move(boundaries), step_.loads, - std::move(surface_pressure_loads), step_.procedure_kind, step_.nlgeom, - step_.static_values[0], step_.static_values[1], step_.static_values[2], - step_.static_values[3], step_.location}); + {model_element_family_ == ElementFamily::kSolid ? step_.name_raw + : "Step-1", + std::move(boundaries), step_.loads, std::move(surface_pressure_loads), + step_.procedure_kind, step_.nlgeom, step_.static_values[0], + step_.static_values[1], step_.static_values[2], step_.static_values[3], + step_.location}); + definition_.steps.back().direct = step_.static_direct; + definition_.steps.back().increment_count = step_.static_increment_count; + definition_.steps.back().maximum_increments = step_.maximum_increments; } const ParsedInput& input_; diff --git a/src/fesa/loads/concentrated_nodal_load.cpp b/src/fesa/loads/concentrated_nodal_load.cpp index 9061895..3f3dd05 100644 --- a/src/fesa/loads/concentrated_nodal_load.cpp +++ b/src/fesa/loads/concentrated_nodal_load.cpp @@ -30,15 +30,18 @@ ConcentratedNodalLoad::ConcentratedNodalLoad( global_components_{global_components}, source_order_{source_order} {} -ConcentratedNodalLoad::ConcentratedNodalLoad(SourceTargetQuery target, - const int source_dof, - const double magnitude, - const std::size_t source_order, - SourceLocation location) +ConcentratedNodalLoad::ConcentratedNodalLoad( + SourceTargetQuery target, const int source_dof, const double magnitude, + const std::size_t source_order, SourceLocation location, + std::optional amplitude_index, std::string amplitude_name_raw, + SourceLocation amplitude_location) : target_{std::move(target)}, source_order_{source_order}, location_{std::move(location)}, - source_dof_{source_dof} { + source_dof_{source_dof}, + amplitude_index_{amplitude_index}, + amplitude_name_raw_{std::move(amplitude_name_raw)}, + amplitude_location_{std::move(amplitude_location)} { if (source_dof >= 1 && source_dof <= static_cast(kDofsPerNode)) { global_components_[static_cast(source_dof - 1)] = magnitude; } @@ -47,6 +50,33 @@ ConcentratedNodalLoad::ConcentratedNodalLoad(SourceTargetQuery target, Result> ConcentratedNodalLoad::ComputeContributions(const LoadContext& context) const { const std::string& identity = target_.target_name_or_label; + const std::size_t components = context.dof_manager.ComponentsPerNode(); + if (context.dof_manager.Profile() == DofProfile::kThermal || + (source_dof_ != 0 && source_dof_ > static_cast(components))) { + return Result>::Failure(LoadFailure( + "invalid-load-dof", location_, identity, + "The source load component is inactive for this physical profile.")); + } + for (std::size_t component = components; component < kDofsPerNode; + ++component) { + if (global_components_[component] != 0.0) { + return Result>::Failure( + LoadFailure("invalid-load-dof", location_, identity, + "A Solid load cannot contain rotational components.")); + } + } + double amplitude = 1.0; + if (amplitude_index_) { + if (!context.domain.IsSolid() || context.domain.Dynamic() == nullptr || + *amplitude_index_ >= context.domain.Amplitudes().size() || + !std::isfinite(context.step_time) || context.step_time < 0.0) { + return Result>::Failure( + LoadFailure("invalid-load-amplitude", location_, identity, + "Amplitude requires a valid Solid dynamic endpoint.")); + } + amplitude = context.domain.Amplitudes()[*amplitude_index_].Evaluate( + context.step_time); + } if (source_dof_ != 0 && (source_dof_ < 1 || source_dof_ > static_cast(kDofsPerNode))) { return Result>::Failure(LoadFailure( @@ -76,7 +106,7 @@ ConcentratedNodalLoad::ComputeContributions(const LoadContext& context) const { std::vector seen(context.domain.Nodes().size(), 0U); std::vector contributions; - contributions.reserve(resolved.Value().size() * kDofsPerNode); + contributions.reserve(resolved.Value().size() * components); for (const auto& target : resolved.Value()) { const EntityIndex node = target.entity_index; if (node >= context.domain.Nodes().size() || seen[node] != 0U) { @@ -86,13 +116,19 @@ ConcentratedNodalLoad::ComputeContributions(const LoadContext& context) const { "identities.")); } seen[node] = 1U; - for (std::size_t component = 0U; component < kDofsPerNode; ++component) { + for (std::size_t component = 0U; component < components; ++component) { + const double value = global_components_[component] * amplitude; + if (!std::isfinite(value)) { + return Result>::Failure( + LoadFailure("nonfinite-load-value", location_, identity, + "Sampled load magnitude must remain finite.")); + } try { contributions.push_back( {source_order_, context.dof_manager.FullDof(node, static_cast(component)), - global_components_[component]}); + value}); } catch (const std::out_of_range&) { return Result>::Failure(LoadFailure( "invalid-load-target", location_, identity, diff --git a/src/fesa/model/domain.cpp b/src/fesa/model/domain.cpp index 5989890..20e49d4 100644 --- a/src/fesa/model/domain.cpp +++ b/src/fesa/model/domain.cpp @@ -1,10 +1,13 @@ #include "fesa/model/domain.h" -#include +#include +#include #include +#include #include #include #include +#include #include #include @@ -48,12 +51,14 @@ Status ValidateThermalTargets( Status ValidateThermalDefinition(const ModelDefinition& definition) { const bool has_structural = !definition.elements.empty() || - !definition.shell_elements.empty(); - const bool has_thermal_data = !definition.thermal_materials.empty() || - !definition.thermal_sections.empty() || - definition.thermal_step.has_value() || - !definition.thermal_initial_conditions.empty() || - !definition.thermal_model_boundaries.empty(); + !definition.shell_elements.empty() || + !definition.solid_elements.empty(); + const bool has_thermal_data = + !definition.thermal_materials.empty() || + !definition.thermal_sections.empty() || + definition.thermal_step.has_value() || + !definition.thermal_initial_conditions.empty() || + !definition.thermal_model_boundaries.empty(); if (has_structural && !definition.thermal_elements.empty()) { return InvalidThermalDefinition({}, "DS4", "Structural and DS4 elements cannot mix."); @@ -161,6 +166,180 @@ Status ValidateThermalDefinition(const ModelDefinition& definition) { return Status::Ok(); } +/// @brief Checks Solid semantic candidates before immutable ownership starts. +Status ValidateSolidDefinition(ModelDefinition& definition) { + const auto fail = [](const SourceLocation& location, + const std::string& identity, + const std::string& message) { + return Status::Failure(FailureCategory::kModel, + {{Severity::kError, "invalid-solid-definition", + location, "SOLID", identity, message}}); + }; + if (definition.solid_elements.empty()) { + if (!definition.solid_sections.empty() || !definition.amplitudes.empty()) + return fail({}, "C3D8", "Solid definitions require a C3D8 model."); + return Status::Ok(); + } + if (!definition.elements.empty() || !definition.shell_elements.empty() || + !definition.thermal_elements.empty() || definition.thermal_step || + !definition.sections.empty() || !definition.shell_sections.empty() || + !definition.element_surfaces.empty()) { + return fail({}, "C3D8", + "Solid models cannot mix element or section families."); + } + if (definition.steps.size() + (definition.frequency_step ? 1U : 0U) + + (definition.dynamic_step ? 1U : 0U) != + 1U) { + return fail({}, "STEP", "Solid models require exactly one procedure."); + } + std::vector active(definition.nodes.size(), false); + std::set> node_ids; + for (const auto& node : definition.nodes) { + if (!node_ids + .emplace(node.source_id.instance_name, node.source_id.source_label) + .second) + return fail(node.location, node.source_id.source_label_text, + "Node identity is duplicated."); + for (const double coordinate : node.coordinates) + if (!std::isfinite(coordinate)) + return fail(node.location, node.source_id.source_label_text, + "Solid coordinates must be finite."); + } + std::set> element_ids; + std::set active_materials; + for (const auto& element : definition.solid_elements) { + if (!element_ids + .emplace(element.source_id.instance_name, + element.source_id.source_label) + .second || + element.section_index >= definition.solid_sections.size() || + element.material_index >= definition.materials.size()) { + return fail( + element.location, element.source_id.source_label_text, + "Solid identity, section, and material must resolve uniquely."); + } + if (definition.solid_sections[element.section_index].MaterialIndex() != + element.material_index) + return fail(element.location, element.source_id.source_label_text, + "Solid material assignment conflicts."); + active_materials.insert(element.material_index); + std::set distinct; + for (const auto node : element.node_indices) { + if (node >= active.size() || !distinct.insert(node).second) + return fail(element.location, element.source_id.source_label_text, + "Solid connectivity requires eight distinct active nodes."); + active[node] = true; + } + } + for (std::size_t node = 0U; node < active.size(); ++node) + if (!active[node]) + return fail(definition.nodes[node].location, + definition.nodes[node].source_id.source_label_text, + "Solid nodes must be connected."); + for (const auto& section : definition.solid_sections) + if (section.MaterialIndex() >= definition.materials.size() || + section.Name().empty()) + return fail(section.Location(), section.Name(), + "Solid material assignment is invalid."); + for (const auto index : active_materials) { + const auto& material = definition.materials[index]; + const double young = material.YoungsModulus(); + const double nu = material.PoissonsRatio(); + const double mu = young / (2.0 * (1.0 + nu)); + const double lambda = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + if (!std::isfinite(young) || !(young > 0.0) || !std::isfinite(nu) || + !(nu > -1.0 && nu < 0.5) || !std::isfinite(mu) || + !std::isfinite(lambda)) + return fail(material.Location(), material.Name(), + "Solid isotropic constants are invalid."); + } + std::set density_materials; + for (const auto& density : definition.densities) { + if (density.material_index >= definition.materials.size() || + !std::isfinite(density.rho) || !(density.rho > 0.0) || + !density_materials.insert(density.material_index).second) + return fail(density.location, "DENSITY", + "Solid density must be finite positive and unique."); + } + if (definition.dynamic_step || definition.frequency_step) + for (const auto material : active_materials) + if (density_materials.count(material) == 0U) + return fail({}, "DENSITY", + "Solid frequency and dynamic materials require density."); + std::set amplitude_names; + for (const auto& amplitude : definition.amplitudes) { + if (!definition.dynamic_step || amplitude.name_raw.empty() || + amplitude.name_canonical.empty() || + !amplitude_names.insert(amplitude.name_canonical).second || + amplitude.points.empty()) + return fail(amplitude.location, amplitude.name_raw, + "Solid amplitude identity is invalid."); + double previous = -1.0; + for (const auto& point : amplitude.points) { + if (!std::isfinite(point[0]) || !std::isfinite(point[1]) || + point[0] < 0.0 || !(point[0] > previous)) + return fail(amplitude.location, amplitude.name_raw, + "Solid amplitude pairs are invalid."); + previous = point[0]; + } + } + const auto validate_boundaries = [&fail](const auto& boundaries, + bool homogeneous) { + for (const auto& boundary : boundaries) { + if (boundary.first_dof < 1 || boundary.last_dof < boundary.first_dof || + boundary.last_dof > 3 || !std::isfinite(boundary.value) || + (homogeneous && boundary.value != 0.0)) + return fail(boundary.location, boundary.target, + "Solid physical boundaries require DOF1..3."); + } + return Status::Ok(); + }; + const auto validate_loads = [&fail, &definition](const auto& loads, + bool dynamic) { + for (const auto& load : loads) { + if (load.dof < 1 || load.dof > 3 || !std::isfinite(load.magnitude) || + (load.amplitude_index && + (!dynamic || *load.amplitude_index >= definition.amplitudes.size()))) + return fail(load.location, load.target, + "Solid load component or amplitude is invalid."); + } + return Status::Ok(); + }; + if (definition.frequency_step) + return validate_boundaries(definition.frequency_step->boundaries, true); + if (definition.dynamic_step) { + const auto& dynamic = *definition.dynamic_step; + if (!std::isfinite(dynamic.initial_increment) || + !(dynamic.initial_increment > 0.0) || + !std::isfinite(dynamic.time_period) || !(dynamic.time_period > 0.0) || + dynamic.increment_count == 0U || + dynamic.increment_count > dynamic.maximum_increments || + dynamic.initial_increment * + static_cast(dynamic.increment_count) != + dynamic.time_period || + dynamic.alpha != -0.05 || dynamic.beta != 0.275625 || + dynamic.gamma != 0.55) + return fail(dynamic.data_location, dynamic.source_step_name_raw, + "Solid dynamic controls are invalid."); + const auto status = validate_boundaries(dynamic.boundaries, true); + if (!status.IsOk()) return status; + return validate_loads(dynamic.loads, true); + } + const auto& step = definition.steps[0]; + if (step.direct && + (!std::isfinite(step.initial_increment) || + !(step.initial_increment > 0.0) || !std::isfinite(step.time_period) || + !(step.time_period > 0.0) || step.increment_count == 0U || + step.increment_count > step.maximum_increments || + step.initial_increment * static_cast(step.increment_count) != + step.time_period)) + return fail(step.location, step.name, + "Solid static fixed schedule is invalid."); + const auto status = validate_boundaries(step.boundaries, false); + if (!status.IsOk()) return status; + return validate_loads(step.loads, false); +} + } // namespace StepDefinition::StepDefinition(StaticStepDefinition definition) @@ -192,12 +371,17 @@ StepDefinition::StepDefinition(StaticStepDefinition definition) auto& load = definition.loads[source_order]; auto owned = std::make_unique( SourceTargetQuery{SourceEntityKind::kNode, "", load.target}, load.dof, - load.magnitude, source_order, std::move(load.location)); + load.magnitude, source_order, std::move(load.location), + load.amplitude_index, std::move(load.amplitude_name_raw), + std::move(load.amplitude_location)); loads_view_.push_back(std::cref(*owned)); concentrated_loads_view_.Add(*owned); loads_.push_back(std::move(owned)); } surface_pressure_loads_ = std::move(definition.surface_pressure_loads); + direct_static_ = definition.direct; + increment_count_ = definition.increment_count; + maximum_increments_ = definition.maximum_increments; } StepDefinition::StepDefinition(FrequencyStepDefinition definition) @@ -227,8 +411,10 @@ StepDefinition::StepDefinition(FrequencyStepDefinition definition) StepDefinition::StepDefinition(DynamicStepDefinition definition) : name_{definition.source_step_name_raw}, - procedure_kind_{StaticProcedureKind::kNonlinearFullNewton}, - nlgeom_{true}, + procedure_kind_{definition.nlgeom + ? StaticProcedureKind::kNonlinearFullNewton + : StaticProcedureKind::kLinearSmallDisplacement}, + nlgeom_{definition.nlgeom}, initial_increment_{definition.initial_increment}, time_period_{definition.time_period}, minimum_increment_{definition.initial_increment}, @@ -255,7 +441,9 @@ StepDefinition::StepDefinition(DynamicStepDefinition definition) auto& load = dynamic_->loads[source_order]; auto owned = std::make_unique( SourceTargetQuery{SourceEntityKind::kNode, "", load.target}, load.dof, - load.magnitude, source_order, std::move(load.location)); + load.magnitude, source_order, std::move(load.location), + load.amplitude_index, std::move(load.amplitude_name_raw), + std::move(load.amplitude_location)); loads_view_.push_back(std::cref(*owned)); concentrated_loads_view_.Add(*owned); loads_.push_back(std::move(owned)); @@ -325,11 +513,21 @@ double StepDefinition::MaximumIncrement() const noexcept { return maximum_increment_; } +bool StepDefinition::IsDirectStatic() const noexcept { return direct_static_; } +std::uint32_t StepDefinition::IncrementCount() const noexcept { + return increment_count_; +} +std::uint32_t StepDefinition::MaximumIncrements() const noexcept { + return maximum_increments_; +} + const SourceLocation& StepDefinition::Location() const noexcept { return location_; } Result Domain::Create(ModelDefinition definition) { + const Status solid_status = ValidateSolidDefinition(definition); + if (!solid_status.IsOk()) return Result::Failure(solid_status); const Status thermal_status = ValidateThermalDefinition(definition); if (!thermal_status.IsOk()) return Result::Failure(thermal_status); return Result::Success(Domain{std::move(definition)}); @@ -354,6 +552,23 @@ const DomainCollectionView& Domain::ShellElements() return shell_elements_view_; } +const DomainCollectionView& Domain::SolidElements() + const noexcept { + return solid_elements_view_; +} +const DomainCollectionView& Domain::SolidSections() + const noexcept { + return solid_sections_view_; +} +bool Domain::IsSolid() const noexcept { + return !solid_elements_view_.Empty() && + solid_elements_view_.Size() == elements_view_.Size(); +} +const std::vector& Domain::Amplitudes() + const noexcept { + return definition_.amplitudes; +} + const DomainCollectionView& Domain::Materials() const noexcept { return materials_view_; } @@ -496,6 +711,14 @@ Domain::Domain(ModelDefinition definition) element_properties_.push_back(std::move(owned)); } definition_.shell_sections.clear(); + const std::size_t solid_property_offset = element_properties_.size(); + for (auto& section : definition_.solid_sections) { + auto owned = std::make_unique(std::move(section)); + properties_view_.Add(*owned); + solid_sections_view_.Add(*owned); + element_properties_.push_back(std::move(owned)); + } + definition_.solid_sections.clear(); element_definitions_.reserve(definition_.elements.size() + definition_.shell_elements.size()); @@ -517,6 +740,16 @@ Domain::Domain(ModelDefinition definition) element_definitions_.push_back(std::move(owned)); } definition_.shell_elements.clear(); + for (auto& element : definition_.solid_elements) { + element.SynchronizeNodeIndices(); + element.SetPropertyIndex(static_cast(solid_property_offset + + element.section_index)); + auto owned = std::make_unique(std::move(element)); + elements_view_.Add(*owned); + solid_elements_view_.Add(*owned); + element_definitions_.push_back(std::move(owned)); + } + definition_.solid_elements.clear(); for (auto& element : definition_.thermal_elements) { const EntityIndex material_index = static_cast( structural_material_count + diff --git a/tests/unit/io/abaqus/hex8_semantics_test.cpp b/tests/unit/io/abaqus/hex8_semantics_test.cpp new file mode 100644 index 0000000..4308888 --- /dev/null +++ b/tests/unit/io/abaqus/hex8_semantics_test.cpp @@ -0,0 +1,449 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/assembly/load_assembler.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/assembly/sparse_assembler.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/io/abaqus/domain_mapper.h" +#include "fesa/io/abaqus/input_reader.h" + +namespace { +std::string Deck(const std::string& procedure = "*Static\n1,1,1,1\n", + const std::string& boundary = "10,1,3\n", + const std::string& extra = "") { + return "*Part,name=P\n*Node\n10,0,0,0\n20,1,0,0\n30,1,1,0\n40,0,1,0\n" + "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n" + "*Element,type=C3D8,elset=E\n99,10,20,30,40,50,60,70,80\n" + "*Solid Section,elset=E,material=M\n,\n*End Part\n" + "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n" + "*End Assembly\n*Material,name=M\n*Elastic\n1000,.3\n" + + extra + "*Step,name=Solid,inc=100\n" + procedure + + (boundary.empty() ? "" : "*Boundary\n" + boundary) + "*End Step\n"; +} +fesa::Result Map(const std::string& text) { + const auto path = + std::filesystem::temp_directory_path() / "fesa-hex8-semantics.inp"; + { + std::ofstream out{path}; + out << text; + } + auto syntax = fesa::AbaqusInputReader{}.Read(path); + std::filesystem::remove(path); + if (!syntax.HasValue()) + return fesa::Result::Failure(syntax.GetStatus()); + return fesa::AbaqusDomainMapper{}.Map(syntax.Value()); +} +TEST(Hex8Semantics, MapsEightSourceNodesAndExactlyThreePhysicalDofs) { + auto domain = Map(Deck()); + ASSERT_TRUE(domain.HasValue()) + << (domain.GetStatus().Diagnostics().empty() + ? "" + : domain.GetStatus().Diagnostics()[0].message); + ASSERT_EQ(domain.Value().Elements().Size(), 1U); + EXPECT_EQ(domain.Value().Elements()[0].SourceElementType(), "C3D8"); + EXPECT_EQ(domain.Value().Elements()[0].SourceId().source_label, 99); + EXPECT_EQ(domain.Value().Elements()[0].NodeIndices().size(), 8U); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + EXPECT_EQ(dofs.Value().FullDofCount(), 24U); + EXPECT_EQ(dofs.Value().ConstrainedDofCount(), 3U); +} +TEST(Hex8Semantics, InactiveZeroRotationsWarnAndDoNotCreateConstraints) { + auto domain = Map(Deck("*Static\n1,1,1,1\n", "10,1,6\n")); + ASSERT_TRUE(domain.HasValue()) + << (domain.GetStatus().Diagnostics().empty() + ? "" + : domain.GetStatus().Diagnostics()[0].message); + EXPECT_EQ(domain.Value().Steps()[0].BoundaryConditions().size(), 3U); + ASSERT_FALSE(domain.Value().Warnings().empty()); + EXPECT_EQ(domain.Value().Warnings().back().keyword, "BOUNDARY"); + EXPECT_EQ(domain.Value().Warnings().back().entity_identity, "10"); + EXPECT_GT(domain.Value().Warnings().back().location.line, 0U); +} +TEST(Hex8Semantics, RejectsNonzeroRotationsAndAllRotationalLoads) { + EXPECT_FALSE(Map(Deck("*Static\n1,1,1,1\n", "10,4,4,1\n")).HasValue()); + EXPECT_FALSE( + Map(Deck("*Static\n1,1,1,1\n", "10,1,3\n*Cload\n20,4,0\n")).HasValue()); +} + +std::string Replace(std::string input, const std::string& from, + const std::string& to) { + const auto position = input.find(from); + if (position == std::string::npos) + throw std::logic_error{"Missing mutation source"}; + input.replace(position, from.size(), to); + return input; +} +TEST(Hex8Semantics, RotationNoOpPreservesPhysicalMapAndLoads) { + auto physical = + Map(Deck("*Static\n1,1,1,1\n", "10,1,3,0\n*Cload\n20,2,17\n")); + auto extra = + Map(Deck("*Static\n1,1,1,1\n", "10,1,6,0\n20,4,6,-0\n*Cload\n20,2,17\n")); + ASSERT_TRUE(physical.HasValue()); + ASSERT_TRUE(extra.HasValue()); + auto pm = fesa::AnalysisModel::Create(physical.Value()); + auto em = fesa::AnalysisModel::Create(extra.Value()); + auto pd = fesa::DofManager::Create(pm.Value()); + auto ed = fesa::DofManager::Create(em.Value()); + ASSERT_TRUE(pd.HasValue()); + ASSERT_TRUE(ed.HasValue()); + EXPECT_EQ(pd.Value().FreeDofs(), ed.Value().FreeDofs()); + EXPECT_EQ(pd.Value().ConstrainedDofs(), ed.Value().ConstrainedDofs()); + EXPECT_EQ(pd.Value().GetSparsePattern().column_indices, + ed.Value().GetSparsePattern().column_indices); + EXPECT_EQ(ed.Value().Profile(), fesa::DofProfile::kSolidStructural); + EXPECT_EQ(ed.Value().ComponentsPerNode(), 3U); + EXPECT_THROW(ed.Value().FullDof(0, fesa::DofComponent::kUrx), + std::out_of_range); + const auto& element = extra.Value().Elements()[0]; + auto scatter = ed.Value().ElementScatter( + {element.SourceId(), + element.NodeIndices(), + {fesa::DofComponent::kUx, fesa::DofComponent::kUy, + fesa::DofComponent::kUz}}); + ASSERT_TRUE(scatter.HasValue()); + EXPECT_EQ(scatter.Value().size(), 24U); + EXPECT_EQ(scatter.Value().front(), 0U); + EXPECT_EQ(scatter.Value().back(), 23U); + auto pf = fesa::LoadAssembler::AssembleFullNodalLoad(pm.Value(), pd.Value()); + auto ef = fesa::LoadAssembler::AssembleFullNodalLoad(em.Value(), ed.Value()); + ASSERT_TRUE(pf.HasValue()); + ASSERT_TRUE(ef.HasValue()); + for (std::size_t i = 0; i < 24U; ++i) EXPECT_EQ(pf.Value()[i], ef.Value()[i]); +} +TEST(Hex8Semantics, RejectsTinyRotationBadTargetsOptionsAndDofRanges) { + for (const std::string row : + {"10,4,6,1e-100\n", "Missing,4,6,0\n", "10,0,3\n", "10,1,7\n"}) + EXPECT_FALSE(Map(Deck("*Static\n1,1,1,1\n", row)).HasValue()) << row; + for (const std::string option : {"type=VELOCITY", "amplitude=A", "op=NEW"}) + EXPECT_FALSE( + Map(Replace(Deck(), "*Boundary", "*Boundary," + option)).HasValue()); + for (const std::string row : {"10,4,6\n", "10,4,6,0\n", "10,4,6,\n"}) { + auto result = Map(Deck("*Static\n1,1,1,1\n", row)); + ASSERT_TRUE(result.HasValue()) << row; + EXPECT_TRUE(result.Value().Steps()[0].BoundaryConditions().empty()); + ASSERT_EQ(result.Value().Warnings().size(), 1U); + EXPECT_EQ(result.Value().Warnings()[0].code, + "inactive-solid-rotation-boundary"); + EXPECT_EQ(result.Value().Warnings()[0].severity, fesa::Severity::kWarning); + } +} +TEST(Hex8Semantics, AllowsEmptyFreeAndNoBoundaryModels) { + std::string all; + for (int node = 10; node <= 80; node += 10) + all += std::to_string(node) + ",1,6\n"; + for (const auto& boundaries : {all, std::string{}}) { + auto domain = Map(Deck("*Static\n1,1,1,1\n", boundaries)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + EXPECT_EQ(dofs.Value().FullDofCount(), 24U); + EXPECT_EQ(dofs.Value().FreeDofCount(), boundaries.empty() ? 24U : 0U); + } +} +TEST(Hex8Semantics, RetainsStaticDensityTargetsAndDirectSchedule) { + auto domain = + Map(Deck("*Static,direct\n.1,1,\n", "10,1,3,.25\n", "*Density\n2\n")); + ASSERT_TRUE(domain.HasValue()) + << (domain.GetStatus().Diagnostics().empty() + ? "" + : domain.GetStatus().Diagnostics()[0].message); + EXPECT_TRUE(domain.Value().IsSolid()); + ASSERT_EQ(domain.Value().SolidSections().Size(), 1U); + EXPECT_EQ(domain.Value().SolidSections()[0].MaterialIndex(), 0U); + EXPECT_EQ(domain.Value().Densities()[0].rho, 2.0); + const auto& step = domain.Value().Steps()[0]; + EXPECT_EQ(step.Name(), "Solid"); + EXPECT_TRUE(step.IsDirectStatic()); + EXPECT_EQ(step.IncrementCount(), 10U); + EXPECT_EQ(step.MaximumIncrements(), 100U); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + EXPECT_EQ(dofs.Value().PrescribedValues()[0], .25); + for (const std::string values : {".3,1\n", "0,1\n", ".1,11\n"}) + EXPECT_FALSE(Map(Deck("*Static,direct\n" + values)).HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "1000,.3", "1000,.5")).HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "1000,.3", "1000,-1")).HasValue()); +} +TEST(Hex8Semantics, SamplesIndependentAmplitudeLoadsAtActualEndpoints) { + const std::string amplitudes = + "*Density\n2\n*Amplitude,name=A\n0,0,.1,1\n10,1\n" + "*Amplitude,name=B,definition=TABULAR,time=STEP " + "TIME,value=RELATIVE\n0,2,1,4\n"; + const auto text = Deck("*Dynamic,direct\n.1,10,\n", + "10,1,6\n*Cload,amplitude=a\n20,1,10\n" + "*Cload,amplitude=B\n20,1,3\n*Cload\n20,1,5\n", + amplitudes); + auto domain = Map(text); + ASSERT_TRUE(domain.HasValue()) + << (domain.GetStatus().Diagnostics().empty() + ? "" + : domain.GetStatus().Diagnostics()[0].message); + ASSERT_EQ(domain.Value().Amplitudes().size(), 2U); + const auto& amplitude = domain.Value().Amplitudes()[0]; + EXPECT_EQ(amplitude.name_raw, "A"); + EXPECT_EQ(amplitude.Evaluate(0), 0); + EXPECT_EQ(amplitude.Evaluate(.05), .5); + EXPECT_EQ(amplitude.Evaluate(.1), 1); + EXPECT_EQ(amplitude.Evaluate(20), 1); + EXPECT_FALSE(domain.Value().Dynamic()->nlgeom); + EXPECT_EQ(domain.Value().Dynamic()->increment_count, 100U); + const auto& loads = domain.Value().Steps()[0].ConcentratedLoads(); + EXPECT_EQ(loads[0].AmplitudeIndex(), std::optional{0U}); + EXPECT_EQ(loads[1].AmplitudeIndex(), std::optional{1U}); + EXPECT_EQ(loads[0].AmplitudeName(), "a"); + EXPECT_GT(loads[0].AmplitudeLocation().line, 0U); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + for (double time : {0., .05, .1, 1., 10.}) { + auto load = fesa::LoadAssembler::AssembleFullNodalLoad(model.Value(), + dofs.Value(), time); + ASSERT_TRUE(load.HasValue()); + EXPECT_DOUBLE_EQ(load.Value()[3], + amplitude.Evaluate(time) * 10 + + domain.Value().Amplitudes()[1].Evaluate(time) * 3 + 5); + } + EXPECT_FALSE(Map(Replace(text, "10,1,6", "10,1,3,1e-100")).HasValue()); + EXPECT_FALSE( + Map(Replace(text, "amplitude=a", "amplitude=missing")).HasValue()); + EXPECT_TRUE( + Map(Replace(text, "name=Solid,inc=100", "name=Solid,nlgeom=YES,inc=100")) + .HasValue()); +} +TEST(Hex8Semantics, RejectsAmplitudeOptionsDataAndCrossFamilyUse) { + const std::string prefix = "*Density\n2\n"; + for (const std::string amplitude : + {"*Amplitude,name=A,time=TOTAL TIME\n0,0\n", + "*Amplitude,name=A,definition=SMOOTH STEP\n0,0\n", + "*Amplitude,name=A,value=ABSOLUTE\n0,0\n", + "*Amplitude,name=A,smooth=0\n0,0\n", "*Amplitude,name=A\n0,0,1\n", + "*Amplitude,name=A\n0,0,0,1\n", "*Amplitude,name=A\n1,0,0,1\n", + "*Amplitude,name=A\n-1,0\n", "*Amplitude,name=A\n0,nan\n", + "*Amplitude,name=A\n", + "*Amplitude,name=A\n0,0\n*Amplitude,name=a\n0,1\n"}) + EXPECT_FALSE( + Map(Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n", prefix + amplitude)) + .HasValue()) + << amplitude; + EXPECT_FALSE( + Map(Deck("*Static\n1,1,1,1\n", "10,1,3\n", "*Amplitude,name=A\n0,0\n")) + .HasValue()); + EXPECT_FALSE(Map(Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n")).HasValue()); +} +TEST(Hex8Semantics, RejectsAliasesConnectivityThicknessAndMixedModels) { + for (const std::string type : {"C3D8R", "C3D8I", "C3D8H", "C3D20"}) + EXPECT_FALSE(Map(Replace(Deck(), "type=C3D8", "type=" + type)).HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "99,10,20,30,40,50,60,70,80", + "99,10,20,30,40,50,60,70,70")) + .HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "99,10,20,30,40,50,60,70,80", + "99,10,20,30,40,50,60,70,999")) + .HasValue()); + EXPECT_FALSE( + Map(Replace(Deck(), "material=M\n,", "material=M\n2")).HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "*Solid Section", + "*Element,type=S4\n100,10,20,30,40\n*Solid Section")) + .HasValue()); +} +TEST(Hex8Semantics, PreservesShellSixComponentProfile) { + auto text = Replace(Deck(), "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n", ""); + text = Replace(text, "*Element,type=C3D8,elset=E\n99,10,20,30,40,50,60,70,80", + "*Element,type=S4\n99,10,20,30,40\n*Elset,elset=E\n99"); + text = Replace(text, "*Solid Section,elset=E,material=M\n,", + "*Shell Section,elset=E,material=M\n.1"); + auto domain = Map(text); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + EXPECT_EQ(dofs.Value().Profile(), fesa::DofProfile::kStructural); + EXPECT_EQ(dofs.Value().ComponentsPerNode(), 6U); + EXPECT_EQ(dofs.Value().FullDofCount(), 24U); + EXPECT_EQ(dofs.Value().FullDof(1, fesa::DofComponent::kUrx), 9U); +} +TEST(Hex8Semantics, MapsAllFiveUnmodifiedRequiredSourceDecks) { + auto root = std::filesystem::path{__FILE__}.parent_path(); + for (int i = 0; i < 4; ++i) root = root.parent_path(); + for (const std::string name : + {"hexa_linearstatic.inp", "hexa_nonlinearstatic.inp", "hexa_freq.inp", + "hexa_dynlinear.inp", "hexa_dynnonlinear.inp"}) { + auto syntax = + fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name); + ASSERT_TRUE(syntax.HasValue()) << name; + auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value()); + ASSERT_TRUE(domain.HasValue()) + << name << ": " + << (domain.GetStatus().Diagnostics().empty() + ? "" + : domain.GetStatus().Diagnostics()[0].message); + EXPECT_EQ(domain.Value().SolidElements().Size(), 300U); + EXPECT_EQ(domain.Value().Nodes().size(), 442U); + EXPECT_EQ(std::count_if(domain.Value().Warnings().begin(), + domain.Value().Warnings().end(), + [](const auto& warning) { + return warning.code == + "inactive-solid-rotation-boundary"; + }), + 3); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()) << name; + EXPECT_EQ(dofs.Value().FullDofCount(), 1326U); + EXPECT_EQ(dofs.Value().ConstrainedDofCount(), 51U); + } +} + +class ConstantMatrixElement final : public fesa::Element { + public: + ConstantMatrixElement(const fesa::ElementDefinition& definition, double entry) + : layout_{definition.SourceId(), + definition.NodeIndices(), + {fesa::DofComponent::kUx, fesa::DofComponent::kUy, + fesa::DofComponent::kUz}}, + entry_{entry} {} + const fesa::ElementDofLayout& DofLayout() const noexcept override { + return layout_; + } + fesa::Result ComputeStiffness() + const override { + fesa::Matrix values{24U, 24U}; + for (std::size_t row = 0U; row < 24U; ++row) + for (std::size_t column = 0U; column < 24U; ++column) + values(row, column) = entry_; + return fesa::Result::Success( + {layout_, std::move(values)}); + } + fesa::Result Recover( + const fesa::Vector&) const override { + return fesa::Result::Failure( + fesa::Status::Failure(fesa::FailureCategory::kModel, + {{fesa::Severity::kError, + "test-no-recovery", + {}, + "", + "", + "Assembly fixture has no recovery."}})); + } + + private: + fesa::ElementDofLayout layout_; + double entry_; +}; +TEST(Hex8Semantics, AssemblesFullSymmetricSolidMatricesInStableSourceOrder) { + const auto text = + Replace(Deck(), "99,10,20,30,40,50,60,70,80\n", + "99,10,20,30,40,50,60,70,80\n100,10,20,30,40,50,60,70,80\n" + "101,10,20,30,40,50,60,70,80\n"); + auto domain = Map(text); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ConstantMatrixElement first{domain.Value().Elements()[0], 1e16}; + ConstantMatrixElement second{domain.Value().Elements()[1], -1e16}; + ConstantMatrixElement third{domain.Value().Elements()[2], 1}; + const fesa::ElementView elements{std::cref(first), std::cref(second), + std::cref(third)}; + fesa::DofManager dofs; + ASSERT_TRUE(dofs.Build(model.Value(), elements).IsOk()); + auto serial = fesa::SparseAssembler::Assemble(elements, dofs, + fesa::SerialParallelFor{}); + auto parallel = + fesa::SparseAssembler::Assemble(elements, dofs, fesa::TbbParallelFor{}); + ASSERT_TRUE(serial.HasValue()); + ASSERT_TRUE(parallel.HasValue()); + EXPECT_EQ(serial.Value().Rows(), 24U); + EXPECT_EQ(serial.Value().Columns(), 24U); + EXPECT_EQ(serial.Value().Values().size(), 576U); + for (double value : serial.Value().Values()) EXPECT_EQ(value, 1.0); + EXPECT_EQ(serial.Value().RowOffsets(), parallel.Value().RowOffsets()); + EXPECT_EQ(serial.Value().ColumnIndices(), parallel.Value().ColumnIndices()); + EXPECT_EQ(serial.Value().Values(), parallel.Value().Values()); +} +TEST(Hex8Semantics, ValidatesOutputNoOpsAndRejectsUnapprovedDynamicChildren) { + for (const std::string output : + {"*Output,field,variable=PRESELECT\n*Node Output\nU,RF\n*Element " + "Output,directions=YES\nS,E\n", + "*Output,history,variable=PRESELECT\n"}) + EXPECT_TRUE( + Map(Replace(Deck(), "*End Step", output + "*End Step")).HasValue()); + EXPECT_FALSE(Map(Replace(Deck(), "*End Step", + "*Output,field,variable=OTHER\n*End Step")) + .HasValue()); + EXPECT_FALSE( + Map(Replace( + Deck(), "*End Step", + "*Output,field\n*Element Output,position=NODES\nS\n*End Step")) + .HasValue()); + auto dynamic = Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n", "*Density\n2\n"); + EXPECT_FALSE( + Map(Replace( + dynamic, "*End Step", + "*Output,field,variable=PRESELECT\n*Node Output\nU\n*End Step")) + .HasValue()); + EXPECT_FALSE(Map(Replace(dynamic, "*End Step", + "*Restart,write,frequency=1\n*End Step")) + .HasValue()); +} + +fesa::ModelDefinition SolidDefinition() { + fesa::ModelDefinition definition; + const std::array, 8> coordinates{{{0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1}, + {0, 1, 1}}}; + for (std::size_t i = 0; i < coordinates.size(); ++i) + definition.nodes.push_back({{"I", static_cast(10 * (i + 1)), + std::to_string(10 * (i + 1))}, + coordinates[i], + {"solid.inp", i + 1}}); + definition.materials.emplace_back("M", 1000, .3, + fesa::SourceLocation{"solid.inp", 20}); + definition.solid_sections.emplace_back(fesa::SourceEntityId{"P", 1, "E"}, "E", + 0, + fesa::SourceLocation{"solid.inp", 21}); + definition.solid_elements.emplace_back( + fesa::SourceEntityId{"I", 99, "99"}, + std::array{{0, 1, 2, 3, 4, 5, 6, 7}}, 0, 0, + fesa::SourceLocation{"solid.inp", 22}); + definition.steps.emplace_back( + "Solid", std::vector{}, + std::vector{}, 1., 1., 1., 1., + fesa::SourceLocation{"solid.inp", 23}); + return definition; +} +TEST(Hex8Semantics, RejectsInvalidProgrammaticSolidCandidates) { + EXPECT_TRUE(fesa::Domain::Create(SolidDefinition()).HasValue()); + auto repeated = SolidDefinition(); + repeated.solid_elements[0].node_indices[7] = 0; + EXPECT_FALSE(fesa::Domain::Create(std::move(repeated)).HasValue()); + auto material = SolidDefinition(); + material.materials[0].poisson_ratio = .5; + EXPECT_FALSE(fesa::Domain::Create(std::move(material)).HasValue()); + auto section = SolidDefinition(); + section.solid_elements[0].section_index = 99; + EXPECT_FALSE(fesa::Domain::Create(std::move(section)).HasValue()); + auto rotational = SolidDefinition(); + rotational.steps[0].loads.push_back({"10", 4, 0, {"solid.inp", 24}}); + EXPECT_FALSE(fesa::Domain::Create(std::move(rotational)).HasValue()); + auto unused = SolidDefinition(); + unused.nodes.push_back({{"I", 90, "90"}, {2, 2, 2}, {"solid.inp", 25}}); + EXPECT_FALSE(fesa::Domain::Create(std::move(unused)).HasValue()); +} +} // namespace