diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc index fa31d29ff..89421f389 100644 --- a/python/py_dof_manager.cc +++ b/python/py_dof_manager.cc @@ -1,224 +1,213 @@ /* * Copyright (©) 2018-2021 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 <http://www.gnu.org/licenses/>. */ /* -------------------------------------------------------------------------- */ #include "py_dof_manager.hh" #include "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include <dof_manager.hh> #include <non_linear_solver.hh> #include <solver_callback.hh> /* -------------------------------------------------------------------------- */ #include <pybind11/operators.h> #include <pybind11/pybind11.h> #include <pybind11/stl.h> /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { class PySolverCallback : public SolverCallback { public: using SolverCallback::SolverCallback; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(MatrixType, SolverCallback, getMatrixType, matrix_id); } /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleLumpedMatrix, matrix_id); } /// callback to assemble the residual (rhs) void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleResidual); } /// callback for the predictor (in case of dynamic simulation) void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, predictor); } /// callback for the corrector (in case of dynamic simulation) void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, corrector); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, beforeSolveStep); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, afterSolveStep, converged); } }; class PyInterceptSolverCallback : public InterceptSolverCallback { public: using InterceptSolverCallback::InterceptSolverCallback; MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(MatrixType, InterceptSolverCallback, getMatrixType, matrix_id); } void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleLumpedMatrix, matrix_id); } void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleResidual); } void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, predictor); } void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, corrector); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, beforeSolveStep); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, afterSolveStep, converged); } }; } // namespace /* -------------------------------------------------------------------------- */ void register_dof_manager(py::module & mod) { py::class_<DOFManager, std::shared_ptr<DOFManager>>(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( "hasMatrix", [](DOFManager & self, const ID & name) -> bool { return self.hasMatrix(name); }, py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("scale_factor") = 1.) .def("assembleToLumpedMatrix", &DOFManager::assembleToLumpedMatrix, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("lumped_mtx"), py::arg("scale_factor") = 1.) .def("assemblePreassembledMatrix", &DOFManager::assemblePreassembledMatrix, py::arg("matrix_id"), py::arg("terms")) .def("zeroResidual", &DOFManager::zeroResidual) - .def("localToGlobalEquationNumber", - &DOFManager::localToGlobalEquationNumber, - py::arg("local")) - .def("globalToLocalEquationNumber", - &DOFManager::globalToLocalEquationNumber, - py::arg("global")) - .def("getLocalEquationsNumbers", - [](DOFManager& self, const ID dof_id) -> decltype(auto) { - return self.getLocalEquationsNumbers(dof_id); }, - py::arg("dof_id"), - py::return_value_policy::reference_internal ) ; py::class_<NonLinearSolver>(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_<SolverCallback, PySolverCallback>(mod, "SolverCallback") .def(py::init_alias<DOFManager &>()) .def("getMatrixType", &SolverCallback::getMatrixType) .def("assembleMatrix", &SolverCallback::assembleMatrix) .def("assembleLumpedMatrix", &SolverCallback::assembleLumpedMatrix) .def("assembleResidual", [](SolverCallback & self) { self.assembleResidual(); }) .def("predictor", &SolverCallback::predictor) .def("corrector", &SolverCallback::corrector) .def("beforeSolveStep", &SolverCallback::beforeSolveStep) .def("afterSolveStep", &SolverCallback::afterSolveStep) .def_property_readonly("dof_manager", &SolverCallback::getSCDOFManager, py::return_value_policy::reference); py::class_<InterceptSolverCallback, SolverCallback, PyInterceptSolverCallback>(mod, "InterceptSolverCallback") .def(py::init_alias<SolverCallback &>()); } } // namespace akantu diff --git a/python/py_solver.cc b/python/py_solver.cc index f0bcc1483..1f16c3acb 100644 --- a/python/py_solver.cc +++ b/python/py_solver.cc @@ -1,118 +1,111 @@ /** * @file py_solver.cc * * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Tue Sep 29 2020 * @date last modification: Sat Mar 06 2021 * * @brief pybind11 interface to Solver and SparseMatrix * * * @section LICENSE * * Copyright (©) 2018-2021 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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "py_solver.hh" #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include <model.hh> #include <non_linear_solver.hh> #include <sparse_matrix_aij.hh> #include <terms_to_assemble.hh> /* -------------------------------------------------------------------------- */ #include <pybind11/operators.h> #include <pybind11/pybind11.h> #include <pybind11/stl.h> /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_solvers(py::module & mod) { py::class_<SparseMatrix>(mod, "SparseMatrix") .def("getMatrixType", &SparseMatrix::getMatrixType) .def("size", &SparseMatrix::size) .def("zero", &SparseMatrix::zero) .def("saveProfile", &SparseMatrix::saveProfile) .def("saveMatrix", &SparseMatrix::saveMatrix) .def( "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); }, "Add entry in the profile") .def( "add", [](SparseMatrix & self, UInt i, UInt j, Real value) { self.add(i, j, value); }, "Add the value to the matrix") .def( "add", [](SparseMatrix & self, SparseMatrix & A, Real alpha) { self.add(A, alpha); }, "Add a matrix to the matrix", py::arg("A"), py::arg("alpha") = 1.) .def("isFinite", &SparseMatrix::isFinite) .def("getRelease", [](const SparseMatrix& self) -> UInt { return self.getRelease(); }) .def("__call__", [](const SparseMatrix & self, UInt i, UInt j) { return self(i, j); }) .def("getRelease", &SparseMatrix::getRelease); py::class_<SparseMatrixAIJ, SparseMatrix>(mod, "SparseMatrixAIJ") .def("getIRN", &SparseMatrixAIJ::getIRN) .def("getJCN", &SparseMatrixAIJ::getJCN) .def("getA", &SparseMatrixAIJ::getA); - py::class_<SolverVector>(mod, "SolverVector") - .def("getValues", - [](SolverVector& self) -> decltype(auto) { - return static_cast<const Array<Real>& >(self); - }, - py::return_value_policy::reference_internal, - "Transform this into a vector, Is not copied.") - ; + py::class_<SolverVector>(mod, "SolverVector"); py::class_<TermsToAssemble::TermToAssemble>(mod, "TermToAssemble") .def(py::init<Int, Int>()) .def(py::self += Real()) .def_property_readonly("i", &TermsToAssemble::TermToAssemble::i) .def_property_readonly("j", &TermsToAssemble::TermToAssemble::j); py::class_<TermsToAssemble>(mod, "TermsToAssemble") .def(py::init<const ID &, const ID &>()) .def("getDOFIdM", &TermsToAssemble::getDOFIdM) .def("getDOFIdN", &TermsToAssemble::getDOFIdN) .def( "__call__", [](TermsToAssemble & self, UInt i, UInt j, Real val) { auto & term = self(i, j); term = val; return term; }, py::arg("i"), py::arg("j"), py::arg("val") = 0., py::return_value_policy::reference); } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index ff933e3b1..4b5b6fd0e 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,406 +1,404 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras <fabian.barras@epfl.ch> * @author Lucas Frerot <lucas.frerot@epfl.ch> * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch> * @author Philip Mueller <philip.paul.mueller@bluemail.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Damien Spielmann <damien.spielmann@epfl.ch> * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Apr 01 2021 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * * @section LICENSE * * Copyright (©) 2010-2021 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 <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #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 <ElementKind kind, class IntegrationOrderFunctor> class IntegratorGauss; template <ElementKind kind> 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<IntegratorGauss, ShapeStructural, _ek_structural>; StructuralMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "structural_mechanics_model"); ~StructuralMechanicsModel() override; /// 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 & matrix_id) const override; /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback to assemble the residual (rhs) void assembleResidual() override; void assembleResidual(const ID & residual_part) override; bool canSplitResidual() const override { return true; } void afterSolveStep(bool converged) override; /// compute kinetic energy Real getKineticEnergy(); /// compute potential energy Real getPotentialEnergy(); /// compute the specified energy Real getEnergy(const ID & energy); - /// Informs the model to reassemble the matrices - void needToReassembleMatrices(); /** * \brief This function computes the an approximation of the lumped mass. * * The mass is computed by looping over all beams and computing their mass. * The mass of a single beam is computed by the (initial) length of the beam, * its cross sectional area and its density. * The beam mass is then equaly distributed among the two nodes. * * For computing the rotational inertia, the function assumes that the mass of * a node is uniformaly distributed inside a disc (2D) or a sphere (3D). The * size of that disc, depends on the volume of the beam. * * Note that the computation of the mass is not unambigius. * The reason for this is, that the units of `StructralMaterial::rho` are not * clear. By default the function assumes that its unit are 'Mass per Volume'. * However, this makes the computed mass different than the consistent mass, * which seams to assume that its units are 'mass per unit length'. * The main difference between thge two are not the values, but that the * first version depends on `StructuralMaterial::A` while the later does not. * By defining the macro `AKANTU_STRUCTURAL_MECHANICS_CONSISTENT_LUMPED_MASS` * the function will compute the mass in a way that is consistent with the * consistent mass matrix. * * \note The lumped mass is not stored inside the DOFManager. * * \param ghost_type Should ghost types be computed. */ void assembleLumpedMassMatrix(); /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple<ID, TimeStepSolverType> getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; static UInt getNbDegreeOfFreedom(ElementType type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// compute the nodal forces void assembleInternalForce(); /// compute the nodal forces for an element type void assembleInternalForce(ElementType type, GhostType gt); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMassMatrix(); protected: /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMassMatrix(GhostType ghost_type); /// computes rho void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template <ElementType type> void assembleStiffnessMatrix(); template <ElementType type> void computeStressOnQuad(); template <ElementType type> void computeTangentModuli(Array<Real> & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr<dumpers::Field> createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr<dumpers::Field> createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr<dumpers::Field> createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array<Real> &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &); /** * Returns a non const reference to the array that stores the lumped * mass. * * The returned array has dimension `N x d` where `N` is the number of nodes * and `d`, is the number of degrees of freedom per node. */ inline Array<Real> & getLumpedMass() { if (this->mass == nullptr) { AKANTU_EXCEPTION("The pointer to the mass was not allocated."); }; return *(this->mass); }; /** * Returns a constant reference to the array that stores the lumped * mass. */ inline const Array<Real> & getLumpedMass() const { if (this->mass == nullptr) { AKANTU_EXCEPTION("The pointer to the mass was not allocated."); }; return *(this->mass); }; /** * Tests if *this has a lumped mass pointer. */ inline bool hasLumpedMass() const { return (this->mass != nullptr); }; 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); /** * \brief This function adds the `StructuralMaterial` material to the list of * materials managed by *this. * * It is important that this function might invalidate all references to * structural materials, that were previously optained by `getMaterial()`. * * \param material The new material. * * \return The ID of the material that was added. * * \note The return type is is new. */ UInt addMaterial(StructuralMaterial & material, const ID & name = ""); const StructuralMaterial & getMaterialByElement(const Element & element) const; /** * \brief Returns the ith material of *this. * \param i The ith material */ const StructuralMaterial & getMaterial(UInt material_index) const; const StructuralMaterial & getMaterial(const ID & name) const; /** * \brief Returns the number of the different materials inside *this. */ UInt getNbMaterials() const { return materials.size(); } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis void computeForcesByGlobalTractionArray(const Array<Real> & traction_global, ElementType type); /// Compute Linear load function set in local axis void computeForcesByLocalTractionArray(const Array<Real> & tractions, ElementType type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array std::unique_ptr<Array<Real>> displacement_rotation; /// velocities array std::unique_ptr<Array<Real>> velocity; /// accelerations array std::unique_ptr<Array<Real>> acceleration; /// forces array std::unique_ptr<Array<Real>> internal_force; /// forces array std::unique_ptr<Array<Real>> external_force; /** * \brief This is the "lumped" mass array. * * It is a bit special, since it is not a one dimensional array, bit it is * actually a matrix. The number of rows equals the number of nodes. The * number of colums equals the number of degrees of freedoms per node. This * layout makes the thing a bit more simple. * * Note that it is only allocated in case, the "Lumped" mode is enabled. */ std::unique_ptr<Array<Real>> mass; /// boundaries array std::unique_ptr<Array<bool>> blocked_dofs; /// stress array ElementTypeMapArray<Real> stress; ElementTypeMapArray<UInt> element_material; // Define sets of beams ElementTypeMapArray<UInt> set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray<Real> 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; bool need_to_reassemble_mass{true}; bool need_to_reassemble_stiffness{true}; /* ------------------------------------------------------------------------ */ std::vector<StructuralMaterial> materials; std::map<std::string, UInt> materials_names_to_id; }; } // namespace akantu #include "structural_mechanics_model_inline_impl.hh" #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ */