diff --git a/python/py_solver.cc b/python/py_solver.cc
index f2c83cfe6..bbcd931c3 100644
--- a/python/py_solver.cc
+++ b/python/py_solver.cc
@@ -1,105 +1,113 @@
/**
* @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("__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");
+ 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.")
+ .def("isDistributed", [](const SolverVector& self) { return self.isDistributed(); })
+ ;
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/solver/solver_vector.hh b/src/solver/solver_vector.hh
index 85fff5ccc..fd43c5ac5 100644
--- a/src/solver/solver_vector.hh
+++ b/src/solver/solver_vector.hh
@@ -1,91 +1,105 @@
/**
* @file solver_vector.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 21 2013
* @date last modification: Tue May 26 2020
*
* @brief Solver vector interface base class
*
*
* @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 "aka_array.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SOLVER_VECTOR_HH_
#define AKANTU_SOLVER_VECTOR_HH_
namespace akantu {
class DOFManager;
}
namespace akantu {
class SolverVector {
public:
SolverVector(DOFManager & dof_manager, const ID & id = "solver_vector")
: id(id), _dof_manager(dof_manager) {}
SolverVector(const SolverVector & vector, const ID & id = "solver_vector")
: id(id), _dof_manager(vector._dof_manager) {}
virtual ~SolverVector() = default;
// resize the vector to the size of the problem
virtual void resize() = 0;
// clear the vector
virtual void set(Real val) = 0;
void zero() { this->set({}); }
virtual operator const Array<Real> &() const = 0;
virtual Int size() const = 0;
virtual Int localSize() const = 0;
virtual SolverVector & operator+(const SolverVector & y) = 0;
virtual SolverVector & operator=(const SolverVector & y) = 0;
UInt & release() { return release_; }
UInt release() const { return release_; }
virtual void printself(std::ostream & stream, int indent = 0) const = 0;
+
+ /**
+ * \brief Returns `true` if `*this` is distributed or not.
+ *
+ * The default implementation returns `false`.
+ *
+ */
+ virtual
+ bool
+ isDistributed()
+ const
+ { return false; }
+
+
protected:
ID id;
/// Underlying dof manager
DOFManager & _dof_manager;
UInt release_{0};
};
inline std::ostream & operator<<(std::ostream & stream, SolverVector & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#endif /* AKANTU_SOLVER_VECTOR_HH_ */
diff --git a/src/solver/solver_vector_default.hh b/src/solver/solver_vector_default.hh
index 060c95668..932b87288 100644
--- a/src/solver/solver_vector_default.hh
+++ b/src/solver/solver_vector_default.hh
@@ -1,144 +1,154 @@
/**
* @file solver_vector_default.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jan 01 2019
* @date last modification: Sat May 23 2020
*
* @brief Solver vector interface to Array
*
*
* @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 "solver_vector.hh"
/* -------------------------------------------------------------------------- */
#include <utility>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SOLVER_VECTOR_DEFAULT_HH_
#define AKANTU_SOLVER_VECTOR_DEFAULT_HH_
namespace akantu {
class DOFManagerDefault;
} // namespace akantu
namespace akantu {
class SolverVectorArray : public SolverVector {
public:
SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id);
SolverVectorArray(const SolverVectorArray & vector, const ID & id);
~SolverVectorArray() override = default;
virtual Array<Real> & getVector() = 0;
virtual const Array<Real> & getVector() const = 0;
void printself(std::ostream & stream, int indent = 0) const override {
std::string space(indent, AKANTU_INDENT);
stream << space << "SolverVectorArray [" << std::endl;
stream << space << " + id: " << id << std::endl;
this->getVector().printself(stream, indent + 1);
stream << space << "]" << std::endl;
}
+
+ using SolverVector::isDistributed;
};
/* -------------------------------------------------------------------------- */
template <class Array_> class SolverVectorArrayTmpl : public SolverVectorArray {
public:
SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, Array_ & vector,
const ID & id = "solver_vector_default")
: SolverVectorArray(dof_manager, id), dof_manager(dof_manager),
vector(vector) {}
template <class A = Array_,
std::enable_if_t<not std::is_reference<A>::value> * = nullptr>
SolverVectorArrayTmpl(DOFManagerDefault & dof_manager,
const ID & id = "solver_vector_default")
: SolverVectorArray(dof_manager, id), dof_manager(dof_manager),
vector(0, 1, id + ":vector") {}
SolverVectorArrayTmpl(const SolverVectorArrayTmpl & vector,
const ID & id = "solver_vector_default")
: SolverVectorArray(vector, id), dof_manager(vector.dof_manager),
vector(vector.vector) {}
operator const Array<Real> &() const override { return getVector(); };
virtual operator Array<Real> &() { return getVector(); };
SolverVector & operator+(const SolverVector & y) override;
SolverVector & operator=(const SolverVector & y) override;
void resize() override {
static_assert(not std::is_const<std::remove_reference_t<Array_>>::value,
"Cannot resize a const Array");
this->vector.resize(this->localSize(), 0.);
++this->release_;
}
void set(Real val) override {
static_assert(not std::is_const<std::remove_reference_t<Array_>>::value,
"Cannot clear a const Array");
this->vector.set(val);
++this->release_;
}
+ virtual
+ bool
+ isDistributed()
+ const
+ override
+ { return false; }
+
+
public:
Array<Real> & getVector() override { return vector; }
const Array<Real> & getVector() const override { return vector; }
Int size() const override;
Int localSize() const override;
virtual Array<Real> & getGlobalVector() { return this->vector; }
virtual void setGlobalVector(const Array<Real> & solution) {
this->vector.copy(solution);
}
protected:
DOFManagerDefault & dof_manager;
Array_ vector;
template <class A> friend class SolverVectorArrayTmpl;
};
/* -------------------------------------------------------------------------- */
using SolverVectorDefault = SolverVectorArrayTmpl<Array<Real>>;
/* -------------------------------------------------------------------------- */
template <class Array>
using SolverVectorDefaultWrap = SolverVectorArrayTmpl<Array &>;
template <class Array>
decltype(auto) make_solver_vector_default_wrap(DOFManagerDefault & dof_manager,
Array & vector) {
return SolverVectorDefaultWrap<Array>(dof_manager, vector);
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
#include "solver_vector_default_tmpl.hh"
/* -------------------------------------------------------------------------- */
#endif /* AKANTU_SOLVER_VECTOR_DEFAULT_HH_ */
diff --git a/src/solver/solver_vector_distributed.hh b/src/solver/solver_vector_distributed.hh
index 26df926f4..e7b988144 100644
--- a/src/solver/solver_vector_distributed.hh
+++ b/src/solver/solver_vector_distributed.hh
@@ -1,59 +1,68 @@
/**
* @file solver_vector_distributed.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 21 2013
* @date last modification: Tue Jan 01 2019
*
* @brief Solver vector interface for distributed arrays
*
*
* @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 "solver_vector_default.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_
#define AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_
namespace akantu {
class SolverVectorDistributed : public SolverVectorDefault {
public:
SolverVectorDistributed(DOFManagerDefault & dof_manager,
const ID & id = "solver_vector_mumps");
SolverVectorDistributed(const SolverVectorDefault & vector,
const ID & id = "solver_vector_mumps");
Array<Real> & getGlobalVector() override;
void setGlobalVector(const Array<Real> & solution) override;
+
+ virtual
+ bool
+ isDistributed()
+ const
+ override
+ { return true; }
+
+
protected:
// full vector in case it needs to be centralized on master
std::unique_ptr<Array<Real>> global_vector;
};
} // namespace akantu
#endif /* AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ */
diff --git a/src/solver/solver_vector_petsc.hh b/src/solver/solver_vector_petsc.hh
index 223d59a98..570220445 100644
--- a/src/solver/solver_vector_petsc.hh
+++ b/src/solver/solver_vector_petsc.hh
@@ -1,208 +1,215 @@
/**
* @file solver_vector_petsc.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jan 01 2019
* @date last modification: Fri Jul 24 2020
*
* @brief Solver vector interface for PETSc
*
*
* @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 "dof_manager_petsc.hh"
#include "solver_vector.hh"
/* -------------------------------------------------------------------------- */
#include <petscvec.h>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SOLVER_VECTOR_PETSC_HH_
#define AKANTU_SOLVER_VECTOR_PETSC_HH_
namespace akantu {
class DOFManagerPETSc;
} // namespace akantu
namespace akantu {
/* -------------------------------------------------------------------------- */
namespace internal {
/* ------------------------------------------------------------------------ */
class PETScVector {
public:
virtual ~PETScVector() = default;
operator Vec &() { return x; }
operator const Vec &() const { return x; }
Int size() const {
PetscInt n;
PETSc_call(VecGetSize, x, &n);
return n;
}
Int local_size() const {
PetscInt n;
PETSc_call(VecGetLocalSize, x, &n);
return n;
}
AKANTU_GET_MACRO_NOT_CONST(Vec, x, auto &);
AKANTU_GET_MACRO(Vec, x, const auto &);
protected:
Vec x{nullptr};
};
} // namespace internal
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
class SolverVectorPETSc : public SolverVector, public internal::PETScVector {
public:
SolverVectorPETSc(DOFManagerPETSc & dof_manager,
const ID & id = "solver_vector_petsc");
SolverVectorPETSc(const SolverVectorPETSc & vector,
const ID & id = "solver_vector_petsc");
SolverVectorPETSc(Vec x, DOFManagerPETSc & dof_manager,
const ID & id = "solver_vector_petsc");
~SolverVectorPETSc() override;
// resize the vector to the size of the problem
void resize() override;
void set(Real val) override;
operator const Array<Real> &() const override;
SolverVector & operator+(const SolverVector & y) override;
SolverVector & operator=(const SolverVector & y) override;
SolverVectorPETSc & operator=(const SolverVectorPETSc & y);
/// get values using processors global indexes
void getValues(const Array<Int> & idx, Array<Real> & values) const;
/// get values using processors local indexes
void getValuesLocal(const Array<Int> & idx, Array<Real> & values) const;
/// adding values to the vector using the global indices
void addValues(const Array<Int> & gidx, const Array<Real> & values,
Real scale_factor = 1.);
/// adding values to the vector using the local indices
void addValuesLocal(const Array<Int> & lidx, const Array<Real> & values,
Real scale_factor = 1.);
Int size() const override { return internal::PETScVector::size(); }
Int localSize() const override { return internal::PETScVector::local_size(); }
void printself(std::ostream & stream, int indent = 0) const override;
+ virtual
+ bool
+ isDistributed()
+ const
+ override
+ { return true; }
+
protected:
void applyModifications();
void updateGhost();
protected:
DOFManagerPETSc & dof_manager;
// used for the conversion operator
Array<Real> cache;
};
/* -------------------------------------------------------------------------- */
namespace internal {
/* ------------------------------------------------------------------------ */
template <class Array> class PETScWrapedVector : public PETScVector {
public:
PETScWrapedVector(Array && array) : array(array) {
PETSc_call(VecCreateSeqWithArray, PETSC_COMM_SELF, 1, array.size(),
array.storage(), &x);
}
~PETScWrapedVector() override { PETSc_call(VecDestroy, &x); }
private:
Array array;
};
/* ------------------------------------------------------------------------ */
template <bool read_only> class PETScLocalVector : public PETScVector {
public:
PETScLocalVector(const Vec & g) : g(g) {
PETSc_call(VecGetLocalVectorRead, g, x);
}
PETScLocalVector(const SolverVectorPETSc & g)
: PETScLocalVector(g.getVec()) {}
~PETScLocalVector() override {
PETSc_call(VecRestoreLocalVectorRead, g, x);
PETSc_call(VecDestroy, &x);
}
private:
const Vec & g;
};
template <> class PETScLocalVector<false> : public PETScVector {
public:
PETScLocalVector(Vec & g) : g(g) {
PETSc_call(VecGetLocalVectorRead, g, x);
}
PETScLocalVector(SolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {}
~PETScLocalVector() override {
PETSc_call(VecRestoreLocalVectorRead, g, x);
PETSc_call(VecDestroy, &x);
}
private:
Vec & g;
};
/* ------------------------------------------------------------------------ */
template <class Array>
decltype(auto) make_petsc_wraped_vector(Array && array) {
return PETScWrapedVector<Array>(std::forward<Array>(array));
}
template <
typename V,
std::enable_if_t<std::is_same<Vec, std::decay_t<V>>::value> * = nullptr>
decltype(auto) make_petsc_local_vector(V && vec) {
constexpr auto read_only = std::is_const<std::remove_reference_t<V>>::value;
return PETScLocalVector<read_only>(vec);
}
template <typename V, std::enable_if_t<std::is_base_of<
SolverVector, std::decay_t<V>>::value> * = nullptr>
decltype(auto) make_petsc_local_vector(V && vec) {
constexpr auto read_only = std::is_const<std::remove_reference_t<V>>::value;
return PETScLocalVector<read_only>(
dynamic_cast<std::conditional_t<read_only, const SolverVectorPETSc,
SolverVectorPETSc> &>(vec));
}
} // namespace internal
} // namespace akantu
#endif /* AKANTU_SOLVER_VECTOR_PETSC_HH_ */