model.hpp

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <span>
#include <tuple>
 
#include "constraints/constraint.hpp"
#include "constraints/constraints.hpp"
#include "copy_nodes_to_state.hpp"
#include "dof_management/assemble_node_freedom_allocation_table.hpp"
#include "dof_management/compute_node_freedom_map_table.hpp"
#include "dof_management/create_constraint_freedom_table.hpp"
#include "dof_management/create_element_freedom_table.hpp"
#include "elements/beams/beams.hpp"
#include "elements/beams/beams_input.hpp"
#include "elements/beams/create_beams.hpp"
#include "elements/elements.hpp"
#include "elements/masses/create_masses.hpp"
#include "elements/masses/masses_input.hpp"
#include "elements/springs/create_springs.hpp"
#include "elements/springs/springs_input.hpp"
#include "mesh_connectivity.hpp"
#include "node.hpp"
#include "solver/solver.hpp"
#include "state/state.hpp"
 
namespace kynema {
 
template <
    typename DeviceType =
        Kokkos::Device<Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space>>
[[nodiscard]] inline Solver<DeviceType> CreateSolver(
    State<DeviceType>& state, Elements<DeviceType>& elements, Constraints<DeviceType>& constraints
) {
    dof::assemble_node_freedom_allocation_table(state, elements, constraints);
    dof::compute_node_freedom_map_table(state);
    dof::create_element_freedom_table(elements, state);
    dof::create_constraint_freedom_table(constraints, state);
 
    return {
        state.ID,
        state.active_dofs,
        state.node_freedom_map_table,
        elements.NumberOfNodesPerElement(),
        elements.NodeStateIndices(),
        constraints.num_dofs,
        constraints.base_active_dofs,
        constraints.target_active_dofs,
        constraints.base_node_freedom_table,
        constraints.target_node_freedom_table,
        constraints.row_range,
    };
}
 
class Model {
public:
    static constexpr size_t InvalidNodeID{0U};
 
    Model() = default;
 
    explicit Model(std::array<double, 3> gravity) : gravity_(gravity) {}
 
    //--------------------------------------------------------------------------
    // Miscellaneous
    //--------------------------------------------------------------------------
 
    void SetGravity(double x, double y, double z) {
        this->gravity_[0] = x;
        this->gravity_[1] = y;
        this->gravity_[2] = z;
    }
 
    //--------------------------------------------------------------------------
    // Nodes
    //--------------------------------------------------------------------------
 
    NodeBuilder AddNode() {
        const auto id = this->nodes_.size();
        this->nodes_.emplace_back(id);
        return NodeBuilder(this->nodes_.back());
    }
 
    [[nodiscard]] const Node& GetNode(size_t id) const { return this->nodes_[id]; }
 
    [[nodiscard]] Node& GetNode(size_t id) { return this->nodes_[id]; }
 
    [[nodiscard]] size_t NumNodes() const { return this->nodes_.size(); }
 
    [[nodiscard]] const std::vector<Node>& GetNodes() const { return this->nodes_; }
 
    //--------------------------------------------------------------------------
    // Beam Elements
    //--------------------------------------------------------------------------
 
    size_t AddBeamElement(
        std::span<const size_t> node_ids, std::span<const BeamSection> sections,
        std::span<const std::array<double, 2>> quadrature,
        const std::array<double, 6>& mu = {0., 0., 0., 0., 0., 0.}
    ) {
        const auto elem_id = this->beam_elements_.size();
        this->beam_elements_.emplace_back(elem_id, node_ids, sections, quadrature, mu);
        this->mesh_connectivity_.AddBeamElementConnectivity(elem_id, node_ids);
        return elem_id;
    }
 
    [[nodiscard]] const BeamElement& GetBeamElement(size_t id) const {
        return this->beam_elements_[id];
    }
 
    [[nodiscard]] BeamElement& GetBeamElement(size_t id) { return this->beam_elements_[id]; }
 
    [[nodiscard]] const std::vector<BeamElement>& GetBeamElements() const {
        return this->beam_elements_;
    }
 
    [[nodiscard]] size_t NumBeamElements() const { return this->beam_elements_.size(); }
 
    [[nodiscard]] BeamsInput CreateBeamsInput() const {
        return {this->beam_elements_, this->gravity_};
    }
 
    template <typename DeviceType>
    [[nodiscard]] Beams<DeviceType> CreateBeams() const {
        return kynema::CreateBeams<DeviceType>(this->CreateBeamsInput(), this->nodes_);
    }
 
    void TranslateBeam(size_t beam_elem_id, std::span<const double, 3> displacement) {
        const auto& beam_elem = this->beam_elements_[beam_elem_id];
        for (const auto& node_id : beam_elem.node_ids) {
            this->GetNode(node_id).Translate(displacement);
        }
    }
 
    void RotateBeamAboutPoint(
        size_t beam_elem_id, std::span<const double, 4> displacement_quaternion,
        std::span<const double, 3> point
    ) {
        const auto& beam_elem = this->beam_elements_[beam_elem_id];
        for (const auto& node_id : beam_elem.node_ids) {
            this->GetNode(node_id).RotateAboutPoint(displacement_quaternion, point);
        }
    }
 
    void SetBeamVelocityAboutPoint(
        size_t beam_elem_id, std::span<const double, 6> velocity, std::span<const double, 3> point
    ) {
        const auto& beam_elem = this->beam_elements_[beam_elem_id];
        for (const auto& node_id : beam_elem.node_ids) {
            this->GetNode(node_id).SetVelocityAboutPoint(velocity, point);
        }
    }
 
    void SetBeamAccelerationAboutPoint(
        size_t beam_elem_id, std::span<const double, 6> acceleration,
        std::span<const double, 3> omega, std::span<const double, 3> point
    ) {
        const auto& beam_elem = this->beam_elements_[beam_elem_id];
        for (const auto& node_id : beam_elem.node_ids) {
            this->GetNode(node_id).SetAccelerationAboutPoint(acceleration, omega, point);
        }
    }
 
    //--------------------------------------------------------------------------
    // Mass Elements
    //--------------------------------------------------------------------------
 
    size_t AddMassElement(const size_t node_id, const std::array<std::array<double, 6>, 6>& mass) {
        const auto elem_id = this->mass_elements_.size();
        this->mass_elements_.emplace_back(elem_id, node_id, mass);
        this->mesh_connectivity_.AddMassElementConnectivity(elem_id, node_id);
        return elem_id;
    }
 
    [[nodiscard]] const MassElement& GetMassElement(size_t id) const {
        return this->mass_elements_[id];
    }
 
    [[nodiscard]] MassElement& GetMassElement(size_t id) { return this->mass_elements_[id]; }
 
    [[nodiscard]] const std::vector<MassElement>& GetMassElements() const {
        return this->mass_elements_;
    }
 
    [[nodiscard]] size_t NumMassElements() const { return this->mass_elements_.size(); }
 
    template <typename DeviceType>
    [[nodiscard]] Masses<DeviceType> CreateMasses() const {
        return kynema::CreateMasses<DeviceType>(
            MassesInput(this->mass_elements_, this->gravity_), this->nodes_
        );
    }
 
    //--------------------------------------------------------------------------
    // Spring Elements
    //--------------------------------------------------------------------------
 
    size_t AddSpringElement(
        const size_t node1_id, const size_t node2_id, const double stiffness,
        const double undeformed_length
    ) {
        const auto elem_id = this->spring_elements_.size();
        this->spring_elements_.emplace_back(
            elem_id, std::array{node1_id, node2_id}, stiffness, undeformed_length
        );
        this->mesh_connectivity_.AddSpringElementConnectivity(
            elem_id, std::array{node1_id, node2_id}
        );
        return elem_id;
    }
 
    [[nodiscard]] const SpringElement& GetSpringElement(size_t id) const {
        return this->spring_elements_[id];
    }
 
    [[nodiscard]] SpringElement& GetSpringElement(size_t id) { return this->spring_elements_[id]; }
 
    [[nodiscard]] size_t NumSpringElements() const { return this->spring_elements_.size(); }
 
    template <typename DeviceType>
    [[nodiscard]] Springs<DeviceType> CreateSprings() const {
        return kynema::CreateSprings<DeviceType>(SpringsInput(this->spring_elements_), this->nodes_);
    }
 
    //--------------------------------------------------------------------------
    // Elements
    //--------------------------------------------------------------------------
 
    template <typename DeviceType>
    [[nodiscard]] Elements<DeviceType> CreateElements() const {
        return {
            this->CreateBeams<DeviceType>(),
            this->CreateMasses<DeviceType>(),
            this->CreateSprings<DeviceType>(),
        };
    }
 
    //--------------------------------------------------------------------------
    // State
    //--------------------------------------------------------------------------
 
    template <typename DeviceType>
    [[nodiscard]] State<DeviceType> CreateState() const {
        auto state = State<DeviceType>(this->nodes_.size());
        model::CopyNodesToState<DeviceType>(state, this->nodes_);
        return state;
    }
 
    //--------------------------------------------------------------------------
    // Constraints
    //--------------------------------------------------------------------------
 
