diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc
index 074a3907b..6633335c7 100644
--- a/python/py_dof_manager.cc
+++ b/python/py_dof_manager.cc
@@ -1,218 +1,229 @@
/**
* Copyright (©) 2021-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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 <time_step_solver.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);
}
};
+ class PyNonLinearSolver : public NonLinearSolver {
+ public:
+ using NonLinearSolver::NonLinearSolver;
+
+ void solve(SolverCallback & callback) override {
+ // NOLINTNEXTLINE
+ PYBIND11_OVERRIDE(void, NonLinearSolver, solve, callback);
+ }
+ };
} // 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);
- py::class_<NonLinearSolver, Parsable>(mod, "NonLinearSolver")
+ py::class_<NonLinearSolver, Parsable, PyNonLinearSolver>(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); });
+ const SolveConvergenceCriteria & val) { self.set(id, val); })
+ .def("solve", &NonLinearSolver::solve);
py::class_<TimeStepSolver>(mod, "TimeStepSolver")
.def("getIntegrationScheme", &TimeStepSolver::getIntegrationScheme);
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/src/common/aka_math.hh b/src/common/aka_math.hh
index b95391d82..3858368d1 100644
--- a/src/common/aka_math.hh
+++ b/src/common/aka_math.hh
@@ -1,150 +1,150 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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_types.hh"
/* -------------------------------------------------------------------------- */
#include <utility>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_AKA_MATH_H_
#define AKANTU_AKA_MATH_H_
namespace akantu {
/* -------------------------------------------------------------------------- */
namespace AKANTU_EXPORT Math {
/// tolerance for functions that need one
extern Real tolerance; // NOLINT
/* ------------------------------------------------------------------------ */
/* Geometry */
/* ------------------------------------------------------------------------ */
/// compute normal a normal to a vector
template <class D1, std::enable_if_t<aka::is_vector_v<D1>> * = nullptr>
inline Vector<Real> normal(const Eigen::MatrixBase<D1> & vec);
template <class D1, aka::enable_if_t<not aka::is_vector_v<D1>> * = nullptr>
inline Vector<Real> normal(const Eigen::MatrixBase<D1> & /*vec*/) {
AKANTU_TO_IMPLEMENT();
}
/// compute normal a normal to a vector
template <class D1, class D2,
std::enable_if_t<aka::are_vectors<D1, D2>::value> * = nullptr>
inline Vector<Real, 3> normal(const Eigen::MatrixBase<D1> & vec1,
const Eigen::MatrixBase<D2> & vec2);
/// compute the tangents to an array of normal vectors
void compute_tangents(const Array<Real> & normals, Array<Real> & tangents);
/// radius of the in-circle of a triangle in 2d space
template <class D1, class D2, class D3>
static inline Real triangle_inradius(const Eigen::MatrixBase<D1> & coord1,
const Eigen::MatrixBase<D2> & coord2,
const Eigen::MatrixBase<D3> & coord3);
/// radius of the in-circle of a tetrahedron
template <class D1, class D2, class D3, class D4>
static inline Real tetrahedron_inradius(const Eigen::MatrixBase<D1> & coord1,
const Eigen::MatrixBase<D2> & coord2,
const Eigen::MatrixBase<D3> & coord3,
const Eigen::MatrixBase<D4> & coord4);
/// volume of a tetrahedron
template <class D1, class D2, class D3, class D4>
static inline Real tetrahedron_volume(const Eigen::MatrixBase<D1> & coord1,
const Eigen::MatrixBase<D2> & coord2,
const Eigen::MatrixBase<D3> & coord3,
const Eigen::MatrixBase<D4> & coord4);
/// compute the barycenter of n points
template <class D1, class D2>
inline void barycenter(const Eigen::MatrixBase<D1> & coord,
Eigen::MatrixBase<D2> & barycenter);
/// test if two scalar are equal within a given tolerance
inline bool are_float_equal(Real x, Real y);
/// test if two vectors are equal within a given tolerance
inline bool are_vector_equal(Int n, Real * x, Real * y);
#ifdef isnan
#error \
"You probably included <math.h> which is incompatible with aka_math please use\
<cmath> or add a \"#undef isnan\" before akantu includes"
#endif
/// test if a real is a NaN
inline bool isnan(Real x);
/// test if the line x and y intersects each other
inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max);
/// test if a is in the range [x_min, x_max]
inline bool is_in_range(Real a, Real x_min, Real x_max);
inline Real getTolerance() { return Math::tolerance; }
inline void setTolerance(Real tol) { Math::tolerance = tol; }
template <Int p, typename T> inline T pow(T x);
template <class T1, class T2,
std::enable_if_t<std::is_integral<T1>::value and
std::is_integral<T2>::value> * = nullptr>
inline Real kronecker(T1 i, T2 j) {
return static_cast<Real>(i == j);
}
/* --------------------------------------------------------------------------
*/
template <typename T> static inline constexpr T pow(T x, int p) {
return p == 0 ? T(1) : (pow(x, p - 1) * x);
}
/// reduce all the values of an array, the summation is done in place and the
/// array is modified
- Real reduce(Array<Real> & array);
+ AKANTU_EXPORT Real reduce(Array<Real> & array);
template <class T> class NewtonRaphson {
public:
NewtonRaphson(Real tolerance, Int max_iteration)
: tolerance(tolerance), max_iteration(max_iteration) {}
template <class Functor> T solve(const Functor & funct, const T & x_0);
private:
Real tolerance;
Int max_iteration;
};
template <class T> struct NewtonRaphsonFunctor {
explicit NewtonRaphsonFunctor(const std::string & name) : name(name) {}
virtual T f(const T & x) const = 0;
virtual T f_prime(const T & x) const = 0;
std::string name;
};
} // namespace Math
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "aka_math_tmpl.hh"
#endif /* AKANTU_AKA_MATH_H_ */
diff --git a/src/model/common/dof_manager/dof_manager_default.cc b/src/model/common/dof_manager/dof_manager_default.cc
index 231bb7dd7..965479386 100644
--- a/src/model/common/dof_manager/dof_manager_default.cc
+++ b/src/model/common/dof_manager/dof_manager_default.cc
@@ -1,471 +1,473 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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 "dof_manager_default.hh"
#include "communicator.hh"
#include "dof_synchronizer.hh"
#include "element_group.hh"
#include "non_linear_solver_default.hh"
#include "periodic_node_synchronizer.hh"
#include "solver_vector_default.hh"
#include "solver_vector_distributed.hh"
#include "sparse_matrix_aij.hh"
#include "time_step_solver_default.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <memory>
#include <numeric>
#include <unordered_map>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
DOFManagerDefault::DOFManagerDefault(const ID & id)
: DOFManager(id), synchronizer(nullptr) {
residual = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":residual"));
solution = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":solution"));
data_cache = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":data_cache"));
}
/* -------------------------------------------------------------------------- */
DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id)
: DOFManager(mesh, id), synchronizer(nullptr) {
if (this->mesh->isDistributed()) {
this->synchronizer = std::make_unique<DOFSynchronizer>(
*this, this->id + ":dof_synchronizer");
residual = std::make_unique<SolverVectorDistributed>(
*this, std::string(id + ":residual"));
solution = std::make_unique<SolverVectorDistributed>(
*this, std::string(id + ":solution"));
data_cache = std::make_unique<SolverVectorDistributed>(
*this, std::string(id + ":data_cache"));
} else {
residual = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":residual"));
solution = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":solution"));
data_cache = std::make_unique<SolverVectorDefault>(
*this, std::string(id + ":data_cache"));
}
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::makeConsistentForPeriodicity(const ID & dof_id,
SolverVector & array) {
auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
if (dof_data.support_type != _dst_nodal) {
return;
}
if (not mesh->isPeriodic()) {
return;
}
this->mesh->getPeriodicNodeSynchronizer()
.reduceSynchronizeWithPBCSlaves<AddOperation>(
aka::as_type<SolverVectorDefault>(array).getVector());
}
/* -------------------------------------------------------------------------- */
template <typename T>
void DOFManagerDefault::assembleToGlobalArray(
const ID & dof_id, const Array<T> & array_to_assemble,
Array<T> & global_array, T scale_factor) {
AKANTU_DEBUG_IN();
auto & dof_data = this->getDOFDataTyped<DOFDataDefault>(dof_id);
AKANTU_DEBUG_ASSERT(dof_data.local_equation_number.size() ==
array_to_assemble.size() *
array_to_assemble.getNbComponent(),
"The array to assemble does not have a correct size."
<< " (" << array_to_assemble.getID() << ")");
if (dof_data.support_type == _dst_nodal and mesh->isPeriodic()) {
for (auto && data :
zip(dof_data.local_equation_number, dof_data.associated_nodes,
make_view(array_to_assemble))) {
auto && equ_num = std::get<0>(data);
// auto && node = std::get<1>(data);
auto && arr = std::get<2>(data);
// Guillaume to Nico:
// This filter of periodic slave should not be.
// Indeed you want to get the contribution even
// from periodic slaves and cumulate to the right
// equation number.
global_array(equ_num) += scale_factor * (arr);
// scale_factor * (arr) * (not this->mesh->isPeriodicSlave(node));
}
} else {
for (auto && data :
zip(dof_data.local_equation_number, make_view(array_to_assemble))) {
auto && equ_num = std::get<0>(data);
auto && arr = std::get<1>(data);
global_array(equ_num) += scale_factor * (arr);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::assembleToGlobalArray(
const ID & dof_id, const Array<Real> & array_to_assemble,
SolverVector & global_array_v, Real scale_factor) {
assembleToGlobalArray(
dof_id, array_to_assemble,
aka::as_type<SolverVectorDefault>(global_array_v).getVector(),
scale_factor);
}
/* -------------------------------------------------------------------------- */
DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id)
: DOFData(dof_id) {}
/* -------------------------------------------------------------------------- */
auto DOFManagerDefault::getNewDOFData(const ID & dof_id)
-> std::unique_ptr<DOFData> {
return std::make_unique<DOFDataDefault>(dof_id);
}
/* -------------------------------------------------------------------------- */
std::tuple<Int, Int, Int>
DOFManagerDefault::registerDOFsInternal(const ID & dof_id,
Array<Real> & dofs_array) {
auto ret = DOFManager::registerDOFsInternal(dof_id, dofs_array);
// update the synchronizer if needed
if (this->synchronizer) {
this->synchronizer->registerDOFs(dof_id);
}
return ret;
}
/* -------------------------------------------------------------------------- */
SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id,
const MatrixType & matrix_type) {
return this->registerSparseMatrix<SparseMatrixAIJ>(*this, id, matrix_type);
}
/* -------------------------------------------------------------------------- */
SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id,
const ID & matrix_to_copy_id) {
return this->registerSparseMatrix<SparseMatrixAIJ>(id, matrix_to_copy_id);
}
/* -------------------------------------------------------------------------- */
SolverVector & DOFManagerDefault::getNewLumpedMatrix(const ID & id) {
return this->registerLumpedMatrix<SolverVectorDefault>(*this, id);
}
/* -------------------------------------------------------------------------- */
SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) {
auto & matrix = DOFManager::getMatrix(id);
return aka::as_type<SparseMatrixAIJ>(matrix);
}
/* -------------------------------------------------------------------------- */
NonLinearSolver &
DOFManagerDefault::getNewNonLinearSolver(const ID & id,
const NonLinearSolverType & type) {
switch (type) {
#if defined(AKANTU_USE_MUMPS)
case NonLinearSolverType::_newton_raphson:
/* FALLTHRU */
/* [[fallthrough]]; un-comment when compiler will get it */
case NonLinearSolverType::_newton_raphson_contact:
case NonLinearSolverType::_newton_raphson_modified: {
return this->registerNonLinearSolver<NonLinearSolverNewtonRaphson>(
*this, id, type);
}
case NonLinearSolverType::_linear: {
return this->registerNonLinearSolver<NonLinearSolverLinear>(*this, id,
type);
}
#endif
case NonLinearSolverType::_lumped: {
return this->registerNonLinearSolver<NonLinearSolverLumped>(*this, id,
type);
}
default:
AKANTU_EXCEPTION("The asked type of non linear solver is not supported by "
"this dof manager");
}
}
/* -------------------------------------------------------------------------- */
TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(
const ID & id, const TimeStepSolverType & type,
NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) {
return this->registerTimeStepSolver<TimeStepSolverDefault>(
*this, id, type, non_linear_solver, solver_callback);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id,
const Array<T> & global_array,
Array<T> & local_array) const {
AKANTU_DEBUG_IN();
const auto & equation_number = this->getLocalEquationsNumbers(dof_id);
auto nb_degree_of_freedoms = equation_number.size();
local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent());
for (auto data : zip(equation_number, make_view(local_array))) {
std::get<1>(data) = global_array(std::get<0>(data));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id,
const SolverVector & global_array,
Array<Real> & local_array) {
getArrayPerDOFs(dof_id,
aka::as_type<SolverVectorDefault>(global_array).getVector(),
local_array);
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::assembleLumpedMatMulVectToResidual(
const ID & dof_id, const ID & A_id, const Array<Real> & x,
Real scale_factor) {
const auto & A =
aka::as_type<SolverVectorArray>(this->getLumpedMatrix(A_id)).getVector();
auto & cache = aka::as_type<SolverVectorArray>(*this->data_cache);
cache.zero();
this->assembleToGlobalArray(dof_id, x, cache.getVector(), scale_factor);
for (auto && data : zip(make_view(A), make_view(cache.getVector()),
make_view(this->getResidualArray()))) {
const auto & A = std::get<0>(data);
const auto & x = std::get<1>(data);
auto & r = std::get<2>(data);
r += A * x;
}
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::assembleElementalMatricesToMatrix(
const ID & matrix_id, const ID & dof_id, const Array<Real> & elementary_mat,
ElementType type, GhostType ghost_type,
const MatrixType & elemental_matrix_type,
const Array<Int> & filter_elements) {
this->addToProfile(matrix_id, dof_id, type, ghost_type);
auto & A = getMatrix(matrix_id);
DOFManager::assembleElementalMatricesToMatrix_(
A, dof_id, elementary_mat, type, ghost_type, elemental_matrix_type,
filter_elements);
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::assemblePreassembledMatrix(
const ID & matrix_id, const TermsToAssemble & terms) {
auto & A = getMatrix(matrix_id);
DOFManager::assemblePreassembledMatrix_(A, terms);
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::assembleMatMulVectToArray(const ID & dof_id,
const ID & A_id,
const Array<Real> & x,
Array<Real> & array,
Real scale_factor) {
if (mesh->isDistributed()) {
DOFManager::assembleMatMulVectToArray_<SolverVectorDistributed>(
dof_id, A_id, x, array, scale_factor);
} else {
DOFManager::assembleMatMulVectToArray_<SolverVectorDefault>(
dof_id, A_id, x, array, scale_factor);
}
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id,
ElementType type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
const auto & dof_data = this->getDOFData(dof_id);
if (dof_data.support_type != _dst_nodal) {
return;
}
auto mat_dof = std::make_pair(matrix_id, dof_id);
auto type_pair = std::make_pair(type, ghost_type);
auto prof_it = this->matrix_profiled_dofs.find(mat_dof);
if (prof_it != this->matrix_profiled_dofs.end() &&
std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) !=
prof_it->second.end()) {
return;
}
auto nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent();
const auto & equation_number = this->getLocalEquationsNumbers(dof_id);
auto & A = this->getMatrix(matrix_id);
A.resize(system_size);
auto size = A.size();
auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto & connectivity = this->mesh->getConnectivity(type, ghost_type);
auto cbegin = connectivity.begin(nb_nodes_per_element);
auto cit = cbegin;
auto nb_elements = connectivity.size();
Int * ge_it = nullptr;
if (dof_data.group_support != "__mesh__") {
const auto & group_elements =
this->mesh->getElementGroup(dof_data.group_support)
.getElements(type, ghost_type);
ge_it = group_elements.data();
nb_elements = group_elements.size();
}
auto size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node;
Vector<Int> element_eq_nb(size_mat);
for (Int e = 0; e < nb_elements; ++e) {
if (ge_it != nullptr) {
cit = cbegin + *ge_it;
}
this->extractElementEquationNumber(
equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb);
std::transform(
element_eq_nb.begin(), element_eq_nb.end(), element_eq_nb.begin(),
[&](auto & local) { return this->localToGlobalEquationNumber(local); });
if (ge_it != nullptr) {
++ge_it;
} else {
++cit;
}
for (Int i = 0; i < size_mat; ++i) {
auto c_irn = element_eq_nb(i);
if (c_irn < size) {
for (Int j = 0; j < size_mat; ++j) {
auto c_jcn = element_eq_nb(j);
if (c_jcn < size) {
A.add(c_irn, c_jcn);
}
}
}
}
}
this->matrix_profiled_dofs[mat_dof].push_back(type_pair);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Array<Real> & DOFManagerDefault::getSolutionArray() {
return dynamic_cast<SolverVectorDefault *>(this->solution.get())->getVector();
}
/* -------------------------------------------------------------------------- */
const Array<Real> & DOFManagerDefault::getResidualArray() const {
return dynamic_cast<SolverVectorDefault *>(this->residual.get())->getVector();
}
/* -------------------------------------------------------------------------- */
Array<Real> & DOFManagerDefault::getResidualArray() {
return dynamic_cast<SolverVectorDefault *>(this->residual.get())->getVector();
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::onNodesAdded(const Array<Idx> & nodes_list,
const NewNodesEvent & event) {
DOFManager::onNodesAdded(nodes_list, event);
if (this->synchronizer) {
this->synchronizer->onNodesAdded(nodes_list);
}
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::resizeGlobalArrays() {
DOFManager::resizeGlobalArrays();
this->global_blocked_dofs.resize(this->local_system_size, 1);
this->previous_global_blocked_dofs.resize(this->local_system_size, 1);
matrix_profiled_dofs.clear();
}
/* -------------------------------------------------------------------------- */
void DOFManagerDefault::updateGlobalBlockedDofs() {
DOFManager::updateGlobalBlockedDofs();
if (this->global_blocked_dofs_release ==
this->previous_global_blocked_dofs_release) {
return;
}
global_blocked_dofs_uint.resize(local_system_size);
global_blocked_dofs_uint.set(false);
for (const auto & dof : global_blocked_dofs) {
global_blocked_dofs_uint[dof] = true;
}
}
/* -------------------------------------------------------------------------- */
Array<bool> & DOFManagerDefault::getBlockedDOFs() {
return global_blocked_dofs_uint;
}
/* -------------------------------------------------------------------------- */
const Array<bool> & DOFManagerDefault::getBlockedDOFs() const {
return global_blocked_dofs_uint;
}
/* -------------------------------------------------------------------------- */
-static bool dof_manager_is_registered =
+const bool dof_manager_is_registered [[maybe_unused]] =
DOFManagerFactory::getInstance().registerAllocator(
"default",
[](Mesh & mesh, const ID & id) -> std::unique_ptr<DOFManager> {
return std::make_unique<DOFManagerDefault>(mesh, id);
});
-static bool dof_manager_is_registered_mumps =
+#if defined(AKANTU_USE_MUMPS)
+const bool dof_manager_is_registered_mumps [[maybe_unused]] =
DOFManagerFactory::getInstance().registerAllocator(
"mumps", [](Mesh & mesh, const ID & id) -> std::unique_ptr<DOFManager> {
return std::make_unique<DOFManagerDefault>(mesh, id);
});
+#endif
} // namespace akantu
diff --git a/src/model/common/non_linear_solver/non_linear_solver.hh b/src/model/common/non_linear_solver/non_linear_solver.hh
index 8c7c5da19..10bb1fd1d 100644
--- a/src/model/common/non_linear_solver/non_linear_solver.hh
+++ b/src/model/common/non_linear_solver/non_linear_solver.hh
@@ -1,104 +1,104 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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_common.hh"
#include "parsable.hh"
/* -------------------------------------------------------------------------- */
#include <set>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_NON_LINEAR_SOLVER_HH_
#define AKANTU_NON_LINEAR_SOLVER_HH_
namespace akantu {
class DOFManager;
class SolverCallback;
} // namespace akantu
namespace akantu {
-class NonLinearSolver : public Parsable {
+class AKANTU_EXPORT NonLinearSolver : public Parsable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
NonLinearSolver(DOFManager & dof_manager,
const NonLinearSolverType & non_linear_solver_type,
const ID & id = "non_linear_solver");
~NonLinearSolver() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// solve the system described by the jacobian matrix, and rhs contained in
/// the dof manager
virtual void solve(SolverCallback & callback) = 0;
/// intercept the call to set for options
template <typename T> void set(const ID & param, T && t) {
if (has_internal_set_param) {
set_param(param, std::to_string(t));
} else {
ParameterRegistry::set(param, t);
}
}
protected:
void checkIfTypeIsSupported();
void assembleResidual(SolverCallback & callback);
/// internal set param for solvers that should intercept the parameters
virtual void set_param(const ID & /*param*/, const std::string & /*value*/) {}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
ID id;
DOFManager & _dof_manager;
/// type of non linear solver
NonLinearSolverType non_linear_solver_type;
/// list of supported non linear solver types
std::set<NonLinearSolverType> supported_type;
/// specifies if the set param should be redirected
bool has_internal_set_param{false};
};
namespace debug {
class NLSNotConvergedException : public Exception {
public:
NLSNotConvergedException(Real threshold, Int niter, Real error)
: Exception("The non linear solver did not converge."),
threshold(threshold), niter(niter), error(error) {}
Real threshold;
Int niter;
Real error;
};
} // namespace debug
} // namespace akantu
#endif /* AKANTU_NON_LINEAR_SOLVER_HH_ */