diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index d04f2059d81a2550a7eccfea190d9a455ae881f4..b9893d29205d357d602e98f32610ea908c033aca 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -7,14 +7,13 @@
* http://www.opengeosys.org/project/license
*/
-#include "BaseLib/Logging.h"
+#include "GeometricBasics.h"
#include <Eigen/Dense>
+#include "BaseLib/Logging.h"
#include "Point3d.h"
-#include "GeometricBasics.h"
-
namespace MathLib
{
double orientation3d(MathLib::Point3d const& p,
@@ -64,31 +63,31 @@ bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a,
MathLib::Point3d const& b, MathLib::Point3d const& c,
MathLib::Point3d const& d, double eps)
{
- double const d0 (MathLib::orientation3d(d,a,b,c));
+ double const d0(MathLib::orientation3d(d, a, b, c));
// if tetrahedron is not coplanar
if (std::abs(d0) > std::numeric_limits<double>::epsilon())
{
- bool const d0_sign (d0>0);
+ bool const d0_sign(d0 > 0);
// if p is on the same side of bcd as a
- double const d1 (MathLib::orientation3d(d, p, b, c));
+ double const d1(MathLib::orientation3d(d, p, b, c));
if (!(d0_sign == (d1 >= 0) || std::abs(d1) < eps))
{
return false;
}
// if p is on the same side of acd as b
- double const d2 (MathLib::orientation3d(d, a, p, c));
+ double const d2(MathLib::orientation3d(d, a, p, c));
if (!(d0_sign == (d2 >= 0) || std::abs(d2) < eps))
{
return false;
}
// if p is on the same side of abd as c
- double const d3 (MathLib::orientation3d(d, a, b, p));
+ double const d3(MathLib::orientation3d(d, a, b, p));
if (!(d0_sign == (d3 >= 0) || std::abs(d3) < eps))
{
return false;
}
// if p is on the same side of abc as d
- double const d4 (MathLib::orientation3d(p, a, b, c));
+ double const d4(MathLib::orientation3d(p, a, b, c));
return d0_sign == (d4 >= 0) || std::abs(d4) < eps;
}
return false;
diff --git a/MathLib/Integration/GaussLegendre.cpp b/MathLib/Integration/GaussLegendre.cpp
index 50d802f321aa69efe4422f4e524f6d3df3dc47cf..bde424f321e1d57e5acdc324194babb3aa40cb71 100644
--- a/MathLib/Integration/GaussLegendre.cpp
+++ b/MathLib/Integration/GaussLegendre.cpp
@@ -14,19 +14,27 @@
namespace MathLib
{
+template <>
+double const GaussLegendre<1>::X[1] = {0.};
+template <>
+double const GaussLegendre<1>::W[1] = {2.};
-template <> double const GaussLegendre<1>::X[1] = {0.};
-template <> double const GaussLegendre<1>::W[1] = {2.};
+template <>
+double const GaussLegendre<2>::X[2] = {0.577350269189626, -0.577350269189626};
+template <>
+double const GaussLegendre<2>::W[2] = {1., 1.};
-template <> double const GaussLegendre<2>::X[2] = {0.577350269189626, -0.577350269189626};
-template <> double const GaussLegendre<2>::W[2] = {1., 1.};
+template <>
+double const GaussLegendre<3>::X[3] = {0.774596669241483, 0.,
+ -0.774596669241483};
+template <>
+double const GaussLegendre<3>::W[3] = {5. / 9, 8. / 9, 5. / 9};
-template <> double const GaussLegendre<3>::X[3] = {0.774596669241483, 0., -0.774596669241483};
-template <> double const GaussLegendre<3>::W[3] = {5./9, 8./9, 5./9};
+template <>
+double const GaussLegendre<4>::X[4] = {-0.861136311594053, -0.339981043584856,
+ 0.339981043584856, 0.861136311594053};
+template <>
+double const GaussLegendre<4>::W[4] = {0.347854845137454, 0.652145154862546,
+ 0.652145154862546, 0.347854845137454};
-template <> double const GaussLegendre<4>::X[4] =
- {-0.861136311594053, -0.339981043584856, 0.339981043584856, 0.861136311594053};
-template <> double const GaussLegendre<4>::W[4] =
- { 0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
-
-} // namespace MathLib
+} // namespace MathLib
diff --git a/MathLib/Integration/GaussLegendrePyramid.cpp b/MathLib/Integration/GaussLegendrePyramid.cpp
index 3fa19342dd03b95685b78428d2437a8ec486f06e..19140b36baca7a0cf8739b0cff2cbe658091e7ac 100644
--- a/MathLib/Integration/GaussLegendrePyramid.cpp
+++ b/MathLib/Integration/GaussLegendrePyramid.cpp
@@ -12,50 +12,45 @@
namespace MathLib
{
-
-template <> const std::array<std::array<double, 3>, GaussLegendrePyramid<1>::NPoints>
-GaussLegendrePyramid<1>::X = {{{{1./3., 1./3., 1./3.}}}};
-template <> double const GaussLegendrePyramid<1>::W[1] = {1.0};
+template <>
+const std::array<std::array<double, 3>, GaussLegendrePyramid<1>::NPoints>
+ GaussLegendrePyramid<1>::X = {{{{1. / 3., 1. / 3., 1. / 3.}}}};
+template <>
+double const GaussLegendrePyramid<1>::W[1] = {1.0};
const std::array<std::array<double, 3>, GaussLegendrePyramid<2>::NPoints>
-GaussLegendrePyramid<2>::X = {{ {{-0.584237394672177188, -0.584237394672177188, -2./3.}},
- {{ 0.584237394672177188, -0.584237394672177188, -2./3.}},
- {{ 0.584237394672177188, 0.584237394672177188, -2./3.}},
- {{-0.584237394672177188, 0.584237394672177188, -2./3.}},
- {{0., 0., 2./5.}}
- }};
-double const GaussLegendrePyramid<2>::W[5] = {81./100., 81./100., 81./100., 81./100., 125./27.};
+ GaussLegendrePyramid<2>::X = {
+ {{{-0.584237394672177188, -0.584237394672177188, -2. / 3.}},
+ {{0.584237394672177188, -0.584237394672177188, -2. / 3.}},
+ {{0.584237394672177188, 0.584237394672177188, -2. / 3.}},
+ {{-0.584237394672177188, 0.584237394672177188, -2. / 3.}},
+ {{0., 0., 2. / 5.}}}};
+double const GaussLegendrePyramid<2>::W[5] = {
+ 81. / 100., 81. / 100., 81. / 100., 81. / 100., 125. / 27.};
const std::array<std::array<double, 3>, GaussLegendrePyramid<3>::NPoints>
-GaussLegendrePyramid<3>::X = {{
- {{-0.673931986207731726, -0.673931986207731726, -0.142857142857142857}},
- {{ 0.673931986207731726, -0.673931986207731726, -0.142857142857142857}},
- {{ 0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
- {{-0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
- {{-0.610639618865075532, 0.0, -0.321428571428571429}},
- {{ 0.610639618865075532, 0.0, -0.321428571428571429}},
- {{ 0.0, -0.610639618865075532, -0.321428571428571429}},
- {{ 0.0, 0.610639618865075532, -0.321428571428571429}},
- {{ 0.0, 0.0, 0.524394036075370072}},
- {{-0.580939660561084423, -0.580939660561084423, -0.830065359477124183}},
- {{ 0.580939660561084423, -0.580939660561084423, -0.830065359477124183}},
- {{ 0.580939660561084423, 0.580939660561084423, -0.830065359477124183}},
- {{-0.580939660561084423, 0.580939660561084423, -0.830065359477124183}}
- }};
+ GaussLegendrePyramid<3>::X = {
+ {{{-0.673931986207731726, -0.673931986207731726,
+ -0.142857142857142857}},
+ {{0.673931986207731726, -0.673931986207731726, -0.142857142857142857}},
+ {{0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
+ {{-0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
+ {{-0.610639618865075532, 0.0, -0.321428571428571429}},
+ {{0.610639618865075532, 0.0, -0.321428571428571429}},
+ {{0.0, -0.610639618865075532, -0.321428571428571429}},
+ {{0.0, 0.610639618865075532, -0.321428571428571429}},
+ {{0.0, 0.0, 0.524394036075370072}},
+ {{-0.580939660561084423, -0.580939660561084423,
+ -0.830065359477124183}},
+ {{0.580939660561084423, -0.580939660561084423, -0.830065359477124183}},
+ {{0.580939660561084423, 0.580939660561084423, -0.830065359477124183}},
+ {{-0.580939660561084423, 0.580939660561084423,
+ -0.830065359477124183}}}};
double const GaussLegendrePyramid<3>::W[13] = {
- 0.515003019323671498,
- 0.515003019323671498,
- 0.515003019323671498,
- 0.515003019323671498,
- 0.257183745242064659,
- 0.257183745242064659,
- 0.257183745242064659,
- 0.257183745242064659,
- 2.474004977113405936,
- 0.419515737191525950,
- 0.419515737191525950,
- 0.419515737191525950,
- 0.419515737191525950
- };
+ 0.515003019323671498, 0.515003019323671498, 0.515003019323671498,
+ 0.515003019323671498, 0.257183745242064659, 0.257183745242064659,
+ 0.257183745242064659, 0.257183745242064659, 2.474004977113405936,
+ 0.419515737191525950, 0.419515737191525950, 0.419515737191525950,
+ 0.419515737191525950};
} // namespace MathLib
diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
index fe8022430283475b14b9242e814b89a8f4c5ba7e..593b251cb16cbe6ccd800e817994150643729c3b 100644
--- a/MathLib/Integration/GaussLegendreTet.cpp
+++ b/MathLib/Integration/GaussLegendreTet.cpp
@@ -19,12 +19,13 @@ template <>
double const GaussLegendreTet<1>::W[1] = {1. / 6.};
const std::array<std::array<double, 3>, GaussLegendreTet<2>::NPoints>
-GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
- {{1./6., 1./6., 1./6.}},
- {{1./2., 1./6., 1./6.}},
- {{1./6., 1./2., 1./6.}},
- {{1./6., 1./6., 1./2.}} }};
-double const GaussLegendreTet<2>::W[5] = {-2./15., 0.075, 0.075, 0.075, 0.075};
+ GaussLegendreTet<2>::X = {{{{1. / 4., 1. / 4., 1. / 4.}},
+ {{1. / 6., 1. / 6., 1. / 6.}},
+ {{1. / 2., 1. / 6., 1. / 6.}},
+ {{1. / 6., 1. / 2., 1. / 6.}},
+ {{1. / 6., 1. / 6., 1. / 2.}}}};
+double const GaussLegendreTet<2>::W[5] = {-2. / 15., 0.075, 0.075, 0.075,
+ 0.075};
static std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints>
initGLTet3X()
diff --git a/MathLib/Integration/GaussLegendreTri.cpp b/MathLib/Integration/GaussLegendreTri.cpp
index 27e1b038474ad5ca9041762cdcc559391354eca2..efee648484e2e420bee89ec20f131836bc634b7f 100644
--- a/MathLib/Integration/GaussLegendreTri.cpp
+++ b/MathLib/Integration/GaussLegendreTri.cpp
@@ -16,23 +16,24 @@
namespace MathLib
{
-
template <>
const std::array<std::array<double, 2>, GaussLegendreTri<1>::NPoints>
GaussLegendreTri<1>::X = {{{{1. / 3., 1. / 3.}}}};
-template <> double const GaussLegendreTri<1>::W[1] = {1.0};
+template <>
+double const GaussLegendreTri<1>::W[1] = {1.0};
const std::array<std::array<double, 2>, GaussLegendreTri<2>::NPoints>
GaussLegendreTri<2>::X = {
{{{1. / 6., 1. / 6.}}, {{2. / 3., 1. / 6.}}, {{1. / 6., 2. / 3.}}}};
-double const GaussLegendreTri<2>::W[3] = {1./3., 1./3., 1./3.};
+double const GaussLegendreTri<2>::W[3] = {1. / 3., 1. / 3., 1. / 3.};
const std::array<std::array<double, 2>, GaussLegendreTri<3>::NPoints>
GaussLegendreTri<3>::X = {{{{1. / 3., 1. / 3.}},
{{1. / 5., 3. / 5.}},
{{1. / 5., 1. / 5.}},
{{3. / 5., 1. / 5.}}}};
-double const GaussLegendreTri<3>::W[4] = {-27./48., 25./48., 25./48., 25./48.};
+double const GaussLegendreTri<3>::W[4] = {-27. / 48., 25. / 48., 25. / 48.,
+ 25. / 48.};
const std::array<std::array<double, 2>, GaussLegendreTri<4>::NPoints>
GaussLegendreTri<4>::X = {{{{1. / 3., 1. / 3.}},
@@ -42,8 +43,12 @@ const std::array<std::array<double, 2>, GaussLegendreTri<4>::NPoints>
{{0.797426985353087, 0.101286507323456}},
{{0.101286507323456, 0.797426985353087}},
{{0.101286507323456, 0.101286507323456}}}};
-double const GaussLegendreTri<4>::W[7] = {
- 0.225, 0.132394152788506, 0.132394152788506, 0.132394152788506,
- 0.125939180544827, 0.125939180544827, 0.125939180544827};
+double const GaussLegendreTri<4>::W[7] = {0.225,
+ 0.132394152788506,
+ 0.132394152788506,
+ 0.132394152788506,
+ 0.125939180544827,
+ 0.125939180544827,
+ 0.125939180544827};
} // namespace MathLib
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index e8d4db27c1b57adacf1af7a52e510ad079332bbd..cdff20a3e61d6ed6b4f1b302e931af44f58ca32e 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
+++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
@@ -11,6 +11,8 @@
* http://www.opengeosys.org/project/license
*
*/
+#include "PiecewiseLinearInterpolation.h"
+
#include <cmath>
#include <limits>
#include <utility>
@@ -18,13 +20,12 @@
#include "BaseLib/Error.h"
#include "BaseLib/quicksort.h"
-#include "PiecewiseLinearInterpolation.h"
-
namespace MathLib
{
PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(
std::vector<double> supporting_points,
- std::vector<double> values_at_supp_pnts,
+ std::vector<double>
+ values_at_supp_pnts,
bool supp_pnts_sorted)
: _supp_pnts(std::move(supporting_points)),
_values_at_supp_pnts(std::move(values_at_supp_pnts))
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 6595a9ac30cfd3b19956a6c4ae7cb5cac999d355..fbcb8cd7a9ee6fee909bb9f02a7201195a1ac534 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -17,21 +17,20 @@
#endif
#ifdef USE_EIGEN_UNSUPPORTED
-#include <Eigen/Sparse>
#include <unsupported/Eigen/src/IterativeSolvers/GMRES.h>
#include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
+
+#include <Eigen/Sparse>
#endif
#include "BaseLib/ConfigTree.h"
-#include "EigenVector.h"
#include "EigenMatrix.h"
#include "EigenTools.h"
-
+#include "EigenVector.h"
#include "MathLib/LinAlg/LinearSolverOptions.h"
namespace MathLib
{
-
// TODO change to LinearSolver
class EigenLinearSolverBase
{
@@ -42,12 +41,12 @@ public:
virtual ~EigenLinearSolverBase() = default;
//! Solves the linear equation system \f$ A x = b \f$ for \f$ x \f$.
- virtual bool solve(Matrix &A, Vector const& b, Vector &x, EigenOption &opt) = 0;
+ virtual bool solve(Matrix& A, Vector const& b, Vector& x,
+ EigenOption& opt) = 0;
};
namespace details
{
-
/// Template class for Eigen direct linear solvers
template <class T_SOLVER>
class EigenDirectLinearSolver final : public EigenLinearSolverBase
@@ -62,13 +61,15 @@ public:
}
_solver.compute(A);
- if(_solver.info()!=Eigen::Success) {
+ if (_solver.info() != Eigen::Success)
+ {
ERR("Failed during Eigen linear solver initialization");
return false;
}
x = _solver.solve(b);
- if(_solver.info()!=Eigen::Success) {
+ if (_solver.info() != Eigen::Success)
+ {
ERR("Failed during Eigen linear solve");
return false;
}
@@ -92,7 +93,8 @@ public:
EigenOption::getPreconName(opt.precon_type));
_solver.setTolerance(opt.error_tolerance);
_solver.setMaxIterations(opt.max_iterations);
- MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setRestart(opt.restart);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setRestart(
+ opt.restart);
if (!A.isCompressed())
{
@@ -100,7 +102,8 @@ public:
}
_solver.compute(A);
- if(_solver.info()!=Eigen::Success) {
+ if (_solver.info() != Eigen::Success)
+ {
ERR("Failed during Eigen linear solver initialization");
return false;
}
@@ -110,7 +113,8 @@ public:
opt.max_iterations);
INFO("\t residual: {:e}\n", _solver.error());
- if(_solver.info()!=Eigen::Success) {
+ if (_solver.info() != Eigen::Success)
+ {
ERR("Failed during Eigen linear solve");
return false;
}
@@ -120,8 +124,7 @@ public:
private:
T_SOLVER _solver;
- void setRestart(int const /*restart*/) {
- }
+ void setRestart(int const /*restart*/) {}
};
/// Specialization for (all) three preconditioners separately
@@ -155,8 +158,8 @@ void EigenIterativeLinearSolver<
template <template <typename, typename> class Solver, typename Precon>
std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
{
- using Slv = EigenIterativeLinearSolver<
- Solver<EigenMatrix::RawMatrixType, Precon>>;
+ using Slv =
+ EigenIterativeLinearSolver<Solver<EigenMatrix::RawMatrixType, Precon>>;
return std::make_unique<Slv>();
}
@@ -164,7 +167,8 @@ template <template <typename, typename> class Solver>
std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
EigenOption::PreconType precon_type)
{
- switch (precon_type) {
+ switch (precon_type)
+ {
case EigenOption::PreconType::NONE:
return createIterativeSolver<Solver,
Eigen::IdentityPreconditioner>();
@@ -173,9 +177,10 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
Solver, Eigen::DiagonalPreconditioner<double>>();
case EigenOption::PreconType::ILUT:
// TODO for this preconditioner further options can be passed.
- // see https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
- return createIterativeSolver<
- Solver, Eigen::IncompleteLUT<double>>();
+ // see
+ // https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
+ return createIterativeSolver<Solver,
+ Eigen::IncompleteLUT<double>>();
default:
OGS_FATAL("Invalid Eigen preconditioner type.");
}
@@ -187,14 +192,18 @@ using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
{
- switch (solver_type) {
- case EigenOption::SolverType::BiCGSTAB: {
+ switch (solver_type)
+ {
+ case EigenOption::SolverType::BiCGSTAB:
+ {
return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
}
- case EigenOption::SolverType::CG: {
+ case EigenOption::SolverType::CG:
+ {
return createIterativeSolver<EigenCGSolver>(precon_type);
}
- case EigenOption::SolverType::GMRES: {
+ case EigenOption::SolverType::GMRES:
+ {
#ifdef USE_EIGEN_UNSUPPORTED
return createIterativeSolver<Eigen::GMRES>(precon_type);
#else
@@ -210,9 +219,8 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
} // namespace details
-EigenLinearSolver::EigenLinearSolver(
- const std::string& /*solver_name*/,
- const BaseLib::ConfigTree* const option)
+EigenLinearSolver::EigenLinearSolver(const std::string& /*solver_name*/,
+ const BaseLib::ConfigTree* const option)
{
using Matrix = EigenMatrix::RawMatrixType;
@@ -223,8 +231,10 @@ EigenLinearSolver::EigenLinearSolver(
// TODO for my taste it is much too unobvious that the default solver type
// currently is SparseLU.
- switch (_option.solver_type) {
- case EigenOption::SolverType::SparseLU: {
+ switch (_option.solver_type)
+ {
+ case EigenOption::SolverType::SparseLU:
+ {
using SolverType =
Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
_solver = std::make_unique<
@@ -237,7 +247,8 @@ EigenLinearSolver::EigenLinearSolver(
_solver = details::createIterativeSolver(_option.solver_type,
_option.precon_type);
return;
- case EigenOption::SolverType::PardisoLU: {
+ case EigenOption::SolverType::PardisoLU:
+ {
#ifdef USE_MKL
using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
_solver.reset(new details::EigenDirectLinearSolver<SolverType>);
@@ -267,27 +278,32 @@ void EigenLinearSolver::setOption(BaseLib::ConfigTree const& option)
if (auto solver_type =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__solver_type}
- ptSolver->getConfigParameterOptional<std::string>("solver_type")) {
+ ptSolver->getConfigParameterOptional<std::string>("solver_type"))
+ {
_option.solver_type = MathLib::EigenOption::getSolverType(*solver_type);
}
if (auto precon_type =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__precon_type}
- ptSolver->getConfigParameterOptional<std::string>("precon_type")) {
+ ptSolver->getConfigParameterOptional<std::string>("precon_type"))
+ {
_option.precon_type = MathLib::EigenOption::getPreconType(*precon_type);
}
if (auto error_tolerance =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__error_tolerance}
- ptSolver->getConfigParameterOptional<double>("error_tolerance")) {
+ ptSolver->getConfigParameterOptional<double>("error_tolerance"))
+ {
_option.error_tolerance = *error_tolerance;
}
if (auto max_iteration_step =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__max_iteration_step}
- ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
+ ptSolver->getConfigParameterOptional<int>("max_iteration_step"))
+ {
_option.max_iterations = *max_iteration_step;
}
if (auto scaling =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__scaling}
- ptSolver->getConfigParameterOptional<bool>("scaling")) {
+ ptSolver->getConfigParameterOptional<bool>("scaling"))
+ {
#ifdef USE_EIGEN_UNSUPPORTED
_option.scaling = *scaling;
#else
@@ -298,7 +314,8 @@ void EigenLinearSolver::setOption(BaseLib::ConfigTree const& option)
}
if (auto restart =
//! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__restart}
- ptSolver->getConfigParameterOptional<int>("restart")) {
+ ptSolver->getConfigParameterOptional<int>("restart"))
+ {
#ifdef USE_EIGEN_UNSUPPORTED
_option.restart = *restart;
#else
@@ -307,10 +324,9 @@ void EigenLinearSolver::setOption(BaseLib::ConfigTree const& option)
"GMRES/GMRES option restart is not available.");
#endif
}
-
}
-bool EigenLinearSolver::solve(EigenMatrix &A, EigenVector& b, EigenVector &x)
+bool EigenLinearSolver::solve(EigenMatrix& A, EigenVector& b, EigenVector& x)
{
INFO("------------------------------------------------------------------");
INFO("*** Eigen solver computation");
diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
index 80cc16cd6993773cbd8217e3bfa047f2475d38f2..dd722f7441b3a37471830b5e1531a92a37a74448 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.cpp
+++ b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -9,11 +9,11 @@
*/
#include "EigenOption.h"
+
#include "BaseLib/Error.h"
namespace MathLib
{
-
EigenOption::EigenOption()
{
solver_type = SolverType::SparseLU;
@@ -26,7 +26,8 @@ EigenOption::EigenOption()
#endif
}
-EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_name)
+EigenOption::SolverType EigenOption::getSolverType(
+ const std::string& solver_name)
{
if (solver_name == "CG")
{
@@ -52,7 +53,8 @@ EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_nam
OGS_FATAL("Unknown Eigen solver type `{:s}'", solver_name);
}
-EigenOption::PreconType EigenOption::getPreconType(const std::string &precon_name)
+EigenOption::PreconType EigenOption::getPreconType(
+ const std::string& precon_name)
{
if (precon_name == "NONE")
{
@@ -72,7 +74,8 @@ EigenOption::PreconType EigenOption::getPreconType(const std::string &precon_nam
std::string EigenOption::getSolverName(SolverType const solver_type)
{
- switch (solver_type) {
+ switch (solver_type)
+ {
case SolverType::CG:
return "CG";
case SolverType::BiCGSTAB:
@@ -89,7 +92,8 @@ std::string EigenOption::getSolverName(SolverType const solver_type)
std::string EigenOption::getPreconName(PreconType const precon_type)
{
- switch (precon_type) {
+ switch (precon_type)
+ {
case PreconType::NONE:
return "NONE";
case PreconType::DIAGONAL:
diff --git a/MathLib/LinAlg/Eigen/EigenTools.cpp b/MathLib/LinAlg/Eigen/EigenTools.cpp
index 97f5b0f9d6d94bf5231e7b0ae03487de814e995b..5985fa16a839f9181e4f4eee83c903b1c7773438 100644
--- a/MathLib/LinAlg/Eigen/EigenTools.cpp
+++ b/MathLib/LinAlg/Eigen/EigenTools.cpp
@@ -10,27 +10,27 @@
#include "EigenTools.h"
-
#include "EigenVector.h"
namespace MathLib
{
-
-void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &/*x*/,
- const std::vector<EigenMatrix::IndexType> &vec_knownX_id,
- const std::vector<double> &vec_knownX_x, double /*penalty_scaling*/)
+void applyKnownSolution(
+ EigenMatrix& A, EigenVector& b, EigenVector& /*x*/,
+ const std::vector<EigenMatrix::IndexType>& vec_knownX_id,
+ const std::vector<double>& vec_knownX_x, double /*penalty_scaling*/)
{
using SpMat = EigenMatrix::RawMatrixType;
static_assert(SpMat::IsRowMajor, "matrix is assumed to be row major!");
- auto &A_eigen = A.getRawMatrix();
- auto &b_eigen = b.getRawVector();
+ auto& A_eigen = A.getRawMatrix();
+ auto& b_eigen = b.getRawVector();
// A_eigen(k, j) = 0.
// set row to zero
for (auto row_id : vec_knownX_id)
{
- for (SpMat::InnerIterator it(A_eigen,row_id); it; ++it) {
+ for (SpMat::InnerIterator it(A_eigen, row_id); it; ++it)
+ {
if (it.col() != decltype(it.col())(row_id))
{
it.valueRef() = 0.0;
@@ -40,7 +40,7 @@ void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &/*x*/,
SpMat AT = A_eigen.transpose();
- for (std::size_t ix=0; ix<vec_knownX_id.size(); ix++)
+ for (std::size_t ix = 0; ix < vec_knownX_id.size(); ix++)
{
SpMat::Index const row_id = vec_knownX_id[ix];
auto const x = vec_knownX_x[ix];
@@ -54,14 +54,17 @@ void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &/*x*/,
continue;
}
- b_eigen[it.col()] -= it.value()*x;
+ b_eigen[it.col()] -= it.value() * x;
it.valueRef() = 0.0;
}
auto& c = AT.coeffRef(row_id, row_id);
- if (c != 0.0) {
+ if (c != 0.0)
+ {
b_eigen[row_id] = x * c;
- } else {
+ }
+ else
+ {
b_eigen[row_id] = x;
c = 1.0;
}
diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
index 05e30cdf91de49df877073a606f7d4f9f6416556..d2495cfd00157117aa1e49f5ce2be267dd058369 100644
--- a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
+++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.cpp
@@ -17,27 +17,26 @@
#include "BaseLib/ConfigTree.h"
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
#include "MathLib/LinAlg/Eigen/EigenVector.h"
+#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
#include "MathLib/LinAlg/Lis/LisMatrix.h"
#include "MathLib/LinAlg/Lis/LisVector.h"
-#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
namespace MathLib
{
EigenLisLinearSolver::EigenLisLinearSolver(
- const std::string /*solver_name*/,
- BaseLib::ConfigTree const* const option)
+ const std::string /*solver_name*/, BaseLib::ConfigTree const* const option)
: _lis_option(option)
{
}
-bool EigenLisLinearSolver::solve(EigenMatrix &A_, EigenVector& b_,
- EigenVector &x_)
+bool EigenLisLinearSolver::solve(EigenMatrix& A_, EigenVector& b_,
+ EigenVector& x_)
{
static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
"Sparse matrix is required to be in row major storage.");
- auto &A = A_.getRawMatrix();
- auto &b = b_.getRawVector();
- auto &x = x_.getRawVector();
+ auto& A = A_.getRawMatrix();
+ auto& b = b_.getRawVector();
+ auto& x = x_.getRawVector();
if (!A.isCompressed())
A.makeCompressed();
@@ -49,14 +48,14 @@ bool EigenLisLinearSolver::solve(EigenMatrix &A_, EigenVector& b_,
LisVector lisb(b.rows(), b.data());
LisVector lisx(x.rows(), x.data());
- LisLinearSolver lissol; // TODO not always create Lis solver here
+ LisLinearSolver lissol; // TODO not always create Lis solver here
lissol.setOption(_lis_option);
bool const status = lissol.solve(lisA, lisb, lisx);
- for (std::size_t i=0; i<lisx.size(); i++)
+ for (std::size_t i = 0; i < lisx.size(); i++)
x[i] = lisx[i];
return status;
}
-} //MathLib
+} // namespace MathLib
diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp
index 6489bff05d9a731fa65aff5abbb6415a899aab36..486601d02e81c3e070d88393392746f0130e5ea0 100644
--- a/MathLib/LinAlg/LinAlg.cpp
+++ b/MathLib/LinAlg/LinAlg.cpp
@@ -15,12 +15,13 @@
// Global PETScMatrix/PETScVector //////////////////////////////////////////
#ifdef USE_PETSC
-#include "MathLib/LinAlg/PETSc/PETScVector.h"
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+#include "MathLib/LinAlg/PETSc/PETScVector.h"
-namespace MathLib { namespace LinAlg
+namespace MathLib
+{
+namespace LinAlg
{
-
// Vector
void setLocalAccessibleVector(PETScVector const& x)
@@ -35,7 +36,8 @@ void set(PETScVector& x, double const a)
void copy(PETScVector const& x, PETScVector& y)
{
- if (!y.getRawVector()) y.shallowCopy(x);
+ if (!y.getRawVector())
+ y.shallowCopy(x);
VecCopy(x.getRawVector(), y.getRawVector());
}
@@ -67,16 +69,16 @@ void axpby(PETScVector& y, double const a, double const b, PETScVector const& x)
// Explicit specialization
// Computes w = x/y componentwise.
-template<>
-void componentwiseDivide(PETScVector& w,
- PETScVector const& x, PETScVector const& y)
+template <>
+void componentwiseDivide(PETScVector& w, PETScVector const& x,
+ PETScVector const& y)
{
VecPointwiseDivide(w.getRawVector(), x.getRawVector(), y.getRawVector());
}
// Explicit specialization
// Computes the Manhattan norm of x
-template<>
+template <>
double norm1(PETScVector const& x)
{
PetscScalar norm = 0.;
@@ -86,7 +88,7 @@ double norm1(PETScVector const& x)
// Explicit specialization
// Computes the Euclidean norm of x
-template<>
+template <>
double norm2(PETScVector const& x)
{
PetscScalar norm = 0.;
@@ -96,7 +98,7 @@ double norm2(PETScVector const& x)
// Explicit specialization
// Computes the Maximum norm of x
-template<>
+template <>
double normMax(PETScVector const& x)
{
PetscScalar norm = 0.;
@@ -104,7 +106,6 @@ double normMax(PETScVector const& x)
return norm;
}
-
// Matrix
void copy(PETScMatrix const& A, PETScMatrix& B)
@@ -123,8 +124,7 @@ void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X)
{
// TODO check sizes
// TODO sparsity pattern, currently they are assumed to be different (slow)
- MatAYPX(Y.getRawMatrix(), a, X.getRawMatrix(),
- DIFFERENT_NONZERO_PATTERN);
+ MatAYPX(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
}
// Y = a*X + Y
@@ -132,11 +132,9 @@ void axpy(PETScMatrix& Y, double const a, PETScMatrix const& X)
{
// TODO check sizes
// TODO sparsity pattern, currently they are assumed to be different (slow)
- MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(),
- DIFFERENT_NONZERO_PATTERN);
+ MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
}
-
// Matrix and Vector
// v3 = A*v1 + v2
@@ -144,18 +142,21 @@ void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y)
{
// TODO check sizes
assert(&x != &y);
- if (!y.getRawVector()) y.shallowCopy(x);
+ if (!y.getRawVector())
+ y.shallowCopy(x);
MatMult(A.getRawMatrix(), x.getRawVector(), y.getRawVector());
}
// v3 = A*v1 + v2
void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
- PETScVector const& v2, PETScVector& v3)
+ PETScVector const& v2, PETScVector& v3)
{
// TODO check sizes
assert(&v1 != &v3);
- if (!v3.getRawVector()) v3.shallowCopy(v1);
- MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(), v3.getRawVector());
+ if (!v3.getRawVector())
+ v3.shallowCopy(v1);
+ MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(),
+ v3.getRawVector());
}
void finalizeAssembly(PETScMatrix& A)
@@ -168,24 +169,23 @@ void finalizeAssembly(PETScVector& x)
x.finalizeAssembly();
}
-}} // namespaces
-
-
+} // namespace LinAlg
+} // namespace MathLib
-// Sparse global EigenMatrix/EigenVector //////////////////////////////////////////
+// Sparse global EigenMatrix/EigenVector
+// //////////////////////////////////////////
#else
-#include "MathLib/LinAlg/Eigen/EigenVector.h"
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
-namespace MathLib { namespace LinAlg
+namespace MathLib
+{
+namespace LinAlg
{
-
// Vector
-void setLocalAccessibleVector(EigenVector const& /*x*/)
-{
-}
+void setLocalAccessibleVector(EigenVector const& /*x*/) {}
void set(EigenVector& x, double const a)
{
@@ -225,17 +225,17 @@ void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
// Explicit specialization
// Computes w = x/y componentwise.
-template<>
-void componentwiseDivide(EigenVector& w,
- EigenVector const& x, EigenVector const& y)
+template <>
+void componentwiseDivide(EigenVector& w, EigenVector const& x,
+ EigenVector const& y)
{
w.getRawVector().noalias() =
- x.getRawVector().cwiseQuotient(y.getRawVector());
+ x.getRawVector().cwiseQuotient(y.getRawVector());
}
// Explicit specialization
// Computes the Manhattan norm of x
-template<>
+template <>
double norm1(EigenVector const& x)
{
return x.getRawVector().lpNorm<1>();
@@ -243,7 +243,7 @@ double norm1(EigenVector const& x)
// Explicit specialization
// Euclidean norm
-template<>
+template <>
double norm2(EigenVector const& x)
{
return x.getRawVector().norm();
@@ -251,13 +251,12 @@ double norm2(EigenVector const& x)
// Explicit specialization
// Computes the Maximum norm of x
-template<>
+template <>
double normMax(EigenVector const& x)
{
return x.getRawVector().lpNorm<Eigen::Infinity>();
}
-
// Matrix
void copy(EigenMatrix const& A, EigenMatrix& B)
@@ -276,17 +275,16 @@ void scale(EigenMatrix& A, double const a)
void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X)
{
// TODO: does that break anything?
- Y.getRawMatrix() = a*Y.getRawMatrix() + X.getRawMatrix();
+ Y.getRawMatrix() = a * Y.getRawMatrix() + X.getRawMatrix();
}
// Y = a*X + Y
void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X)
{
// TODO: does that break anything?
- Y.getRawMatrix() = a*X.getRawMatrix() + Y.getRawMatrix();
+ Y.getRawMatrix() = a * X.getRawMatrix() + Y.getRawMatrix();
}
-
// Matrix and Vector
// v3 = A*v1 + v2
@@ -297,11 +295,13 @@ void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y)
}
// v3 = A*v1 + v2
-void matMultAdd(EigenMatrix const& A, EigenVector const& v1, EigenVector const& v2, EigenVector& v3)
+void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
+ EigenVector const& v2, EigenVector& v3)
{
assert(&v1 != &v3);
// TODO: does that break anything?
- v3.getRawVector() = v2.getRawVector() + A.getRawMatrix()*v1.getRawVector();
+ v3.getRawVector() =
+ v2.getRawVector() + A.getRawMatrix() * v1.getRawVector();
}
void finalizeAssembly(EigenMatrix& x)
@@ -311,8 +311,8 @@ void finalizeAssembly(EigenMatrix& x)
void finalizeAssembly(EigenVector& /*x*/) {}
-} // namespace LinAlg
+} // namespace LinAlg
-} // namespace MathLib
+} // namespace MathLib
#endif
diff --git a/MathLib/LinAlg/LinAlgEnums.cpp b/MathLib/LinAlg/LinAlgEnums.cpp
index 282433ecb579d76edf2bd36fc5eeb7ff86cfa0cc..4a4331583d7f29e37b6b4398fd6ac1bd0f4910e0 100644
--- a/MathLib/LinAlg/LinAlgEnums.cpp
+++ b/MathLib/LinAlg/LinAlgEnums.cpp
@@ -14,19 +14,22 @@
namespace MathLib
{
-
std::string convertVecNormTypeToString(VecNormType normType)
{
switch (normType)
{
- case VecNormType::NORM1: return "NORM1";
- case VecNormType::NORM2: return "NORM2";
- case VecNormType::INFINITY_N: return "INFINITY_N";
- default: return "INVALID";
+ case VecNormType::NORM1:
+ return "NORM1";
+ case VecNormType::NORM2:
+ return "NORM2";
+ case VecNormType::INFINITY_N:
+ return "INFINITY_N";
+ default:
+ return "INVALID";
}
}
-VecNormType convertStringToVecNormType(const std::string &str)
+VecNormType convertStringToVecNormType(const std::string& str)
{
if (str == "NORM1")
{
@@ -43,4 +46,4 @@ VecNormType convertStringToVecNormType(const std::string &str)
return VecNormType::INVALID;
}
-} // end namespace MathLib
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/LinearSolverOptions.cpp b/MathLib/LinAlg/LinearSolverOptions.cpp
index fc1e1b9a6341680b066f90d28a9113a9f112536f..2bb3381ed8743572390e4514f4a27e6c6936bb9b 100644
--- a/MathLib/LinAlg/LinearSolverOptions.cpp
+++ b/MathLib/LinAlg/LinearSolverOptions.cpp
@@ -5,19 +5,15 @@
//! Configuration tag names of all known linear solvers for their
//! configuration in the project file.
//! Add your tag name here when you add a new solver.
-static
-std::set<std::string>
-known_linear_solvers { "eigen", "lis", "petsc" };
+static std::set<std::string> known_linear_solvers{"eigen", "lis", "petsc"};
namespace MathLib
{
-
-
-void
-ignoreOtherLinearSolvers(const BaseLib::ConfigTree &config,
- const std::string &solver_name)
+void ignoreOtherLinearSolvers(const BaseLib::ConfigTree& config,
+ const std::string& solver_name)
{
- for (auto const& s : known_linear_solvers) {
+ for (auto const& s : known_linear_solvers)
+ {
if (s != solver_name)
{
config.ignoreConfigParameter(s);
diff --git a/MathLib/LinAlg/Lis/LisLinearSolver.cpp b/MathLib/LinAlg/Lis/LisLinearSolver.cpp
index c8a68d174b7ad192090ed279f1d1a6effa963a49..ff0548db2c6a27c319c2004567f3446342ca4001 100644
--- a/MathLib/LinAlg/Lis/LisLinearSolver.cpp
+++ b/MathLib/LinAlg/Lis/LisLinearSolver.cpp
@@ -16,25 +16,21 @@
#include <omp.h>
#endif
-#include "LisLinearSolver.h"
-
#include "BaseLib/Logging.h"
-
#include "LisCheck.h"
+#include "LisLinearSolver.h"
#include "LisMatrix.h"
#include "LisVector.h"
namespace MathLib
{
-
-LisLinearSolver::LisLinearSolver(
- const std::string /*solver_name*/,
- const BaseLib::ConfigTree* const option)
-: _lis_option(option)
+LisLinearSolver::LisLinearSolver(const std::string /*solver_name*/,
+ const BaseLib::ConfigTree* const option)
+ : _lis_option(option)
{
}
-bool LisLinearSolver::solve(LisMatrix &A, LisVector &b, LisVector &x)
+bool LisLinearSolver::solve(LisMatrix& A, LisVector& b, LisVector& x)
{
finalizeMatrixAssembly(A);
@@ -47,8 +43,8 @@ bool LisLinearSolver::solve(LisMatrix &A, LisVector &b, LisVector &x)
if (!checkLisError(ierr))
return false;
- lis_solver_set_option(
- const_cast<char*>(_lis_option._option_string.c_str()), solver);
+ lis_solver_set_option(const_cast<char*>(_lis_option._option_string.c_str()),
+ solver);
#ifdef _OPENMP
INFO("-> number of threads: {:i}", (int)omp_get_max_threads());
#endif
@@ -69,7 +65,8 @@ bool LisLinearSolver::solve(LisMatrix &A, LisVector &b, LisVector &x)
// solve
INFO("-> solve");
- ierr = lis_solve(A.getRawMatrix(), b.getRawVector(), x.getRawVector(), solver);
+ ierr =
+ lis_solve(A.getRawMatrix(), b.getRawVector(), x.getRawVector(), solver);
if (!checkLisError(ierr))
return false;
@@ -97,8 +94,8 @@ bool LisLinearSolver::solve(LisMatrix &A, LisVector &b, LisVector &x)
}
{
double time, itime, ptime, p_ctime, p_itime;
- ierr = lis_solver_get_timeex(solver, &time, &itime,
- &ptime, &p_ctime, &p_itime);
+ ierr = lis_solver_get_timeex(solver, &time, &itime, &ptime, &p_ctime,
+ &p_itime);
if (!checkLisError(ierr))
return false;
INFO("-> time total (s): {:g}", time);
@@ -117,4 +114,4 @@ bool LisLinearSolver::solve(LisMatrix &A, LisVector &b, LisVector &x)
return linear_solver_status == LIS_SUCCESS;
}
-} //MathLib
+} // namespace MathLib
diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
index e8c927aad2cfdc749db5e6302999ee3e719e728b..a1af8ba42158febbe482c02607298d2decd8d877 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.cpp
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -17,16 +17,17 @@
#include <cmath>
#include <cstdlib>
-
#include "BaseLib/Error.h"
-#include "LisVector.h"
#include "LisCheck.h"
+#include "LisVector.h"
namespace MathLib
{
-
LisMatrix::LisMatrix(std::size_t n_rows, MatrixType mat_type)
- : _n_rows(n_rows), _mat_type(mat_type), _is_assembled(false), _use_external_arrays(false)
+ : _n_rows(n_rows),
+ _mat_type(mat_type),
+ _is_assembled(false),
+ _use_external_arrays(false)
{
int ierr = lis_matrix_create(0, &_AA);
checkLisError(ierr);
@@ -37,8 +38,8 @@ LisMatrix::LisMatrix(std::size_t n_rows, MatrixType mat_type)
checkLisError(ierr);
}
-LisMatrix::LisMatrix(std::size_t n_rows, int nnz, IndexType *row_ptr,
- IndexType *col_idx, double *data)
+LisMatrix::LisMatrix(std::size_t n_rows, int nnz, IndexType* row_ptr,
+ IndexType* col_idx, double* data)
: _n_rows(n_rows),
_mat_type(MatrixType::CRS),
_is_assembled(false),
@@ -74,8 +75,8 @@ LisMatrix::~LisMatrix()
void LisMatrix::setZero()
{
- // A matrix has to be destroyed and created again because Lis doesn't provide a
- // function to set matrix entries to zero
+ // A matrix has to be destroyed and created again because Lis doesn't
+ // provide a function to set matrix entries to zero
int ierr = lis_matrix_destroy(_AA);
checkLisError(ierr);
ierr = lis_matrix_create(0, &_AA);
@@ -90,9 +91,10 @@ void LisMatrix::setZero()
int LisMatrix::setValue(IndexType rowId, IndexType colId, double v)
{
- if (v == 0.0) return 0;
+ if (v == 0.0)
+ return 0;
lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
- if (rowId==colId)
+ if (rowId == colId)
lis_vector_set_value(LIS_INS_VALUE, rowId, v, _diag);
_is_assembled = false;
return 0;
@@ -100,15 +102,16 @@ int LisMatrix::setValue(IndexType rowId, IndexType colId, double v)
int LisMatrix::add(IndexType rowId, IndexType colId, double v)
{
- if (v == 0.0) return 0;
+ if (v == 0.0)
+ return 0;
lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
- if (rowId==colId)
+ if (rowId == colId)
lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _diag);
_is_assembled = false;
return 0;
}
-void LisMatrix::write(const std::string &filename) const
+void LisMatrix::write(const std::string& filename) const
{
if (!_is_assembled)
{
@@ -117,12 +120,14 @@ void LisMatrix::write(const std::string &filename) const
lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str()));
}
-bool finalizeMatrixAssembly(LisMatrix &mat)
+bool finalizeMatrixAssembly(LisMatrix& mat)
{
- LIS_MATRIX &A = mat.getRawMatrix();
+ LIS_MATRIX& A = mat.getRawMatrix();
- if (!mat.isAssembled()) {
- int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
+ if (!mat.isAssembled())
+ {
+ int ierr =
+ lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
checkLisError(ierr);
ierr = lis_matrix_assemble(A);
checkLisError(ierr);
@@ -131,5 +136,4 @@ bool finalizeMatrixAssembly(LisMatrix &mat)
return true;
}
-
-} //MathLib
+} // namespace MathLib
diff --git a/MathLib/LinAlg/Lis/LisVector.cpp b/MathLib/LinAlg/Lis/LisVector.cpp
index 37f758b388f1d4a3b452b6c3fb9524bcebdc5c2a..f5b68251bc9b47dc355896a9a4f9e5375bb20c7b 100644
--- a/MathLib/LinAlg/Lis/LisVector.cpp
+++ b/MathLib/LinAlg/Lis/LisVector.cpp
@@ -12,9 +12,10 @@
*
*/
+#include "LisVector.h"
+
#include <cassert>
-#include "LisVector.h"
#include "LisCheck.h"
namespace MathLib
@@ -59,4 +60,4 @@ void LisVector::write(const std::string& filename) const
lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename));
}
-} // MathLib
+} // namespace MathLib
diff --git a/MathLib/LinAlg/MatrixVectorTraits.cpp b/MathLib/LinAlg/MatrixVectorTraits.cpp
index ff52b18cccd83d2a3b2fc248f082fdb1c6aee6e2..4d6a334b6c9a2e729848fdaa33793ee2c1f83e29 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.cpp
+++ b/MathLib/LinAlg/MatrixVectorTraits.cpp
@@ -9,30 +9,26 @@
*/
#include "MatrixVectorTraits.h"
+
#include "MatrixSpecifications.h"
#ifdef USE_PETSC
namespace MathLib
{
-
-std::unique_ptr<PETScMatrix>
-MatrixVectorTraits<PETScMatrix>::
-newInstance()
+std::unique_ptr<PETScMatrix> MatrixVectorTraits<PETScMatrix>::newInstance()
{
return std::make_unique<PETScMatrix>();
}
-std::unique_ptr<PETScMatrix>
-MatrixVectorTraits<PETScMatrix>::
-newInstance(PETScMatrix const& A)
+std::unique_ptr<PETScMatrix> MatrixVectorTraits<PETScMatrix>::newInstance(
+ PETScMatrix const& A)
{
return std::make_unique<PETScMatrix>(A);
}
-std::unique_ptr<PETScMatrix>
-MatrixVectorTraits<PETScMatrix>::
-newInstance(MatrixSpecifications const& spec)
+std::unique_ptr<PETScMatrix> MatrixVectorTraits<PETScMatrix>::newInstance(
+ MatrixSpecifications const& spec)
{
auto const nrows = spec.nrows;
auto const ncols = spec.ncols;
@@ -54,30 +50,29 @@ newInstance(MatrixSpecifications const& spec)
return std::make_unique<PETScMatrix>(nrows, ncols);
}
-std::unique_ptr<PETScVector>
-MatrixVectorTraits<PETScVector>::
-newInstance()
+std::unique_ptr<PETScVector> MatrixVectorTraits<PETScVector>::newInstance()
{
return std::make_unique<PETScVector>();
}
-std::unique_ptr<PETScVector>
-MatrixVectorTraits<PETScVector>::
-newInstance(PETScVector const& x)
+std::unique_ptr<PETScVector> MatrixVectorTraits<PETScVector>::newInstance(
+ PETScVector const& x)
{
return std::make_unique<PETScVector>(x);
}
-std::unique_ptr<PETScVector>
-MatrixVectorTraits<PETScVector>::
-newInstance(MatrixSpecifications const& spec)
+std::unique_ptr<PETScVector> MatrixVectorTraits<PETScVector>::newInstance(
+ MatrixSpecifications const& spec)
{
auto const is_global_size = false;
- if (spec.ghost_indices != nullptr) {
+ if (spec.ghost_indices != nullptr)
+ {
return std::make_unique<PETScVector>(spec.nrows, *spec.ghost_indices,
is_global_size);
- } else {
+ }
+ else
+ {
return std::make_unique<PETScVector>(spec.nrows, is_global_size);
}
}
@@ -89,31 +84,25 @@ std::unique_ptr<PETScVector> MatrixVectorTraits<PETScVector>::newInstance(
return std::make_unique<PETScVector>(length, is_global_size);
}
-} // namespace MathLib
-
+} // namespace MathLib
#else
namespace MathLib
{
-
-std::unique_ptr<EigenMatrix>
-MatrixVectorTraits<EigenMatrix>::
-newInstance()
+std::unique_ptr<EigenMatrix> MatrixVectorTraits<EigenMatrix>::newInstance()
{
return std::make_unique<EigenMatrix>(0, 0); // TODO default constructor
}
-std::unique_ptr<EigenMatrix>
-MatrixVectorTraits<EigenMatrix>::
-newInstance(EigenMatrix const& A)
+std::unique_ptr<EigenMatrix> MatrixVectorTraits<EigenMatrix>::newInstance(
+ EigenMatrix const& A)
{
return std::make_unique<EigenMatrix>(A);
}
-std::unique_ptr<EigenMatrix>
-MatrixVectorTraits<EigenMatrix>::
-newInstance(MatrixSpecifications const& spec)
+std::unique_ptr<EigenMatrix> MatrixVectorTraits<EigenMatrix>::newInstance(
+ MatrixSpecifications const& spec)
{
auto A = std::make_unique<EigenMatrix>(spec.nrows);
@@ -125,23 +114,19 @@ newInstance(MatrixSpecifications const& spec)
return A;
}
-std::unique_ptr<EigenVector>
-MatrixVectorTraits<EigenVector>::
-newInstance()
+std::unique_ptr<EigenVector> MatrixVectorTraits<EigenVector>::newInstance()
{
return std::make_unique<EigenVector>();
}
-std::unique_ptr<EigenVector>
-MatrixVectorTraits<EigenVector>::
-newInstance(EigenVector const& x)
+std::unique_ptr<EigenVector> MatrixVectorTraits<EigenVector>::newInstance(
+ EigenVector const& x)
{
return std::make_unique<EigenVector>(x);
}
-std::unique_ptr<EigenVector>
-MatrixVectorTraits<EigenVector>::
-newInstance(MatrixSpecifications const& spec)
+std::unique_ptr<EigenVector> MatrixVectorTraits<EigenVector>::newInstance(
+ MatrixSpecifications const& spec)
{
return std::make_unique<EigenVector>(spec.nrows);
}
@@ -151,6 +136,6 @@ std::unique_ptr<EigenVector> MatrixVectorTraits<EigenVector>::newInstance(
{
return std::make_unique<EigenVector>(length);
}
-} // namespace MathLib
+} // namespace MathLib
#endif
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index 2445daa1d456ac7b2cb91332b73e6dea7ba0dab1..80fe3b42588a3bb80c49b1a88404c4f2c0ea8c16 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -155,4 +155,4 @@ bool finalizeMatrixAssembly(PETScMatrix& mat, const MatAssemblyType asm_type)
return true;
}
-} // end of namespace
+} // namespace MathLib
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index 34f908d3150afb96389b6154b5241f95ddbfe69b..4de8178a9dd9d0c8c40d46fe4ca51336e9425462 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -323,4 +323,4 @@ void finalizeVectorAssembly(PETScVector& vec)
vec.finalizeAssembly();
}
-} // end of namespace
+} // namespace MathLib
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.cpp b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
index 18985bf423958adccafdf34d4382d6d0591aea8b..2678cbcb5d319bc4ff3e84e058d927eaab44d367 100644
--- a/MathLib/LinAlg/UnifiedMatrixSetters.cpp
+++ b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
@@ -8,21 +8,22 @@
*
*/
-#include <cassert>
#include "UnifiedMatrixSetters.h"
+#include <cassert>
+
#ifdef USE_PETSC
// Global PETScMatrix/PETScVector //////////////////////////////////////////
#include <numeric>
-#include "MathLib/LinAlg/PETSc/PETScVector.h"
+
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+#include "MathLib/LinAlg/PETSc/PETScVector.h"
namespace MathLib
{
-void setVector(PETScVector& v,
- std::initializer_list<double> values)
+void setVector(PETScVector& v, std::initializer_list<double> values)
{
std::vector<double> const vals(values);
std::vector<PETScVector::IndexType> idcs(vals.size());
@@ -31,14 +32,14 @@ void setVector(PETScVector& v,
v.set(idcs, vals);
}
-void setVector(PETScVector& v, MatrixVectorTraits<PETScVector>::Index const index,
+void setVector(PETScVector& v,
+ MatrixVectorTraits<PETScVector>::Index const index,
double const value)
{
- v.set(index, value); // TODO handle negative indices
+ v.set(index, value); // TODO handle negative indices
}
-void setMatrix(PETScMatrix& m,
- std::initializer_list<double> values)
+void setMatrix(PETScMatrix& m, std::initializer_list<double> values)
{
m.setZero();
addToMatrix(m, values);
@@ -61,26 +62,29 @@ void setMatrix(PETScMatrix& m, Eigen::MatrixXd const& tmp)
std::iota(col_idcs.begin(), col_idcs.end(), 0);
// PETSc wants row-major
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp_ = tmp;
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+ tmp_ = tmp;
m.add(row_idcs, col_idcs, tmp_);
}
-void addToMatrix(PETScMatrix& m,
- std::initializer_list<double> values)
+void addToMatrix(PETScMatrix& m, std::initializer_list<double> values)
{
using IndexType = PETScMatrix::IndexType;
auto const rows = m.getNumberOfRows();
auto const cols = m.getNumberOfColumns();
- assert((IndexType) values.size() == rows*cols);
+ assert((IndexType)values.size() == rows * cols);
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(rows, cols);
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(
+ rows, cols);
auto it = values.begin();
- for (IndexType r=0; r<rows; ++r) {
- for (IndexType c=0; c<cols; ++c) {
+ for (IndexType r = 0; r < rows; ++r)
+ {
+ for (IndexType c = 0; c < cols; ++c)
+ {
tmp(r, c) = *(it++);
}
}
@@ -94,21 +98,19 @@ void addToMatrix(PETScMatrix& m,
m.add(row_idcs, col_idcs, tmp);
}
-} // namespace MathLib
-
+} // namespace MathLib
#else
-// Sparse global EigenMatrix/EigenVector //////////////////////////////////////////
+// Sparse global EigenMatrix/EigenVector
+// //////////////////////////////////////////
-#include "MathLib/LinAlg/Eigen/EigenVector.h"
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
namespace MathLib
{
-
-void setVector(EigenVector& v_,
- std::initializer_list<double> values)
+void setVector(EigenVector& v_, std::initializer_list<double> values)
{
auto& v(v_.getRawVector());
assert((std::size_t)v.size() == values.size());
@@ -119,25 +121,26 @@ void setVector(EigenVector& v_,
}
}
-void setVector(EigenVector& v, MatrixVectorTraits<EigenVector>::Index const index,
+void setVector(EigenVector& v,
+ MatrixVectorTraits<EigenVector>::Index const index,
double const value)
{
v.getRawVector()[index] = value;
}
-
-void setMatrix(EigenMatrix& m,
- std::initializer_list<double> values)
+void setMatrix(EigenMatrix& m, std::initializer_list<double> values)
{
auto const rows = m.getNumberOfRows();
auto const cols = m.getNumberOfColumns();
- assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows*cols);
+ assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows * cols);
Eigen::MatrixXd tmp(rows, cols);
auto it = values.begin();
- for (GlobalIndexType r=0; r<rows; ++r) {
- for (GlobalIndexType c=0; c<cols; ++c) {
+ for (GlobalIndexType r = 0; r < rows; ++r)
+ {
+ for (GlobalIndexType c = 0; c < cols; ++c)
+ {
tmp(r, c) = *(it++);
}
}
@@ -150,18 +153,20 @@ void setMatrix(EigenMatrix& m, Eigen::MatrixXd const& tmp)
m.getRawMatrix() = tmp.sparseView();
}
-void addToMatrix(EigenMatrix& m,
- std::initializer_list<double> values)
+void addToMatrix(EigenMatrix& m, std::initializer_list<double> values)
{
auto const rows = m.getNumberOfRows();
auto const cols = m.getNumberOfColumns();
- assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows*cols);
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(rows, cols);
+ assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows * cols);
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(
+ rows, cols);
auto it = values.begin();
- for (GlobalIndexType r=0; r<rows; ++r) {
- for (GlobalIndexType c=0; c<cols; ++c) {
+ for (GlobalIndexType r = 0; r < rows; ++r)
+ {
+ for (GlobalIndexType c = 0; c < cols; ++c)
+ {
tmp(r, c) = *(it++);
}
}
@@ -169,6 +174,6 @@ void addToMatrix(EigenMatrix& m,
m.getRawMatrix() += tmp.sparseView();
}
-} // namespace MathLib
+} // namespace MathLib
#endif
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index cb5b3d6f3c16e2018d3ac5491d87b57fef0cab4d..503445b3426579157a8e75ad6001e86781549451 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -16,7 +16,6 @@
namespace MathLib
{
-
double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
Point3d const& pb, double& lambda, double& d0)
{
diff --git a/MathLib/ODE/CVodeSolver.cpp b/MathLib/ODE/CVodeSolver.cpp
index c9a885f6ce50b9ba1d9b3f6326d2182a154c14a1..266b2e5d194eb5325ae136a74d1fcdcd9ef659ef 100644
--- a/MathLib/ODE/CVodeSolver.cpp
+++ b/MathLib/ODE/CVodeSolver.cpp
@@ -9,17 +9,17 @@
*/
#include "CVodeSolver.h"
-#include <cassert>
-#include "BaseLib/Logging.h"
-
#include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
-#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
#include <cvode/cvode_dense.h> /* prototype for CVDense */
+#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
#include <sundials/sundials_types.h> /* definition of type realtype */
+#include <cassert>
+
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Error.h"
+#include "BaseLib/Logging.h"
//! \addtogroup ExternalODESolverInterface
//! @{
@@ -132,8 +132,9 @@ CVodeSolverImpl::CVodeSolverImpl(const BaseLib::ConfigTree& config,
const unsigned num_equations)
{
if (auto const param =
- //! \ogs_file_param{ode_solver__CVODE__linear_multistep_method}
- config.getConfigParameterOptional<std::string>("linear_multistep_method"))
+ //! \ogs_file_param{ode_solver__CVODE__linear_multistep_method}
+ config.getConfigParameterOptional<std::string>(
+ "linear_multistep_method"))
{
DBUG("setting linear multistep method (config: {:s})", param->c_str());
@@ -152,8 +153,9 @@ CVodeSolverImpl::CVodeSolverImpl(const BaseLib::ConfigTree& config,
}
if (auto const param =
- //! \ogs_file_param{ode_solver__CVODE__nonlinear_solver_iteration}
- config.getConfigParameterOptional<std::string>("nonlinear_solver_iteration"))
+ //! \ogs_file_param{ode_solver__CVODE__nonlinear_solver_iteration}
+ config.getConfigParameterOptional<std::string>(
+ "nonlinear_solver_iteration"))
{
DBUG("setting nonlinear solver iteration (config: {:s})",
param->c_str());
@@ -186,8 +188,7 @@ CVodeSolverImpl::CVodeSolverImpl(const BaseLib::ConfigTree& config,
}
auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
- void* function_handles) -> int
- {
+ void* function_handles) -> int {
bool successful =
static_cast<detail::FunctionHandles*>(function_handles)
->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
@@ -257,8 +258,7 @@ void CVodeSolverImpl::preSolve()
const N_Vector ydot, const DlsMat jac,
void* function_handles, N_Vector /*tmp1*/,
N_Vector /*tmp2*/, N_Vector /*tmp3*/
- ) -> int
- {
+ ) -> int {
(void)N; // prevent warnings during non-debug build
auto* fh = static_cast<detail::FunctionHandles*>(function_handles);
assert(N == fh->getNumberOfEquations());
diff --git a/MathLib/Point3d.cpp b/MathLib/Point3d.cpp
index 13873692de5549bc3541143519de685ae19533d6..9e131e186a6cb0da9020cc3c0795601f78f6f36e 100644
--- a/MathLib/Point3d.cpp
+++ b/MathLib/Point3d.cpp
@@ -12,5 +12,5 @@
namespace MathLib
{
-extern const Point3d ORIGIN{ {{0.0, 0.0, 0.0}} };
+extern const Point3d ORIGIN{{{0.0, 0.0, 0.0}}};
}