    size_t AddFixedBC(const size_t node_id) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::FixedBC, std::array{InvalidNodeID, node_id}
        );
        this->mesh_connectivity_.AddConstraintConnectivity(id, std::vector<size_t>{node_id});
        return id;
    }
 
    size_t AddPrescribedBC(
        const size_t node_id,
        const std::array<double, 7>& initial_displacement = {0., 0., 0., 1., 0., 0., 0.}
    ) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::PrescribedBC, std::array{InvalidNodeID, node_id},
            std::array{0., 0., 0.}, initial_displacement
        );
        this->mesh_connectivity_.AddConstraintConnectivity(id, std::vector<size_t>{node_id});
        return id;
    }
 
    size_t AddRigidJointConstraint(std::span<const size_t, 2> node_ids) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::RigidJoint, std::array{node_ids[0], node_ids[1]}
        );
        this->mesh_connectivity_.AddConstraintConnectivity(
            id, std::vector<size_t>{node_ids[0], node_ids[1]}
        );
        return id;
    }
 
    size_t AddRevoluteJointConstraint(
        std::span<const size_t, 2> node_ids, std::span<const double, 3> axis, double* torque
    ) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::RevoluteJoint, std::array{node_ids[0], node_ids[1]},
            std::array{axis[0], axis[1], axis[2]}, std::array{0., 0., 0., 1., 0., 0., 0.}, torque
        );
        this->mesh_connectivity_.AddConstraintConnectivity(
            id, std::vector<size_t>{node_ids[0], node_ids[1]}
        );
        return id;
    }
 
    size_t AddRotationControl(
        std::span<const size_t, 2> node_ids, std::span<const double, 3> axis, double* control
    ) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::RotationControl, std::array{node_ids[0], node_ids[1]},
            std::array{axis[0], axis[1], axis[2]}, std::array{0., 0., 0., 1., 0., 0., 0.}, control
        );
        this->mesh_connectivity_.AddConstraintConnectivity(
            id, std::vector<size_t>{node_ids[0], node_ids[1]}
        );
        return id;
    }
 
    size_t AddFixedBC3DOFs(const size_t node_id) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::FixedBC3DOFs, std::array{InvalidNodeID, node_id}
        );
        this->mesh_connectivity_.AddConstraintConnectivity(id, std::vector<size_t>{node_id});
        return id;
    }
 
    size_t AddPrescribedBC3DOFs(const size_t node_id) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::PrescribedBC3DOFs, std::array{InvalidNodeID, node_id}
        );
        this->mesh_connectivity_.AddConstraintConnectivity(id, std::vector<size_t>{node_id});
        return id;
    }
 
    size_t AddRigidJoint6DOFsTo3DOFs(std::span<const size_t, 2> node_ids) {
        const auto id = this->constraints_.size();
        this->constraints_.emplace_back(
            id, constraints::ConstraintType::RigidJoint6DOFsTo3DOFs,
            std::array{node_ids[0], node_ids[1]}
        );
        this->mesh_connectivity_.AddConstraintConnectivity(
            id, std::vector<size_t>{node_ids[0], node_ids[1]}
        );
        return id;
    }
 
    [[nodiscard]] size_t NumConstraints() const { return this->constraints_.size(); }
 
    template <typename DeviceType>
    [[nodiscard]] Constraints<DeviceType> CreateConstraints() const {
        return Constraints<DeviceType>(this->constraints_, this->nodes_);
    }
 
    template <
        typename DeviceType = Kokkos::Device<
            Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space>>
    [[nodiscard]] std::tuple<State<DeviceType>, Elements<DeviceType>, Constraints<DeviceType>>
    CreateSystem() const {
        return {
            this->CreateState<DeviceType>(), this->CreateElements<DeviceType>(),
            this->CreateConstraints<DeviceType>()
        };
    }
 
    template <
        typename DeviceType = Kokkos::Device<
            Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space>>
    [[nodiscard]] std::tuple<
        State<DeviceType>, Elements<DeviceType>, Constraints<DeviceType>, Solver<DeviceType>>
    CreateSystemWithSolver() const {
        auto [state, elements, constraints] = this->CreateSystem<DeviceType>();
        auto solver = CreateSolver<DeviceType>(state, elements, constraints);
        return {state, elements, constraints, solver};
    }
 
    //--------------------------------------------------------------------------
    // Mesh Connectivity
    //--------------------------------------------------------------------------
 
    [[nodiscard]] const model::MeshConnectivity& GetMeshConnectivity() const {
        return this->mesh_connectivity_;
    }
 
    [[nodiscard]] model::MeshConnectivity& GetMeshConnectivity() { return this->mesh_connectivity_; }
 
    void ExportMeshConnectivityToYAML(const std::string& filename = "mesh_connectivity.yaml") const {
        this->mesh_connectivity_.ExportToYAML(filename);
    }
 
private:
    std::array<double, 3> gravity_ = {0., 0., 0.};      //< Gravity components
    std::vector<Node> nodes_;                           //< Nodes in the model
    std::vector<BeamElement> beam_elements_;            //< Beam elements in the model
    std::vector<MassElement> mass_elements_;            //< Mass elements in the model
    std::vector<SpringElement> spring_elements_;        //< Spring elements in the model
    std::vector<constraints::Constraint> constraints_;  //< Constraints in the model
    model::MeshConnectivity
        mesh_connectivity_;  //< Mesh connectivity tracking element-node relationships
};
 
}  // namespace kynema