diff --git a/cmake/modules/AmdisMacros.cmake b/cmake/modules/AmdisMacros.cmake
index 3ee268187aff02a90ce4e3ae91ba3fb0e56b64a6..0862c33eb4d36324aa3ec47feedf42145640ac77 100644
--- a/cmake/modules/AmdisMacros.cmake
+++ b/cmake/modules/AmdisMacros.cmake
@@ -1,21 +1,29 @@
-# File for module specific CMake tests.
-
include(AMDiSCXXFeatures)
# some additional packages and flags
-find_package(MTL REQUIRED
- PATHS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4)
+find_package(MTL
+ PATHS ${MTL_DIR}
+ HINTS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4
+)
+
+if (MTL_FOUND)
+ set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
+ set(MTL_COMPILE_DEFINITIONS "")
+ foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
+ list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
+ endforeach()
-set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
-set(MTL_COMPILE_DEFINITIONS "")
-foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
- list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
-endforeach()
+ find_package(SuiteSparse QUIET)
+ if (SuiteSparse_FOUND)
+ list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
+ endif (SuiteSparse_FOUND)
+endif (MTL_FOUND)
-if(HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
- list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
-endif()
+set(HAVE_MTL MTL_FOUND)
-dune_register_package_flags(
- COMPILE_DEFINITIONS ${MTL_COMPILE_DEFINITIONS}
- INCLUDE_DIRS ${MTL_INCLUDE_DIRS})
+if (MTL_FOUND)
+ list(APPEND MTL_COMPILE_DEFINITIONS "ENABLE_MTL=1")
+ dune_register_package_flags(
+ COMPILE_DEFINITIONS ${MTL_COMPILE_DEFINITIONS}
+ INCLUDE_DIRS ${MTL_INCLUDE_DIRS})
+endif (MTL_FOUND)
\ No newline at end of file
diff --git a/config.h.cmake b/config.h.cmake
index e53392bc7ccd643d244e463b52e0df06dcb2fa1b..23a127dba08ec64ebac14df6d1a68429d9f5bbfe 100644
--- a/config.h.cmake
+++ b/config.h.cmake
@@ -1,4 +1,4 @@
-/* begin dune-amdis
+/* begin amdis
put the definitions for config.h specific to
your project here. Everything above will be
overwritten
@@ -40,12 +40,13 @@
/* Define to the revision of amdis */
#define AMDIS_VERSION_REVISION @AMDIS_VERSION_REVISION@
+/* Define to ENABLE_UMFPACK if the MTL library is available */
+#cmakedefine HAVE_MTL ENABLE_MTL
/* some detected compiler features may be used in AMDiS */
#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
-
#cmakedefine AMDIS_HAS_CXX_CONSTEXPR_IF 1
-/* end dune-amdis
+/* end amdis
Everything below here will be overwritten
*/
diff --git a/examples/convection_diffusion.cc b/examples/convection_diffusion.cc
index 40871822f4fe7a9483e404fe2ab8bb8f38dc6124..14bea3a42cbf8d4e46213281d25ec3d3236cbbb2 100644
--- a/examples/convection_diffusion.cc
+++ b/examples/convection_diffusion.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
#include <iostream>
#include <amdis/AMDiS.hpp>
@@ -26,13 +26,13 @@ int main(int argc, char** argv)
// -div(A*grad(u)) + div(b*u) + c*u = f
auto opCD = convectionDiffusion(/*A=*/1.0, /*b=*/0.0, /*c=*/1.0, /*f=*/1.0);
- prob.addMatrixOperator(opCD, _0, _0);
- prob.addVectorOperator(opCD, _0);
+ prob.addMatrixOperator(opCD, 0, 0);
+ prob.addVectorOperator(opCD, 0);
// set boundary condition
auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
auto dbcValues = [](auto const& x){ return 0.0; }; // set value
- prob.addDirichletBC(predicate, _0, _0, dbcValues);
+ prob.addDirichletBC(predicate, 0, 0, dbcValues);
AdaptInfo adaptInfo("adapt");
prob.buildAfterCoarsen(adaptInfo, Flag(0));
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 396737c3643a17662486482fe66a5fc45f540577..fc7830e8fee7e0b9b46ec87198a53e30593fb627 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
#include <iostream>
#include <fmt/core.h>
@@ -15,7 +15,7 @@ using namespace AMDiS;
using namespace Dune::Indices;
// 1 component with polynomial degree 1
-using Param = YaspGridBasis<AMDIS_DIM, 1>;
+using Param = YaspGridBasis<AMDIS_DIM, 2>;
using ElliptProblem = ProblemStat<Param>;
int main(int argc, char** argv)
diff --git a/examples/heat.cc b/examples/heat.cc
index d232264159f3ca35290fbaf364d9b3c06b032400..c2fa5e08d1c7d8afbd2fa611393bed74afa48822 100644
--- a/examples/heat.cc
+++ b/examples/heat.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
#include <iostream>
#include <amdis/AMDiS.hpp>
diff --git a/examples/init/ellipt.dat.2d b/examples/init/ellipt.dat.2d
index f12565729e95e06f7902cf81744b29450c3445e1..0ef4d635178cbd825b31cbd818eab17a036d5c72 100644
--- a/examples/init/ellipt.dat.2d
+++ b/examples/init/ellipt.dat.2d
@@ -3,10 +3,10 @@ elliptMesh->global refinements: 0
ellipt->mesh: elliptMesh
ellipt->solver->name: bicgstab
-ellipt->solver->max iteration: 1000
-ellipt->solver->relative tolerance: 1.e-8
-ellipt->solver->info: 10
-ellipt->solver->left precon: diag
+ellipt->solver->max iteration: 10000
+ellipt->solver->relative tolerance: 1.e-10
+ellipt->solver->info: 1
+ellipt->solver->left precon: ilu
ellipt->solver->right precon: no
ellipt->output[0]->filename: ellipt.2d
diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 4fcf941ba26aa0d16f10dc0ef26832dc01a881d7..762b369deba45fa66ffe002b519cb6bebf03a7ca 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <iostream>
#include <ctime>
#include <cmath>
@@ -41,8 +45,8 @@ int main(int argc, char** argv)
AdaptInfo adaptInfo("adapt");
// tree-paths for components
- auto _v = 0_c;
- auto _p = 1_c;
+ auto _v = Dune::Indices::_0;
+ auto _p = Dune::Indices::_1;
// <1/tau * u, v>
auto opTime = makeOperator(tag::testvec_trialvec{}, density);
diff --git a/examples/stokes0.cc b/examples/stokes0.cc
index 7118684a3780769425f17ceac524c369e6335ecb..40a8c425474d611be7ca3d5845f55270f2055eab 100644
--- a/examples/stokes0.cc
+++ b/examples/stokes0.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <iostream>
#include <ctime>
#include <cmath>
@@ -30,8 +34,8 @@ int main(int argc, char** argv)
Parameters::get("stokes->viscosity", viscosity);
// tree-paths for components
- auto _v = 0_c;
- auto _p = 1_c;
+ auto _v = Dune::Indices::_0;
+ auto _p = Dune::Indices::_1;
// <viscosity*grad(u_i), grad(v_i)>
for (std::size_t i = 0; i < DOW; ++i) {
@@ -76,10 +80,12 @@ int main(int argc, char** argv)
// assemble and solve system
prob.buildAfterCoarsen(adaptInfo, Flag(0));
+#ifdef DEBUG_MTL
// write matrix to file
mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
- out << prob.getSystemMatrix().getMatrix();
- std::cout << prob.getSystemMatrix().getMatrix() << '\n';
+ out << prob.getSystemMatrix().matrix();
+ std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
prob.solve(adaptInfo);
diff --git a/examples/stokes1.cc b/examples/stokes1.cc
index e0effccebad0af27d55163287295bf5a13eaa598..bdb95a3b5c11de137c2ba754850a7c708ad38fc0 100644
--- a/examples/stokes1.cc
+++ b/examples/stokes1.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <iostream>
#include <ctime>
#include <cmath>
@@ -30,8 +34,8 @@ int main(int argc, char** argv)
Parameters::get("stokes->viscosity", viscosity);
// tree-paths for components
- auto _v = 0_c;
- auto _p = 1_c;
+ auto _v = Dune::Indices::_0;
+ auto _p = Dune::Indices::_1;
// <viscosity*grad(u_i), grad(v_i)>
for (std::size_t i = 0; i < DOW; ++i) {
@@ -77,10 +81,12 @@ int main(int argc, char** argv)
// assemble and solve system
prob.buildAfterCoarsen(adaptInfo, Flag(0));
+#ifdef DEBUG_MTL
// write matrix to file
mtl::io::matrix_market_ostream out("matrix_stokes1.mtx");
- out << prob.getSystemMatrix().getMatrix();
- std::cout << prob.getSystemMatrix().getMatrix() << '\n';
+ out << prob.getSystemMatrix().matrix();
+ std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
prob.solve(adaptInfo);
diff --git a/examples/stokes3.cc b/examples/stokes3.cc
index 54c75d0da556897759920a782f3cf0b029875902..44ccdd6a4db128ab582a7907dd4296ad7d7172c8 100644
--- a/examples/stokes3.cc
+++ b/examples/stokes3.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <iostream>
#include <ctime>
#include <cmath>
@@ -24,8 +28,8 @@ int main(int argc, char** argv)
Parameters::get("stokes->viscosity", viscosity);
// tree-paths for components
- auto _v = 0_c;
- auto _p = 1_c;
+ auto _v = Dune::Indices::_0;
+ auto _p = Dune::Indices::_1;
auto opStokes = makeOperator(tag::stokes{}, viscosity);
prob.addMatrixOperator(opStokes, treepath(), treepath());
diff --git a/examples/vecellipt.cc b/examples/vecellipt.cc
index e87f870d172ba675e8d7a49dbefce7e77d6cfe7a..a4dfb40735a29230dbad1e02bedb43c9e4cf8927 100644
--- a/examples/vecellipt.cc
+++ b/examples/vecellipt.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <iostream>
#include <amdis/AMDiS.hpp>
@@ -38,13 +42,15 @@ int main(int argc, char** argv)
prob.buildAfterCoarsen(adaptInfo, Flag(0));
+#ifdef DEBUG_MTL
// write matrix to file
if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 4) {
mtl::io::matrix_market_ostream out("matrix.mtx");
- out << prob.getSystemMatrix().getMatrix();
+ out << prob.getSystemMatrix().matrix();
- std::cout << prob.getSystemMatrix().getMatrix() << '\n';
+ std::cout << prob.getSystemMatrix().matrix() << '\n';
}
+#endif
prob.solve(adaptInfo);
prob.writeFiles(adaptInfo, true);
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index ef99fb52f1cd519f54b83db8a0e34dfd3c7c219c..e487c99f6b2814ea3195cd9de0ca5da46fa843f1 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -3,14 +3,9 @@
#include <memory>
#include <tuple>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-
#include <amdis/DirichletBC.hpp>
-#include <amdis/LinearAlgebra.hpp>
#include <amdis/LocalAssemblerList.hpp>
#include <amdis/common/Mpl.hpp>
-#include <amdis/common/TypeDefs.hpp>
namespace AMDiS
{
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index f073f40c9c90ddb077b4ad20bc14ac223cc126c3..2854493e03077196120f629a32ecb0fccbca0308 100644
--- a/src/amdis/Assembler.inc.hpp
+++ b/src/amdis/Assembler.inc.hpp
@@ -1,5 +1,7 @@
#pragma once
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <amdis/utility/TreePath.hpp>
#include <amdis/utility/Visitor.hpp>
@@ -24,14 +26,14 @@ void Assembler<Traits>::assemble(
// 2. create a local matrix and vector
std::size_t localSize = localView.maxSize();
- Impl::ElementMatrix elementMatrix(localSize, localSize);
- Impl::ElementVector elementVector(localSize);
+ Dune::DynamicMatrix<double> elementMatrix(localSize, localSize);
+ Dune::DynamicVector<double> elementVector(localSize);
// 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(globalBasis_.gridView()))
{
- set_to_zero(elementMatrix);
- set_to_zero(elementVector);
+ elementMatrix = 0;
+ elementVector = 0;
localView.bind(element);
auto geometry = element.geometry();
@@ -75,24 +77,9 @@ void Assembler<Traits>::assemble(
});
});
- // add element-matrix to system-matrix
- for (std::size_t i = 0; i < localView.size(); ++i) {
- auto const row = localView.index(i);
- for (std::size_t j = 0; j < localView.size(); ++j) {
- if (std::abs(elementMatrix(i,j)) > threshold<double>) {
- auto const col = localView.index(j);
- matrix(row,col) += elementMatrix(i,j);
- }
- }
- }
-
- // add element-vector to system-vector
- for (std::size_t i = 0; i < localView.size(); ++i) {
- if (std::abs(elementVector[i]) > threshold<double>) {
- auto const idx = localView.index(i);
- rhs[idx] += elementVector[i];
- }
- }
+ // add element-matrix to system-matrix and element-vector to rhs
+ matrix.insert(localView, localView, elementMatrix);
+ rhs.insert(localView, elementVector);
// unbind all operators
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
@@ -203,7 +190,7 @@ std::size_t Assembler<Traits>::finishMatrixVector(
});
});
- return matrix.getNnz();
+ return matrix.nnz();
}
} // end namespace AMDiS
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index 34a281e1c992a40192133583403bac3484a1fa91..30960d343be66ef19a143304054ef5847f032d0c 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -11,12 +11,14 @@
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ValueCategory.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
+#include <amdis/linear_algebra/HierarchicWrapper.hpp>
#include <amdis/utility/RangeType.hpp>
#include <amdis/utility/TreeData.hpp>
namespace AMDiS
{
+ struct BoundaryType { int b; };
+
/// Implements a boundary condition of Dirichlet-type.
/**
* By calling the methods \ref init() and \ref finish before and after
@@ -82,13 +84,13 @@ namespace AMDiS
using Dune::Functions::interpolate;
Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
[&](auto id) {
- auto rhsWrapper = wrapper(rhs.getVector());
+ auto rhsWrapper = hierarchicVectorWrapper(rhs.vector());
interpolate(id(rowBasis), rhsWrapper, values_, dirichletNodes_);
});
Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
[&](auto id) {
- auto solutionWrapper = wrapper(solution.getVector());
+ auto solutionWrapper = hierarchicVectorWrapper(solution.vector());
interpolate(id(colBasis), solutionWrapper, values_, dirichletNodes_);
});
diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp
index 8f4fe3750c62984c61803f26431f50ebc0104730..c2536f44d9a2c2b181694a47395f5d1151e6aa88 100644
--- a/src/amdis/DirichletBC.inc.hpp
+++ b/src/amdis/DirichletBC.inc.hpp
@@ -3,9 +3,6 @@
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <amdis/LinearAlgebra.hpp>
-#include <amdis/linear_algebra/HierarchicWrapper.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
namespace AMDiS
{
diff --git a/src/amdis/FileWriter.hpp b/src/amdis/FileWriter.hpp
index 50a40c1d7d71969f84d9fe48f62b1a4b336e2a49..7e0686d3fef0c5a45713fc661cf81416406ab04d 100644
--- a/src/amdis/FileWriter.hpp
+++ b/src/amdis/FileWriter.hpp
@@ -43,14 +43,13 @@ namespace AMDiS
constexpr std::size_t VTKFieldSize = Size<Range>;
- template <class Traits, class TreePath>
+ template <class GlobalBasis, class RangeType, class TreePath>
class FileWriter
: public FileWriterInterface
{
private: // typedefs and static constants
- using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView;
- using Vector = DOFVectorConstView<Traits,TreePath>;
+ using Vector = DOFVectorConstView<GlobalBasis,RangeType,TreePath>;
using Range = typename Vector::Range;
/// Dimension of the mesh
@@ -123,12 +122,12 @@ namespace AMDiS
};
- template <class Traits, class TreePath>
- std::shared_ptr<FileWriter<Traits,TreePath>> makeFileWriterPtr(
- std::string baseName,
- DOFVectorConstView<Traits,TreePath> const& dofvector)
+ template <class GlobalBasis, class Range, class TreePath>
+ std::shared_ptr<FileWriter<GlobalBasis,Range,TreePath>>
+ makeFileWriterPtr(std::string baseName,
+ DOFVectorConstView<GlobalBasis,Range,TreePath> const& dofvector)
{
- return std::make_shared<FileWriter<Traits,TreePath>>(baseName, dofvector);
+ return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, dofvector);
}
} // end namespace AMDiS
diff --git a/src/amdis/GridFunctionOperator.hpp b/src/amdis/GridFunctionOperator.hpp
index 318506fa12473ce9c561258b8a9f781d4cbb2975..167892ffa6c2c05ef08667453d05cb9ed154e319 100644
--- a/src/amdis/GridFunctionOperator.hpp
+++ b/src/amdis/GridFunctionOperator.hpp
@@ -6,6 +6,7 @@
#include <amdis/GridFunctions.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/QuadratureFactory.hpp>
+#include <amdis/common/Transposed.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/utility/FiniteElementType.hpp>
@@ -172,7 +173,7 @@ namespace AMDiS
ColNode const& colNode,
ElementMatrix& elementMatrix)
{
- auto elementMatrixTransposed = trans(elementMatrix);
+ auto elementMatrixTransposed = transposed(elementMatrix);
transposedOp_.getElementMatrix(
context, colNode, rowNode, elementMatrixTransposed);
}
diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp
index 569c7457b973b15f23645ab52a57761317b6c9d5..d0ec769f15b0ce3d8c9b72a40afc026009218e0b 100644
--- a/src/amdis/LinearAlgebra.hpp
+++ b/src/amdis/LinearAlgebra.hpp
@@ -1,31 +1,16 @@
#pragma once
-#include <amdis/linear_algebra/LinearSolverInterface.hpp>
+#include <amdis/linear_algebra/LinearSolver.hpp>
#include <amdis/linear_algebra/SolverInfo.hpp>
-#if defined(AMDIS_BACKEND_ISTL)
-
-#include <amdis/linear_algebra/istl/SystemVector.hpp>
-#include <amdis/linear_algebra/istl/SystemMatrix.hpp>
-#include <amdis/linear_algebra/istl/LinearSolver.hpp>
-
-#elif defined(AMDIS_BACKEND_MTL)
-
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-#include <amdis/linear_algebra/mtl/LinearSolver.hpp>
+#if HAVE_MTL
+#include <amdis/linear_algebra/mtl/DOFVector.hpp>
+#include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
-#include <amdis/linear_algebra/mtl/BITL_Solver.hpp>
-
-#elif defined(AMDIS_BACKEND_PETSC)
-
-#include <amdis/linear_algebra/petsc/SystemVector.hpp>
-#include <amdis/linear_algebra/petsc/SystemMatrix.hpp>
-#include <amdis/linear_algebra/petsc/LinearSolver.hpp>
-
+#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
#else
-
-#error "Unknown linear algebra backend!. Set corresponding variable \
- AMDIS_BACKEND_ISTL, AMDIS_BACKEND_MTL or AMDIS_BACKEND_PETSC."
-
+#include <amdis/linear_algebra/istl/DOFVector.hpp>
+#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
#endif
\ No newline at end of file
diff --git a/src/amdis/LocalAssemblerBase.hpp b/src/amdis/LocalAssemblerBase.hpp
index f6320729f7ddfc14d37557fa8635e03908948379..2504c42c1edb863fa8992202e055829d2b0a6c73 100644
--- a/src/amdis/LocalAssemblerBase.hpp
+++ b/src/amdis/LocalAssemblerBase.hpp
@@ -2,9 +2,10 @@
#include <type_traits>
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
#include <amdis/ContextGeometry.hpp>
-// #include <amdis/common/ConceptsBase.hpp>
-#include <amdis/common/TypeDefs.hpp>
namespace AMDiS
{
@@ -22,10 +23,13 @@ namespace AMDiS
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
+ using ElementMatrix = Dune::DynamicMatrix<double>; // TODO: choose correct value_type
+ using ElementVector = Dune::DynamicVector<double>;
+
/// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
using ElementMatrixVector = std::conditional_t<
- (sizeof...(Nodes)==1), Impl::ElementVector, std::conditional_t<
- (sizeof...(Nodes)==2), Impl::ElementMatrix, void>>;
+ (sizeof...(Nodes)==1), ElementVector, std::conditional_t<
+ (sizeof...(Nodes)==2), ElementMatrix, void>>;
public:
/// Virtual destructor
diff --git a/src/amdis/ProblemInstat.hpp b/src/amdis/ProblemInstat.hpp
index 2a67680284ac7f7943224caaec4a2c9b13a12aae..a61b9b3b55b5e88736fb289a74e5b8b4817c6c83 100644
--- a/src/amdis/ProblemInstat.hpp
+++ b/src/amdis/ProblemInstat.hpp
@@ -60,6 +60,8 @@ namespace AMDiS
/// Returns \ref oldSolution.
std::unique_ptr<SystemVector> getOldSolutionVector() const
{
+ test_exit_dbg(oldSolution,
+ "OldSolution need to be created. Call initialize with INIT_UH_OLD.");
return *oldSolution;
}
diff --git a/src/amdis/ProblemInstat.inc.hpp b/src/amdis/ProblemInstat.inc.hpp
index e6ee9c7fbad04dfd169026268226e5d647afe9fa..5e93f0e889caf94fec46c4bd1c18675e98b71f5b 100644
--- a/src/amdis/ProblemInstat.inc.hpp
+++ b/src/amdis/ProblemInstat.inc.hpp
@@ -45,7 +45,7 @@ void ProblemInstat<Traits>::createUhOld()
if (oldSolution)
warning("oldSolution already created\n");
else // create oldSolution
- oldSolution.reset(new SystemVector(problemStat.globalBasis(), name + "_uOld"));
+ oldSolution.reset(new SystemVector(problemStat.globalBasis()));
}
@@ -53,7 +53,7 @@ template <class Traits>
void ProblemInstat<Traits>::initTimestep(AdaptInfo&)
{
if (oldSolution)
- oldSolution->copy(problemStat.getSolutionVector());
+ *oldSolution = problemStat.getSolutionVector();
}
} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 9e36f2ca8ddaa34689db11e5d47400ac691a8971..0a53c64db5594422c3182544185357dde60ab29a 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -27,7 +27,6 @@
#include <amdis/StandardProblemIteration.hpp>
#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/TypeDefs.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/GridFunctions.hpp>
@@ -68,7 +67,7 @@ namespace AMDiS
static constexpr int dow = Grid::dimensionworld;
using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
- using SystemVector = DOFVector<Traits, double>;
+ using SystemVector = DOFVector<GlobalBasis, double>;
using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
@@ -275,9 +274,9 @@ namespace AMDiS
void createMatricesAndVectors()
{
- systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
- solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
- rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
+ systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
+ solution_ = std::make_shared<SystemVector>(*globalBasis_);
+ rhs_ = std::make_shared<SystemVector>(*globalBasis_);
AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
{
@@ -286,8 +285,6 @@ namespace AMDiS
for (std::size_t j = 0; j < estimates_[i].size(); j++)
estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
});
-
- rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
}
void createSolver()
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 92a7b2155a095f204eea5f93351a73925bec6fc7..833b31c882d8ba0e9d665f100a2d564af54da3a0 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -281,7 +281,7 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solution_->compress();
- linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
+ linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
solverInfo);
if (solverInfo.getInfo() > 0) {
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index 03492346bd8b0e507998ca66b8be429a9b1f9a94..7b4b4ac25b8cc8d50d016821c0f0c2bfb9565e00 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -18,8 +18,8 @@ install(FILES
ScalarTypes.hpp
Size.hpp
Tags.hpp
+ Transposed.hpp
TupleUtility.hpp
- TypeDefs.hpp
Utility.hpp
ValueCategory.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/common)
diff --git a/src/amdis/common/Math.hpp b/src/amdis/common/Math.hpp
index 555d1585c0cbb0cb97bd83e727f7afa681845834..9897ad0c7898ec30295d730648b2ad83a7d420af 100644
--- a/src/amdis/common/Math.hpp
+++ b/src/amdis/common/Math.hpp
@@ -118,7 +118,7 @@ namespace AMDiS
}
template <class T>
- constexpr T threshold = T(1.e-18L); //Math::sqr(std::numeric_limits<T>::epsilon());
+ constexpr T threshold = std::numeric_limits<T>::epsilon();
/// Calculates factorial of i
diff --git a/src/amdis/common/Transposed.hpp b/src/amdis/common/Transposed.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4778e95f956c1ee1635e7d7ab18a5271a47faed5
--- /dev/null
+++ b/src/amdis/common/Transposed.hpp
@@ -0,0 +1,79 @@
+#pragma once
+
+#include <type_traits>
+
+namespace AMDiS
+{
+ /// \brief The transposed view onto a matrix
+ template <class Matrix>
+ class TransposedMatrix
+ {
+ using RawMatrix = std::decay_t<Matrix>;
+
+ public:
+ using size_type = typename RawMatrix::size_type;
+ using value_type = typename RawMatrix::value_type;
+
+ private:
+ struct ConstRowProxy
+ {
+ RawMatrix const* mat;
+ size_type row;
+
+ value_type const& operator[](size_type col) const
+ {
+ return (*mat)[col][row];
+ }
+ };
+
+ struct MutableRowProxy
+ {
+ RawMatrix* mat;
+ size_type row;
+
+ value_type& operator[](size_type col)
+ {
+ return (*mat)[col][row];
+ }
+ };
+
+ public:
+ template <class M>
+ TransposedMatrix(M&& matrix)
+ : matrix_(std::forward<M>(matrix))
+ {}
+
+ ConstRowProxy operator[](size_type row) const
+ {
+ return ConstRowProxy{&matrix_, row};
+ }
+
+ template <class M = Matrix,
+ std::enable_if_t<not std::is_const<M>::value, int> = 0>
+ MutableRowProxy operator[](size_type row)
+ {
+ return MutableRowProxy{&matrix_, row};
+ }
+
+ size_type N() const
+ {
+ return matrix_.M();
+ }
+
+ size_type M() const
+ {
+ return matrix_.N();
+ }
+
+ private:
+ Matrix& matrix_;
+ };
+
+ template <class Matrix>
+ auto transposed(Matrix&& matrix)
+ {
+ using M = std::remove_reference_t<Matrix>;
+ return TransposedMatrix<M>{std::forward<Matrix>(matrix)};
+ }
+
+} // end namespace AMDiS
diff --git a/src/amdis/common/TypeDefs.hpp b/src/amdis/common/TypeDefs.hpp
deleted file mode 100644
index 735d9cb6a50b7022c0ef0b414b8ed850a5f194ba..0000000000000000000000000000000000000000
--- a/src/amdis/common/TypeDefs.hpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#pragma once
-
-#include <type_traits>
-
-#include <amdis/linear_algebra/Mtl.hpp>
-
-namespace AMDiS
-{
- namespace Impl
- {
- using ElementVector = mtl::dense_vector<double>;
- using ElementMatrix = mtl::dense2D<double>;
-
- using SystemVector = mtl::dense_vector<double>;
- using SystemMatrix = mtl::compressed2D<double>;
-
- } // end namespace Impl
-
- struct BoundaryType { int b; };
-
-} // end namespace AMDiS
diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp
index 0cfccfc24e91911cdb18f4ec96b69cda4d50ca34..86b52e7e7771030973eeae3b3a83932e9672cd47 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -20,13 +20,12 @@ namespace AMDiS
* @{
**/
- template <class Traits, class TreePathType>
+ template <class GlobalBasisType, class RangeType, class TreePathType>
class DOFVectorConstView
{
public:
- using GlobalBasis = typename Traits::GlobalBasis;
+ using GlobalBasis = GlobalBasisType;
using TreePath = TreePathType;
- using Vector = DOFVector<Traits>;
using Tree = typename GlobalBasis::LocalView::Tree;
using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
@@ -37,6 +36,8 @@ namespace AMDiS
using Domain = typename EntitySet::GlobalCoordinate;
using Range = RangeType_t<SubTree>;
+ static_assert(std::is_arithmetic<RangeType>::value, "");
+ // Don't know how to determine Range with non-trivial RangeType
using RawSignature = typename Dune::Functions::SignatureTraits<Range(Domain)>::RawSignature;
using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
@@ -201,11 +202,11 @@ namespace AMDiS
public:
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
- DOFVectorConstView(DOFVector<Traits> const& dofVector, TreePath const& treePath)
+ DOFVectorConstView(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath const& treePath)
: dofVector_(&dofVector)
, treePath_(treePath)
- , entitySet_(dofVector.getFeSpace().gridView())
- , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.getFeSpace(), treePath))
+ , entitySet_(dofVector.basis().gridView())
+ , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath))
{}
/// Evaluate the view on this DOFVector in global coordinates
@@ -230,7 +231,7 @@ namespace AMDiS
/// Return global basis
GlobalBasis const& basis() const
{
- return dofVector_->getFeSpace();
+ return dofVector_->basis();
}
/// Return treePath associated with this view
@@ -240,13 +241,13 @@ namespace AMDiS
}
/// Return const coefficient vector
- DOFVector<Traits> const& coefficients() const
+ DOFVector<GlobalBasis,RangeType> const& coefficients() const
{
return *dofVector_;
}
protected:
- DOFVector<Traits> const* dofVector_;
+ DOFVector<GlobalBasis,RangeType> const* dofVector_;
TreePath const treePath_;
EntitySet entitySet_;
@@ -255,18 +256,18 @@ namespace AMDiS
// A mutable version of DOFVectorView
- template <class Traits, class TreePathType>
+ template <class GlobalBasisType, class RangeType, class TreePathType>
class DOFVectorMutableView
- : public DOFVectorConstView<Traits, TreePathType>
+ : public DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>
{
- using Super = DOFVectorConstView<Traits, TreePathType>;
+ using Super = DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>;
- using GlobalBasis = typename Traits::GlobalBasis;
+ using GlobalBasis = GlobalBasisType;
using TreePath = TreePathType;
public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
- DOFVectorMutableView(DOFVector<Traits>& dofVector, TreePath const& treePath)
+ DOFVectorMutableView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath)
: Super(dofVector, treePath)
, mutableDofVector_(&dofVector)
{}
@@ -281,11 +282,11 @@ namespace AMDiS
auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
- DOFVector<Traits> tmp(basis, "tmp");
+ DOFVector<GlobalBasis,RangeType> tmp(*mutableDofVector_);
Dune::Functions::interpolate(basis, treePath, tmp, std::forward<decltype(gridFct)>(gridFct));
// move data from temporary vector into stored DOFVector
- mutableDofVector_->getVector() = std::move(tmp.getVector());
+ mutableDofVector_->vector() = std::move(tmp.vector());
return *this;
}
@@ -297,13 +298,13 @@ namespace AMDiS
/// Return the mutable DOFVector
- DOFVector<Traits>& coefficients() { return *mutableDofVector_; }
+ DOFVector<GlobalBasis,RangeType>& coefficients() { return *mutableDofVector_; }
/// Return the const DOFVector
using Super::coefficients;
protected:
- DOFVector<Traits>* mutableDofVector_;
+ DOFVector<GlobalBasis,RangeType>* mutableDofVector_;
};
/** @} **/
@@ -311,34 +312,34 @@ namespace AMDiS
#ifndef DOXYGEN
// A Generator for a const \ref DOFVectorView.
- template <class Traits, class TreePath>
- auto makeDOFVectorView(DOFVector<Traits> const& dofVector, TreePath const& treePath)
+ template <class GlobalBasis, class RangeType, class TreePath>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
{
- return DOFVectorConstView<Traits, TreePath>{dofVector, treePath};
+ return DOFVectorConstView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
}
// A Generator for a mutable \ref DOFVectorView.
- template <class Traits, class TreePath>
- auto makeDOFVectorView(DOFVector<Traits>& dofVector, TreePath const& treePath)
+ template <class GlobalBasis, class RangeType, class TreePath>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector, TreePath const& treePath)
{
- return DOFVectorMutableView<Traits, TreePath>{dofVector, treePath};
+ return DOFVectorMutableView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
}
// A Generator for a const \ref DOFVectorView.
- template <class Traits>
- auto makeDOFVectorView(DOFVector<Traits> const& dofVector)
+ template <class GlobalBasis, class RangeType>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorConstView<Traits, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorConstView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
}
// A Generator for a mutable \ref DOFVectorView.
- template <class Traits>
- auto makeDOFVectorView(DOFVector<Traits>& dofVector)
+ template <class GlobalBasis, class RangeType>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorMutableView<Traits, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorMutableView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
}
#endif
diff --git a/src/amdis/gridfunctions/DOFVectorView.inc.hpp b/src/amdis/gridfunctions/DOFVectorView.inc.hpp
index a97dc1a47b1922b0db421f80b19468858a4f2076..139a1ab4c98883747d4882eee5bac08975a663ea 100644
--- a/src/amdis/gridfunctions/DOFVectorView.inc.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.inc.hpp
@@ -4,8 +4,8 @@
namespace AMDiS {
-template <class Traits, class TreePath>
-typename DOFVectorConstView<Traits, TreePath>::Range DOFVectorConstView<Traits, TreePath>::
+template <class GB, class RT, class TP>
+typename DOFVectorConstView<GB,RT,TP>::Range DOFVectorConstView<GB,RT,TP>::
LocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
@@ -48,8 +48,8 @@ LocalFunction::operator()(LocalDomain const& x) const
}
-template <class Traits, class TreePath>
-typename DOFVectorConstView<Traits, TreePath>::DerivativeRange DOFVectorConstView<Traits, TreePath>::
+template <class GB, class RT, class TP>
+typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
GradientLocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
diff --git a/src/amdis/linear_algebra/CMakeLists.txt b/src/amdis/linear_algebra/CMakeLists.txt
index 200338d51155f07071e20d69efbbafe1ab409853..ff79554f15af37891a18244704545418050e4a33 100644
--- a/src/amdis/linear_algebra/CMakeLists.txt
+++ b/src/amdis/linear_algebra/CMakeLists.txt
@@ -1,10 +1,12 @@
#install headers
install(FILES
+ Common.hpp
+ DOFMatrixBase.hpp
+ DOFVectorBase.hpp
+ DOFVectorInterface.hpp
HierarchicWrapper.hpp
- LinearAlgebraBase.hpp
LinearSolverInterface.hpp
- Mtl.hpp
PreconditionerInterface.hpp
RunnerInterface.hpp
SolverInfo.hpp
diff --git a/src/amdis/linear_algebra/Common.hpp b/src/amdis/linear_algebra/Common.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7bf16d923776bf233d6114d4706986d9929d9b3
--- /dev/null
+++ b/src/amdis/linear_algebra/Common.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <cstddef>
+
+namespace AMDiS
+{
+ template <class T>
+ struct Triplet
+ {
+ std::size_t row, cols;
+ T value;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/DOFMatrixBase.hpp b/src/amdis/linear_algebra/DOFMatrixBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9080a167cf70a842aaef3dbad509b17db985e851
--- /dev/null
+++ b/src/amdis/linear_algebra/DOFMatrixBase.hpp
@@ -0,0 +1,122 @@
+#pragma once
+
+#include <cmath>
+#include <amdis/common/Math.hpp>
+
+namespace AMDiS
+{
+ template <class RowBasisType, class ColBasisType, class Backend>
+ class DOFMatrixBase
+ {
+ public:
+ /// The type of the finite element space / basis of the row
+ using RowBasis = RowBasisType;
+
+ /// The type of the finite element space / basis of the column
+ using ColBasis = ColBasisType;
+
+ /// The index/size - type
+ using size_type = typename RowBasis::size_type;
+
+ using BaseMatrix = typename Backend::BaseMatrix;
+
+ public:
+ /// Constructor. Constructs new BaseVector.
+ DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : rowBasis_(&rowBasis)
+ , colBasis_(&colBasis)
+ {}
+
+ /// Return the row-basis \ref rowBasis of the matrix
+ RowBasis const& rowBasis() const
+ {
+ return *rowBasis_;
+ }
+
+ /// Return the col-basis \ref colBasis of the matrix
+ ColBasis const& colBasis() const
+ {
+ return *colBasis_;
+ }
+
+ /// Return the data-matrix
+ BaseMatrix const& matrix() const
+ {
+ return backend_.matrix();
+ }
+
+ /// Return the data-matrix
+ BaseMatrix& matrix()
+ {
+ return backend_.matrix();
+ }
+
+ /// Return the size of the \ref rowBasis_
+ size_type rows() const
+ {
+ return rowBasis_->dimension();
+ }
+
+ /// Return the size of the \ref colBasis_
+ size_type cols() const
+ {
+ return colBasis_->dimension();
+ }
+
+ /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern
+ /// If \p setToZero is true, the matrix is set to 0
+ void init(bool setToZero)
+ {
+ backend_.init(*rowBasis_, *colBasis_, setToZero);
+ }
+
+ /// Finish the matrix insertion, e.g. cleanup or final insertion
+ void finish()
+ {
+ backend_.finish();
+ }
+
+ /// Insert a block of values into the matrix (add to existing values)
+ template <class ElementMatrix>
+ void insert(typename RowBasis::LocalView const& rowLocalView,
+ typename ColBasis::LocalView const& colLocalView,
+ ElementMatrix const& elementMatrix)
+ {
+ for (size_type i = 0; i < rowLocalView.size(); ++i) {
+ size_type const row = flatMultiIndex(rowLocalView.index(i));
+ for (size_type j = 0; j < colLocalView.size(); ++j) {
+ if (std::abs(elementMatrix[i][j]) > threshold<double>) {
+ size_type const col = flatMultiIndex(colLocalView.index(j));
+ backend_.insert(row, col, elementMatrix[i][j]);
+ }
+ }
+ }
+ }
+
+ /// Insert a single value into the matrix (add to existing value)
+ template <class RowIndex, class ColIndex>
+ void insert(RowIndex row, ColIndex col, typename Backend::value_type const& value)
+ {
+ backend_.insert(flatMultiIndex(row), flatMultiIndex(col), value);
+ }
+
+ /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
+ /// a one on the diagonal of the matrix.
+ auto applyDirichletBC(std::vector<char> const& dirichletNodes)
+ {
+ return backend_.applyDirichletBC(dirichletNodes);
+ }
+
+ size_type nnz() const
+ {
+ return backend_.nnz();
+ }
+
+ protected:
+ RowBasis const* rowBasis_;
+ ColBasis const* colBasis_;
+
+ Backend backend_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/DOFVectorBase.hpp b/src/amdis/linear_algebra/DOFVectorBase.hpp
index 2b3902eed018517804ee296c1ba3f36288269f73..33fa8ea8277796806c2c4dbebd00de08fff58c0c 100644
--- a/src/amdis/linear_algebra/DOFVectorBase.hpp
+++ b/src/amdis/linear_algebra/DOFVectorBase.hpp
@@ -1,30 +1,124 @@
#pragma once
+#include <cmath>
+
+#include <dune/functions/functionspacebases/sizeinfo.hh>
+
+#include <amdis/common/Math.hpp>
+#include <amdis/common/ScalarTypes.hpp>
+#include <amdis/linear_algebra/DOFVectorInterface.hpp>
+#include <amdis/utility/MultiIndex.hpp>
+
namespace AMDiS
{
- template <class GlobalBasisType>
+ /// The basic container that stores a base vector and a corresponding basis
+ template <class BasisType, class Backend>
class DOFVectorBase
+ : public DOFVectorInterface
{
+ using Self = DOFVectorBase;
+
+ public:
+ /// The type of the \ref basis
+ using Basis = BasisType;
+
+ /// The index/size - type
+ using size_type = typename BasisType::size_type;
+
+ /// The type of the elements of the DOFVector
+ using value_type = typename Backend::value_type;
+
+ using BaseVector = typename Backend::BaseVector;
+
public:
- using GlobalBasis = GlobalBasisType;
- using SizeType = std::size_t;
+ /// Constructor. Constructs new BaseVector.
+ DOFVectorBase(BasisType const& basis)
+ : basis_(&basis)
+ {
+ compress();
+ }
- DOFVectorBase(GlobalBasisType const& basis)
- : basis_(basis)
- {}
+ /// Return the basis \ref basis_ associated to the vector
+ Basis const& basis() const
+ {
+ return *basis_;
+ }
- SizeType size() const
+ /// Return the data-vector
+ BaseVector const& vector() const
{
- return basis_.dimension();
+ return backend_.vector();
}
- GlobalBasis const& basis() const
+ /// Return the data-vector
+ BaseVector& vector()
{
- return GlobalBasis;
+ return backend_.vector();
+ }
+
+ /// Return the size of the \ref basis
+ size_type size() const
+ {
+ return basis_->dimension();
+ }
+
+ void resize(Dune::Functions::SizeInfo<Basis> const& s)
+ {
+ backend_.resize(size_type(s));
+ }
+
+ /// Resize the \ref vector to the size of the \ref basis.
+ virtual void compress() override
+ {
+ if (backend_.size() != size()) {
+ backend_.resize(size());
+ vector() = 0;
+ }
+ }
+
+ /// Access the entry \p idx of the \ref vector with read-access.
+ template <class Index>
+ auto const& operator[](Index idx) const
+ {
+ size_type i = flatMultiIndex(idx);
+ return backend_[i];
+ }
+
+ /// Access the entry \p idx of the \ref vector with write-access.
+ template <class Index>
+ auto& operator[](Index idx)
+ {
+ size_type i = flatMultiIndex(idx);
+ return backend_[i];
+ }
+
+ /// Insert a block of values into the matrix (add to existing values)
+ template <class ElementVector>
+ void insert(typename Basis::LocalView const& localView,
+ ElementVector const& elementVector)
+ {
+ for (size_type i = 0; i < localView.size(); ++i) {
+ if (std::abs(elementVector[i]) > threshold<double>) {
+ size_type const idx = flatMultiIndex(localView.index(i));
+ backend_[idx] += elementVector[i];
+ }
+ }
+ }
+
+ /// Sets each DOFVector to the scalar \p value.
+ template <class Scalar,
+ std::enable_if_t<Concepts::Arithmetic<Scalar>, int> = 0>
+ Self& operator=(Scalar value)
+ {
+ vector() = value;
+ return *this;
}
private:
- GlobalBasis const& basis_;
+ /// The finite element space / basis associated with the data vector
+ Basis const* basis_;
+
+ Backend backend_;
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/DOFVectorInterface.hpp b/src/amdis/linear_algebra/DOFVectorInterface.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f2b5aa84301a305ff5682f01582793410d633132
--- /dev/null
+++ b/src/amdis/linear_algebra/DOFVectorInterface.hpp
@@ -0,0 +1,15 @@
+#pragma once
+
+namespace AMDiS
+{
+ class DOFVectorInterface
+ {
+ public:
+ /// Virtual destructor
+ ~DOFVectorInterface() = default;
+
+ /// Change dimension of DOFVector to dimension of basis
+ virtual void compress() = 0;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/HierarchicWrapper.hpp b/src/amdis/linear_algebra/HierarchicWrapper.hpp
index 94617e6de355bc2fbd9b7bb8ea3953fad3abb998..79f1f4c1b576b65bcffe2787fd35cf0deb536a6e 100644
--- a/src/amdis/linear_algebra/HierarchicWrapper.hpp
+++ b/src/amdis/linear_algebra/HierarchicWrapper.hpp
@@ -1,5 +1,9 @@
#pragma once
+#if HAVE_MTL
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#endif
+
#include <dune/functions/functionspacebases/sizeinfo.hh>
#include <amdis/utility/MultiIndex.hpp>
@@ -85,11 +89,13 @@ namespace AMDiS
vec.change_dim(std::size_t(s));
}
+#if HAVE_MTL
template <class... Params>
- std::size_t sizeImpl(mtl::dense_vector<Params...> const& v) const
+ std::size_t sizeImpl(mtl::vec::dense_vector<Params...> const& v) const
{
return mtl::size(v);
}
+#endif
template <class V>
using HasSizeMember = decltype(std::declval<V>().size());
diff --git a/src/amdis/linear_algebra/LinearAlgebraBase.hpp b/src/amdis/linear_algebra/LinearAlgebraBase.hpp
deleted file mode 100644
index 41c3d1601381e12200481f63d4ea7c23ac205a3c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/LinearAlgebraBase.hpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#pragma once
-
-namespace AMDiS
-{
- namespace Impl
- {
- // Type-Trait needs to be specialized by linear algebra libraries
- template <class Vector>
- struct BaseVector { using type = Vector; };
-
- template <class Matrix>
- struct BaseMatrix { using type = Matrix; };
- }
-
- template <class Vector>
- using BaseVector_t = typename Impl::BaseVector<Vector>::type;
-
- template <class Matrix>
- using BaseMatrix_t = typename Impl::BaseMatrix<Matrix>::type;
-
- template <class T>
- struct Triplet
- {
- std::size_t row;
- std::size_t col;
- T value;
- };
-
- namespace Impl
- {
- template <class Vector>
- class BaseWrapper
- {
- public:
- using value_type = typename Vector::value_type;
- using size_type = typename Vector::size_type;
-
- explicit BaseWrapper(Vector& vec) : vec(vec) {}
-
- void resize(size_type s)
- {
- vec.resize(s);
- }
-
- size_type size() const
- {
- return vec.size();
- }
-
- value_type& operator[](size_type i) { return vec[i]; }
- value_type const& operator[](size_type i) const { return vec[i]; }
-
- protected:
- Vector& vec;
- };
-
- } // end namespace Impl
-
- template <class Vector>
- Impl::BaseWrapper<std::remove_reference_t<Vector>> wrapper(Vector&& vec)
- {
- return {std::forward<Vector>(vec)};
- }
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/LinearSolver.hpp b/src/amdis/linear_algebra/LinearSolver.hpp
similarity index 57%
rename from src/amdis/linear_algebra/mtl/LinearSolver.hpp
rename to src/amdis/linear_algebra/LinearSolver.hpp
index b8a1f70e1a1929b999d6f9abd35c5d603ec14751..b132ef31ba780e5591000c50bf304a5848b29837 100644
--- a/src/amdis/linear_algebra/mtl/LinearSolver.hpp
+++ b/src/amdis/linear_algebra/LinearSolver.hpp
@@ -8,25 +8,24 @@
#include <amdis/CreatorInterface.hpp>
#include <amdis/Output.hpp>
#include <amdis/linear_algebra/LinearSolverInterface.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
namespace AMDiS
{
/** \ingroup Solver
*
* \brief
- * Wrapper class for various MTL4 solvers. These solvers
+ * Wrapper class for various solvers. These solvers
* are parametrized by MatrixType and VectorType. The different
* algorithms, like krylov subspace methods or other external
- * solvers where MTL4 provides an interface, can be assigned
+ * solvers where the backend provides an interface, can be assigned
* by different Runner objects.
**/
- template <class Matrix, class Vector, class Runner>
+ template <class Matrix, class VectorX, class VectorB, class Runner>
class LinearSolver
- : public LinearSolverInterface<Matrix, Vector, Vector>
+ : public LinearSolverInterface<Matrix, VectorX, VectorB>
{
using Self = LinearSolver;
- using Super = LinearSolverInterface<Matrix, Vector, Vector>;
+ using Super = LinearSolverInterface<Matrix, VectorX, VectorB>;
using RunnerBase = typename Super::RunnerBase;
public:
@@ -42,44 +41,39 @@ namespace AMDiS
public:
/// Constructor
explicit LinearSolver(std::string prefix)
- : runner(std::make_shared<Runner>(prefix))
+ : runner_(std::make_shared<Runner>(prefix))
{}
- /// Implements \ref LinearSolverInterface::getRunner()
- virtual std::shared_ptr<RunnerBase> getRunner() override
+ /// Implements \ref LinearSolverInterface::runner()
+ virtual std::shared_ptr<RunnerBase> runner() override
{
- return runner;
+ return runner_;
}
private:
/// Implements \ref LinearSolverInterface::solveSystemImpl()
- virtual void solveImpl(Matrix const& A, Vector& x, Vector const& b,
+ virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo) override
{
Dune::Timer t;
if (solverInfo.doCreateMatrixData()) {
// init matrix or wrap block-matrix or ...
- runner->init(A);
+ runner_->init(A);
}
- // wrappers copy the data back on destruction
- auto X = blockWrapper(x);
- auto B = blockWrapper(b);
+ if (solverInfo.getInfo() > 0)
+ msg("fill matrix needed ", t.elapsed(), " seconds");
- if (solverInfo.getInfo() > 0) {
- msg("fill MTL4 matrix needed ", t.elapsed(), " seconds");
- }
-
- int error = runner->solve(A, X.getVector(), B.getVector(), solverInfo);
+ int error = runner_->solve(A, x, b, solverInfo);
solverInfo.setError(error);
if (!solverInfo.doStoreMatrixData())
- runner->exit();
+ runner_->exit();
}
private:
/// redirect the implementation to a runner. Implements a \ref RunnerInterface
- std::shared_ptr<Runner> runner;
+ std::shared_ptr<Runner> runner_;
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/LinearSolverInterface.hpp b/src/amdis/linear_algebra/LinearSolverInterface.hpp
index 5d1534f9f522db80170bdd3ee36bea85d25d56cd..54228ef66b9004fd42e0c2c9d1b8907f24cf0d4f 100644
--- a/src/amdis/linear_algebra/LinearSolverInterface.hpp
+++ b/src/amdis/linear_algebra/LinearSolverInterface.hpp
@@ -12,7 +12,6 @@
*/
#include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
#include <amdis/linear_algebra/PreconditionerInterface.hpp>
#include <amdis/linear_algebra/RunnerInterface.hpp>
@@ -25,12 +24,11 @@ namespace AMDiS
class LinearSolverInterface
{
protected:
- using RunnerBase = RunnerInterface<Matrix, BaseVector_t<VectorX>,
- BaseVector_t<VectorB>>;
+ using RunnerBase = RunnerInterface<Matrix, VectorX, VectorB>;
public:
/// Destructor.
- virtual ~LinearSolverInterface() {}
+ virtual ~LinearSolverInterface() = default;
/// Public method to call in order to solve a linear system Ax = b.
/**
@@ -45,12 +43,11 @@ namespace AMDiS
void solve(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo)
{
- assert( num_cols(A) == size(x) && num_rows(A) == size(b) );
solveImpl(A, x, b, solverInfo);
}
// return the runner/worker corresponding to this solver (optional)
- virtual std::shared_ptr<RunnerBase> getRunner() { return {}; };
+ virtual std::shared_ptr<RunnerBase> runner() { return {}; };
private:
/// main methods that all solvers must implement
diff --git a/src/amdis/linear_algebra/Mtl.hpp b/src/amdis/linear_algebra/Mtl.hpp
deleted file mode 100644
index 5d89ea427fd21f3cd57d2926100c19586478b532..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/Mtl.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#pragma once
-
-#include <dune/common/ftraits.hh>
-
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
-namespace Dune
-{
- template <class Elt, class Parameters>
- struct FieldTraits<mtl::compressed2D<Elt,Parameters>>
- {
- using field_type = typename FieldTraits<Elt>::field_type;
- using real_type = typename FieldTraits<Elt>::real_type;
- };
-
- template <class Value, class Parameters>
- struct FieldTraits<mtl::dense2D<Value,Parameters>>
- {
- using field_type = typename FieldTraits<Value>::field_type;
- using real_type = typename FieldTraits<Value>::real_type;
- };
-
- template <class Value, class Parameters>
- struct FieldTraits<mtl::dense_vector<Value,Parameters>>
- {
- using field_type = typename FieldTraits<Value>::field_type;
- using real_type = typename FieldTraits<Value>::real_type;
- };
-
-} // end namespace Dune
diff --git a/src/amdis/linear_algebra/PreconditionerInterface.hpp b/src/amdis/linear_algebra/PreconditionerInterface.hpp
index 692a35250d904b2b5e63bf71f261ffdceaee705f..055ee01b886b6ee3e600583b811fc00c869c6265 100644
--- a/src/amdis/linear_algebra/PreconditionerInterface.hpp
+++ b/src/amdis/linear_algebra/PreconditionerInterface.hpp
@@ -7,7 +7,7 @@ namespace AMDiS
struct PreconditionerInterface
{
/// Virtual destructor.
- virtual ~PreconditionerInterface() {}
+ virtual ~PreconditionerInterface() = default;
/// Is called a the beginning of a solution procedure
virtual void init(Matrix const& fullMatrix) = 0;
diff --git a/src/amdis/linear_algebra/RunnerInterface.hpp b/src/amdis/linear_algebra/RunnerInterface.hpp
index 0c81227ce650ab7b3a18929dc7c0e9480bc5aa89..ae7fc03da6e7f36ff8de4cd2a474f6ad93e1b889 100644
--- a/src/amdis/linear_algebra/RunnerInterface.hpp
+++ b/src/amdis/linear_algebra/RunnerInterface.hpp
@@ -14,7 +14,7 @@ namespace AMDiS
public:
/// virtual destructor
- virtual ~RunnerInterface() {}
+ virtual ~RunnerInterface() = default;
/// Is called at the beginning of a solution procedure
virtual void init(Matrix const& A) = 0;
@@ -34,8 +34,8 @@ namespace AMDiS
return 0;
}
- virtual std::shared_ptr<PreconBase> getLeftPrecon() { return {}; }
- virtual std::shared_ptr<PreconBase> getRightPrecon() { return {}; }
+ virtual std::shared_ptr<PreconBase> leftPrecon() { return {}; }
+ virtual std::shared_ptr<PreconBase> rightPrecon() { return {}; }
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/CMakeLists.txt b/src/amdis/linear_algebra/istl/CMakeLists.txt
index cd5ee97f5f7fbaf3fa22c39f74a54b6eeb33f9ce..f1e73c05ad38ceff6076ef8155c52301e49c5f6b 100644
--- a/src/amdis/linear_algebra/istl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/istl/CMakeLists.txt
@@ -4,14 +4,11 @@
#)
install(FILES
+ DirectRunner.hpp
DOFMatrix.hpp
DOFVector.hpp
+ Fwd.hpp
ISTL_Preconditioner.hpp
ISTL_Solver.hpp
ISTLRunner.hpp
- LinearSolver.hpp
- Preconditioner.hpp
- SystemMatrix.hpp
- SystemVector.hpp
- UmfpackRunner.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/istl)
diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
index e256a05eecca3bf3460ae2b4a39a102311be6147..61837bf6083cdb3830c306245e1edcff9880edfd 100644
--- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <list>
#include <string>
#include <memory>
@@ -7,112 +8,73 @@
#include <dune/istl/matrixindexset.hh>
#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
namespace AMDiS
{
- template <class RowFeSpaceType, class ColFeSpaceType,
- class ValueType = Dune::FieldMatrix<double,1,1>>
- class DOFMatrix
+ template <class T, class = void>
+ struct BlockMatrixType
+ {
+ using type = Dune::FieldMatrix<T,1,1>;
+ };
+
+ template <class T>
+ struct BlockMatrixType<T, typename T::field_type>
+ {
+ using type = T;
+ };
+
+ template <class ValueType>
+ class IstlMatrix
{
public:
- using RowFeSpace = RowFeSpaceType;
- using ColFeSpace = ColFeSpaceType;
- using BaseMatrix = Dune::BCRSMatrix<ValueType>;
+ /// The type of the elements of the DOFMatrix
+ using value_type = typename BlockMatrixType<ValueType>::type;
- using size_type = typename RowFeSpaceType::size_type;
- using field_type = typename ValueType::field_type;
+ /// The matrix type of the underlying base matrix
+ using BaseMatrix = Dune::BCRSMatrix<value_type>;
- using value_type = ValueType;
+ /// The index/size - type
+ using size_type = typename BaseMatrix::size_type;
+ public:
/// Constructor. Constructs new BaseVector.
- DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
- : rowFeSpace(rowFeSpace)
- , colFeSpace(colFeSpace)
- , name(name)
- , matrix(ClonablePtr<BaseMatrix>::make())
- {}
-
- /// Constructor. Takes pointer of data-vector.
- DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name,
- BaseMatrix& matrix_ref)
- : rowFeSpace(rowFeSpace)
- , colFeSpace(colFeSpace)
- , name(name)
- , matrix(matrix_ref)
- {}
-
- /// Return the row-basis \ref rowFeSpace of the matrix
- RowFeSpace const& getRowFeSpace() const
- {
- return rowFeSpace;
- }
-
- /// Return the col-basis \ref colFeSpace of the matrix
- ColFeSpace const& getColFeSpace() const
- {
- return colFeSpace;
- }
+ IstlMatrix() = default;
/// Return the data-vector \ref vector
- BaseMatrix const& getMatrix() const
+ BaseMatrix const& matrix() const
{
- return *matrix;
+ return matrix_;
}
/// Return the data-vector \ref vector
- BaseMatrix& getMatrix()
- {
- return *matrix;
- }
-
- /// Return the size of the \ref feSpace
- size_type N() const
- {
- return rowFeSpace.size();
- }
-
- size_type M() const
- {
- return colFeSpace.size();
- }
-
- /// Return the \ref name of this vector
- std::string getName() const
+ BaseMatrix& matrix()
{
- return name;
+ return matrix_;
}
- /// Access the entry \p i of the \ref vector with read-access.
- value_type const& operator()(size_type r, size_type c) const
+ /// Insert a single value into the matrix (add to existing value)
+ void insert(size_type r, size_type c, value_type const& value)
{
- test_exit_dbg( initialized, "Occupation pattern not initialized!");
- test_exit_dbg( r < N() && c < M() ,
- "Indices out of range [0,", N(), ")x[0,", M(), ")" );
- return (*matrix)[r][c];
- }
-
- /// Access the entry \p i of the \ref vector with write-access.
- value_type& operator()(size_type r, size_type c)
- {
- test_exit_dbg( initialized, "Occupation pattern not initialized!");
- test_exit_dbg( r < N() && c < M() ,
- "Indices out of range [0,", N(), ")x[0,", M(), ")" );
- return (*matrix)[r][c];
+ test_exit_dbg( initialized_, "Occupation pattern not initialized!");
+ test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
+ "Indices out of range [0,", matrix_.N(), ")x[0,", matrix_.M(), ")" );
+ matrix_[r][c] += value;
}
/// create occupation pattern and apply it to the matrix
- void init()
+ template <class RowBasis, class ColBasis>
+ void init(RowBasis const& rowBasis, ColBasis const& colBasis,
+ bool prepareForInsertion)
{
- occupationPattern.resize(rowFeSpace.size(), colFeSpace.size());
- auto meshView = rowFeSpace.gridView();
-
- // A loop over all elements of the grid
- auto rowLocalView = rowFeSpace.localView();
- auto colLocalView = colFeSpace.localView();
+ auto occupationPattern = Dune::MatrixIndexSet{rowBasis.dimension(), colBasis.dimension()};
- for (const auto& element : elements(meshView)) {
+ auto rowLocalView = rowBasis.localView();
+ auto colLocalView = colBasis.localView();
+ for (const auto& element : elements(rowBasis.gridView())) {
rowLocalView.bind(element);
colLocalView.bind(element);
@@ -126,50 +88,44 @@ namespace AMDiS
}
}
}
- occupationPattern.exportIdx(*matrix);
+ occupationPattern.exportIdx(matrix_);
- initialized = true;
+ initialized_ = true;
}
void finish()
{
- initialized = false;
+ initialized_ = false;
}
- std::size_t getNnz()
+ std::size_t nnz() const
{
- return matrix->nonzeroes();
+ return matrix_.nonzeroes();
}
- auto clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
+ auto applyDirichletBC(std::vector<char> const& dirichletNodes)
{
// loop over the matrix rows
- for (std::size_t i = 0; i < matrix->N(); ++i) {
+ for (std::size_t i = 0; i < matrix_.N(); ++i) {
if (dirichletNodes[i]) {
- auto cIt = (*matrix)[i].begin();
- auto cEndIt = (*matrix)[i].end();
+ auto cIt = matrix_[i].begin();
+ auto cEndIt = matrix_[i].end();
// loop over nonzero matrix entries in current row
for (; cIt != cEndIt; ++cIt)
- *cIt = (apply && i == cIt.index() ? 1 : 0);
+ *cIt = (i == cIt.index() ? 1 : 0);
}
}
- std::vector<std::map<size_t,value_type>> columns; // symmetric dbc not yet implemented
- return columns;
+ return std::list<Triplet<value_type>>{};
}
private:
- RowFeSpace const& rowFeSpace;
- ColFeSpace const& colFeSpace;
- std::string name;
+ BaseMatrix matrix_;
- ClonablePtr<BaseMatrix> matrix;
- Dune::MatrixIndexSet occupationPattern;
-
- bool initialized = false;
-
- // friend class declarations
- template <class, class, class> friend class SystemMatrix;
+ bool initialized_ = false;
};
+ template <class RowBasisType, class ColBasisType, class ValueType = double>
+ using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, IstlMatrix<ValueType>>;
+
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DOFVector.hpp b/src/amdis/linear_algebra/istl/DOFVector.hpp
index 77fa9d66c4cae779646ccf1a69a5638b19abbff4..a1234090740cb542c806b6e13b54bf5817ff1c2a 100644
--- a/src/amdis/linear_algebra/istl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/istl/DOFVector.hpp
@@ -1,124 +1,97 @@
#pragma once
-#include <string>
-#include <memory>
-
#include <dune/istl/bvector.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
-
#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
namespace AMDiS
{
- template <class FeSpaceType, class ValueType = Dune::FieldVector<double,1>>
- class DOFVector
- {
- using Self = DOFVector;
-
- public:
- using FeSpace = FeSpaceType;
- using BaseVector = Dune::BlockVector<ValueType>;
-
- using size_type = typename FeSpace::size_type;
- using field_type = typename ValueType::field_type;
-
- using value_type = ValueType;
-
- /// Constructor. Constructs new BaseVector.
- DOFVector(FeSpace const& feSpace, std::string name)
- : feSpace(feSpace)
- , name(name)
- , vector(ClonablePtr<BaseVector>::make())
- {
- compress();
- }
-
- /// Constructor. Takes pointer of data-vector.
- DOFVector(FeSpace const& feSpace, std::string name,
- BaseVector& vector_ref)
- : feSpace(feSpace)
- , name(name)
- , vector(vector_ref)
- {}
-
- /// Return the basis \ref feSpace of the vector
- FeSpace const& getFeSpace() const
- {
- return feSpace;
- }
-
- /// Return the data-vector \ref vector
- BaseVector const& getVector() const
+ template <class T, class = void>
+ struct BlockVectorType
+ {
+ using type = Dune::FieldVector<T,1>;
+ };
+
+ template <class T>
+ struct BlockVectorType<T, typename T::field_type>
+ {
+ using type = T;
+ };
+
+ template <class ValueType>
+ class IstlVector
{
- return *vector;
- }
+ public:
+ /// The type of the elements of the DOFVector
+ using value_type = typename BlockVectorType<ValueType>::type;
- /// Return the data-vector \ref vector
- BaseVector& getVector()
- {
- return *vector;
- }
+ /// The vector type of the underlying base vector
+ using BaseVector = Dune::BlockVector<value_type>;
- /// Return the size of the \ref feSpace
- size_type getSize() const
- {
- return feSpace.size();
- }
+ /// The index/size - type
+ using size_type = typename BaseVector::size_type;
- /// Return the \ref name of this vector
- std::string getName() const
- {
- return name;
- }
+ public:
+ /// Constructor. Constructs new BaseVector.
+ IstlVector() = default;
- /// Resize the \ref vector to the size of the \ref feSpace.
- void compress()
- {
- vector->resize(getSize() / value_type::dimension);
- }
+ /// Return the data-vector \ref vector
+ BaseVector const& vector() const
+ {
+ return vector_;
+ }
+ /// Return the data-vector \ref vector
+ BaseVector& vector()
+ {
+ return vector_;
+ }
- /// Access the entry \p i of the \ref vector with read-access.
- value_type const& operator[](size_type i) const
- {
- test_exit_dbg(i < vector->size() ,
- "Index ", i, " out of range [0,", vector->size(), ")" );
- return (*vector)[i];
- }
- /// Access the entry \p i of the \ref vector with write-access.
- value_type& operator[](size_type i)
- {
- test_exit_dbg(i < vector->size() ,
- "Index ", i, " out of range [0,", vector->size(), ")" );
- return (*vector)[i];
- }
-
- /// interpolate a function \p f to the basis \ref feSpace and store the
- /// coefficients in \ref vector.
- template <class F>
- void interpol(F const& f)
- {
- Dune::Functions::interpolate(feSpace, *vector, f);
- }
+ /// Return the current size of the \ref vector_
+ size_type size() const
+ {
+ return vector_.size();
+ }
- /// Calls the copy assignment operator of the BaseVector \ref vector
- void copy(Self const& that)
- {
- *vector = that.getVector();
- }
+ /// Resize the \ref vector_ to the size \p s
+ void resize(size_type s)
+ {
+ vector_.resize(s);
+ }
- private:
- FeSpace const& feSpace;
- std::string name;
- ClonablePtr<BaseVector> vector;
+ /// Access the entry \p i of the \ref vector with read-access.
+ value_type const& operator[](size_type i) const
+ {
+ test_exit_dbg(i < vector_.size(),
+ "Index ", i, " out of range [0,", vector_.size(), ")" );
+ return vector_[i];
+ }
- // friend class declarations
- template <class, class>
- friend class SystemVector;
- };
+ /// Access the entry \p i of the \ref vector with write-access.
+ value_type& operator[](size_type i)
+ {
+ test_exit_dbg(i < vector_.size(),
+ "Index ", i, " out of range [0,", vector_.size(), ")" );
+ return vector_[i];
+ }
+
+ private:
+ BaseVector vector_;
+ };
+
+
+ template <class BasisType, class VectorType = double>
+ using DOFVector = DOFVectorBase<BasisType, IstlVector<VectorType>>;
+
+ /// Constructor a dofvector from given basis and name
+ template <class ValueType = double, class Basis>
+ DOFVector<Basis, ValueType>
+ makeDOFVector(Basis const& basis)
+ {
+ return {basis};
+ }
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DirectRunner.hpp b/src/amdis/linear_algebra/istl/DirectRunner.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..187abadfa5c0767810abf3147a4c28d8ef1e0733
--- /dev/null
+++ b/src/amdis/linear_algebra/istl/DirectRunner.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <algorithm>
+#include <string>
+
+#include <amdis/Output.hpp>
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \ingroup Solver
+ * \class AMDiS::DirectRunner
+ * \brief \implements RunnerInterface for the (external) direct solvers
+ */
+ template <class Solver, class Matrix, class VectorX, class VectorB>
+ class DirectRunner
+ : public RunnerInterface<Matrix, VectorX, VectorB>
+ {
+ protected:
+ using Super = RunnerInterface<Matrix, VectorX, VectorB>;
+ using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
+
+ public:
+ /// Constructor.
+ DirectRunner(std::string prefix)
+ : solverCreator_(prefix)
+ {}
+
+ /// Implementes \ref RunnerInterface::init()
+ virtual void init(Matrix const& A) override
+ {
+ solver_ = solverCreator_.create(A);
+ }
+
+ /// Implementes \ref RunnerInterface::exit()
+ virtual void exit() override
+ {
+ solver_.reset();
+ }
+
+ /// Implementes \ref RunnerInterface::solve()
+ virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
+ SolverInfo& solverInfo) override
+ {
+ // storing some statistics
+ Dune::InverseOperatorResult statistics;
+
+ // solve the linear system
+ VectorB _b = b;
+ solver_->apply(x, _b, statistics);
+
+ solverInfo.setRelResidual(statistics.reduction);
+
+ return 0;
+ }
+
+ private:
+ ISTLSolverCreator<Solver> solverCreator_;
+ std::shared_ptr<BaseSolver> solver_;
+ };
+}
diff --git a/src/amdis/linear_algebra/istl/Fwd.hpp b/src/amdis/linear_algebra/istl/Fwd.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2ad22c295e1f6de8e8926c94be5b7bf36fc7f9c
--- /dev/null
+++ b/src/amdis/linear_algebra/istl/Fwd.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+namespace AMDiS
+{
+ template <class ISTLSolver, bool specialized = true>
+ struct ISTLSolverCreator;
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTLRunner.hpp b/src/amdis/linear_algebra/istl/ISTLRunner.hpp
index 46d7790fed6db10f62cb510ea65b0e24f200e266..99c4927e876456f1bfffc9c2586467c8ec6fb69e 100644
--- a/src/amdis/linear_algebra/istl/ISTLRunner.hpp
+++ b/src/amdis/linear_algebra/istl/ISTLRunner.hpp
@@ -1,56 +1,45 @@
#pragma once
-#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
-#include <amdis/linear_algebra/istl/Preconditioner.hpp>
+#include <dune/istl/preconditioner.hh>
+
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+#include <amdis/linear_algebra/istl/Fwd.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
namespace AMDiS
{
- template <class Matrix, class VectorX, class VectorB>
+ template <class ISTLSolver, class Matrix, class VectorX, class VectorB>
class ISTLRunner
- : public RunnerInterface<Matrix, VectorX, VectorB>
+ : public RunnerInterface<Matrix, VectorX, VectorB>
{
using LinOperator = Dune::MatrixAdapter<Matrix, VectorX, VectorB>;
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
-
- using PreconRunner = ISTLPreconditioner<Matrix, VectorX, VectorB>;
- using PreconBase = PreconditionerInterface<Matrix, VectorX, VectorB>;
+ using BasePrecon = Dune::Preconditioner<VectorX, VectorB>;
+ using ISTLPrecon = ISTLPreconInterface<Matrix, VectorX, VectorB>;
public:
- ISTLRunner(std::string solverName, std::string prefix)
- : solverName(solverName)
- , prefix(prefix)
- , solverCreator(prefix)
+ ISTLRunner(std::string prefix)
+ : solverCreator_(prefix)
{
- // get creator string for preconditioner
- std::string preconName = "no";
- Parameters::get(prefix + "->left precon", preconName);
- if (preconName.empty() || preconName == "no")
- Parameters::get(prefix + "->right precon", preconName);
- if (preconName.empty() || preconName == "no")
- Parameters::get(prefix + "->precon->name", preconName);
-
- precon = std::make_shared<PreconRunner>(preconName, prefix + "->precon");
+ initPrecon(prefix);
}
/// Implementes \ref RunnerInterface::init()
virtual void init(Matrix const& A) override
{
- precon->init(A);
-
- Matrix& _A = const_cast<Matrix&>(A); // Unschoen!!!
- linOperator = std::make_shared<LinOperator>(_A);
- solver = solverCreator.create(solverName, *linOperator, *precon);
+ precon_ = preconCreator_->create(A);
+ linOperator_ = std::make_shared<LinOperator>(A);
+ solver_ = solverCreator_.create(*linOperator_, *precon_);
}
/// Implementes \ref RunnerInterface::exit()
virtual void exit() override
{
- solver.reset();
- linOperator.reset();
-
- precon->exit();
- precon.reset();
+ solver_.reset();
+ linOperator_.reset();
+ precon_.reset();
}
/// Implementes \ref RunnerInterface::solve()
@@ -62,26 +51,42 @@ namespace AMDiS
// solve the linear system
VectorB _b = b;
- solver->apply(x, _b, statistics);
+ solver_->apply(x, _b, statistics);
solverInfo.setRelResidual(statistics.reduction);
+
+ return statistics.converged ? 0 : 1;
}
- /// Implementes \ref RunnerInterface::getLeftPrecon()
- virtual std::shared_ptr<PreconBase> getLeftPrecon() override
+ protected:
+ void initPrecon(std::string const& prefix)
{
- return preconAdapter;
+ // get creator string for preconditioner
+ std::string preconName = "no";
+ std::string initFileStr = "";
+ for (std::string postfix : {"left precon", "right precon", "precon->name"}) {
+ Parameters::get(prefix + "->" + postfix, preconName);
+ if (!preconName.empty() && preconName != "no") {
+ initFileStr = prefix + "->" + postfix;
+ break;
+ }
+ }
+
+ if (preconName.empty() || preconName == "no")
+ preconName = "default";
+
+ auto* creator = named(CreatorMap<ISTLPrecon>::getCreator(preconName, initFileStr));
+ preconCreator_ = creator->create(initFileStr);
+ assert(preconCreator_);
}
private:
- std::string solverName;
- std::string prefix;
-
- ISTLSolverCreator<Matrix, VectorX, VectorB> solverCreator;
+ ISTLSolverCreator<ISTLSolver> solverCreator_;
+ std::shared_ptr<ISTLPrecon> preconCreator_;
- std::shared_ptr<PreconRunner> precon;
- std::shared_ptr<LinOperator> linOperator;
- std::shared_ptr<BaseSolver> solver;
+ std::shared_ptr<BasePrecon> precon_;
+ std::shared_ptr<LinOperator> linOperator_;
+ std::shared_ptr<BaseSolver> solver_;
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
index 750b562b2e714b7a5d9720962017a65decbe23be..4a653fa796949ba28b3e3511da3628f4aed5cd1f 100644
--- a/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
+++ b/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
@@ -2,71 +2,116 @@
#include <dune/istl/preconditioners.hh>
+#include <amdis/CreatorInterface.hpp>
+
namespace AMDiS
{
+ template <class Matrix, class VectorX, class VectorB>
+ struct ISTLPreconInterface
+ {
+ using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+ virtual std::unique_ptr<PreconBase> create(Matrix const& matrix) const = 0;
+ };
+
+ template <class> struct Type {};
// creators for preconditioners, depending on matrix-type
// ---------------------------------------------------------------------------
- template <class Matrix, class VectorX, class VectorB = VectorX>
- class ISTLPreconCreator;
-
- template <class... MatrixRows, class VectorX, class VectorB>
- class ISTLPreconCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
+ template <class Precon, class Matrix>
+ struct ISTLPrecon
+ : public ISTLPreconInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
{
- using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
+ using VectorX = typename Precon::domain_type;
+ using VectorB = typename Precon::range_type;
+ using Super = ISTLPreconInterface<Matrix, VectorX, VectorB>;
+ using Self = ISTLPrecon;
- public:
- ISTLPreconCreator(std::string prefix)
- : prefix(prefix)
- {}
+ struct Creator : CreatorInterfaceName<Super>
+ {
+ virtual std::shared_ptr<Super> create(std::string prefix) override
+ {
+ return std::make_shared<Self>(prefix);
+ }
+ };
- template <class Matrix>
- std::unique_ptr<BasePreconditioner> create(std::string /*name*/, Matrix const& A)
+ ISTLPrecon(std::string prefix)
{
- double w = 1.0;
- Parameters::get(prefix + "->relaxation", w);
+ Parameters::get(prefix + "->relaxation", w_);
+ }
- warning("Only Richardson preconditioner available for MultiTypeBlockMatrix!");
- return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
+ using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+ virtual std::unique_ptr<PreconBase> create(Matrix const& A) const override
+ {
+ return createImpl(A, Type<Precon>{});
}
private:
- std::string prefix;
+ template <class P>
+ std::unique_ptr<P> createImpl(Matrix const& A, Type<P>) const
+ {
+ return std::make_unique<P>(A, 1, w_);
+ }
+
+ std::unique_ptr<Dune::SeqILU0<Matrix, VectorX, VectorB>>
+ createImpl(Matrix const& A, Type<Dune::SeqILU0<Matrix, VectorX, VectorB>>) const
+ {
+ return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w_);
+ }
+
+ std::unique_ptr<Dune::Richardson<VectorX, VectorB>>
+ createImpl(Matrix const& /*A*/, Type<Dune::Richardson<VectorX, VectorB>>) const
+ {
+ return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w_);
+ }
+
+ double w_ = 1.0;
};
- template <class Block, class Alloc, class VectorX, class VectorB>
- class ISTLPreconCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
+
+ /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
+ /**
+ * Adds creators for full-matrix aware solvers.
+ * - *cg*: conjugate gradient method, \see Dune::CGSolver
+ * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
+ * - *minres*: Minimal residul method, \see Dune::MINRESSolver
+ * - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
+ * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
+ **/
+ template <class Matrix, class VectorX, class VectorB>
+ class DefaultCreators<ISTLPreconInterface<Matrix, VectorX, VectorB>>
{
- using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
+ template <class Precon>
+ using PreconCreator
+ = typename ISTLPrecon<Precon, Matrix>::Creator;
- public:
- ISTLPreconCreator(std::string prefix)
- : prefix(prefix)
- {}
+ using Map = CreatorMap<ISTLPreconInterface<Matrix, VectorX, VectorB>>;
- template <class Matrix>
- std::unique_ptr<BasePreconditioner> create(std::string name, Matrix const& A)
+ public:
+ static void init()
{
- double w = 1.0;
- Parameters::get(prefix + "->relaxation", w);
-
- if (name == "diag" || name == "jacobi")
- return std::make_unique<Dune::SeqJac<Matrix, VectorX, VectorB>>(A, 1, w);
- else if (name == "gs" || name == "gauss_seidel")
- return std::make_unique<Dune::SeqGS<Matrix, VectorX, VectorB>>(A, 1, w);
- else if (name == "sor")
- return std::make_unique<Dune::SeqSOR<Matrix, VectorX, VectorB>>(A, 1, w);
- else if (name == "ssor")
- return std::make_unique<Dune::SeqSSOR<Matrix, VectorX, VectorB>>(A, 1, w);
- else if (name == "ilu" || name == "ilu0")
- return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w);
- else
- return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
- }
+ auto jacobi = new PreconCreator<Dune::SeqJac<Matrix, VectorX, VectorB>>;
+ Map::addCreator("diag", jacobi);
+ Map::addCreator("jacobi", jacobi);
- private:
- std::string prefix;
+ auto gs = new PreconCreator<Dune::SeqGS<Matrix, VectorX, VectorB>>;
+ Map::addCreator("gs", gs);
+ Map::addCreator("gauss_seidel", gs);
+
+ auto sor = new PreconCreator<Dune::SeqSOR<Matrix, VectorX, VectorB>>;
+ Map::addCreator("sor", sor);
+
+ auto ssor = new PreconCreator<Dune::SeqSSOR<Matrix, VectorX, VectorB>>;
+ Map::addCreator("ssor", ssor);
+
+ auto ilu = new PreconCreator<Dune::SeqILU0<Matrix, VectorX, VectorB>>;
+ Map::addCreator("ilu", ilu);
+ Map::addCreator("ilu0", ilu);
+
+ auto richardson = new PreconCreator<Dune::Richardson<VectorX, VectorB>>;
+ Map::addCreator("richardson", richardson);
+ Map::addCreator("default", richardson);
+ }
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTL_Solver.hpp b/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
index 558677dbdf152e31750adcbc4efefe9604381d15..e32d52fe3c726a57a951e40ad99819cc59d542a6 100644
--- a/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
+++ b/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
@@ -1,106 +1,196 @@
#pragma once
+#include <memory>
+
#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>
+#include <dune/istl/superlu.hh>
+
+#include <amdis/CreatorMap.hpp>
+#include <amdis/Initfile.hpp>
+
+#include <amdis/linear_algebra/istl/Fwd.hpp>
+#include <amdis/linear_algebra/istl/DirectRunner.hpp>
+#include <amdis/linear_algebra/istl/ISTLRunner.hpp>
namespace AMDiS
{
- // creators for solvers, depending on matrix-type
- // ---------------------------------------------------------------------------
+ template <class ISTLSolver, bool specialized>
+ struct ISTLSolverCreator
+ {
+ ISTLSolverCreator(std::string prefix)
+ {
+ Parameters::get(prefix + "->info", info_);
+ Parameters::get(prefix + "->max iteration", maxIter_);
+ Parameters::get(prefix + "->relative tolerance", rTol_);
+ }
- template <class Matrix, class VectorX, class VectorB = VectorX>
- class ISTLSolverCreator;
+ template <class LinOperator, class Precon>
+ std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
+ {
+ return std::make_unique<ISTLSolver>(A, P, rTol_, maxIter_, info_);
+ }
+
+ int info_ = 2;
+ std::size_t maxIter_ = 500;
+ double rTol_ = 1.e-6;
+ };
- template <class... MatrixRows, class VectorX, class VectorB>
- class ISTLSolverCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
+ template <class VectorX>
+ struct ISTLSolverCreator<Dune::RestartedGMResSolver<VectorX>, true>
+ : public ISTLSolverCreator<Dune::RestartedGMResSolver<VectorX>, false>
{
- using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
+ using ISTLSolver = Dune::RestartedGMResSolver<VectorX>;
+ using Super = ISTLSolverCreator<ISTLSolver, false>;
- public:
ISTLSolverCreator(std::string prefix)
- : prefix(prefix)
- , info(2)
+ : Super(prefix)
{
- Parameters::get(prefix + "->info", info);
+ Parameters::get(prefix + "->restart", restart_);
}
- template <class Matrix, class Precon>
- std::unique_ptr<BaseSolver> create(std::string name, Matrix& A, Precon& P)
+ template <class LinOperator, class Precon>
+ std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
{
- std::size_t maxIter = 500;
- Parameters::get(prefix + "->max iteration", maxIter);
-
- double rTol = 1.e-6;
- Parameters::get(prefix + "->reduction", rTol);
-
- if (name == "cg")
- return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "bicgstab")
- return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "minres")
- return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "gmres")
- {
- int restart = 30;
- Parameters::get(prefix + "->restart", restart);
- return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
- }
- else
- error_exit("Unknown solver name!");
+ return std::make_unique<ISTLSolver>(A, P, this->rTol_, restart_, this->maxIter_, this->info_);
}
- private:
- std::string prefix;
- int info;
+ int restart_ = 30;
};
-
- template <class Block, class Alloc, class VectorX, class VectorB>
- class ISTLSolverCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
+ template <class VectorX>
+ struct ISTLSolverCreator<Dune::GeneralizedPCGSolver<VectorX>, true>
+ : public ISTLSolverCreator<Dune::GeneralizedPCGSolver<VectorX>, false>
{
- using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
- using Matrix = Dune::BCRSMatrix<Block,Alloc>;
+ using ISTLSolver = Dune::GeneralizedPCGSolver<VectorX>;
+ using Super = ISTLSolverCreator<ISTLSolver, false>;
- public:
ISTLSolverCreator(std::string prefix)
- : prefix(prefix)
- , info(2)
+ : Super(prefix)
{
- Parameters::get(prefix + "->info", info);
+ Parameters::get(prefix + "->restart", restart_);
}
template <class LinOperator, class Precon>
- std::unique_ptr<BaseSolver> create(std::string name, LinOperator& A, Precon& P)
+ std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
{
- std::size_t maxIter = 500;
- Parameters::get(prefix + "->max iteration", maxIter);
-
- double rTol = 1.e-6;
- Parameters::get(prefix + "->relative tolerance", rTol);
-
- if (name == "cg")
- return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "bicgstab")
- return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "minres")
- return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
- else if (name == "gmres")
- {
- int restart = 30;
- Parameters::get(prefix + "->restart", restart);
- return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
- }
+ return std::make_unique<ISTLSolver>(A, P, this->rTol_, this->maxIter_, this->info_, restart_);
+ }
+
+ int restart_ = 30;
+ };
+
+
#if HAVE_SUITESPARSE_UMFPACK
- else if (name == "umfpack" || name == "direct")
- return std::make_unique<Dune::UMFPack<Matrix>>(A, info);
+ template <class Matrix>
+ struct ISTLSolverCreator<Dune::UMFPack<Matrix>, true>
+ {
+ using Solver = Dune::UMFPack<Matrix>;
+ ISTLSolverCreator(std::string prefix)
+ {
+ Parameters::get(prefix + "->info", verbose_);
+ }
+
+ std::unique_ptr<Solver> create(Matrix const& A) const
+ {
+ return std::make_unique<Solver>(A, verbose_);
+ }
+
+ private:
+ int verbose_ = 2;
+ };
#endif
- else
- error_exit("Unknown solver name!");
+
+#if HAVE_SUPERLU
+ template <class Matrix>
+ struct ISTLSolverCreator<Dune::SuperLU<Matrix>, true>
+ {
+ using Solver = Dune::SuperLU<Matrix>;
+ ISTLSolverCreator(std::string prefix)
+ {
+ int info = 2;
+ Parameters::get(prefix + "->info", info);
+ verbose_ = (info != 0);
+
+ Parameters::get(prefix + "->reuse vector", reuseVector_);
+ }
+
+ std::unique_ptr<Solver> create(Matrix const& A) const
+ {
+ return std::make_unique<Solver>(A, verbose_, reuseVector_);
}
private:
- std::string prefix;
- int info;
+ bool verbose_ = true;
+ bool reuseVector_ = true;
+ };
+#endif
+
+ /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
+ /**
+ * Adds creators for full-matrix aware solvers.
+ * - *cg*: conjugate gradient method, \see Dune::CGSolver
+ * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
+ * - *minres*: Minimal residul method, \see Dune::MINRESSolver
+ * - *gmres*: Generalized minimal residual method, \see Dune::RestartedGMResSolver
+ * - *fcg*: Generalized preconditioned conjugate gradient solver, \see Dune::GeneralizedPCGSolver
+ * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
+ * - *superlu*: external SuperLU solver, \see Dune::SuperLU
+ **/
+ template <class Block, class Alloc, class VectorX, class VectorB>
+ class DefaultCreators<LinearSolverInterface<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>>
+ {
+ using Matrix = Dune::BCRSMatrix<Block,Alloc>;
+ using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
+
+ template <class ISTLSolver>
+ using SolverCreator
+ = typename LinearSolver<Matrix, VectorX, VectorB,
+ ISTLRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+
+ template <class ISTLSolver>
+ using DirectSolverCreator
+ = typename LinearSolver<Matrix, VectorX, VectorB,
+ DirectRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+
+ using Map = CreatorMap<SolverBase>;
+
+ public:
+ static void init()
+ {
+ auto cg = new SolverCreator<Dune::CGSolver<VectorX>>;
+ Map::addCreator("cg", cg);
+
+ auto bicgstab = new SolverCreator<Dune::BiCGSTABSolver<VectorX>>;
+ Map::addCreator("bicgstab", bicgstab);
+ Map::addCreator("bcgs", bicgstab);
+
+ auto minres = new SolverCreator<Dune::MINRESSolver<VectorX>>;
+ Map::addCreator("minres", minres);
+
+ auto gmres = new SolverCreator<Dune::RestartedGMResSolver<VectorX>>;
+ Map::addCreator("gmres", gmres);
+
+ // Generalized preconditioned conjugate gradient solver.
+ auto fcg = new SolverCreator<Dune::GeneralizedPCGSolver<VectorX>>;
+ Map::addCreator("fcg", fcg);
+
+#if HAVE_SUITESPARSE_UMFPACK
+ auto umfpack = new DirectSolverCreator<Dune::UMFPack<Matrix>>;
+ Map::addCreator("umfpack", umfpack);
+#endif
+
+#if HAVE_SUPERLU
+ auto superlu = new DirectSolverCreator<Dune::SuperLU<Matrix>>;
+ Map::addCreator("superlu", superlu);
+#endif
+
+#if HAVE_SUITESPARSE_UMFPACK
+ Map::addCreator("direct", umfpack);
+#elif HAVE_SUPERLU
+ Map::addCreator("direct", superlu);
+#endif
+ }
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/LinearSolver.hpp b/src/amdis/linear_algebra/istl/LinearSolver.hpp
deleted file mode 100644
index 7799eaa29ae2c1d1407bf049d3f2d6d7f87043c8..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/LinearSolver.hpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#pragma once
-
-#include <amdis/linear_algebra/LinearSolverInterface.hpp>
-#include <amdis/linear_algebra/istl/ISTLRunner.hpp>
-
-namespace AMDiS
-{
- // A general solver class for istl solvers
- // ---------------------------------------------------------------------------
-
- template <class Matrix, class VectorX, class VectorB, class Runner>
- class LinearSolver
- : public LinearSolverInterface<Matrix, VectorX, VectorB>
- {
- using Super = LinearSolverInterface<Matrix, Vector, Vector>;
- using RunnerBase = typename Super::RunnerBase;
-
- public:
- /// Constructor
- LinearSolver(std::string name, std::string prefix)
- : runner(std::make_shared<Runner>(name, prefix))
- {}
-
-
- private:
- /// Implements \ref LinearSolverInterface::solveImpl()
- virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
- SolverInfo& solverInfo) override
- {
- if (solverInfo.doCreateMatrixData()) {
- // init matrix or wrap block-matrix or ...
- runner->init(A);
- }
-
- // solve the linear system
- int error = runner->solve(A, x, b, solverInfo);
- solverInfo.setError(error);
-
- if (!solverInfo.doStoreMatrixData())
- runner->exit();
- }
-
- private:
- /// redirect the implementation to a runner. Implements a \ref RunnerInterface
- std::shared_ptr<Runner> runner;
- };
-
-
- template <class Matrix, class VectorX, class VectorB>
- class SolverCreator<Matrix, VectorX, VectorB>
- {
- using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
-
- template <class Runner>
- using Solver = LinearSolver<Matrix, VectorX, VectorB, Runner>;
-
- /// Instantiate a new linear solver
- static std::unique_ptr<SolverBase> create(std::string name, std::string prefix)
- {
- using ISTL = ISTLRunner<Matrix, VectorX, VectorB>;
-
- if (name == "cg" || name == "bicgstab" || name == "minres" || name == "gmres")
- return std::make_unique<Solver<ISTL>>(name, prefix);
- else
- error_exit("Unknown solver name!"); // TODO: add direct solver interface here
- }
- };
-
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/Preconditioner.hpp b/src/amdis/linear_algebra/istl/Preconditioner.hpp
deleted file mode 100644
index 6c331d53e789198dee9dd064479ef48aca5dd97c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/Preconditioner.hpp
+++ /dev/null
@@ -1,80 +0,0 @@
-#pragma once
-
-#include <type_traits>
-
-#include <dune/istl/preconditioner.hh>
-#include <dune/istl/solvercategory.hh>
-
-#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
-
-namespace AMDiS {
-
-namespace detail
-{
- template <class T, class = typename T::category >
- std::integral_constant<bool, true> HasCategoryImpl(int);
-
- template <class T>
- std::integral_constant<bool, false> HasCategoryImpl(...);
-
- template <class T>
- using HasCategory = decltype(HasCategoryImpl<T>(int{}));
-
-} // end namespace detail
-
-
-template <class Matrix, class VectorX, class VectorB>
-class ISTLPreconditioner
- : public Dune::Preconditioner<VectorX, VectorB>
- , public PreconditionerInterface<Matrix, VectorX, VectorB>
-{
-public:
- enum {category=Dune::SolverCategory::sequential};
-
- ISTLPreconditioner(std::string preconName, std::string prefix)
- : preconName(preconName)
- , prefix(prefix)
- , preconCreator(prefix)
- {}
-
- /// Implementation of \ref Dune::Preconditioner::pre()
- virtual void pre(VectorX& x, VectorB& b) override { runner.pre(x, b); }
-
- /// Implementation of \ref Dune::Preconditioner::apply()
- virtual void apply(VectorX& v, const VectorB& d) override { runner.apply(v, d); }
-
- /// Implementation of \ref Dune::Preconditioner::post()
- virtual void post(VectorX& x) override { runner.post(x); }
-
- // -------------------------------------------------------------------------
-
- /// Implementation of \ref PreconditionerInterface::init()
- virtual void init(Matrix const& A) override
- {
- Matrix& _A = const_cast<Matrix&>(A); // Unschoen!!!
-
- runner = preconCreator.create(preconName, _A);
- };
-
- /// Implementation of \ref PreconditionerInterface::exit()
- virtual void exit() override
- {
- runner.reset();
- };
-
- /// Implementation of \ref PreconditionerInterface::solve()
- virtual void solve(VectorX const& x, VectorX& y) const override
- {
- runner.apply(v, d);
- }
-
-private:
- std::string preconName;
- std::string prefix;
-
- ISTLPreconCreator<Matrix, VectorX, VectorB> preconCreator;
-
- std::shared_ptr<PreconRunner> runner;
-};
-
-} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/linear_algebra/istl/SystemMatrix.cpp b/src/amdis/linear_algebra/istl/SystemMatrix.cpp
deleted file mode 100644
index 26d4d975b2d79c524632b40f846bbb0c2b84c0cc..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemMatrix.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMMATRIX
-#include <amdis/linear_algebra/istl/SystemMatrix.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMMATRIX
-
-
-namespace AMDiS
-{
- // explicit template instatiation
- template class SystemMatrix<typename TestTraits<2,1>::FeSpaces>;
- template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemMatrix.hpp b/src/amdis/linear_algebra/istl/SystemMatrix.hpp
deleted file mode 100644
index fa3e50b219658eef03b2ea288132132eba7c9f9c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemMatrix.hpp
+++ /dev/null
@@ -1,318 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <dune/istl/multitypeblockmatrix.hh>
-#include <dune/istl/multitypeblockvector.hh>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
- namespace Impl
- {
- // DOFMatrices = std::tuple< std::tuple< DOFMatrix<RowFeSpace, ColFeSpace, ValueType>... >... >
- template <class DOFMatrices>
- struct BuildDOFMatrices
- {
- template <std::size_t R>
- using DOFMatrixRow = std::tuple_element_t<R, DOFMatrices>;
-
- template <std::size_t R, std::size_t C>
- using DOFMatrixRowCol = std::tuple_element_t<C, DOFMatrixRow<R>>;
-
- // add arg to repeated constructor argument list
- template <std::size_t... R, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
- static DOFMatrices make(Indices<R...>,
- RowFeSpaces&& rowFeSpaces,
- ColFeSpaces&& colFeSpace,
- MultiMatrix&& multiMatrix)
- {
- return DOFMatrices{make_row(
- index_<R>,
- MakeSeq_t<std::tuple_size<std::decay_t<ColFeSpaces>>::value>(),
- std::get<R>(std::forward<RowFeSpaces>(rowFeSpaces)),
- std::forward<ColFeSpaces>(colFeSpace),
- std::get<R>(std::forward<MultiMatrix>(multiMatrix))
- )...};
- }
-
- template <std::size_t R, std::size_t... C, class RowFeSpace, class ColFeSpaces, class MultiVector>
- static DOFMatrixRow<R> make_row(index_t<R>,
- Indices<C...>,
- RowFeSpace&& rowFeSpace,
- ColFeSpaces&& colFeSpace,
- MultiVector&& multiVector)
- {
- return DOFMatrixRow<R>{DOFMatrixRowCol<R,C>(
- std::forward<RowFeSpace>(rowFeSpace),
- std::get<C>(std::forward<ColFeSpaces>(colFeSpace)),
- "matrix[" + std::to_string(R) + "][" + std::to_string(C) + "]",
- std::get<C>(std::forward<MultiVector>(multiVector))
- )...};
- }
- };
-
- template <class DOFMatrices, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
- DOFMatrices buildDOFMatrices(RowFeSpaces&& rowFeSpaces,
- ColFeSpaces&& colFeSpace,
- MultiMatrix&& multiMatrix)
- {
- return BuildDOFMatrices<DOFMatrices>::make(
- MakeSeq_t<std::tuple_size<std::decay_t<RowFeSpaces>>::value>(), // traverse rows first
- std::forward<RowFeSpaces>(rowFeSpaces),
- std::forward<ColFeSpaces>(colFeSpace),
- std::forward<MultiMatrix>(multiMatrix));
- }
-
- } // end namespace Impl
-
-
-
-
- /// \brief Container that repesents multiple data-Vectors of different value types.
- /**
- * Represents a \ref Dune::MultiTypeBlockVector + a tuple of corresponding
- * feSpaces. This I'th \ref Dune::BlockVector compbined with the I'th FeSpace
- * builds a \ref DOFVector and can be return by the \ref operator[].
- *
- * By default, the ValueTypes are all \ref Dune::FieldMatrix of 1x1 double entry.
- **/
- template <class FeSpaces,
- class FieldType = double,
- class Dimensions = Repeat_t<std::tuple_size<FeSpaces>::value, index_t<1>> >
- class SystemMatrix
- {
- template <class RowFeSpace, class RowDim>
- struct MultiMatrixRowGenerator
- {
- template <class ColDim>
- using ValueTypeGenerator = Dune::FieldMatrix<FieldType, RowDim::value, ColDim::value>;
- using ValueTypes = MakeTuple_t<ValueTypeGenerator, Dimensions>;
-
- template <class ColFeSpace, class VT>
- using DOFMatrixGenerator = DOFMatrix<RowFeSpace, ColFeSpace, VT>;
- using DOFMatrices = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, ValueTypes>;
-
- template <class VT>
- using BaseMatrixGenerator = Dune::BCRSMatrix<VT, std::allocator<VT>>;
- using BaseMatrices = MakeTuple_t<BaseMatrixGenerator, ValueTypes>;
-
- using MultiVector = ExpandTuple_t<Dune::MultiTypeBlockVector, BaseMatrices>;
- };
-
- template <class RowFeSpace, class RowDim>
- using MultiVectorGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::MultiVector;
-
- template <class RowFeSpace, class RowDim>
- using ValueTypeGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::ValueTypes;
-
- template <class RowFeSpace, class RowDim>
- using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::DOFMatrices;
-
- public:
- using MultiVectors = MakeTuple2_t<MultiVectorGenerator, FeSpaces, Dimensions>;
- using ValueTypes = MakeTuple2_t<ValueTypeGenerator, FeSpaces, Dimensions>;
- using DOFMatrices = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, Dimensions>;
- using MultiMatrix = ExpandTuple_t<Dune::MultiTypeBlockMatrix, MultiVectors>;
-
- static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
-
- AMDIS_STATIC_ASSERT( nComponents > 0 );
- AMDIS_STATIC_ASSERT( nComponents == std::tuple_size<Dimensions>::value );
-
- /// Constructor that constructs new base-vectors
- SystemMatrix(FeSpaces const& feSpaces)
- : feSpaces(feSpaces)
- , matrix{}
- , dofmatrices(Impl::buildDOFMatrices<DOFMatrices>(feSpaces, feSpaces, matrix))
- {}
-
-
- /// Return the number of vector entries
- static constexpr std::size_t N()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- static constexpr std::size_t M()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t R, std::size_t C>
- auto& getMatrix(const index_t<R> = {}, const index_t<C> = {})
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(matrix));
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t R, std::size_t C>
- auto const& getMatrix(const index_t<R> = {}, const index_t<C> = {}) const
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(matrix));
- }
-
- /// Return the underlying multi vector
- MultiMatrix& getMatrix()
- {
- return matrix;
- }
-
- MultiMatrix const& getMatrix() const
- {
- return matrix;
- }
-
- /// Return the I'th finite element space
- template <std::size_t I = 0>
- auto& getRowFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < N(), "Index out of range [0,N)");
- return std::get<I>(feSpaces);
- }
-
- template <std::size_t I = 0>
- auto const& getColFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < M(), "Index out of range [0,M)");
- return std::get<I>(feSpaces);
- }
-
- template <std::size_t R = 0, std::size_t C = 0>
- auto& operator()(const index_t<R> = {}, const index_t<C> = {})
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(dofmatrices));
- }
-
- /// Create a DOFVector with references to the unerlaying data
- template <std::size_t R = 0, std::size_t C = 0>
- auto const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(dofmatrices));
- }
-
- private:
- FeSpaces const& feSpaces; ///< a tuple of feSpaces
- MultiMatrix matrix; ///< a tuple of data-vectors
- DOFMatrices dofmatrices;
- };
-
-
- // specialization for 1 component
- // -------------------------------------------------------------------------
-
- template <class FeSpace,
- class FieldType,
- int dim>
- class SystemMatrix<std::tuple<FeSpace>, FieldType, std::tuple<index_t<dim>>>
- {
- using ValueType = Dune::FieldMatrix<FieldType, dim, dim>;
- using DOFMatrixType = DOFMatrix<FeSpace, FeSpace, ValueType>;
-
- public:
- using MultiMatrix = typename DOFMatrixType::BaseMatrix;
-
- static constexpr std::size_t nComponent = 1;
-
-
- /// Constructor that constructs new base-vectors
- template <class FeSpaces>
- SystemMatrix(FeSpaces const& feSpaces)
- : dofmatrix(std::get<0>(feSpaces), std::get<0>(feSpaces), "")
- {}
-
-
- /// Return the number of vector entries
- static constexpr std::size_t N()
- {
- return 1;
- }
-
- static constexpr std::size_t M()
- {
- return 1;
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t R = 0, std::size_t C = 0>
- MultiMatrix& getMatrix(const index_t<R> = {}, const index_t<C> = {})
- {
- static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
- return dofmatrix.getMatrix();
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t R = 0, std::size_t C = 0>
- MultiMatrix const& getMatrix(const index_t<R> = {}, const index_t<C> = {}) const
- {
- static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
- return dofmatrix.getMatrix();
- }
-
- /// Return the underlying multi vector
- MultiMatrix const& getMatrix() const
- {
- return dofmatrix.getMatrix();
- }
-
- MultiMatrix& getMatrix()
- {
- return dofmatrix.getMatrix();
- }
-
- /// Return the I'th finite element space
- template <std::size_t I = 0>
- FeSpace const& getRowFeSpace(const index_t<I> = {}) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofmatrix.getRowFeSpace();
- }
-
- template <std::size_t I = 0>
- FeSpace const& getColFeSpace(const index_t<I> = {}) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofmatrix.getColFeSpace();
- }
-
- template <std::size_t R = 0, std::size_t C = 0>
- DOFMatrixType& operator()(const index_t<R> = {}, const index_t<C> = {})
- {
- static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
- return dofmatrix;
- }
-
- /// Create a DOFVector with references to the unerlaying data
- template <std::size_t R = 0, std::size_t C = 0>
- DOFMatrixType const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
- {
- static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
- return dofmatrix;
- }
-
- private:
- DOFMatrixType dofmatrix;
- };
-
-#ifndef AMDIS_NO_EXTERN_SYSTEMMATRIX
- // explicit instantiation in SystemMatrix.cpp
- extern template class SystemMatrix<typename TestTraits<2,1>::FeSpaces>;
- extern template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-#endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemVector.cpp b/src/amdis/linear_algebra/istl/SystemVector.cpp
deleted file mode 100644
index fc04f99bd3495ed19e973913814e6a57c97732c0..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemVector.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMVECTOR
-#include <amdis/linear_algebra/istl/SystemVector.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMVECTOR
-
-
-namespace AMDiS
-{
- // explicit template instatiation
- template class SystemVector<typename TestTraits<2,1>::FeSpaces>;
- template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemVector.hpp b/src/amdis/linear_algebra/istl/SystemVector.hpp
deleted file mode 100644
index 7dfcbb2bad000c96ee0a7fd6ee8e5daf012a08e4..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemVector.hpp
+++ /dev/null
@@ -1,383 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <dune/istl/multitypeblockvector.hh>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-
-#include <amdis/linear_algebra/istl/DOFVector.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
-
- namespace Impl
- {
- // DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
- template <class DOFVectors>
- struct BuildDOFVectors
- {
- template <std::size_t I>
- using DOFVectorIdx = std::tuple_element_t<I, DOFVectors>;
-
- template <std::size_t... I, class FeSpaces, class MultiVector>
- static DOFVectors make(Indices<I...>,
- FeSpaces&& feSpaces,
- std::vector<std::string> const& names,
- MultiVector&& multiVector)
- {
- return DOFVectors{DOFVectorIdx<I>(
- std::get<I>(std::forward<FeSpaces>(feSpaces)),
- names[I],
- std::get<I>(std::forward<MultiVector>(multiVector))
- )...};
- }
- };
-
- template <class DOFVectors, class FeSpaces, class MultiVector>
- DOFVectors buildDOFVectors(FeSpaces&& feSpaces,
- std::vector<std::string> const& names,
- MultiVector&& multiVector)
- {
- return BuildDOFVectors<DOFVectors>::make(
- MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
- std::forward<FeSpaces>(feSpaces),
- names,
- std::forward<MultiVector>(multiVector));
- }
-
- } // end namespace Impl
-
-
- /// \brief Container that repesents multiple data-Vectors of different value types.
- /**
- * Represents a \ref Dune::MultiTypeBlockVector + a tuple of corresponding
- * feSpaces. This I'th \ref Dune::BlockVector compbined with the I'th FeSpace
- * builds a \ref DOFVector and can be return by the \ref operator[].
- *
- * By default, the ValueTypes are all \ref Dune::FieldVector of 1 double entry.
- **/
- template <class FeSpaces,
- class ValueTypes = Repeat_t<std::tuple_size<FeSpaces>::value, Dune::FieldVector<double,1>> >
- class SystemVector
- {
- using Self = SystemVector;
-
- public:
- template <class T> using BlockVectorGenerator = Dune::BlockVector<T, std::allocator<T>>;
-
- using BaseVectors = MakeTuple_t<BlockVectorGenerator, ValueTypes>;
- using MultiVector = ExpandTuple_t<Dune::MultiTypeBlockVector, BaseVectors>;
-
- using DOFVectors = MakeTuple2_t<DOFVector, FeSpaces, ValueTypes>;
-
- static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
-
- AMDIS_STATIC_ASSERT( nComponents > 0 );
- AMDIS_STATIC_ASSERT( nComponents == std::tuple_size<ValueTypes>::value );
-
- /// Constructor that constructs new base-vectors
- SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
- : feSpaces(feSpaces)
- , names(names)
- , vector{}
- , dofvectors(Impl::buildDOFVectors<DOFVectors>(feSpaces, names, vector))
- {
- compress();
- }
-
-
- /// Return the number of vector entries
- static constexpr std::size_t size()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- /// Return number of elements
- int count() const
- {
- return int(size());
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t I>
- auto& getVector(const index_t<I> = {})
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(vector);
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t I>
- auto const& getVector(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(vector);
- }
-
- /// Return the underlying multi vector
- MultiVector& getVector()
- {
- return vector;
- }
-
- MultiVector const& getVector() const
- {
- return vector;
- }
-
- /// Return the I'th finite element space
- template <std::size_t I = 0>
- auto const& getFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(feSpaces);
- }
-
- /// Return the name of the i'th DOFVector
- std::string getName(std::size_t i) const
- {
- return names[i];
- }
-
- template <std::size_t I>
- auto& operator[](const index_t<I>)
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- template <std::size_t I>
- auto const& operator[](const index_t<I>) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- template <std::size_t I>
- auto& getDOFVector(const index_t<I> = {})
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- template <std::size_t I>
- auto const& getDOFVector(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- /// Resize the I'th \ref vector to the sizes of the I'th \ref feSpaces.
- template <std::size_t I = 0>
- void compress(const index_t<I> = {})
- {
- std::get<I>(dofvectors).compress();
- }
-
- /// Resize the \ref vectors to the sizes of the \ref feSpaces.
- void compress()
- {
- forEach(range_<0, size()>, [this](const auto I) {
- std::get<I>(dofvectors).compress();
- });
- }
-
- void copy(Self const& that)
- {
- forEach(range_<0, size()>, [this, &that](const auto I) {
- std::get<I>(dofvectors).copy(that[I]);
- });
- }
-
- template <class Scalar>
- std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
- operator*=(Scalar s)
- {
- forEach(range_<0, size()>, [this, s](const auto I) {
- std::get<I>(vector) *= s;
- });
- return *this;
- }
-
- template <class Scalar>
- std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
- operator=(Scalar s)
- {
- forEach(range_<0, size()>, [this, s](const auto I) {
- std::get<I>(vector) = s;
- });
- return *this;
- }
-
- private:
- FeSpaces const& feSpaces; ///< a tuple of feSpaces
- std::vector<std::string> names;
-
- MultiVector vector; ///< a tuple of data-vectors
- DOFVectors dofvectors; ///< a tuple of dofvectors referencing \ref vector
- };
-
-
- // specialization for the case of only 1 component.
-
- template <class FeSpace, class ValueType>
- class SystemVector<std::tuple<FeSpace>, std::tuple<ValueType>>
- {
- using Self = SystemVector;
-
- public:
- using MultiVector = Dune::BlockVector<ValueType>;
- using DOFVectors = DOFVector<FeSpace, ValueType>;
-
- static constexpr std::size_t nComponent = 1;
-
- /// Constructor that constructs new base-vectors
- template <class FeSpaces>
- SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
- : dofvector(std::get<0>(feSpaces), names.front())
- {
- compress();
- }
-
-
- /// Return the number of vector entries
- static constexpr std::size_t size()
- {
- return 1;
- }
-
- /// Return number of elements
- int count() const
- {
- return 1;
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t I = 0>
- MultiVector& getVector(const index_t<I> = {})
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector.getVector();
- }
-
- /// Return a shared pointer to the i'th underlaying base vector
- template <std::size_t I = 0>
- MultiVector const& getVector(const index_t<I> = {}) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector.getVector();
- }
-
- /// Return the underlying multi vector
- MultiVector& getVector()
- {
- return dofvector.getVector();
- }
-
- MultiVector const& getVector() const
- {
- return dofvector.getVector();
- }
-
- /// Return the I'th finite element space
- template <std::size_t I = 0>
- auto const& getFeSpace(const index_t<I> = {}) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector.getFeSpace();
- }
-
- /// Return the name of the i'th DOFVector
- std::string getName(std::size_t i = 0) const
- {
- test_exit(i == 0, "Index out of range [0,1)");
- return dofvector.getName();
- }
-
- template <std::size_t I>
- auto& operator[](const index_t<I>)
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector;
- }
-
- template <std::size_t I>
- auto const& operator[](const index_t<I>) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector;
- }
-
- template <std::size_t I>
- auto& getDOFVector(const index_t<I> = {})
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector;
- }
-
- template <std::size_t I>
- auto const& getDOFVector(const index_t<I> = {}) const
- {
- static_assert(I == 0, "Index out of range [0,1)");
- return dofvector;
- }
-
- /// Resize the I'th \ref vector to the sizes of the I'th \ref feSpaces.
- template <std::size_t I = 0>
- void compress(const index_t<I> = {})
- {
- static_assert(I == 0, "Index out of range [0,1)");
- dofvector.compress();
- }
-
- /// Resize the \ref vectors to the sizes of the \ref feSpaces.
- void compress()
- {
- dofvector.compress();
- }
-
- void copy(Self const& that)
- {
- dofvector.copy(that[_0]);
- }
-
- template <class Scalar>
- std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
- operator*=(Scalar s)
- {
- dofvector.getVector() *= s;
- return *this;
- }
-
- template <class Scalar>
- std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
- operator=(Scalar s)
- {
- dofvector.getVector() = s;
- return *this;
- }
-
- private:
- DOFVectors dofvector; ///< a tuple of sofvectors referencing \ref vector
- };
-
-
-
-
-#ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
- // explicit instantiation in SystemVector.cpp
- extern template class SystemVector<typename TestTraits<2,1>::FeSpaces>;
- extern template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-#endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp b/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
deleted file mode 100644
index 8a5739689d56b9406fc696096879fb6a6ff1dc84..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#pragma once
-
-#include "itl/block_diagonal.hpp"
-
-#include <amdis/CreatorMap.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
-
-namespace AMDiS
-{
-
- template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
- class DefaultCreators<PreconditionerInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>>
- {
- using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
- using PreconBase = PreconditionerInterface<Matrix, Vector>;
-
- template <template<class> class ITLPrecon>
- using PreconCreator
- = typename Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>::Creator;
-
- using Map = CreatorMap<PreconBase>;
-
- public:
- static void init()
- {
- auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
- Map::addCreator("diag", pc_diag);
- Map::addCreator("jacobi", pc_diag);
-
- auto pc_id = new PreconCreator<IdentityPreconditioner>;
- Map::addCreator("identity", pc_id);
- Map::addCreator("no", pc_id);
-
- Map::addCreator("default", pc_id);
- }
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/BITL_Solver.hpp b/src/amdis/linear_algebra/mtl/BITL_Solver.hpp
deleted file mode 100644
index c011bdc769633a0d49c109d654763dc56efd6b41..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BITL_Solver.hpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#pragma once
-
-#include <amdis/CreatorMap.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
-
-namespace AMDiS
-{
- /// Adds default creators for linear solvers based on `BlockMTLMatrix`.
- /**
- * Adds creators for block matrix aware solvers.
- * - *cg*: conjugate gradient method, \see cg_solver_type
- * - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
- * - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
- * - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
- * - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
- * - *gmres*: Generalized minimal residula method, \see gmres_type
- * - *fgmres*: Flexible GMRES, \see fgmres_type
- * - *minres*: Minimal residul method, \see minres_solver_type
- * - *gcr*: Generalized conjugate residual method, \see gcr_type
- * - *umfpack*: external UMFPACK solver, \see UmfpackRunner
- **/
- template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
- class DefaultCreators<LinearSolverInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector, Vector>>
- {
- using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
- using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
-
- template <class ITLSolver>
- using SolverCreator
- = typename LinearSolver<Matrix, Vector,
- KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
-
-#ifdef HAVE_UMFPACK
- using UmfpackSolverCreator
- = typename LinearSolver<Matrix, Vector,
- UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
-#endif
-
- using Map = CreatorMap<SolverBase>;
-
- public:
- static void init()
- {
- auto cg = new SolverCreator<cg_solver_type>;
- Map::addCreator("cg", cg);
-
- auto cgs = new SolverCreator<cgs_solver_type>;
- Map::addCreator("cgs", cgs);
-
- auto bcgs = new SolverCreator<bicgstab_type>;
- Map::addCreator("bicgstab", bcgs);
- Map::addCreator("bcgs", bcgs);
-
- auto tfqmr = new SolverCreator<tfqmr_solver_type>;
- Map::addCreator("tfqmr", tfqmr);
-
- auto bcgsl = new SolverCreator<bicgstab_ell_type>;
- Map::addCreator("bicgstab_ell", bcgsl);
- Map::addCreator("bcgsl", bcgsl);
-
- auto gmres = new SolverCreator<gmres_type>;
- Map::addCreator("gmres", gmres);
-
- auto fgmres = new SolverCreator<fgmres_type>;
- Map::addCreator("fgmres", fgmres);
-
- auto minres = new SolverCreator<minres_solver_type>;
- Map::addCreator("minres", minres);
-
- auto gcr = new SolverCreator<gcr_type>;
- Map::addCreator("gcr", gcr);
-
-#ifdef HAVE_UMFPACK
- auto umfpack = new UmfpackSolverCreator;
- Map::addCreator("umfpack", umfpack);
- Map::addCreator("direct", umfpack);
-#endif
-
- Map::addCreator("default", gmres);
- }
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp b/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
deleted file mode 100644
index ef661b363706530a4b5672b1a9b2f808f91295ca..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
+++ /dev/null
@@ -1,208 +0,0 @@
-/** \file BlockMTLMatrix.hpp */
-
-#pragma once
-
-#include <array>
-
-#include <boost/numeric/mtl/matrices.hpp>
-
-#include <amdis/common/Utility.hpp>
-#include <amdis/common/Literals.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-
-namespace AMDiS
-{
- /// A wrapper for AMDiS::SolverMatrix to be used in MTL/ITL solvers
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- class BlockMTLMatrix
- : public std::array<std::array<MTLMatrix, _M>, _N>
- {
- using Self = BlockMTLMatrix;
-
- public:
- /// The index/size - type
- using size_type = typename MTLMatrix::size_type;
-
- /// The type of the elements of the MTLMatrix
- using value_type = typename MTLMatrix::value_type;
-
- /// The underlying mtl matrix type
- using BaseMatrix = MTLMatrix;
-
- public:
- /// Return the (R,C)'th matrix block
- template <std::size_t R, std::size_t C>
- auto& operator()(const index_t<R>, const index_t<C>)
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(*this));
- }
-
- /// Return the (R,C)'th matrix block
- template <std::size_t R, std::size_t C>
- auto const& operator()(const index_t<R>, const index_t<C>) const
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(*this));
- }
-
- /// Return the number of row blocks
- static constexpr std::size_t N() { return _N; }
-
- /// Return the number of column blocks
- static constexpr std::size_t M() { return _M; }
-
- /// perform blockwise multiplication A*b -> x
- template <class VectorIn, class VectorOut, class Assign>
- void mult(VectorIn const& b, VectorOut& x, Assign) const
- {
- // create iranges to access array blocks
- std::array<mtl::irange, _N> r_rows;
- std::array<mtl::irange, _M> r_cols;
- getRanges(r_rows, r_cols);
-
- forEach(range_<0, _N>, [&](const auto _i)
- {
- bool first = true;
-
- // a reference to the i'th block of x
- VectorOut x_i(x[r_rows[_i]]);
- forEach(range_<0, _M>, [&](const auto _j)
- {
- auto const& A_ij = this->operator()(_i, _j);
- if (num_rows(A_ij) > 0 && A_ij.nnz() > 0) {
- // a reference to the j'th block of b
- const VectorIn b_j(b[r_cols[_j]]);
-
- if (first) {
- Assign::first_update(x_i, A_ij * b_j);
- first = false;
- }
- else {
- Assign::update(x_i, A_ij * b_j);
- }
- }
- });
- });
- }
-
- /// A Multiplication operator returns a multiplication-expresssion.
- /// Calls \ref mult internally.
- template <class VectorIn>
- mtl::vec::mat_cvec_multiplier<Self, VectorIn>
- operator*(VectorIn const& v) const
- {
- return {*this, v};
- }
-
- /// Fill an array of irange corresponding to the row-sizes, used
- /// to access sub-vectors
- void getRowRanges(std::array<mtl::irange, _N>& r_rows) const
- {
- std::size_t start = 0;
- forEach(range_<0, _N>, [&](const auto _r) {
- std::size_t finish = start + num_rows((*this)(_r, 0_c));
- r_rows[_r].set(start, finish);
- start = finish;
- });
- }
-
- /// Fill an array of irange corresponding to the column-sizes, used
- /// to access sub-vectors
- void getColRanges(std::array<mtl::irange, _M>& r_cols) const
- {
- std::size_t start = 0;
- forEach(range_<0, _M>, [&](const auto _c) {
- std::size_t finish = start + num_cols((*this)(0_c, _c));
- r_cols[_c].set(start, finish);
- start = finish;
- });
- }
-
- /// Fill two arrays of irange corresponding to row and column sizes.
- /// \see getRowRanges() and \see getColRanges()
- void getRanges(std::array<mtl::irange, _N>& r_rows,
- std::array<mtl::irange, _M>& r_cols) const
- {
- getRowRanges(r_rows);
- getColRanges(r_cols);
- }
- };
-
-
- namespace Impl
- {
- /// Specialization of Impl::MTLMatrix from \file LinearAlgebraBase.hpp
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- struct BaseMatrix<BlockMTLMatrix<MTLMatrix, _N, _M>>
- {
- using type = MTLMatrix;
- };
- }
-
- /// Return the number of overall rows of a BlockMTLMatrix
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- inline std::size_t num_rows(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
- {
- std::size_t nRows = 0;
- forEach(range_<0, _N>, [&](const auto _r) {
- nRows += num_rows(A(_r, 0_c));
- });
- return nRows;
- }
-
- /// Return the number of overall columns of a BlockMTLMatrix
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- inline std::size_t num_cols(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
- {
- std::size_t nCols = 0;
- forEach(range_<0, _M>, [&](const auto _c) {
- nCols += num_cols(A(0_c, _c));
- });
- return nCols;
- }
-
- /// Return the size, i.e. rows*columns of a BlockMTLMatrix
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- inline std::size_t size(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
- {
- return num_rows(A) * num_cols(A);
- }
-
- /// Nullify a BlockMTLMatrix, i.e. nullify each block.
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- inline void set_to_zero(BlockMTLMatrix<MTLMatrix, _N, _M>& A)
- {
- forEach(range_<0, _N>, [&](const auto _r) {
- forEach(range_<0, _M>, [&](const auto _c) {
- set_to_zero(A(_r,_c));
- });
- });
- }
-
-} // end namespace AMDiS
-
-
-/// \cond HIDDEN_SYMBOLS
-namespace mtl
-{
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- struct Collection<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
- {
- using value_type = typename MTLMatrix::value_type;
- using size_type = typename MTLMatrix::size_type;
- };
-
- namespace ashape
- {
- template <class MTLMatrix, std::size_t _N, std::size_t _M>
- struct ashape_aux<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
- {
- using type = nonscal;
- };
-
- } // end namespace ashape
-
-} // end namespace mtl
-/// \endcond
diff --git a/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp b/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp
deleted file mode 100644
index 05d93d238c0798722d1682b17689493d009c1671..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#pragma once
-
-#include <array>
-#include <type_traits>
-
-#include <boost/numeric/mtl/utility/irange.hpp>
-
-#include <amdis/common/Mpl.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
-
-namespace AMDiS
-{
- /// A simple block-vector class
- template <class MTLVector, std::size_t N>
- class BlockMTLVector
- : public std::array<MTLVector, N>
- {
- public:
- /// The type of the blocks
- using BaseVector = MTLVector;
-
- /// The index/size - type
- using size_type = typename MTLVector::size_type;
-
- /// The type of the elements of the MTLVector
- using value_type = typename MTLVector::value_type;
- };
-
-
- namespace Impl
- {
- /// Specialization of Impl::BaseVector from \file LinearAlgebraBase.hpp
- template <class MTLVector, std::size_t N>
- struct BaseVector<BlockMTLVector<MTLVector, N>>
- {
- using type = MTLVector;
- };
- }
-
-
- /// Return the number of overall rows of a BlockMTLVector
- template <class MTLVector, std::size_t N>
- inline std::size_t num_rows(BlockMTLVector<MTLVector, N> const& vec)
- {
- std::size_t nRows = 0;
- forEach(range_<0, N>, [&](const auto _i) {
- nRows += num_rows(std::get<_i>(vec));
- });
- return nRows;
- }
-
- /// Return the number of overall columns of a BlockMTLVector
- template <class MTLVector, std::size_t N>
- inline std::size_t num_cols(BlockMTLVector<MTLVector, N> const& vec)
- {
- return 1;
- }
-
- /// Return the size, i.e. rows*columns of a BlockMTLVector
- template <class MTLVector, std::size_t N>
- inline std::size_t size(BlockMTLVector<MTLVector, N> const& vec)
- {
- return num_rows(vec);
- }
-
- /// Nullify a BlockMTLVector, i.e. nullify each block.
- template <class MTLVector, std::size_t N>
- inline void set_to_zero(BlockMTLVector<MTLVector, N>& vec)
- {
- forEach(range_<0, N>, [&](const auto _i) {
- set_to_zero(std::get<_i>(vec));
- });
- }
-
-
- /// A wrapper, that creates a contiguos vector corresponding to a block-vector
- /// and copy the value on construction and eventually back on destruction, if
- /// required.
- template <class BlockVector, class Vector, class NonConstBlockVector>
- class BlockVectorWrapper;
-
- // specialization for BlockMTLVector
- template <class BlockVector, class Vector, class MTLVector, std::size_t N>
- class BlockVectorWrapper<BlockVector, Vector, BlockMTLVector<MTLVector, N>>
- {
- static_assert( std::is_same< std::remove_const_t<BlockVector>, BlockMTLVector<MTLVector, N> >::value,
- "This specialization is for BlockMTLVectors only.");
-
- public:
- explicit BlockVectorWrapper(BlockVector& blockVector,
- bool copyBack = !std::is_const<BlockVector>::value)
- : blockVector(blockVector)
- , vector(num_rows(blockVector))
- , copyBack(copyBack)
- {
- assignTo();
- }
-
- ~BlockVectorWrapper()
- {
- if (copyBack)
- assignFrom(bool_<!std::is_const<BlockVector>::value>);
- }
-
- /// Return a reference to the block-vector
- BlockVector const& getBlockVector() const { return blockVector; }
-
- /// Return a reference to the contiguose-vector
- Vector& getVector() { return vector; }
- Vector const& getVector() const { return vector; }
-
- private:
- /// copy from block-vector to vector
- void assignTo()
- {
- std::size_t start = 0;
- for (std::size_t r = 0; r < N; ++r) {
- std::size_t finish = start + num_rows(blockVector[r]);
- mtl::irange range(start, finish);
-
- vector[range] = blockVector[r];
- start = finish;
- }
- }
-
- /// do not copy vector to block-vector since this is const
- template <bool Copy>
- std::enable_if_t<!Copy> assignFrom(bool_t<Copy>) { /* Do nothing */ }
-
- /// copy from vector to block-vector
- template <bool Copy>
- std::enable_if_t<Copy> assignFrom(bool_t<Copy>)
- {
- std::size_t start = 0;
- for (std::size_t r = 0; r < N; ++r) {
- std::size_t finish = start + num_rows(blockVector[r]);
- mtl::irange range(start, finish);
-
- blockVector[r] = vector[range];
- start = finish;
- }
- }
-
- private:
- BlockVector& blockVector;
- Vector vector;
-
- /// Copy data back to block-vector on destruction
- bool copyBack;
- };
-
-
- // Specialization for non-block-vectors
- template <class BlockVector, class Vector>
- class BlockVectorWrapper<BlockVector, Vector, Vector>
- {
- static_assert( std::is_same< std::remove_const_t<BlockVector>, Vector >::value,
- "This specialization is for contiguose vectors only.");
- public:
- explicit BlockVectorWrapper(BlockVector& vector, bool = false)
- : vector(vector)
- {}
-
- /// Return a reference to the vector
- BlockVector const& getBlockVector() const { return vector; }
- BlockVector& getVector() { return vector; }
- BlockVector const& getVector() const { return vector; }
-
- private:
- BlockVector& vector;
- };
-
- template <class BlockVector, class Vector>
- using BlockVectorWrapper_t
- = BlockVectorWrapper<BlockVector, Vector, std::remove_const_t<BlockVector>>;
-
-
- template <class BlockVector, class Vector = MTLDenseVector<typename BlockVector::value_type>>
- auto blockWrapper(BlockVector& bvec)
- {
- return BlockVectorWrapper_t<BlockVector, Vector>{bvec};
- }
-
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/CMakeLists.txt b/src/amdis/linear_algebra/mtl/CMakeLists.txt
index 0e6d6a6706c12d6bc4d862d2dfb09b1ff4649065..e295f73f7767c0b95987825dd993f89720bc30f5 100644
--- a/src/amdis/linear_algebra/mtl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/mtl/CMakeLists.txt
@@ -1,25 +1,10 @@
-#dune_library_add_sources(amdis SOURCES
-# SystemMatrix.cpp
-# SystemVector.cpp
-#)
-
install(FILES
- BITL_Preconditioner.hpp
- BITL_Solver.hpp
- BlockMTLMatrix.hpp
- BlockMTLVector.hpp
- Copy.hpp
DOFMatrix.hpp
DOFVector.hpp
ITL_Preconditioner.hpp
ITL_Solver.hpp
KrylovRunner.hpp
- LinearSolver.hpp
- Mapper.hpp
- MTLDenseVector.hpp
Preconditioner.hpp
- SystemMatrix.hpp
- SystemVector.hpp
UmfpackRunner.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/mtl)
diff --git a/src/amdis/linear_algebra/mtl/Copy.hpp b/src/amdis/linear_algebra/mtl/Copy.hpp
deleted file mode 100644
index 1724b311ac95756a312d1c5bf8d5f3b348e8e842..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/Copy.hpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#pragma once
-
-#include <algorithm>
-
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-
-#include <amdis/linear_algebra/mtl/Mapper.hpp>
-
-namespace AMDiS
-{
- // requires MatrixType to have an inserter
- template <class Value, std::size_t _N, std::size_t _M, class TargetType>
- void copy(BlockMTLMatrix<Value, _N, _M> const& source, TargetType& target)
- {
- using Mapper = std::conditional_t<(_N == _M), BlockMapper, RectangularMapper>;
- using value_type = typename TargetType::value_type;
- using BaseInserter = mtl::mat::inserter<TargetType, mtl::operations::update_plus<value_type>>;
- using MappedInserter = typename Mapper::template Inserter<BaseInserter>::type;
-
- std::size_t nnz = 0;
- for (std::size_t rb = 0; rb < _N; ++rb)
- for (std::size_t cb = 0; cb < _M; ++cb)
- nnz += source[rb][cb].nnz();
-
- {
- target.change_dim(num_rows(source), num_cols(source));
- set_to_zero(target);
-
- Mapper mapper(source);
- MappedInserter ins(target, mapper, int(1.2 * nnz / target.dim1()));
-
- for (std::size_t rb = 0; rb < _N; ++rb) {
- mapper.setRow(rb);
- for (std::size_t cb = 0; cb < _M; ++cb) {
- mapper.setCol(cb);
- ins << source[rb][cb];
- }
- }
- }
- }
-
- template <class FeSpaces, class Value, class TargetType>
- void copy(SystemMatrix<FeSpaces, Value> const& source, TargetType& target)
- {
- copy(source.getMatrix(), target);
- }
-
-} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index dba705f9dba0809107058fb1ea8e3e12dd71ccb3..5d3824c10d697d0d442f7714904e57f3c0b6a69c 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -1,7 +1,6 @@
#pragma once
#include <list>
-#include <map>
#include <memory>
#include <string>
#include <vector>
@@ -11,148 +10,74 @@
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>
-#include <dune/common/std/optional.hh>
-#include <dune/common/reservedvector.hh>
-#include <dune/functions/functionspacebases/flatmultiindex.hh>
-
#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
namespace AMDiS
{
/// \brief The basic container that stores a base matrix and a corresponding
/// row/column feSpace.
- template <class RowFeSpaceType,
- class ColFeSpaceType,
- class ValueType = double>
- class DOFMatrix
+ template <class ValueType>
+ class MtlMatrix
{
public:
- /// The type of the finite element space / basis of the row
- using RowFeSpace = RowFeSpaceType;
-
- /// The type of the finite element space / basis of the column
- using ColFeSpace = ColFeSpaceType;
-
/// The matrix type of the underlying base matrix
using BaseMatrix = mtl::compressed2D<ValueType>;
- /// The index/size - type
- using size_type = typename BaseMatrix::size_type;
-
/// The type of the elements of the DOFMatrix
using value_type = ValueType;
- /// The underlying field-type (typically the same as \ref value_type)
- using field_type = typename BaseMatrix::value_type;
+ /// The index/size - type
+ using size_type = typename BaseMatrix::size_type;
/// The type of the inserter used when filling the matrix. NOTE: need not be public
using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
public:
/// Constructor. Constructs new BaseMatrix.
- DOFMatrix(RowFeSpace const& rowFeSpace,
- ColFeSpace const& colFeSpace,
- std::string name)
- : rowFeSpace(rowFeSpace)
- , colFeSpace(colFeSpace)
- , name(name)
- , matrix(ClonablePtr<BaseMatrix>::make())
- {}
-
- /// Constructor. Takes pointer of data-matrix.
- DOFMatrix(RowFeSpace const& rowFeSpace,
- ColFeSpace const& colFeSpace,
- std::string name,
- BaseMatrix& matrix_ref)
- : rowFeSpace(rowFeSpace)
- , colFeSpace(colFeSpace)
- , name(name)
- , matrix(matrix_ref)
- {}
-
- /// Return the row-basis \ref rowFeSpace of the matrix
- RowFeSpace const& getRowFeSpace() const
- {
- return rowFeSpace;
- }
-
- /// Return the col-basis \ref colFeSpace of the matrix
- ColFeSpace const& getColFeSpace() const
- {
- return colFeSpace;
- }
+ MtlMatrix() = default;
/// Return a reference to the data-matrix \ref matrix
- BaseMatrix& getMatrix()
+ BaseMatrix& matrix()
{
- assert( !insertionMode );
- return *matrix;
+ assert( !inserter_ );
+ return matrix_;
}
/// Return a reference to the data-matrix \ref matrix
- BaseMatrix const& getMatrix() const
+ BaseMatrix const& matrix() const
{
- assert( !insertionMode );
- return *matrix;
+ assert( !inserter_ );
+ return matrix_;
}
- /// Return the size of the row \ref feSpace
- size_type N() const
- {
- return rowFeSpace.size();
- }
-
- /// Return the size of the column \ref feSpace
- size_type M() const
- {
- return colFeSpace.size();
- }
-
- /// Return the \ref name of this matrix
- std::string getName() const
- {
- return name;
- }
-
-
- using FlatIndex = Dune::Functions::FlatMultiIndex<size_type>;
- using BlockIndex = Dune::ReservedVector<size_type,1>;
/// \brief Returns an update-proxy of the inserter, to inserte/update a value at
/// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
/// be created by \ref init and must be closed by \ref finish after insertion.
- auto operator()(FlatIndex row, FlatIndex col)
+ void insert(size_type r, size_type c, value_type const& value)
{
- size_type r = row[0], c = col[0];
- test_exit_dbg( insertionMode, "Inserter not initilized!");
- test_exit_dbg( r < N() && c < M() ,
- "Indices out of range [0,", N(), ")x[0,", M(), ")" );
- return (*inserter)[r][c];
+ test_exit_dbg(inserter_, "Inserter not initilized!");
+ test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
+ "Indices out of range [0,", num_rows(matrix_), ")x[0,", num_cols(matrix_), ")" );
+ (*inserter_)[r][c] += value;
}
- auto operator()(BlockIndex row, BlockIndex col)
- {
- size_type r = row[0], c = col[0];
- test_exit_dbg( insertionMode, "Inserter not initilized!");
- test_exit_dbg( r < N() && c < M() ,
- "Indices out of range [0,", N(), ")x[0,", M(), ")" );
- return (*inserter)[r][c];
- }
/// Create inserter. Assumes that no inserter is currently active on this matrix.
- void init(bool prepareForInsertion)
+ template <class RowBasis, class ColBasis>
+ void init(RowBasis const& rowBasis, ColBasis const& colBasis,
+ bool prepareForInsertion)
{
- test_exit(!insertionMode && !inserter, "Matrix already in insertion mode!");
+ test_exit(!inserter_, "Matrix already in insertion mode!");
calculateSlotSize();
- matrix->change_dim(rowFeSpace.dimension(), colFeSpace.dimension());
- test_exit(num_rows(*matrix) == rowFeSpace.dimension() && num_cols(*matrix) == colFeSpace.dimension(),
- "Wrong dimensions in matrix");
+ matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
if (prepareForInsertion) {
- set_to_zero(*matrix);
- inserter = new Inserter(*matrix, slot_size);
- insertionMode = true;
+ set_to_zero(matrix_);
+ inserter_ = new Inserter(matrix_, slotSize_);
}
}
@@ -160,45 +85,27 @@ namespace AMDiS
/// final construction of the matrix.
void finish()
{
- delete inserter;
- inserter = NULL;
- insertionMode = false;
+ delete inserter_;
+ inserter_ = nullptr;
}
/// Return the number of nonzeros in the matrix
- std::size_t getNnz() const
+ std::size_t nnz() const
{
- return matrix->nnz();
+ return matrix_.nnz();
}
/// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
/// a one on the diagonal of the matrix.
auto applyDirichletBC(std::vector<char> const& dirichletNodes)
{
- std::list<Triplet<value_type>> columns;
-
// Define the property maps
- auto row = mtl::mat::row_map(*matrix);
- auto col = mtl::mat::col_map(*matrix);
- auto value = mtl::mat::value_map(*matrix);
+ auto row = mtl::mat::row_map(matrix_);
+ auto col = mtl::mat::col_map(matrix_);
+ auto value = mtl::mat::value_map(matrix_);
// iterate over the matrix
-#if 0
- for (auto r : mtl::major_of(*matrix)) { // rows or columns
- for (auto i : mtl::nz_of(r)) { // non-zeros within
- if (dirichletNodes[row(i)]) {
- // set identity row
- value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
- }
- else if (dirichletNodes[col(i)]) {
- columns.push_back({row(i), col(i), value(i)});
- value(i, value_type(0));
- }
- }
- }
-#endif
-
- for (auto r : mtl::rows_of(*matrix)) { // rows or columns
+ for (auto r : mtl::rows_of(matrix_)) { // rows of the matrix
if (dirichletNodes[r.value()]) {
for (auto i : mtl::nz_of(r)) { // non-zeros within
// set identity row
@@ -207,107 +114,36 @@ namespace AMDiS
}
}
- return columns;
- }
-
-#if 0
- void applyPeriodicBC(std::vector<char> const& periodicNodes,
- std::map<size_t, size_t> const& association, bool applyRow, bool applyCol)
- {
- bool apply = applyRow && applyCol;
-
- // Define the property maps
- auto row = mtl::mat::row_map(*matrix);
- auto col = mtl::mat::col_map(*matrix);
- auto value = mtl::mat::value_map(*matrix);
-
- std::vector<std::map<size_t, std::list<value_type>>> row_values(num_cols(*matrix));
- std::vector<std::map<size_t, std::list<value_type>>> col_values(num_rows(*matrix));
- std::vector<char> dualNodes(periodicNodes.size(), 0);
-
- // iterate over the matrix to collect the values and erase the source-row/col
- for (auto r : mtl::major_of(*matrix)) { // rows or columns
- for (auto i : mtl::nz_of(r)) { // non-zeros within
- if (applyRow && periodicNodes[row(i)]) {
- row_values[col(i)][association[row(i)]].push_back(value(i));
- value(i, value_type(0));
- dualNodes[association[row(i)]] = 1;
- } else if (applyCol && periodicNodes[col(i)]) {
- col_values[row(i)][association[col(i)]].push_back(value(i));
- value(i, value_type(0));
- }
- }
- }
-
- // TODO: use inserter for assignment of values!!!
- // iterate over the matrix to assign the values
- for (auto r : mtl::major_of(*matrix)) { // rows or columns
- for (auto i : mtl::nz_of(r)) { // non-zeros within
- if (applyRow && dualNodes[row(i)]) {
- value(i, std::accumulate(row_values[col(i)][row(i)].begin(),
- row_values[col(i)][row(i)].end(),
- value(i)) );
- }
- if (applyCol && dualNodes[col(i)]) {
- value(i, std::accumulate(col_values[row(i)][col(i)].begin(),
- col_values[row(i)][col(i)].end(),
- value(i)) );
- }
- }
- }
-
- // assign values 1, -1 to diagonal and associated entries
- if (apply) {
- Inserter ins(*matrix);
- for (auto const& a : association) {
- ins[a.first][a.first] << value_type( 1);
- ins[a.second][a.second] << value_type( 1);
- ins[a.first][a.second] << value_type(-1);
- ins[a.second][a.first] << value_type(-1);
- }
- }
+ return std::list<Triplet<value_type>>{};
}
-#endif
protected:
// Estimates the slot-size used in the inserter to reserve enough space per row.
void calculateSlotSize()
{
- slot_size = 0;
+ slotSize_ = 0;
- if (num_rows(*matrix) != 0)
- slot_size = int(double(matrix->nnz()) / num_rows(*matrix) * 1.2);
- if (slot_size < MIN_NNZ_PER_ROW)
- slot_size = MIN_NNZ_PER_ROW;
+ if (num_rows(matrix_) != 0)
+ slotSize_ = int(double(matrix_.nnz()) / num_rows(matrix_) * 1.2);
+ if (slotSize_ < MIN_NNZ_PER_ROW)
+ slotSize_ = MIN_NNZ_PER_ROW;
}
private:
/// The minimal number of nonzeros per row
static constexpr int MIN_NNZ_PER_ROW = 50;
- /// The finite element space / basis of the row
- RowFeSpace const& rowFeSpace;
-
- /// The finite element space / basis of the column
- ColFeSpace const& colFeSpace;
-
- /// The name of the DOFMatrix
- std::string name;
-
/// The data-matrix (can hold a new BaseMatrix or a pointer to external data
- ClonablePtr<BaseMatrix> matrix;
+ BaseMatrix matrix_;
/// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
- Inserter* inserter = NULL;
-
- /// A flag that indicates whether the matrix is in insertion mode
- bool insertionMode = false;
+ Inserter* inserter_ = nullptr;
/// The size of the slots used to insert values per row
- int slot_size = MIN_NNZ_PER_ROW;
-
- // friend class declarations
- template <class, class> friend class SystemMatrix;
+ int slotSize_ = MIN_NNZ_PER_ROW;
};
+ template <class RowBasisType, class ColBasisType, class ValueType = double>
+ using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, MtlMatrix<ValueType>>;
+
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp
index 6cfb293aac137b7c371a73642be289c0f64e54cb..01842741a5ffae416f7fe4866204bab847b3072a 100644
--- a/src/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -1,179 +1,87 @@
#pragma once
-#include <algorithm>
-#include <memory>
-#include <string>
-
-#include <dune/common/reservedvector.hh>
-#include <dune/functions/functionspacebases/flatmultiindex.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
-#include <amdis/common/ScalarTypes.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
-#include <amdis/utility/MultiIndex.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
namespace AMDiS
{
- /// The basic container that stores a base vector and a corresponding feSpace
- template <class Traits, class ValueType = double>
- class DOFVector
+ /// The basic container that stores a base vector and a corresponding basis
+ template <class ValueType>
+ class MtlVector
{
- using Self = DOFVector;
-
public:
- /// The type of the \ref feSpace
- using GlobalBasis = typename Traits::GlobalBasis;
-
/// The type of the base vector
- using BaseVector = MTLDenseVector<ValueType>;
+ using BaseVector = mtl::vec::dense_vector<ValueType>;
/// The index/size - type
- using size_type = typename GlobalBasis::size_type;
+ using size_type = typename BaseVector::size_type;
/// The type of the elements of the DOFVector
using value_type = ValueType;
- /// The underlying field-type (typically the same as \ref value_type)
- using field_type = typename BaseVector::value_type;
-
+ public:
/// Constructor. Constructs new BaseVector.
- DOFVector(GlobalBasis const& feSpace, std::string name)
- : feSpace(feSpace)
- , name(name)
- , vector(ClonablePtr<BaseVector>::make())
- {
- compress();
- }
-
- /// Constructor. Takes reference to existing BaseVector
- DOFVector(GlobalBasis const& feSpace, std::string name,
- BaseVector& vector_ref)
- : feSpace(feSpace)
- , name(name)
- , vector(vector_ref)
- {}
-
- /// Return the basis \ref feSpace of the vector
- GlobalBasis const& getFeSpace() const
- {
- return feSpace;
- }
-
- /// Return the data-vector \ref vector
- BaseVector const& getVector() const
- {
- return *vector;
- }
-
- /// Return the data-vector \ref vector
- BaseVector& getVector()
- {
- return *vector;
- }
+ MtlVector() = default;
- /// Return the size of the \ref feSpace
- size_type getSize() const
+ /// Return the data-vector \ref vector_
+ BaseVector const& vector() const
{
- return feSpace.dimension();
+ return vector_;
}
- /// Return the \ref name of this vector
- std::string getName() const
+ /// Return the data-vector \ref vector_
+ BaseVector& vector()
{
- return name;
+ return vector_;
}
- /// Resize the \ref vector to the size of the \ref feSpace.
- void compress()
+ /// Return the current size of the \ref vector_
+ size_type size() const
{
- if (num_rows(*vector) != getSize()) {
- resize(getSize());
- *vector = ValueType(0);
- }
+ return mtl::vec::size(vector_);
}
- template <class SizeInfo>
- void resize(SizeInfo s)
+ /// Resize the \ref vector_ to the size \p s
+ void resize(size_type s)
{
- vector->change_dim(size_type(s));
+ vector_.change_dim(s);
}
/// Access the entry \p i of the \ref vector with read-access.
- template <class Index>
- value_type const& operator[](Index idx) const
+ value_type const& operator[](size_type i) const
{
- auto const i = flatMultiIndex(idx);
- test_exit_dbg( i < num_rows(*vector) ,
- "Index ", i, " out of range [0,", num_rows(*vector), ")" );
- return (*vector)[i];
+ test_exit_dbg(i < mtl::vec::size(vector_),
+ "Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
+ return vector_[i];
}
/// Access the entry \p i of the \ref vector with write-access.
- template <class Index>
- value_type& operator[](Index idx)
+ value_type& operator[](size_type i)
{
- auto const i = flatMultiIndex(idx);
- test_exit_dbg( i < num_rows(*vector) ,
- "Index ", i, " out of range [0,", num_rows(*vector), ")" );
- return (*vector)[i];
- }
-
- /// \brief interpolate a function \p f to the basis \ref feSpace and store the
- /// coefficients in \ref vector.
- template <class F>
- void interpol(F&& f)
- {
- Dune::Functions::interpolate(feSpace, *this, std::forward<F>(f));
- }
-
- /// Scale each DOFVector by the factor \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator*=(Scalar s)
- {
- (*vector) *= s;
- return *this;
- }
-
- /// Sets each DOFVector to the scalar \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator=(Scalar s)
- {
- (*vector) = s;
- return *this;
- }
-
- /// Calls the copy assignment operator of the BaseVector \ref vector
- void copy(Self const& that)
- {
- *vector = that.getVector();
+ test_exit_dbg(i < mtl::vec::size(vector_),
+ "Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
+ return vector_[i];
}
private:
- /// The finite element space / basis associated with the data vector
- GlobalBasis const& feSpace;
-
- /// The name of the DOFVector (used in file-writing)
- std::string name;
-
/// The data-vector (can hold a new BaseVector or a pointer to external data
- ClonablePtr<BaseVector> vector;
-
- // friend class declarations
- template <class, class> friend class SystemVector;
+ BaseVector vector_;
};
- /// Constructor a dofvector from given feSpace and name
- template <class Traits, class ValueType = double>
- DOFVector<Traits, ValueType>
- makeDOFVector(typename Traits::GlobalBasis const& basis, std::string name)
+ template <class BasisType, class VectorType = double>
+ using DOFVector = DOFVectorBase<BasisType, MtlVector<VectorType>>;
+
+
+ /// Constructor a dofvector from given basis and name
+ template <class ValueType = double, class Basis>
+ DOFVector<Basis, ValueType>
+ makeDOFVector(Basis const& basis)
{
- return {basis, name};
+ return {basis};
}
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp b/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
index 2e8ac6ef19f5ed4dbd640d5628e66fb90fdcb9d7..42e25d04fb451bbd41363d9d318ed7b2aed2b4c9 100644
--- a/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
+++ b/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Preconditioner.h */
-
#pragma once
// MTL4 headers
diff --git a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
index c60b372cd5f30d18eb71b83b38cb37adc187c174..258e37500039cb67f5c59f617abb1f111c465aef 100644
--- a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
@@ -1,21 +1,9 @@
-/** \file ITL_Solver.h */
-
#pragma once
#include <string>
// MTL4 includes
#include <boost/numeric/itl/itl.hpp>
-// #include <boost/numeric/itl/krylov/bicg.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab_2.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab.hpp>
-// #include <boost/numeric/itl/krylov/cg.hpp>
-// #include <boost/numeric/itl/krylov/cgs.hpp>
-// #include <boost/numeric/itl/krylov/gmres.hpp>
-// #include <boost/numeric/itl/krylov/idr_s.hpp>
-// #include <boost/numeric/itl/krylov/qmr.hpp>
-// #include <boost/numeric/itl/krylov/tfqmr.hpp>
// more solvers defined in AMDiS
#include <amdis/linear_algebra/mtl/itl/minres.hpp>
@@ -28,7 +16,7 @@
#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
-#include <amdis/linear_algebra/mtl/LinearSolver.hpp>
+#include <amdis/linear_algebra/LinearSolver.hpp>
#include <amdis/linear_algebra/mtl/KrylovRunner.hpp>
#include <amdis/linear_algebra/mtl/UmfpackRunner.hpp>
@@ -405,13 +393,13 @@ namespace AMDiS
template <class ITLSolver>
using SolverCreator
- = typename LinearSolver<Matrix, Vector,
- KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
+ = typename LinearSolver<Matrix, Vector, Vector,
+ KrylovRunner<ITLSolver, Matrix, Vector>>::Creator;
#ifdef HAVE_UMFPACK
using UmfpackSolverCreator
- = typename LinearSolver<Matrix, Vector,
- UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
+ = typename LinearSolver<Matrix, Vector, Vector,
+ UmfpackRunner<Matrix, Vector>>::Creator;
#endif
using Map = CreatorMap<SolverBase>;
diff --git a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
index 45b9ea2a1475e777a9cf7a671b098cf4e34443ad..2927909c207baba7066ab66cba64476561ef54c2 100644
--- a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Runner.h */
-
#pragma once
#include <string>
@@ -14,7 +12,6 @@
#include <amdis/linear_algebra/SolverInfo.hpp>
#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
-#include <amdis/linear_algebra/mtl/BITL_Preconditioner.hpp>
namespace AMDiS
{
@@ -22,10 +19,11 @@ namespace AMDiS
*
* \brief
* Wrapper class for different MTL4 itl-solvers. These solvers
- * are parametrized by Matrix- and Vector.
+ * are parametrized by Matrix and Vector.
**/
template <class ITLSolver, class Matrix, class Vector>
- class KrylovRunner : public RunnerInterface<Matrix, Vector>
+ class KrylovRunner
+ : public RunnerInterface<Matrix, Vector>
{
using Super = RunnerInterface<Matrix, Vector>;
using PreconBase = typename Super::PreconBase;
@@ -33,71 +31,71 @@ namespace AMDiS
public:
/// Constructor.
explicit KrylovRunner(std::string prefix)
- : itlSolver(prefix)
+ : itlSolver_(prefix)
{
initPrecon(prefix);
- Parameters::get(prefix + "->max iteration", maxIter);
- Parameters::get(prefix + "->print cycle", printCycle);
+ Parameters::get(prefix + "->max iteration", maxIter_);
+ Parameters::get(prefix + "->print cycle", printCycle_);
}
- /// Implements \ref RunnerInterface::getLeftPrecon(), Returns \ref L
- virtual std::shared_ptr<PreconBase> getLeftPrecon() override
+ /// Implements \ref RunnerInterface::lLeftPrecon(), Returns \ref L_
+ virtual std::shared_ptr<PreconBase> leftPrecon() override
{
- return L;
+ return L_;
}
- /// Implements \ref RunnerInterface::getRightPrecon(), Returns \ref R
- virtual std::shared_ptr<PreconBase> getRightPrecon() override
+ /// Implements \ref RunnerInterface::rightPrecon(), Returns \ref R_
+ virtual std::shared_ptr<PreconBase> rightPrecon() override
{
- return R;
+ return R_;
}
- /// Set a new left preconditioner \ref L
+ /// Set a new left preconditioner \ref L_
void setLeftPrecon(std::shared_ptr<PreconBase> const& precon)
{
- L = precon;
+ L_ = precon;
}
- /// Set a new right preconditioner \ref R
+ /// Set a new right preconditioner \ref R_
void setRightPrecon(std::shared_ptr<PreconBase> const& precon)
{
- R = precon;
+ R_ = precon;
}
/// Implementation of \ref RunnerInterface::init()
virtual void init(Matrix const& A) override
{
- L->init(A);
- R->init(A);
+ L_->init(A);
+ R_->init(A);
}
/// Implementation of \ref RunnerInterface::exit()
virtual void exit() override
{
- L->exit();
- R->exit();
+ L_->exit();
+ R_->exit();
}
/// Implementation of \ref RunnerInterface::solve()
virtual int solve(Matrix const& A, Vector& x, Vector const& b,
SolverInfo& solverInfo) override
{
- test_exit(bool(L), "There is no left preconditioner");
- test_exit(bool(R), "There is no right preconditioner");
+ test_exit(bool(L_), "There is no left preconditioner");
+ test_exit(bool(R_), "There is no right preconditioner");
auto r = Vector(b);
if (two_norm(x) != 0)
r -= A * x;
// print information about the solution process
- itl::cyclic_iteration<double> iter(r, maxIter, solverInfo.getRelTolerance(),
- solverInfo.getAbsTolerance(),
- printCycle);
+ itl::cyclic_iteration<double> iter(r, maxIter_, solverInfo.getRelTolerance(),
+ solverInfo.getAbsTolerance(),
+ printCycle_);
iter.set_quite(solverInfo.getInfo() == 0);
iter.suppress_resume(solverInfo.getInfo() == 0);
- int error = itlSolver(A, x, b, *L, *R, iter);
+ int error = itlSolver_(A, x, b, *L_, *R_, iter);
solverInfo.setAbsResidual(iter.resid());
solverInfo.setRelResidual(iter.relresid());
@@ -109,32 +107,31 @@ namespace AMDiS
void initPrecon(std::string prefix)
{
// Creator for the left preconditioner
- std::string leftPrecon = "", rightPrecon = "";
- Parameters::get(prefix + "->left precon", leftPrecon);
- Parameters::get(prefix + "->right precon", rightPrecon);
+ std::string leftPreconName = "", rightPreconName = "";
+ Parameters::get(prefix + "->left precon", leftPreconName);
+ Parameters::get(prefix + "->right precon", rightPreconName);
auto leftCreator
- = CreatorMap<PreconBase>::getCreator(leftPrecon, prefix + "->left precon");
+ = CreatorMap<PreconBase>::getCreator(leftPreconName, prefix + "->left precon");
auto rightCreator
- = CreatorMap<PreconBase>::getCreator(rightPrecon, prefix + "->right precon");
+ = CreatorMap<PreconBase>::getCreator(rightPreconName, prefix + "->right precon");
- L = leftCreator->create();
- R = rightCreator->create();
+ L_ = leftCreator->create();
+ R_ = rightCreator->create();
}
private:
/// The itl solver object that performes the actual solve
- ITLSolver itlSolver;
+ ITLSolver itlSolver_;
/// Pointer to the left and right preconditioner
- std::shared_ptr<PreconBase> L, R;
+ std::shared_ptr<PreconBase> L_, R_;
/// The maximal number of iterations
- std::size_t maxIter = 1000;
+ std::size_t maxIter_ = 1000;
/// Print solver residuals every \ref printCycle 's iteration
- std::size_t printCycle = 100;
-
+ std::size_t printCycle_ = 100;
};
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp b/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
deleted file mode 100644
index 4e08e4698dca62edaafda5c74647266decfc4911..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
+++ /dev/null
@@ -1,50 +0,0 @@
-#pragma once
-
-
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-
-namespace AMDiS
-{
- template <class Value>
- using MTLDenseVector = mtl::vec::dense_vector<Value>;
-
- namespace Impl
- {
- template <class Vector>
- class MTLWrapper
- : public BaseWrapper<Vector>
- {
- using Super = BaseWrapper<Vector>;
-
- public:
- using value_type = typename Super::value_type;
- using size_type = typename Super::size_type;
-
- template <class Vector_>
- MTLWrapper(Vector_&& vec)
- : Super(std::forward<Vector_>(vec))
- {}
-
- void resize(size_type s)
- {
- this->vec.change_dim(s);
- }
-
- size_type size() const
- {
- return num_rows(this->vec);
- }
- };
-
- } // end namespace Impl
-
- /// wrap the mtl dense_vector class to provide resize and size member-functions
- template <class Value>
- Impl::MTLWrapper<MTLDenseVector<Value>> wrapper(MTLDenseVector<Value>& vec)
- {
- return {vec};
- }
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/Mapper.hpp b/src/amdis/linear_algebra/mtl/Mapper.hpp
deleted file mode 100644
index 6249f0949554a1e13a49edcac61b74855e868df0..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/Mapper.hpp
+++ /dev/null
@@ -1,241 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include <boost/numeric/mtl/matrix/mapped_inserter.hpp>
-
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-
-namespace AMDiS
-{
-
- /**
- * \brief BaseClass for all mapper types
- *
- * A mapper assigned a local row/column index of a cell of
- * a block matrix to a global matrix index.
- * Call \ref setRow and \ref setCol first, to select a block
- * in the block matrix. Then get with \ref row, respective \ref col,
- * the global matrix index to the assigned local index.
- **/
- template <class Derived>
- struct MapperBase
- {
- using size_type = std::size_t;
-
- template <class BaseInserter>
- struct Inserter {
- using type = mtl::mat::mapped_inserter<BaseInserter, Derived>;
- };
-
- /// set the current block row
- void setRow(size_type rb) { self().setRow(rb); }
-
- /// set the current block columns
- void setCol(size_type cb) { self().setCol(cb); }
-
- /// return global matrix row, for local row in the current matrix block
- size_type row(size_type r) const { return self().row(r); }
-
- /// return global matrix column, for local column in the current matrix block
- size_type col(size_type c) const { return self().col(c); }
-
- /// return overall number of rows
- size_type getNumRows() const { return self().getNumRows(); }
-
- /// return overall number of columns
- size_type getNumCols() const { return self().getNumCols(); }
-
- size_type getNumRows(size_type comp) const { return self().getNumRows(comp); }
- size_type getNumCols(size_type comp) const { return self().getNumCols(comp); }
-
- /// return number of components/blocks
- size_type getNumComponents() const { return self().getNumComponents(); }
-
- size_type row(size_type rb, size_type cb, size_type r)
- {
- setRow(rb);
- setCol(cb);
- return row(r);
- }
-
- size_type col(size_type rb, size_type cb, size_type c)
- {
- setRow(rb);
- setCol(cb);
- return col(c);
- }
-
- Derived& self() { return static_cast<Derived&>(*this); }
- Derived const& self() const { return static_cast<Derived const&>(*this); }
- };
-
-
- /// Mapper implementation for non-parallel block matrices
- struct BlockMapper : public MapperBase<BlockMapper>
- {
- /// Default constructor
- BlockMapper() = default;
- BlockMapper(BlockMapper const& other) = default;
- BlockMapper& operator=(BlockMapper const& other) = default;
-
- /// Constructor for BlockMTLMatrix
- template <class Value, std::size_t _N>
- explicit BlockMapper(BlockMTLMatrix<Value, _N, _N> const& mat)
- : nComp(_N)
- , rowOffset(0), colOffset(0)
- , nrow(num_rows(mat)), ncol(nrow)
- , sizes(nComp)
- {
- for (std::size_t i = 0; i < _N; ++i)
- sizes[i] = num_rows(mat[i][i]);
- }
-
- /// Constructor for SystemMatrix
- template <class FeSpaces, class Value>
- explicit BlockMapper(SystemMatrix<FeSpaces, Value> const& mat)
- : BlockMapper(mat.getMatrix())
- {}
-
- /// Constructor for single DOFMatrix
- template <class FeSpace, class Value>
- explicit BlockMapper(DOFMatrix<FeSpace, Value> const& mat)
- : nComp(1)
- , rowOffset(0), colOffset(0)
- , nrow(0), ncol(0)
- , sizes(nComp)
- {
- sizes[0] = mat.N();
- nrow += sizes[0];
- ncol = nrow;
- }
-
- /// Constructor for system with equal components
- BlockMapper(size_type nComp, size_type nDOFperComp)
- : nComp(nComp),
- rowOffset(0), colOffset(0),
- nrow(nComp*nDOFperComp), ncol(nrow),
- sizes(nComp)
- {
- for (size_type i = 0; i < nComp; ++i)
- sizes[i] = nDOFperComp;
- }
-
- /// calculate row offset for row component \param r
- void setRow(size_type rb)
- {
- test_exit_dbg(rb <= sizes.size(), "row nr out of range!");
- rowOffset = sum(rb);
- }
-
- /// calculate column offset for col component \param c
- void setCol(size_type cb)
- {
- test_exit_dbg(cb <= sizes.size(), "column nr out of range!");
- colOffset = sum(cb);
- }
-
- size_type row(size_type r) const { return r + rowOffset; }
- size_type col(size_type c) const { return c + colOffset; }
-
- size_type getNumRows() const { return nrow; }
- size_type getNumCols() const { return ncol; }
-
- size_type getNumRows(size_type comp) const { return sizes[comp]; }
- size_type getNumCols(size_type comp) const { return sizes[comp]; }
-
- size_type getNumComponents() const { return nComp; }
-
- private: // methods
-
- ///compute the sum of sizes from [0, end)
- size_type sum(size_type end) const
- {
- size_type ret = 0;
- for (size_type i = 0; i < end; ++i)
- ret += sizes[i];
- return ret;
- }
-
- private: // variables
-
- size_type nComp = 0;
- size_type rowOffset = 0;
- size_type colOffset = 0;
- size_type nrow = 0;
- size_type ncol = 0;
-
- std::vector<size_type> sizes;
- };
-
-
- /// Mapper implementation for non-parallel rectangular block matrices
- struct RectangularMapper : public MapperBase<RectangularMapper>
- {
- /// Constructor for block-matrices
- template <class Value, std::size_t _N, std::size_t _M>
- explicit RectangularMapper(BlockMTLMatrix<Value, _N, _M> const& mat)
- : nRowComp(_N), nColComp(_M)
- , rowOffset(0), colOffset(0)
- , nrow(num_rows(mat)), ncol(num_cols(mat))
- , sizes_rows(nRowComp), sizes_cols(nColComp)
- {
- for (size_type i = 0; i < _N; ++i)
- sizes_rows[i] = num_rows(mat[i][0]);
- for (size_type j = 0; j < _M; ++j)
- sizes_cols[j] = num_cols(mat[0][j]);
- }
-
- /// calculate row offset for row component \param r
- void setRow(size_type rb)
- {
- test_exit_dbg(rb <= sizes_rows.size(), "row nr out of range!");
- rowOffset = sum(rb, sizes_rows);
- }
-
- /// calculate column offset for col component \param c
- void setCol(size_type cb)
- {
- test_exit_dbg(cb <= sizes_cols.size(), "column nr out of range!");
- colOffset = sum(cb, sizes_cols);
- }
-
- size_type row(size_type r) const { return r + rowOffset; }
- size_type col(size_type c) const { return c + colOffset; }
-
- size_type getNumRows() const { return nrow; }
- size_type getNumCols() const { return ncol; }
-
- size_type getNumRows(size_type comp) const { return sizes_rows[comp]; }
- size_type getNumCols(size_type comp) const { return sizes_cols[comp]; }
-
- size_type getNumComponents() const { return nRowComp; }
- size_type getNumRowComponents() const { return nRowComp; }
- size_type getNumColComponents() const { return nColComp; }
-
- private: // methods
-
- ///compute the sum of sizes from [0, end)
- size_type sum(size_type end, std::vector<size_type> const& sizes) const
- {
- size_type ret = 0;
- for (size_type i = 0; i < end; ++i)
- ret += sizes[i];
- return ret;
- }
-
- private: // variables
-
- size_type nRowComp = 0;
- size_type nColComp = 0;
- size_type rowOffset = 0;
- size_type colOffset = 0;
- size_type nrow = 0;
- size_type ncol = 0;
-
- std::vector<size_type> sizes_rows;
- std::vector<size_type> sizes_cols;
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/Preconditioner.hpp b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
index 822d1221841e73676e74b282f74571fa1460fc3b..cc610a82cacf5a498ca9c67bcf55de638e4296bc 100644
--- a/src/amdis/linear_algebra/mtl/Preconditioner.hpp
+++ b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Preconditioner.h */
-
#pragma once
#include <amdis/linear_algebra/PreconditionerInterface.hpp>
@@ -12,7 +10,8 @@ namespace AMDiS
* \brief Wrapper for using ITL preconditioners in AMDiS.
*/
template <class Matrix, class Vector, class PreconRunner>
- class Preconditioner : public PreconditionerInterface<Matrix, Vector>
+ class Preconditioner
+ : public PreconditionerInterface<Matrix, Vector>
{
using Self = Preconditioner;
using Super = PreconditionerInterface<Matrix, Vector>;
@@ -31,31 +30,31 @@ namespace AMDiS
/// Implementation of \ref PreconditionerBase::init()
virtual void init(Matrix const& fullMatrix) override
{
- precon.reset(new PreconRunner(fullMatrix));
+ precon_.reset(new PreconRunner(fullMatrix));
}
/// Implementation of \ref PreconditionerInterface::exit()
virtual void exit() override
{
- precon.reset();
+ precon_.reset();
}
/// Implementation of \ref PreconditionerBase::solve()
virtual void solve(Vector const& vin, Vector& vout) const override
{
- test_exit_dbg(bool(precon), "No preconditioner initialized!");
- precon->solve(vin, vout);
+ test_exit_dbg(bool(precon_), "No preconditioner initialized!");
+ precon_->solve(vin, vout);
}
/// Implementation of \ref PreconditionerBase::adjoint_solve()
virtual void adjoint_solve(Vector const& vin, Vector& vout) const override
{
- test_exit_dbg(bool(precon), "No preconditioner initialized!");
- precon->adjoint_solve(vin, vout);
+ test_exit_dbg(bool(precon_), "No preconditioner initialized!");
+ precon_->adjoint_solve(vin, vout);
}
private:
- std::shared_ptr<PreconRunner> precon;
+ std::shared_ptr<PreconRunner> precon_;
};
} // namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemMatrix.cpp b/src/amdis/linear_algebra/mtl/SystemMatrix.cpp
deleted file mode 100644
index 6182e50796fac5ccbea79435d9e3c444df2f91f3..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemMatrix.cpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMMATRIX
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMMATRIX
-
-
-namespace AMDiS
-{
- // explicit template instatiation
- // template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemMatrix.hpp b/src/amdis/linear_algebra/mtl/SystemMatrix.hpp
deleted file mode 100644
index a674fd16a44fc8afa275cfd2a00faa4ebaeea3de..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemMatrix.hpp
+++ /dev/null
@@ -1,207 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-#include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
- /// \cond HIDDEN_SYMBOLS
- namespace Impl
- {
- // DOFMatrices = std::tuple< std::tuple< DOFMatrix<RowFeSpace, ColFeSpace, ValueType>... >... >
- template <class DOFMatrices>
- struct BuildDOFMatrices
- {
- template <std::size_t R>
- using DOFMatrixRow = std::tuple_element_t<R, DOFMatrices>;
-
- template <std::size_t R, std::size_t C>
- using DOFMatrixRowCol = std::tuple_element_t<C, DOFMatrixRow<R>>;
-
- // add arg to repeated constructor argument list
- template <std::size_t... R, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
- static DOFMatrices make(Indices<R...>,
- RowFeSpaces&& rowFeSpaces,
- ColFeSpaces&& colFeSpace,
- MultiMatrix&& multiMatrix)
- {
- return DOFMatrices{make_row(
- index_<R>,
- MakeSeq_t<std::tuple_size<std::decay_t<ColFeSpaces>>::value>(),
- std::get<R>(std::forward<RowFeSpaces>(rowFeSpaces)),
- std::forward<ColFeSpaces>(colFeSpace),
- std::forward<MultiMatrix>(multiMatrix)
- )...};
- }
-
- template <std::size_t R, std::size_t... C, class RowFeSpace, class ColFeSpaces, class MultiMatrix>
- static DOFMatrixRow<R> make_row(index_t<R>,
- Indices<C...>,
- RowFeSpace&& rowFeSpace,
- ColFeSpaces&& colFeSpace,
- MultiMatrix&& multiMatrix)
- {
- return DOFMatrixRow<R>{DOFMatrixRowCol<R,C>(
- std::forward<RowFeSpace>(rowFeSpace),
- std::get<C>(std::forward<ColFeSpaces>(colFeSpace)),
- "matrix[" + std::to_string(R) + "][" + std::to_string(C) + "]",
- multiMatrix(index_<R>, index_<C>)
- )...};
- }
- };
-
- // construction method to construct a tuple of tuples of DOFMatrices
- template <class DOFMatrices, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
- DOFMatrices buildDOFMatrices(RowFeSpaces&& rowFeSpaces,
- ColFeSpaces&& colFeSpace,
- MultiMatrix&& multiMatrix)
- {
- return BuildDOFMatrices<DOFMatrices>::make(
- MakeSeq_t<std::tuple_size<std::decay_t<RowFeSpaces>>::value>(), // traverse rows first
- std::forward<RowFeSpaces>(rowFeSpaces),
- std::forward<ColFeSpaces>(colFeSpace),
- std::forward<MultiMatrix>(multiMatrix));
- }
-
- } // end namespace Impl
- /// \endcond
-
-
- /// \brief Container that repesents multiple data-matrices.
- /**
- * Represents a fixed-size matrix of \ref mtl::compressed2D + a tuple of corresponding
- * feSpaces. This (R,C)'th matrix combined with the (R,C)'th FeSpace
- * builds a \ref DOFMatrix and can be return by the \ref operator().
- **/
- template <class FeSpaces, class ValueType = double>
- class SystemMatrix
- {
- template <class RowFeSpace>
- struct MultiMatrixRowGenerator
- {
- template <class ColFeSpace>
- using DOFMatrixGenerator = DOFMatrix<RowFeSpace, ColFeSpace, ValueType>;
- using DOFMatrices = MakeTuple_t<DOFMatrixGenerator, FeSpaces>;
- };
-
- template <class RowFeSpace>
- using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace>::DOFMatrices;
-
- public:
- /// The number of blocks per row/column
- static constexpr std::size_t nComponent = std::tuple_size<FeSpaces>::value;
- AMDIS_STATIC_ASSERT( nComponent > 0 );
-
- /// The type of the matrix in each block
- using BaseMatrix = mtl::compressed2D<ValueType>;
-
- /// The type of the matrix of DOFMatrices
- using DOFMatrices = MakeTuple_t<DOFMatrixGenerator, FeSpaces>;
-
- /// The type of the matrix of base-matrices
- using MultiMatrix = BlockMTLMatrix<BaseMatrix, nComponent, nComponent>;
-
-
- /// Constructor that constructs new base-matrices
- explicit SystemMatrix(FeSpaces const& feSpaces)
- : feSpaces(feSpaces)
- , matrix()
- , dofmatrices(Impl::buildDOFMatrices<DOFMatrices>(feSpaces, feSpaces, matrix))
- {}
-
-
- /// Return the number of blocks per row
- static constexpr std::size_t N()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- /// Return the number of blocks per column
- static constexpr std::size_t M()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- /// Return the (R,C)'th underlaying base matrix
- template <std::size_t R, std::size_t C>
- auto& getMatrix(const index_t<R> _r = {}, const index_t<C> _c = {})
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return matrix(_r, _c);
- }
-
- /// Return the (R,C)'th underlaying base matrix
- template <std::size_t R, std::size_t C>
- auto const& getMatrix(const index_t<R> _r = {}, const index_t<C> _c = {}) const
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return matrix(_r, _c);
- }
-
- /// Return the underlying multi matrix
- MultiMatrix& getMatrix() { return matrix; }
- MultiMatrix const& getMatrix() const { return matrix; }
-
- /// Return the I'th row finite element space
- template <std::size_t I = 0>
- auto& getRowFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < N(), "Index out of range [0,N)");
- return std::get<I>(feSpaces);
- }
-
- /// Return the I'th column finite element space
- template <std::size_t I = 0>
- auto const& getColFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < M(), "Index out of range [0,M)");
- return std::get<I>(feSpaces);
- }
-
- /// Return the (R,C)'th DOFMatrix
- template <std::size_t R = 0, std::size_t C = 0>
- auto& operator()(const index_t<R> = {}, const index_t<C> = {})
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(dofmatrices));
- }
-
- /// Return the (R,C)'th DOFMatrix
- template <std::size_t R = 0, std::size_t C = 0>
- auto const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
- {
- static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
- return std::get<C>(std::get<R>(dofmatrices));
- }
-
- private:
- /// a tuple of feSpaces for rows and columns of the matrix
- FeSpaces const& feSpaces;
-
- /// a tuple of data-matrices
- MultiMatrix matrix;
-
- /// a tuple of DOFMatrices
- DOFMatrices dofmatrices;
- };
-
-// #ifndef AMDIS_NO_EXTERN_SYSTEMMATRIX
-// // explicit instantiation in SystemMatrix.cpp
-// extern template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-// #endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemVector.cpp b/src/amdis/linear_algebra/mtl/SystemVector.cpp
deleted file mode 100644
index 6fe23e982458095ec98ca59a2564cafb541cc638..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemVector.cpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMVECTOR
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMVECTOR
-
-
-namespace AMDiS
-{
- // explicit template instatiation
- // template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemVector.hpp b/src/amdis/linear_algebra/mtl/SystemVector.hpp
deleted file mode 100644
index 93123b79b26a328fcd345c29a51fab599c31e567..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemVector.hpp
+++ /dev/null
@@ -1,263 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <memory>
-#include <string>
-#include <tuple>
-#include <vector>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/ScalarTypes.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
-#include <amdis/linear_algebra/mtl/DOFVector.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
- namespace Impl
- {
- // DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
- template <class DOFVectors>
- struct BuildDOFVectors
- {
- template <std::size_t I>
- using DOFVectorIdx = std::tuple_element_t<I, DOFVectors>;
-
- template <std::size_t... I, class FeSpaces, class MultiVector>
- static DOFVectors make(Indices<I...>,
- FeSpaces&& feSpaces,
- std::vector<std::string> const& names,
- MultiVector&& multiVector)
- {
- return DOFVectors{DOFVectorIdx<I>(
- std::get<I>(std::forward<FeSpaces>(feSpaces)),
- names[I],
- multiVector[I]
- )...};
- }
- };
-
- // construction method to construct a tuple of DOFVectors
- template <class DOFVectors, class FeSpaces, class MultiVector>
- DOFVectors buildDOFVectors(FeSpaces&& feSpaces,
- std::vector<std::string> const& names,
- MultiVector&& multiVector)
- {
- return BuildDOFVectors<DOFVectors>::make(
- MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
- std::forward<FeSpaces>(feSpaces),
- names,
- std::forward<MultiVector>(multiVector));
- }
-
- } // end namespace Impl
-
-
- /// \brief Container that repesents multiple data-Vectors of different value types.
- /**
- * Represents a std::array of \ref mtl::dense_vector + a tuple of corresponding
- * feSpaces. This I'th matrix combined with the I'th FeSpace
- * builds a \ref DOFVector and can be return by the \ref operator().
- **/
- template <class FeSpaces, class ValueType = double>
- class SystemVector
- {
- using Self = SystemVector;
-
- template <class FeSpace>
- using DOFVectorGenerator = DOFVector<FeSpace, ValueType>;
-
- public:
- /// The number of blocks
- static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
- AMDIS_STATIC_ASSERT( nComponents > 0 );
-
- /// The type of the block-vector
- using BaseVector = typename DOFVectorGenerator<std::tuple_element_t<0, FeSpaces>>::BaseVector;
- using MultiVector = BlockMTLVector<BaseVector, nComponents>;
-
- /// The type of the vector of DOFVectors
- using DOFVectors = MakeTuple_t<DOFVectorGenerator, FeSpaces>;
-
- /// Constructor that constructs new base-vectors
- SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
- : feSpaces(feSpaces)
- , names(names)
- , vector{}
- , dofvectors(Impl::buildDOFVectors<DOFVectors>(feSpaces, names, vector))
- {
- compress();
- }
-
-
- /// Return the number of blocks
- static constexpr std::size_t size()
- {
- return std::tuple_size<FeSpaces>::value;
- }
-
- /// Return the number of blocks
- int count() const
- {
- return int(size());
- }
-
- /// Return a shared pointer to the I'th underlaying base vector
- template <std::size_t I>
- auto& getVector(const index_t<I> = {})
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return vector[I];
- }
-
- /// Return a shared pointer to the I'th underlaying base vector
- template <std::size_t I>
- auto const& getVector(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return vector[I];
- }
-
- /// Return the underlying multi vector
- MultiVector& getVector() { return vector; }
- MultiVector const& getVector() const { return vector; }
-
- /// Return the I'th finite element space
- template <std::size_t I = 0>
- auto const& getFeSpace(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(feSpaces);
- }
-
- /// Return the name of the i'th DOFVector
- std::string getName(std::size_t i) const
- {
- return names[i];
- }
-
- /// Return the I'th DOFVector
- template <std::size_t I>
- auto& operator[](const index_t<I>)
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- /// Return the I'th DOFVector
- template <std::size_t I>
- auto const& operator[](const index_t<I>) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- /// Return the I'th DOFVector
- template <std::size_t I>
- auto& getDOFVector(const index_t<I> = {})
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- /// Return the I'th DOFVector
- template <std::size_t I>
- auto const& getDOFVector(const index_t<I> = {}) const
- {
- static_assert(I < size(), "Index out of range [0,SIZE)" );
- return std::get<I>(dofvectors);
- }
-
- /// Resize the I'th \ref vector to the sizes of the I'th \ref feSpaces.
- template <std::size_t I = 0>
- void compress(const index_t<I> = {})
- {
- std::get<I>(dofvectors).compress();
- }
-
- /// Resize the \ref vectors to the sizes of the \ref feSpaces.
- void compress()
- {
- forEach(range_<0, size()>, [this](const auto I) {
- std::get<I>(dofvectors).compress();
- });
- }
-
- /// Copy the Systemvector \p that element-wise
- void copy(Self const& that)
- {
- forEach(range_<0, size()>, [this, &that](const auto I) {
- std::get<I>(dofvectors).copy(that[I]);
- });
- }
-
- /// Scale each DOFVector by the factor \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator*=(Scalar s)
- {
- forEach(range_<0, size()>, [this, s](const auto I) {
- vector[I] *= s;
- });
- return *this;
- }
-
- /// Sets each DOFVector to the scalar \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator=(Scalar s)
- {
- forEach(range_<0, size()>, [this, s](const auto I) {
- vector[I] = s;
- });
- return *this;
- }
-
- private:
- /// a tuple of feSpaces
- FeSpaces const& feSpaces;
-
- /// The names of the DOFVectors
- std::vector<std::string> names;
-
- /// a tuple of base vectors, i.e. mtl::dense_vectors
- MultiVector vector;
-
- /// a tuple of dofvectors referencing \ref vector
- DOFVectors dofvectors;
- };
-
-
- /// Construct a systemvector from given feSpace-tuple and names-vector
- template <class FeSpacesType, class ValueType = double>
- SystemVector<FeSpacesType, ValueType>
- makeSystemVector(FeSpacesType const& feSpaces, std::vector<std::string> names)
- {
- return {feSpaces, names};
- }
-
- /// Construct a systemvector from given feSpace-tuple. Create vector of names
- /// from base_name, i.e. {base_name_0, base_name_1, ...}
- template <class FeSpaces, class ValueType = double>
- SystemVector<FeSpaces, ValueType>
- makeSystemVector(FeSpaces const& feSpaces, std::string base_name)
- {
- std::vector<std::string> names;
- for (std::size_t i = 0; i < std::tuple_size<FeSpaces>::value; ++i)
- names.push_back(base_name + "_" + std::to_string(i));
- return {feSpaces, names};
- }
-
-
-// #ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
-// // explicit instantiation in SystemVector.cpp
-// extern template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-// #endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp b/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
index fbca866f1b4cef556ec22cb702df3c28f21f8eb8..26b0cb0805edca1fa5c349f0268be75264b56a38 100644
--- a/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
@@ -1,23 +1,17 @@
#pragma once
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
#ifdef HAVE_UMFPACK
#include <algorithm>
-#include <iostream>
#include <string>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
+#include <amdis/Output.hpp>
#include <amdis/linear_algebra/RunnerInterface.hpp>
#include <amdis/linear_algebra/SolverInfo.hpp>
-#include <amdis/linear_algebra/mtl/Copy.hpp>
-
namespace AMDiS
{
/**
@@ -40,16 +34,16 @@ namespace AMDiS
using Super = RunnerInterface<Matrix, Vector>;
using PreconBase = typename Super::PreconBase;
- using FullMatrix = BaseMatrix_t<Matrix>;
+ using FullMatrix = Matrix;
using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
public:
/// Constructor. Reads UMFPACK parameters from initfile
UmfpackRunnerBase(std::string prefix)
{
- Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
- Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
- Parameters::get(prefix + "->alloc init", alloc_init);
+ Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ?
+ Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
+ Parameters::get(prefix + "->alloc init", allocInit_);
}
/// Destructor
@@ -60,11 +54,11 @@ namespace AMDiS
SolverInfo& solverInfo) override
{
AMDIS_FUNCNAME("Umfpack_Runner::solve()");
- test_exit(bool(solver), "The umfpack solver was not initialized\n");
+ test_exit(bool(solver_), "The umfpack solver was not initialized\n");
int code = 0;
try {
- code = (*solver)(x, b);
+ code = (*solver_)(x, b);
} catch (mtl::mat::umfpack::error& e) {
error_exit("UMFPACK_ERROR(solve, ", e.code, ") = ", e.what());
}
@@ -79,47 +73,17 @@ namespace AMDiS
return code;
}
- /// Implementation of \ref RunnerInterface::adjoint_solve()
+ /// Implementation of \ref RunnerInterface::exit()
virtual void exit() override {}
protected:
- std::shared_ptr<SolverType> solver;
+ std::shared_ptr<SolverType> solver_;
- int store_symbolic = 0;
- int symmetric_strategy = 0;
- double alloc_init = 0.7;
+ int storeSymbolic_ = 0;
+ int symmetricStrategy_ = 0;
+ double allocInit_ = 0.7;
};
- // specialization for block-matrices
- template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
- class UmfpackRunner<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
- : public UmfpackRunnerBase<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
- {
- using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
- using Super = UmfpackRunnerBase<Matrix, Vector>;
-
- using SolverType = typename Super::SolverType;
-
- public:
- UmfpackRunner(std::string prefix)
- : Super(prefix)
- {}
-
- /// Implementation of \ref RunnerInterface::init()
- virtual void init(Matrix const& A) override
- {
- copy(A, fullMatrix);
-
- try {
- Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
- } catch (mtl::mat::umfpack::error const& e) {
- error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
- }
- }
-
- private:
- typename Super::FullMatrix fullMatrix;
- };
// specialization for full-matrices
template <class T, class Param, class Vector>
@@ -140,7 +104,7 @@ namespace AMDiS
virtual void init(FullMatrix const& fullMatrix) override
{
try {
- Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
+ Super::solver_.reset(new SolverType(fullMatrix, Super::symmetricStrategy_, Super::allocInit_));
} catch (mtl::mat::umfpack::error const& e) {
error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
}
diff --git a/src/amdis/utility/MultiIndex.hpp b/src/amdis/utility/MultiIndex.hpp
index 196f12d70829ba7f3abbbfb24f72ce625e9e02b4..424e70aedd8da13f50bd1f5b780525d1e733b37a 100644
--- a/src/amdis/utility/MultiIndex.hpp
+++ b/src/amdis/utility/MultiIndex.hpp
@@ -1,17 +1,23 @@
#pragma once
#include <dune/common/reservedvector.hh>
-#include <dune/functions/common/indexaccess.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
-#include <amdis/common/Mpl.hpp>
-
namespace AMDiS
{
- std::size_t flatMultiIndex(std::size_t idx) { return idx; }
+ std::size_t flatMultiIndex(std::size_t idx)
+ {
+ return idx;
+ }
- std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx) { return idx[0]; }
+ std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx)
+ {
+ return idx[0];
+ }
- std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx) { return idx[0]; }
+ std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx)
+ {
+ return idx[0];
+ }
} // end namespace AMDiS
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b5eb3aaba4cf03f3d7803773612b3456642963a3..afb4a3ccb701712555aa9c954ca56d6c346907e1 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -8,6 +8,9 @@ dune_add_test(SOURCES ClonablePtrTest.cpp
dune_add_test(SOURCES ConceptsTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES DOFVectorTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES ExpressionsTest.cpp
LINK_LIBRARIES amdis
CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d")
diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1706a5b90537cb71dbd8bb2ec952e1fc5f5285e
--- /dev/null
+++ b/test/DOFVectorTest.cpp
@@ -0,0 +1,56 @@
+#include <dune/common/filledarray.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <amdis/LinearAlgebra.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+using namespace Dune::Functions::BasisFactory;
+
+template <class B, class T>
+void test_dofvector(B const& basis, DOFVector<B,T>& vec)
+{
+ AMDIS_TEST( vec.size() == basis.dimension() );
+ vec.compress();
+
+ vec = T(0);
+ vec[0] = T(1);
+ auto m0 = std::min_element(std::begin(vec.vector()), std::end(vec.vector()));
+ auto m1 = std::max_element(std::begin(vec.vector()), std::end(vec.vector()));
+ AMDIS_TEST( *m0 == T(0) );
+ AMDIS_TEST( *m1 == T(1) );
+}
+
+int main()
+{
+ // create grid
+ Dune::FieldVector<double, 2> L; L = 1.0;
+ auto s = Dune::filledArray<2>(1);
+ Dune::YaspGrid<2> grid(L, s);
+
+ // create basis
+ auto basis = makeBasis(grid.leafGridView(),
+ composite(power<2>(lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic()));
+
+ using Basis = decltype(basis);
+ DOFVector<Basis> vec1(basis);
+ test_dofvector(basis, vec1);
+
+ DOFVector<Basis, float> vec2(basis);
+ test_dofvector(basis, vec2);
+
+ auto vec3 = makeDOFVector(basis);
+ test_dofvector(basis, vec3);
+
+ auto vec4 = makeDOFVector<float>(basis);
+ test_dofvector(basis, vec4);
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+ DOFVector vec5(basis);
+ test_dofvector(basis, vec5);
+#endif
+}
\ No newline at end of file