diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 55e56f6d6..a0fecce27 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,434 +1,441 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Oct 15 2015 * * @brief Model implementation for StucturalMechanics elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_elemental_field.hh" #include "dumper_paraview.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.cc" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), time_step(NAN), f_m2a(1.0), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_DEBUG_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("paraview_all", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { // <<<< This is the SolidMechanicsModel implementation for future ref >>>> // material_index.initialize(mesh, _element_kind = _ek_not_defined, // _default_value = UInt(-1), _with_nb_element = true); // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, // _with_nb_element = true); // Model::initFullImpl(options); // // initialize pbc // if (this->pbc_pair.size() != 0) // this->initPBC(); // // initialize the materials // if (this->parser.getLastParsedFile() != "") { // this->instantiateMaterials(); // } // this->initMaterials(); // this->initBC(*this, *displacement, *displacement_increment, *external_force); // <<<< END >>>> Model::initFullImpl(options); } +/* -------------------------------------------------------------------------- */ + +void StructuralMechanicsModel::initFEEngineBoundary() { + /// TODO: this function should not be reimplemented + /// we're just avoiding a call to Model::initFEEngineBoundary() +} + /* -------------------------------------------------------------------------- */ // void StructuralMechanicsModel::setTimeStep(Real time_step) { // this->time_step = time_step; // #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); // #endif // } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force_momentum, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force_momentum, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(velocity, spatial_dimension, "velocity"); this->allocNodalField(acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) { computeRotationMatrix(type); element_material.alloc(mesh.getNbElement(type), 1, type); } getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); getDOFManager().getMatrix("K").clear(); for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); auto & K = getDOFManager().getMatrix("K"); internal_force_momentum->clear(); K.matVecMul(*displacement_rotation, *internal_force_momentum); *internal_force_momentum *= -1.; getDOFManager().assembleToResidual("displacement", *external_force_momentum, 1.); getDOFManager().assembleToResidual("displacement", *internal_force_momentum, 1.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if (!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom * nb_nodes_per_element * nb_degree_of_freedom * nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Array rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); auto T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { auto & T = *T_it; auto & R = *R_it; for (UInt k = 0; k < nb_nodes_per_element; ++k) { for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0, padding_flag); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "force") { return mesh.createStridedNodalField(external_force_momentum, group_name, n, 0, padding_flag); } if (field_name == "momentum") { return mesh.createStridedNodalField(external_force_momentum, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force_momentum, group_name, n, 0, padding_flag); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField(internal_force_momentum, group_name, nb_degree_of_freedom - n, n, padding_flag); } return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool, const ElementKind & kind, const std::string &) { dumper::Field * field = NULL; if (field_name == "element_index_by_material") field = mesh.createElementalField( field_name, group_name, this->spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") assembleStiffnessMatrix(); } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { this->getDOFManager().assembleToResidual("displacement", *this->external_force_momentum, 1); } /* -------------------------------------------------------------------------- */ /* Virtual methods from MeshEventHandler */ /* -------------------------------------------------------------------------- */ /// function to implement to react on akantu::NewNodesEvent void StructuralMechanicsModel::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent & /*event*/) {} /// function to implement to react on akantu::RemovedNodesEvent void StructuralMechanicsModel::onNodesRemoved( const Array & /*nodes_list*/, const Array & /*new_numbering*/, const RemovedNodesEvent & /*event*/) {} /// function to implement to react on akantu::NewElementsEvent void StructuralMechanicsModel::onElementsAdded( const Array & /*elements_list*/, const NewElementsEvent & /*event*/) {} /// function to implement to react on akantu::RemovedElementsEvent void StructuralMechanicsModel::onElementsRemoved( const Array & /*elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const RemovedElementsEvent & /*event*/) {} /// function to implement to react on akantu::ChangedElementsEvent void StructuralMechanicsModel::onElementsChanged( const Array & /*old_elements_list*/, const Array & /*new_elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const ChangedElementsEvent & /*event*/) {} /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID( const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_static: { options.non_linear_solver_type = _nls_linear; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index bef47cb68..8ec1f8739 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,329 +1,334 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Jan 21 2016 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_named_argument.hh" #include "boundary_condition.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeStructural; } // namespace akantu namespace akantu { struct StructuralMaterial { Real E{0}; Real A{1}; Real I{0}; Real Iz{0}; Real Iy{0}; Real GJ{0}; Real rho{0}; Real t{0}; Real nu{0}; }; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; StructuralMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "structural_mechanics_model", const MemoryID & memory_id = 0); virtual ~StructuralMechanicsModel(); + /// Init full model void initFullImpl(const ModelOptions & options) override; + /// Init boundary FEEngine + void initFEEngineBoundary() override; + + /* ------------------------------------------------------------------------ */ /* Virtual methods from SolverCallback */ /* ------------------------------------------------------------------------ */ /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; /// callback to assemble the residual (rhs) void assembleResidual() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from MeshEventHandler */ /* ------------------------------------------------------------------------ */ /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); /// update the residual vector void updateResidual(); void computeRotationMatrix(const ElementType & type); protected: /// compute Rotation Matrices template void computeRotationMatrix(__attribute__((unused)) Array & rotations) {} /* ------------------------------------------------------------------------ */ /* Mass (structural_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template void assembleStiffnessMatrix(); template void assembleMass(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const ElementKind & kind, const std::string & fe_engine_id = ""); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step // void setTimeStep(Real time_step, const ID & solver_id = "") override; /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force_momentum, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force_momentum, Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); void addMaterial(StructuralMaterial & material) { materials.push_back(material); } /** * @brief set the StructuralMechanicsModel::increment_flag to on, the * activate the * update of the StructuralMechanicsModel::increment vector */ void setIncrementFlagOn(); const StructuralMaterial & getMaterial(const Element & element) const { return materials[element_material(element)]; } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis template void computeForcesByGlobalTractionArray(const Array & tractions); /// Compute Linear load function set in local axis template void computeForcesByLocalTractionArray(const Array & tractions); /// compute force vector from a function(x,y,momentum) that describe stresses // template // void computeForcesFromFunction(BoundaryFunction in_function, // BoundaryFunctionType function_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement_rotation{nullptr}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// forces array Array * internal_force_momentum{nullptr}; /// forces array Array * external_force_momentum{nullptr}; /// lumped mass array Array * mass{nullptr}; /// boundaries array Array * blocked_dofs{nullptr}; /// stress array ElementTypeMapArray stress; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; // /// analysis method check the list in akantu::AnalysisMethod // AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; /* ------------------------------------------------------------------------ */ std::vector materials; }; } // namespace akantu #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ */