feat: add solver core skeleton

This commit is contained in:
NINI
2026-06-12 02:25:07 +09:00
parent 4e7fd1087d
commit cbd1a6c5d7
46 changed files with 1911 additions and 19 deletions
+34
View File
@@ -0,0 +1,34 @@
#pragma once
namespace fesa::analysis {
class Analysis {
public:
virtual ~Analysis() = default;
void run()
{
initialize();
build_analysis_model();
build_dof_map();
build_sparse_pattern();
assemble();
apply_boundary_conditions();
solve();
update_state();
write_results();
}
protected:
virtual void initialize() {}
virtual void build_analysis_model() {}
virtual void build_dof_map() {}
virtual void build_sparse_pattern() {}
virtual void assemble() {}
virtual void apply_boundary_conditions() {}
virtual void solve() {}
virtual void update_state() {}
virtual void write_results() {}
};
} // namespace fesa::analysis
+73
View File
@@ -0,0 +1,73 @@
#pragma once
#include <fesa/model/domain.hpp>
#include <stdexcept>
#include <vector>
namespace fesa::analysis {
class AnalysisModel {
public:
AnalysisModel(const model::Domain& domain, core::StepId step_id)
: domain_(domain), step_(domain.find_step(step_id))
{
if (step_ == nullptr) {
throw std::invalid_argument("analysis step not found");
}
for (const auto& element : domain_.elements()) {
active_elements_.push_back(&element);
}
for (const auto& boundary_condition : step_->boundary_conditions()) {
active_boundary_conditions_.push_back(&boundary_condition);
}
for (const auto& load : step_->loads()) {
active_loads_.push_back(&load);
}
}
const model::Domain& domain() const
{
return domain_;
}
const model::AnalysisStep& step() const
{
return *step_;
}
const std::vector<const model::Element*>& active_elements() const
{
return active_elements_;
}
const std::vector<const model::BoundaryCondition*>& active_boundary_conditions() const
{
return active_boundary_conditions_;
}
const std::vector<const model::Load*>& active_loads() const
{
return active_loads_;
}
const model::Property* property_for(const model::Element& element) const
{
return domain_.find_property(element.property_id());
}
const model::Material* material_for(const model::Property& property) const
{
return domain_.find_material(property.material_id());
}
private:
const model::Domain& domain_;
const model::AnalysisStep* step_;
std::vector<const model::Element*> active_elements_;
std::vector<const model::BoundaryCondition*> active_boundary_conditions_;
std::vector<const model::Load*> active_loads_;
};
} // namespace fesa::analysis
+146
View File
@@ -0,0 +1,146 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <stdexcept>
#include <utility>
#include <vector>
namespace fesa::analysis {
struct IterationState {
double time = 0.0;
int increment = 0;
int iteration = 0;
};
class AnalysisState {
public:
explicit AnalysisState(int total_dof_count)
: displacement_(vector_of(total_dof_count)),
velocity_(vector_of(total_dof_count)),
acceleration_(vector_of(total_dof_count)),
temperature_(vector_of(total_dof_count)),
external_force_(vector_of(total_dof_count)),
internal_force_(vector_of(total_dof_count)),
residual_(vector_of(total_dof_count))
{
}
const std::vector<double>& displacement() const
{
return displacement_;
}
const std::vector<double>& velocity() const
{
return velocity_;
}
const std::vector<double>& acceleration() const
{
return acceleration_;
}
const std::vector<double>& temperature() const
{
return temperature_;
}
const std::vector<double>& external_force() const
{
return external_force_;
}
const std::vector<double>& internal_force() const
{
return internal_force_;
}
const std::vector<double>& residual() const
{
return residual_;
}
void set_displacement(std::vector<double> values)
{
assign_same_size(displacement_, std::move(values));
}
void set_external_force(std::vector<double> values)
{
assign_same_size(external_force_, std::move(values));
}
void set_internal_force(std::vector<double> values)
{
assign_same_size(internal_force_, std::move(values));
}
void update_residual()
{
for (std::size_t index = 0; index < residual_.size(); ++index) {
residual_[index] = external_force_[index] - internal_force_[index];
}
}
IterationState& iteration_state()
{
return iteration_state_;
}
const IterationState& iteration_state() const
{
return iteration_state_;
}
void set_element_state(core::ElementId element_id, std::vector<double> state)
{
for (auto& entry : element_states_) {
if (entry.first.value == element_id.value) {
entry.second = std::move(state);
return;
}
}
element_states_.push_back({element_id, std::move(state)});
}
const std::vector<double>* element_state(core::ElementId element_id) const
{
for (const auto& entry : element_states_) {
if (entry.first.value == element_id.value) {
return &entry.second;
}
}
return nullptr;
}
private:
static std::vector<double> vector_of(int size)
{
if (size < 0) {
throw std::invalid_argument("negative dof count");
}
return std::vector<double>(static_cast<std::size_t>(size), 0.0);
}
static void assign_same_size(std::vector<double>& target, std::vector<double> values)
{
if (target.size() != values.size()) {
throw std::invalid_argument("vector size mismatch");
}
target = std::move(values);
}
std::vector<double> displacement_;
std::vector<double> velocity_;
std::vector<double> acceleration_;
std::vector<double> temperature_;
std::vector<double> external_force_;
std::vector<double> internal_force_;
std::vector<double> residual_;
IterationState iteration_state_;
std::vector<std::pair<core::ElementId, std::vector<double>>> element_states_;
};
} // namespace fesa::analysis
@@ -0,0 +1,66 @@
#pragma once
#include <fesa/analysis/analysis.hpp>
#include <fesa/analysis/analysis_model.hpp>
#include <fesa/analysis/analysis_state.hpp>
#include <fesa/fem/dof_manager.hpp>
#include <memory>
namespace fesa::analysis {
class LinearStaticAnalysis : public Analysis {
public:
LinearStaticAnalysis(const model::Domain& domain, core::StepId step_id)
: domain_(domain), step_id_(step_id)
{
}
const AnalysisModel* analysis_model() const
{
return analysis_model_.get();
}
const AnalysisState* state() const
{
return state_.get();
}
protected:
void build_analysis_model() override
{
analysis_model_ = std::make_unique<AnalysisModel>(domain_, step_id_);
}
void build_dof_map() override
{
dof_manager_ = std::make_unique<fem::DofManager>();
for (const auto* element : analysis_model_->active_elements()) {
for (const auto node_id : element->node_ids()) {
dof_manager_->define_node_dofs(node_id, {
model::DofComponent::ux,
model::DofComponent::uy,
model::DofComponent::uz
});
}
}
for (const auto* boundary_condition : analysis_model_->active_boundary_conditions()) {
dof_manager_->apply_boundary_condition(*boundary_condition);
}
dof_manager_->number_equations();
}
void update_state() override
{
state_ = std::make_unique<AnalysisState>(dof_manager_->total_dof_count());
}
private:
const model::Domain& domain_;
core::StepId step_id_;
std::unique_ptr<AnalysisModel> analysis_model_;
std::unique_ptr<fem::DofManager> dof_manager_;
std::unique_ptr<AnalysisState> state_;
};
} // namespace fesa::analysis
+19
View File
@@ -0,0 +1,19 @@
#pragma once
#include <string>
namespace fesa::core {
enum class Severity {
info,
warning,
error
};
struct Diagnostic {
Severity severity;
std::string code;
std::string message;
};
} // namespace fesa::core
+25
View File
@@ -0,0 +1,25 @@
#pragma once
namespace fesa::core {
struct NodeId {
int value;
};
struct ElementId {
int value;
};
struct MaterialId {
int value;
};
struct PropertyId {
int value;
};
struct StepId {
int value;
};
} // namespace fesa::core
+43
View File
@@ -0,0 +1,43 @@
#pragma once
#include <fesa/core/diagnostic.hpp>
#include <utility>
#include <vector>
namespace fesa::core {
class Status {
public:
static Status ok()
{
return Status{};
}
static Status failure(Diagnostic diagnostic)
{
Status status;
status.add(std::move(diagnostic));
return status;
}
bool is_ok() const
{
return diagnostics_.empty();
}
const std::vector<Diagnostic>& diagnostics() const
{
return diagnostics_;
}
void add(Diagnostic diagnostic)
{
diagnostics_.push_back(std::move(diagnostic));
}
private:
std::vector<Diagnostic> diagnostics_;
};
} // namespace fesa::core
+18
View File
@@ -0,0 +1,18 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <fesa/model/boundary_condition.hpp>
namespace fesa::fem {
struct DofKey {
core::NodeId node_id;
model::DofComponent component;
};
inline bool operator==(const DofKey& lhs, const DofKey& rhs)
{
return lhs.node_id.value == rhs.node_id.value && lhs.component == rhs.component;
}
} // namespace fesa::fem
+151
View File
@@ -0,0 +1,151 @@
#pragma once
#include <fesa/fem/dof_key.hpp>
#include <algorithm>
#include <optional>
#include <stdexcept>
#include <utility>
#include <vector>
namespace fesa::fem {
class DofManager {
public:
void define_node_dofs(core::NodeId node_id, std::vector<model::DofComponent> components)
{
for (const auto component : components) {
DofKey key{node_id, component};
if (find_record(key) == records_.end()) {
records_.push_back({key});
}
}
}
void apply_boundary_condition(const model::BoundaryCondition& bc)
{
auto record = find_record({bc.node_id(), bc.component()});
if (record == records_.end()) {
throw std::invalid_argument("boundary condition references undefined dof");
}
record->constrained = true;
}
void number_equations()
{
std::sort(records_.begin(), records_.end(), [](const Record& lhs, const Record& rhs) {
if (lhs.key.node_id.value != rhs.key.node_id.value) {
return lhs.key.node_id.value < rhs.key.node_id.value;
}
return static_cast<int>(lhs.key.component) < static_cast<int>(rhs.key.component);
});
int free_id = 0;
for (int equation_id = 0; equation_id < static_cast<int>(records_.size()); ++equation_id) {
records_[equation_id].equation_id = equation_id;
if (records_[equation_id].constrained) {
records_[equation_id].free_equation_id = std::nullopt;
} else {
records_[equation_id].free_equation_id = free_id++;
}
}
sparse_pattern_.clear();
for (int row = 0; row < free_id; ++row) {
for (int column = 0; column < free_id; ++column) {
sparse_pattern_.push_back({row, column});
}
}
}
int total_dof_count() const
{
return static_cast<int>(records_.size());
}
int free_dof_count() const
{
return static_cast<int>(
std::count_if(records_.begin(), records_.end(), [](const Record& record) {
return !record.constrained;
})
);
}
int constrained_dof_count() const
{
return total_dof_count() - free_dof_count();
}
bool is_constrained(DofKey key) const
{
return require_record(key).constrained;
}
int equation_id(DofKey key) const
{
return require_record(key).equation_id;
}
std::optional<int> free_equation_id(DofKey key) const
{
return require_record(key).free_equation_id;
}
std::vector<double> expand_free_vector(const std::vector<double>& free_values) const
{
if (free_values.size() != static_cast<std::size_t>(free_dof_count())) {
throw std::invalid_argument("free vector size does not match dof manager");
}
std::vector<double> full(records_.size(), 0.0);
for (const auto& record : records_) {
if (record.free_equation_id.has_value()) {
full[static_cast<std::size_t>(record.equation_id)] =
free_values[static_cast<std::size_t>(*record.free_equation_id)];
}
}
return full;
}
const std::vector<std::pair<int, int>>& sparse_pattern() const
{
return sparse_pattern_;
}
private:
struct Record {
DofKey key;
bool constrained = false;
int equation_id = -1;
std::optional<int> free_equation_id;
};
std::vector<Record>::iterator find_record(DofKey key)
{
return std::find_if(records_.begin(), records_.end(), [key](const Record& record) {
return record.key == key;
});
}
std::vector<Record>::const_iterator find_record(DofKey key) const
{
return std::find_if(records_.begin(), records_.end(), [key](const Record& record) {
return record.key == key;
});
}
const Record& require_record(DofKey key) const
{
const auto record = find_record(key);
if (record == records_.end()) {
throw std::invalid_argument("dof is not defined");
}
return *record;
}
std::vector<Record> records_;
std::vector<std::pair<int, int>> sparse_pattern_;
};
} // namespace fesa::fem
+56
View File
@@ -0,0 +1,56 @@
#pragma once
#include <fesa/model/boundary_condition.hpp>
#include <fesa/model/load.hpp>
#include <string>
#include <utility>
#include <vector>
namespace fesa::model {
class AnalysisStep {
public:
AnalysisStep(core::StepId id, std::string name)
: id_(id), name_(std::move(name))
{
}
core::StepId id() const
{
return id_;
}
const std::string& name() const
{
return name_;
}
void add_boundary_condition(BoundaryCondition bc)
{
boundary_conditions_.push_back(std::move(bc));
}
void add_load(Load load)
{
loads_.push_back(std::move(load));
}
const std::vector<BoundaryCondition>& boundary_conditions() const
{
return boundary_conditions_;
}
const std::vector<Load>& loads() const
{
return loads_;
}
private:
core::StepId id_;
std::string name_;
std::vector<BoundaryCondition> boundary_conditions_;
std::vector<Load> loads_;
};
} // namespace fesa::model
+45
View File
@@ -0,0 +1,45 @@
#pragma once
#include <fesa/core/ids.hpp>
namespace fesa::model {
enum class DofComponent {
ux,
uy,
uz,
rx,
ry,
rz,
temperature
};
class BoundaryCondition {
public:
BoundaryCondition(core::NodeId node_id, DofComponent component, double value)
: node_id_(node_id), component_(component), value_(value)
{
}
core::NodeId node_id() const
{
return node_id_;
}
DofComponent component() const
{
return component_;
}
double value() const
{
return value_;
}
private:
core::NodeId node_id_;
DofComponent component_;
double value_;
};
} // namespace fesa::model
+124
View File
@@ -0,0 +1,124 @@
#pragma once
#include <fesa/model/analysis_step.hpp>
#include <fesa/model/element.hpp>
#include <fesa/model/material.hpp>
#include <fesa/model/node.hpp>
#include <fesa/model/property.hpp>
#include <utility>
#include <vector>
namespace fesa::model {
class Domain {
public:
void add_node(Node node)
{
nodes_.push_back(std::move(node));
}
void add_element(Element element)
{
elements_.push_back(std::move(element));
}
void add_material(Material material)
{
materials_.push_back(std::move(material));
}
void add_property(Property property)
{
properties_.push_back(std::move(property));
}
void add_step(AnalysisStep step)
{
steps_.push_back(std::move(step));
}
const std::vector<Node>& nodes() const
{
return nodes_;
}
const std::vector<Element>& elements() const
{
return elements_;
}
const std::vector<Material>& materials() const
{
return materials_;
}
const std::vector<Property>& properties() const
{
return properties_;
}
const std::vector<AnalysisStep>& steps() const
{
return steps_;
}
const Node* find_node(core::NodeId id) const
{
for (const auto& node : nodes_) {
if (node.id().value == id.value) {
return &node;
}
}
return nullptr;
}
const Element* find_element(core::ElementId id) const
{
for (const auto& element : elements_) {
if (element.id().value == id.value) {
return &element;
}
}
return nullptr;
}
const Material* find_material(core::MaterialId id) const
{
for (const auto& material : materials_) {
if (material.id().value == id.value) {
return &material;
}
}
return nullptr;
}
const Property* find_property(core::PropertyId id) const
{
for (const auto& property : properties_) {
if (property.id().value == id.value) {
return &property;
}
}
return nullptr;
}
const AnalysisStep* find_step(core::StepId id) const
{
for (const auto& step : steps_) {
if (step.id().value == id.value) {
return &step;
}
}
return nullptr;
}
private:
std::vector<Node> nodes_;
std::vector<Element> elements_;
std::vector<Material> materials_;
std::vector<Property> properties_;
std::vector<AnalysisStep> steps_;
};
} // namespace fesa::model
+56
View File
@@ -0,0 +1,56 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <utility>
#include <vector>
namespace fesa::model {
enum class ElementTopology {
truss2,
bar2,
unknown
};
class Element {
public:
Element(core::ElementId id,
ElementTopology topology,
std::vector<core::NodeId> node_ids,
core::PropertyId property_id)
: id_(id),
topology_(topology),
node_ids_(std::move(node_ids)),
property_id_(property_id)
{
}
core::ElementId id() const
{
return id_;
}
ElementTopology topology() const
{
return topology_;
}
const std::vector<core::NodeId>& node_ids() const
{
return node_ids_;
}
core::PropertyId property_id() const
{
return property_id_;
}
private:
core::ElementId id_;
ElementTopology topology_;
std::vector<core::NodeId> node_ids_;
core::PropertyId property_id_;
};
} // namespace fesa::model
+35
View File
@@ -0,0 +1,35 @@
#pragma once
#include <fesa/model/boundary_condition.hpp>
namespace fesa::model {
class Load {
public:
Load(core::NodeId node_id, DofComponent component, double value)
: node_id_(node_id), component_(component), value_(value)
{
}
core::NodeId node_id() const
{
return node_id_;
}
DofComponent component() const
{
return component_;
}
double value() const
{
return value_;
}
private:
core::NodeId node_id_;
DofComponent component_;
double value_;
};
} // namespace fesa::model
+32
View File
@@ -0,0 +1,32 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <string>
#include <utility>
namespace fesa::model {
class Material {
public:
Material(core::MaterialId id, std::string name)
: id_(id), name_(std::move(name))
{
}
core::MaterialId id() const
{
return id_;
}
const std::string& name() const
{
return name_;
}
private:
core::MaterialId id_;
std::string name_;
};
} // namespace fesa::model
+31
View File
@@ -0,0 +1,31 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <array>
namespace fesa::model {
class Node {
public:
Node(core::NodeId id, std::array<double, 3> coordinates)
: id_(id), coordinates_(coordinates)
{
}
core::NodeId id() const
{
return id_;
}
const std::array<double, 3>& coordinates() const
{
return coordinates_;
}
private:
core::NodeId id_;
std::array<double, 3> coordinates_;
};
} // namespace fesa::model
+38
View File
@@ -0,0 +1,38 @@
#pragma once
#include <fesa/core/ids.hpp>
#include <string>
#include <utility>
namespace fesa::model {
class Property {
public:
Property(core::PropertyId id, std::string name, core::MaterialId material_id)
: id_(id), name_(std::move(name)), material_id_(material_id)
{
}
core::PropertyId id() const
{
return id_;
}
const std::string& name() const
{
return name_;
}
core::MaterialId material_id() const
{
return material_id_;
}
private:
core::PropertyId id_;
std::string name_;
core::MaterialId material_id_;
};
} // namespace fesa::model
+108
View File
@@ -0,0 +1,108 @@
#pragma once
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace fesa::results {
enum class FieldLocation {
nodal,
element,
integration_point
};
struct FieldOutput {
std::string name;
FieldLocation location;
std::vector<std::string> components;
std::vector<int> entity_ids;
std::vector<double> values;
};
struct HistoryOutput {
std::string name;
std::vector<double> time;
std::vector<double> values;
};
class ResultFrame {
public:
ResultFrame(int frame_id, double time)
: frame_id_(frame_id), time_(time)
{
}
int frame_id() const
{
return frame_id_;
}
double time() const
{
return time_;
}
void add_field_output(FieldOutput output)
{
if (output.components.empty()) {
throw std::invalid_argument("field output must have components");
}
if (output.entity_ids.size() * output.components.size() != output.values.size()) {
throw std::invalid_argument("field output values do not match row shape");
}
field_outputs_.push_back(std::move(output));
}
void add_history_output(HistoryOutput output)
{
history_outputs_.push_back(std::move(output));
}
const std::vector<FieldOutput>& field_outputs() const
{
return field_outputs_;
}
const std::vector<HistoryOutput>& history_outputs() const
{
return history_outputs_;
}
private:
int frame_id_;
double time_;
std::vector<FieldOutput> field_outputs_;
std::vector<HistoryOutput> history_outputs_;
};
class ResultStep {
public:
explicit ResultStep(std::string name)
: name_(std::move(name))
{
}
const std::string& name() const
{
return name_;
}
ResultFrame& add_frame(int frame_id, double time)
{
frames_.push_back(ResultFrame{frame_id, time});
return frames_.back();
}
const std::vector<ResultFrame>& frames() const
{
return frames_;
}
private:
std::string name_;
std::vector<ResultFrame> frames_;
};
} // namespace fesa::results