diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc
index a9cb7c2fa..1a073da25 100644
--- a/python/py_dof_manager.cc
+++ b/python/py_dof_manager.cc
@@ -1,217 +1,216 @@
/*
* 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 .
*/
/* -------------------------------------------------------------------------- */
#include "py_dof_manager.hh"
#include "py_aka_array.hh"
#include "py_akantu_pybind11_compatibility.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
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_>(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("zeroResidual", &DOFManager::zeroResidual);
py::class_(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_(mod, "TimeStepSolver")
.def("getIntegrationScheme", &TimeStepSolver::getIntegrationScheme);
py::class_(mod, "SolverCallback")
.def(py::init_alias())
.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_(mod, "InterceptSolverCallback")
.def(py::init_alias());
}
} // namespace akantu
diff --git a/python/py_solver.cc b/python/py_solver.cc
index c1239854a..d5e20670d 100644
--- a/python/py_solver.cc
+++ b/python/py_solver.cc
@@ -1,119 +1,119 @@
/**
* @file py_solver.cc
*
* @author Nicolas Richart
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "py_solver.hh"
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
void register_solvers(py::module & mod) {
py::class_(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(); })
+ [](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_(mod, "SparseMatrixAIJ")
.def("getIRN", &SparseMatrixAIJ::getIRN)
.def("getJCN", &SparseMatrixAIJ::getJCN)
.def("getA", &SparseMatrixAIJ::getA);
py::class_(mod, "SolverVector")
- .def("getValues",
- [](SolverVector& self) -> decltype(auto) {
- return static_cast& >(self);
- },
- py::return_value_policy::reference_internal,
- "Transform this into a vector, Is not copied.")
- .def("isDistributed", [](const SolverVector& self) { return self.isDistributed(); })
- ;
+ .def(
+ "getValues",
+ [](SolverVector & self) -> decltype(auto) {
+ return static_cast &>(self);
+ },
+ py::return_value_policy::reference_internal,
+ "Transform this into a vector, Is not copied.")
+ .def("isDistributed",
+ [](const SolverVector & self) { return self.isDistributed(); });
py::class_(mod, "TermToAssemble")
.def(py::init())
.def(py::self += Real())
.def_property_readonly("i", &TermsToAssemble::TermToAssemble::i)
.def_property_readonly("j", &TermsToAssemble::TermToAssemble::j);
py::class_(mod, "TermsToAssemble")
.def(py::init())
.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/python/py_structural_mechanics_model.cc b/python/py_structural_mechanics_model.cc
index 79801516f..366953881 100644
--- a/python/py_structural_mechanics_model.cc
+++ b/python/py_structural_mechanics_model.cc
@@ -1,177 +1,179 @@
/**
* @file py_structural_mechanics_model.cc
*
* @author Philip Mueller
* @author Mohit Pundir
* @author Nicolas Richart
*
* @date creation: Wed Feb 03 2021
* @date last modification: Thu Apr 01 2021
*
* @brief pybind11 interface to StructuralMechanicsModel
*
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
#define def_deprecated(func_name, mesg) \
def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); })
#define def_function_nocopy(func_name) \
def( \
#func_name, \
[](StructuralMechanicsModel & self) -> decltype(auto) { \
return self.func_name(); \
}, \
py::return_value_policy::reference)
#define def_function_(func_name) \
def(#func_name, [](StructuralMechanicsModel & self) -> decltype(auto) { \
return self.func_name(); \
})
#define def_plainmember(M) def_readwrite(#M, &StructuralMaterial::M)
/* -------------------------------------------------------------------------- */
void register_structural_mechanics_model(pybind11::module & mod) {
/* First we have to register the material class
* The wrapper aims to mimic the behaviour of the real material.
*/
py::class_(mod, "StructuralMaterial")
.def(py::init<>())
.def(py::init())
.def_plainmember(E)
.def_plainmember(A)
.def_plainmember(I)
.def_plainmember(Iz)
.def_plainmember(Iy)
.def_plainmember(GJ)
.def_plainmember(rho)
.def_plainmember(t)
.def_plainmember(nu);
/* Now we create the structural model wrapper
* Note that this is basically a port from the solid mechanic part.
*/
py::class_(mod, "StructuralMechanicsModel")
.def(py::init(), py::arg("mesh"),
py::arg("spatial_dimension") = _all_dimensions,
py::arg("id") = "structural_mechanics_model")
.def(
"initFull",
[](StructuralMechanicsModel & self,
const AnalysisMethod & analysis_method) -> void {
self.initFull(_analysis_method = analysis_method);
},
py::arg("_analysis_method"))
.def("initFull",
[](StructuralMechanicsModel & self) -> void { self.initFull(); })
.def_function_nocopy(getExternalForce)
.def_function_nocopy(getDisplacement)
.def_function_nocopy(getInternalForce)
.def_function_nocopy(getVelocity)
.def_function_nocopy(getAcceleration)
.def_function_nocopy(getInternalForce)
.def_function_nocopy(getBlockedDOFs)
.def_function_nocopy(getMesh)
.def("setTimeStep", &StructuralMechanicsModel::setTimeStep,
py::arg("time_step"), py::arg("solver_id") = "")
.def(
"getElementMaterial",
[](StructuralMechanicsModel & self, const ElementType & type,
GhostType ghost_type) -> decltype(auto) {
return self.getElementMaterial(type, ghost_type);
},
"This function returns the map that maps elements to materials.",
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getMaterialByElement",
[](StructuralMechanicsModel & self, Element element)
-> decltype(auto) { return self.getMaterialByElement(element); },
"This function returns the `StructuralMaterial` instance that is "
"associated with element `element`.",
py::arg("element"), py::return_value_policy::reference)
.def(
"addMaterial",
[](StructuralMechanicsModel & self, StructuralMaterial & mat,
const ID & name) -> UInt { return self.addMaterial(mat, name); },
"This function adds the `StructuralMaterial` `mat` to `self`."
" The function returns the ID of the new material.",
py::arg("mat"), py::arg("name") = "")
.def(
"getMaterial",
- [](const StructuralMechanicsModel & self, UInt material_index)
- -> StructuralMaterial { return self.getMaterial(material_index); },
+ [](const StructuralMechanicsModel & self,
+ UInt material_index) -> StructuralMaterial {
+ return self.getMaterial(material_index);
+ },
"This function returns the `i`th material of `self`."
" The function returns a copy of the material.",
py::arg("material_index"), py::return_value_policy::copy)
.def(
"getMaterial",
[](const StructuralMechanicsModel & self, const ID & name)
-> StructuralMaterial { return self.getMaterial(name); },
"This function returns the `i`th material of `self`."
" The function returns a copy of the material.",
py::arg("material_index"), py::return_value_policy::copy)
.def(
"getNbMaterials",
[](StructuralMechanicsModel & self) { return self.getNbMaterials(); },
"Returns the number of different materials inside `self`.")
.def("getKineticEnergy", &StructuralMechanicsModel::getKineticEnergy,
"Compute kinetic energy")
.def("getPotentialEnergy", &StructuralMechanicsModel::getPotentialEnergy,
"Compute potential energy")
.def("getEnergy", &StructuralMechanicsModel::getEnergy,
"Compute the specified energy")
.def(
"getLumpedMass",
[](const StructuralMechanicsModel & self) -> decltype(auto) {
return self.getLumpedMass();
},
py::return_value_policy::reference_internal)
.def(
"getMass",
[](const StructuralMechanicsModel & self) -> decltype(auto) {
return self.getMass();
},
py::return_value_policy::reference_internal)
.def("assembleLumpedMassMatrix",
&StructuralMechanicsModel::assembleLumpedMassMatrix,
"Assembles the lumped mass matrix")
.def("hasLumpedMass", &StructuralMechanicsModel::hasLumpedMass);
}
} // namespace akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index 6c01e0cb1..04e05dc5e 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,777 +1,775 @@
/**
* @file structural_mechanics_model.cc
*
* @author Fabian Barras
* @author Lucas Frerot
* @author Sébastien Hartmann
* @author Nicolas Richart
* @author Damien Spielmann
*
* @date creation: Fri Jul 15 2011
* @date last modification: Mon Mar 15 2021
*
* @brief Model implementation for Structural Mechanics elements
*
*
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "structural_mechanics_model.hh"
#include "dof_manager.hh"
#include "integrator_gauss.hh"
#include "mesh.hh"
#include "shape_structural.hh"
#include "sparse_matrix.hh"
#include "time_step_solver.hh"
/* -------------------------------------------------------------------------- */
#include "dumpable_inline_impl.hh"
#include "dumper_elemental_field.hh"
#include "dumper_internal_material_field.hh"
#include "dumper_iohelper_paraview.hh"
#include "group_manager_inline_impl.hh"
/* -------------------------------------------------------------------------- */
#include "structural_element_bernoulli_beam_2.hh"
#include "structural_element_bernoulli_beam_3.hh"
#include "structural_element_kirchhoff_shell.hh"
/* -------------------------------------------------------------------------- */
//#include "structural_mechanics_model_inline_impl.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) {
UInt ndof = 0;
#define GET_(type) ndof = ElementClass::getNbDegreeOfFreedom()
AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural);
#undef GET_
return ndof;
}
/* -------------------------------------------------------------------------- */
StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim,
const ID & id)
: Model(mesh, ModelType::_structural_mechanics_model, dim, id), f_m2a(1.0),
stress("stress", id), element_material("element_material", id),
set_ID("beam sets", 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_TO_IMPLEMENT();
}
this->mesh.registerDumper("structural_mechanics_model", id,
true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural);
this->initDOFManager();
this->dumper_default_element_kind = _ek_structural;
mesh.getElementalData("extra_normal")
.initialize(mesh, _element_kind = _ek_structural,
_nb_component = spatial_dimension, _with_nb_element = true,
_default_value = 0.);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
StructuralMechanicsModel::~StructuralMechanicsModel() = default;
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) {
Model::initFullImpl(options);
// Initializing stresses
ElementTypeMap stress_components;
/// TODO this is ugly af, maybe add a function to FEEngine
for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
_element_kind = _ek_structural)) {
UInt nb_components = 0;
// Getting number of components for each element type
#define GET_(type) nb_components = ElementClass::getNbStressComponents()
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_);
#undef GET_
stress_components(nb_components, type);
}
stress.initialize(
getFEEngine(), _spatial_dimension = _all_dimensions,
_element_kind = _ek_structural,
_nb_component = [&stress_components](ElementType type,
GhostType /*unused*/) -> UInt {
return stress_components(type);
});
}
/* -------------------------------------------------------------------------- */
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,
const ID & solver_id) {
Model::setTimeStep(time_step, solver_id);
this->mesh.getDumper().setTimeStep(time_step);
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initSolver(
TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) {
AKANTU_DEBUG_IN();
this->allocNodalField(displacement_rotation, nb_degree_of_freedom,
"displacement");
this->allocNodalField(external_force, nb_degree_of_freedom, "external_force");
this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force");
this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs");
auto & dof_manager = this->getDOFManager();
if (not dof_manager.hasDOFs("displacement")) {
dof_manager.registerDOFs("displacement", *displacement_rotation,
_dst_nodal);
dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
}
if (time_step_solver_type == TimeStepSolverType::_dynamic ||
time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
this->allocNodalField(velocity, nb_degree_of_freedom, "velocity");
this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration");
if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
dof_manager.registerDOFsDerivative("displacement", 2,
*this->acceleration);
}
/* Only allocate the mass if the "lumped" mode is ennabled.
* Also it is not a 1D array, but has an element for every
* DOF, which are most of the time equal, but makes handling
* some operations a bit simpler. */
if (time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
this->allocateLumpedMassArray();
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::initModel() {
element_material.initialize(mesh, _element_kind = _ek_structural,
_default_value = 0, _with_nb_element = true);
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
if (not need_to_reassemble_stiffness) {
return;
}
if (not getDOFManager().hasMatrix("K")) {
getDOFManager().getNewMatrix("K", getMatrixType("K"));
}
this->getDOFManager().zeroMatrix("K");
for (const 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
}
need_to_reassemble_stiffness = false;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::computeStresses() {
AKANTU_DEBUG_IN();
for (const 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();
}
/* -------------------------------------------------------------------------- */
bool StructuralMechanicsModel::allocateLumpedMassArray() {
if (this->mass != nullptr) // Already allocated, so nothing to do.
{
return true;
};
// now allocate it
this->allocNodalField(this->mass, this->nb_degree_of_freedom, "lumped_mass");
return true;
};
/* -------------------------------------------------------------------------- */
std::shared_ptr 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.get();
return mesh.createNodalField(uint_nodal_fields[field_name], group_name);
}
/* -------------------------------------------------------------------------- */
std::shared_ptr
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;
}
UInt padding_size = 0;
if (padding_flag) {
padding_size = 3;
}
if (field_name == "displacement") {
return mesh.createStridedNodalField(displacement_rotation.get(), group_name,
n, 0, padding_size);
}
if (field_name == "velocity") {
return mesh.createStridedNodalField(velocity.get(), group_name, n, 0,
padding_size);
}
if (field_name == "acceleration") {
return mesh.createStridedNodalField(acceleration.get(), group_name, n, 0,
padding_size);
}
if (field_name == "rotation") {
return mesh.createStridedNodalField(displacement_rotation.get(), group_name,
nb_degree_of_freedom - n, n,
padding_size);
}
if (field_name == "force") {
return mesh.createStridedNodalField(external_force.get(), group_name, n, 0,
padding_size);
}
if (field_name == "external_force") {
return mesh.createStridedNodalField(external_force.get(), group_name, n, 0,
padding_size);
}
if (field_name == "momentum") {
return mesh.createStridedNodalField(external_force.get(), group_name,
nb_degree_of_freedom - n, n,
padding_size);
}
if (field_name == "internal_force") {
return mesh.createStridedNodalField(internal_force.get(), group_name, n, 0,
padding_size);
}
if (field_name == "internal_momentum") {
return mesh.createStridedNodalField(internal_force.get(), group_name,
nb_degree_of_freedom - n, n,
padding_size);
}
if (field_name == "mass") {
AKANTU_DEBUG_ASSERT(this->mass.get() != nullptr,
"The lumped mass matrix was not allocated.");
return mesh.createStridedNodalField(this->mass.get(), group_name, n, 0,
padding_size);
}
return nullptr;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr StructuralMechanicsModel::createElementalField(
const std::string & field_name, const std::string & group_name,
bool /*unused*/, UInt spatial_dimension, ElementKind kind) {
std::shared_ptr field;
if (field_name == "element_index_by_material") {
field = mesh.createElementalField(
field_name, group_name, spatial_dimension, kind);
}
if (field_name == "stress") {
ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(stress);
field = mesh.createElementalField(
stress, group_name, this->spatial_dimension, kind, nb_data_per_elem);
}
return field;
}
/* -------------------------------------------------------------------------- */
/* Virtual methods from SolverCallback */
/* -------------------------------------------------------------------------- */
/// get the type of matrix needed
MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) const {
return _symmetric;
}
/// callback to assemble a Matrix
void StructuralMechanicsModel::assembleMatrix(const ID & id) {
if (id == "K") {
assembleStiffnessMatrix();
} else if (id == "M") {
assembleMassMatrix();
}
}
/// callback to assemble a lumped Matrix
void StructuralMechanicsModel::assembleLumpedMatrix(const ID & id) {
if ("M" == id) {
this->assembleLumpedMassMatrix();
}
return;
}
/// callback to assemble the residual StructuralMechanicsModel::(rhs)
void StructuralMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
auto & dof_manager = getDOFManager();
assembleInternalForce();
// Ensures that the matrix are assembled.
if (dof_manager.hasMatrix("K")) {
this->assembleMatrix("K");
}
if (dof_manager.hasMatrix("M")) {
this->assembleMatrix("M");
}
if (dof_manager.hasLumpedMatrix("M")) {
this->assembleLumpedMassMatrix();
}
/* This is essentially a summing up of forces
* first the external forces are counted for and then stored inside the
* residual.
*/
dof_manager.assembleToResidual("displacement", *external_force, 1);
dof_manager.assembleToResidual("displacement", *internal_force, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
auto & dof_manager = this->getDOFManager();
if ("external" == residual_part) {
- dof_manager.assembleToResidual("displacement",
- *this->external_force, 1);
+ dof_manager.assembleToResidual("displacement", *this->external_force, 1);
AKANTU_DEBUG_OUT();
return;
}
if ("internal" == residual_part) {
this->assembleInternalForce();
- dof_manager.assembleToResidual("displacement",
- *this->internal_force, 1);
+ dof_manager.assembleToResidual("displacement", *this->internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
AKANTU_CUSTOM_EXCEPTION(
debug::SolverCallbackResidualPartUnknown(residual_part));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* 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", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
case _explicit_lumped_mass: { // Taken from the solid mechanics part
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: { // Taken from the solid mechanics part
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
default:
std::cout << "UNKOWN." << std::endl;
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* ------------------------------------------------------------------------ */
ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: { // Taken from the solid mechanic
// part
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type =
NonLinearSolverType::_newton_raphson; // _linear;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
#if 1
case TimeStepSolverType::_dynamic: { // Copied from solid
if (this->method == _explicit_consistent_mass) {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
} else {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_trapezoidal_rule_2;
options.solution_type["displacement"] = IntegrationScheme::_displacement;
}
break;
}
#else
case TimeStepSolverType::_dynamic: { // Original
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_trapezoidal_rule_2;
options.solution_type["displacement"] = IntegrationScheme::_displacement;
break;
}
#endif
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
}
return options;
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::assembleInternalForce() {
internal_force->zero();
computeStresses();
for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
_element_kind = _ek_structural)) {
assembleInternalForce(type, _not_ghost);
// assembleInternalForce(type, _ghost);
}
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::assembleInternalForce(ElementType type,
GhostType gt) {
auto & fem = getFEEngine();
auto & sigma = stress(type, gt);
auto ndof = getNbDegreeOfFreedom(type);
auto nb_nodes = mesh.getNbNodesPerElement(type);
auto ndof_per_elem = ndof * nb_nodes;
Array BtSigma(fem.getNbIntegrationPoints(type) *
mesh.getNbElement(type),
ndof_per_elem, "BtSigma");
fem.computeBtD(sigma, BtSigma, type, gt);
Array intBtSigma(0, ndof_per_elem, "intBtSigma");
fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt);
getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force,
type, gt, -1.);
}
/* -------------------------------------------------------------------------- */
Real StructuralMechanicsModel::getKineticEnergy() {
const UInt nb_nodes = mesh.getNbNodes();
const UInt nb_degree_of_freedom = this->nb_degree_of_freedom;
Real ekin = 0.; // used to sum up energy (is divided by two at the very end)
- //if mass matrix was not assembled, assemble it now
+ // if mass matrix was not assembled, assemble it now
this->assembleMassMatrix();
if (this->getDOFManager().hasLumpedMatrix("M")) {
/* This code computes the kinetic energy for the case when the mass is
* lumped. It is based on the solid mechanic equivalent.
*/
AKANTU_DEBUG_ASSERT(this->mass != nullptr,
"The lumped mass is not allocated.");
if (this->need_to_reassemble_lumpedMass) {
this->assembleLumpedMatrix("M");
}
/* Iterating over all nodes.
* Important the velocity and mass also contains the rotational parts.
* However, they can be handled in an uniform way. */
for (auto && data :
zip(arange(nb_nodes), make_view(*this->velocity, nb_degree_of_freedom),
make_view(*this->mass, nb_degree_of_freedom))) {
const UInt n = std::get<0>(data); // This is the ID of the current node
if (not mesh.isLocalOrMasterNode(
n)) // Only handle the node if it belongs to us.
{
continue;
}
const auto & v =
std::get<1>(data); // Get the velocity and mass of that node.
const auto & m = std::get<2>(data);
Real mv2 = 0.; // Contribution of this node.
for (UInt i = 0; i < nb_degree_of_freedom; ++i) {
/* In the solid mechanics part, only masses that are above a certain
* value are considered.
* However, the structural part, does not do this. */
const Real v_ = v(i);
const Real m_ = m(i);
mv2 += v_ * v_ * m_;
} // end for(i): going through the components
ekin += mv2; // add continution
} // end for(n): iterating through all nodes
} else if (this->getDOFManager().hasMatrix("M")) {
/* Handle the case where no lumped mass is there.
* This is basically the original code.
*/
if (this->need_to_reassemble_mass) {
this->assembleMassMatrix();
}
Array Mv(nb_nodes, nb_degree_of_freedom);
this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
*this->velocity, Mv);
for (auto && data :
zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom),
make_view(*this->velocity, nb_degree_of_freedom))) {
if (mesh.isLocalOrMasterNode(std::get<0>(
data))) // only consider the node if we are blonging to it
{
ekin += std::get<2>(data).dot(std::get<1>(data));
}
}
} else {
/* This is the case where no mass is present, for whatever reason, such as
* the static case. We handle it specially be returning directly zero.
* However, by doing that there will not be a syncronizing event as in the
* other cases. Which is faster, but could be a problem in case the user
* expects this.
*
* Another not is, that the solid mechanics part, would generate an error in
* this clause. But, since the original implementation of the structural
* part, did not do that, I, Philip, decided to refrain from that. However,
* it is an option that should be considered.
*/
return 0.;
}
// Sum up across the comunicator
mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
return ekin / 2.; // finally divide the energy by two
}
/* -------------------------------------------------------------------------- */
Real StructuralMechanicsModel::getPotentialEnergy() {
Real epot = 0.;
UInt nb_nodes = mesh.getNbNodes();
- //if stiffness matrix is not assembled, do it
+ // if stiffness matrix is not assembled, do it
this->assembleStiffnessMatrix();
Array Ku(nb_nodes, nb_degree_of_freedom);
this->getDOFManager().assembleMatMulVectToArray(
"displacement", "K", *this->displacement_rotation, Ku);
for (auto && data :
zip(arange(nb_nodes), make_view(Ku, nb_degree_of_freedom),
make_view(*this->displacement_rotation, nb_degree_of_freedom))) {
epot += std::get<2>(data).dot(std::get<1>(data)) *
static_cast(mesh.isLocalOrMasterNode(std::get<0>(data)));
}
mesh.getCommunicator().allReduce(epot, SynchronizerOperation::_sum);
return epot / 2.;
}
/* -------------------------------------------------------------------------- */
Real StructuralMechanicsModel::getEnergy(const ID & energy) {
if (energy == "kinetic") {
return getKineticEnergy();
}
if (energy == "potential") {
return getPotentialEnergy();
}
return 0;
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::computeForcesByLocalTractionArray(
const Array & tractions, ElementType type) {
AKANTU_DEBUG_IN();
auto nb_element = getFEEngine().getMesh().getNbElement(type);
auto nb_nodes_per_element =
getFEEngine().getMesh().getNbNodesPerElement(type);
auto nb_quad = getFEEngine().getNbIntegrationPoints(type);
// check dimension match
AKANTU_DEBUG_ASSERT(
Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(),
"element type dimension does not match the dimension of boundaries : "
<< getFEEngine().getElementDimension()
<< " != " << Mesh::getSpatialDimension(type));
// check size of the vector
AKANTU_DEBUG_ASSERT(
tractions.size() == nb_quad * nb_element,
"the size of the vector should be the total number of quadrature points");
// check number of components
AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom,
"the number of components should be the spatial "
"dimension of the problem");
Array Ntbs(nb_element * nb_quad,
nb_degree_of_freedom * nb_nodes_per_element);
auto & fem = getFEEngine();
fem.computeNtb(tractions, Ntbs, type);
// allocate the vector that will contain the integrated values
auto name = id + std::to_string(type) + ":integral_boundary";
Array int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element,
name);
// do the integration
getFEEngine().integrate(Ntbs, int_funct,
nb_degree_of_freedom * nb_nodes_per_element, type);
// assemble the result into force vector
getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force,
type, _not_ghost, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::computeForcesByGlobalTractionArray(
const Array & traction_global, ElementType type) {
AKANTU_DEBUG_IN();
UInt nb_element = mesh.getNbElement(type);
UInt nb_quad = getFEEngine().getNbIntegrationPoints(type);
Array traction_local(nb_element * nb_quad, nb_degree_of_freedom,
id + ":structuralmechanics:imposed_linear_load");
auto R_it = getFEEngineClass()
.getShapeFunctions()
.getRotations(type)
.begin(nb_degree_of_freedom, nb_degree_of_freedom);
auto Te_it = traction_global.begin(nb_degree_of_freedom);
auto te_it = traction_local.begin(nb_degree_of_freedom);
for (UInt e = 0; e < nb_element; ++e, ++R_it) {
for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) {
// turn the traction in the local referential
te_it->template mul(*R_it, *Te_it);
}
}
computeForcesByLocalTractionArray(traction_local, type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StructuralMechanicsModel::afterSolveStep(bool converged) {
if (converged) {
assembleInternalForce();
}
}
} // namespace akantu