Commit 7749968f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files
parents adf8a07d 951bc381
...@@ -34,25 +34,26 @@ DecGrid<GridBase> grid(*gridBase); ...@@ -34,25 +34,26 @@ DecGrid<GridBase> grid(*gridBase);
This grid adapter provides iterators for incidence relations, like edges of a vertex or faces of an edge. This grid adapter provides iterators for incidence relations, like edges of a vertex or faces of an edge.
Example: Iterate over all edges of the grid and store the coefficients of a Grad-Div laplacian Example: Iterate over all edges of the grid (grid-view) and store the coefficients of a Grad-Div laplacian
in a sparse matrix: in a sparse matrix:
```c++ ```c++
DOFMatrix<double> A(grid.size(1), grid.size(1)); // sparse matrix with num_rows=#edges auto gv = grid.leafGridView();
DOFMatrix<double> A(gv.size(1), gv.size(1)); // sparse matrix with num_rows=#edges
auto const& I = grid.indexSet(); auto const& I = gv.indexSet();
{ auto ins = A.inserter(); { auto ins = A.inserter();
// iterator over all edges in the grid // iterator over all edges in the grid
for (auto const& e : edges(grid.gridView())) { for (auto const& e : edges(gv)) {
// iterate over all vertices incident to edge e // iterate over all vertices incident to edge e
for (auto const& v : vertices(e)) { for (auto const& v : vertices(e)) {
auto factor1 = v.sign(e) * grid.dual_volume(v); auto factor1 = v.sign(e) * gv.dual_volume(v);
// iterate over all edges incident to vertex v // iterate over all edges incident to vertex v
for (auto const& e_ : edges(v)) { for (auto const& e_ : edges(v)) {
auto factor2 = v.sign(e_) * grid.dual_volume(e_) / grid.volume(e_); auto factor2 = v.sign(e_) * gv.dual_volume(e_) / gv.volume(e_);
ins(I.index(e), I.index(e_)) << factor1 * factor2; ins(I.index(e), I.index(e_)) << factor1 * factor2;
} }
......
...@@ -107,7 +107,7 @@ BRIEF_MEMBER_DESC = YES ...@@ -107,7 +107,7 @@ BRIEF_MEMBER_DESC = YES
# brief descriptions will be completely suppressed. # brief descriptions will be completely suppressed.
# The default value is: YES. # The default value is: YES.
REPEAT_BRIEF = YES REPEAT_BRIEF = NO
# This tag implements a quasi-intelligent brief description abbreviator that is # This tag implements a quasi-intelligent brief description abbreviator that is
# used to form the text in various listings. Each string in this list, if found # used to form the text in various listings. Each string in this list, if found
...@@ -756,14 +756,22 @@ WARN_LOGFILE = ...@@ -756,14 +756,22 @@ WARN_LOGFILE =
# Note: If this tag is empty the current directory is searched. # Note: If this tag is empty the current directory is searched.
INPUT = @doxy_main_page@ \ INPUT = @doxy_main_page@ \
@PROJECT_SOURCE_DIR@/src \ @PROJECT_SOURCE_DIR@/dune/dec \
@PROJECT_SOURCE_DIR@/src/common \ @PROJECT_SOURCE_DIR@/dune/dec/common \
@PROJECT_SOURCE_DIR@/src/linear_algebra \ @PROJECT_SOURCE_DIR@/dune/dec/fixmatvec \
@PROJECT_SOURCE_DIR@/src/linear_algebra/interface \ @PROJECT_SOURCE_DIR@/dune/dec/halfedgegrid \
@PROJECT_SOURCE_DIR@/src/linear_algebra/krylov \ @PROJECT_SOURCE_DIR@/dune/dec/linear_algebra \
@PROJECT_SOURCE_DIR@/src/operators \ @PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/block \
@PROJECT_SOURCE_DIR@/src/simplegrid \ @PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/interface \
@PROJECT_SOURCE_DIR@/src/utility @PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/iteration \
@PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/krylov \
@PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/multigrid \
@PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/preconditioner \
@PROJECT_SOURCE_DIR@/dune/dec/linear_algebra/smoother \
@PROJECT_SOURCE_DIR@/dune/dec/operators \
@PROJECT_SOURCE_DIR@/dune/dec/ranges \
@PROJECT_SOURCE_DIR@/dune/dec/simplegrid \
@PROJECT_SOURCE_DIR@/dune/dec/utility
# This tag can be used to specify the character encoding of the source files # This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
...@@ -1562,7 +1570,7 @@ EXTRA_SEARCH_MAPPINGS = ...@@ -1562,7 +1570,7 @@ EXTRA_SEARCH_MAPPINGS =
# If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output. # If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output.
# The default value is: YES. # The default value is: YES.
GENERATE_LATEX = YES GENERATE_LATEX = NO
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a # The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a
# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of # relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
...@@ -2117,7 +2125,7 @@ CLASS_GRAPH = YES ...@@ -2117,7 +2125,7 @@ CLASS_GRAPH = YES
# The default value is: YES. # The default value is: YES.
# This tag requires that the tag HAVE_DOT is set to YES. # This tag requires that the tag HAVE_DOT is set to YES.
COLLABORATION_GRAPH = YES COLLABORATION_GRAPH = NO
# If the GROUP_GRAPHS tag is set to YES then doxygen will generate a graph for # If the GROUP_GRAPHS tag is set to YES then doxygen will generate a graph for
# groups, showing the direct groups dependencies. # groups, showing the direct groups dependencies.
......
...@@ -16,10 +16,20 @@ namespace Dec ...@@ -16,10 +16,20 @@ namespace Dec
struct adjoint {}; struct adjoint {};
} }
/// Interpolation matrices to transfer between two grid levels
template <int dim, class Strategy = tag::full_weighting> template <int dim, class Strategy = tag::full_weighting>
class GridTransfer class GridTransfer
{ {
template <class, class> friend class DefaultTransfer;
public: public:
/// Constructor, stores copy of `mat1` and `mat2`.
/**
* \param mat1 Transfer-matrix from coarse to fine
* \param mat2 Transfer-matrix from fine to coarse
* \param map Index mapping: map[coarse_idx] == fine_idx
**/
template <class Matrix1, class Matrix2, class Map> template <class Matrix1, class Matrix2, class Map>
GridTransfer(Matrix1&& mat1, Matrix2&& mat2, Map&& map) GridTransfer(Matrix1&& mat1, Matrix2&& mat2, Map&& map)
: mat1_(std::forward<Matrix1>(mat1)) : mat1_(std::forward<Matrix1>(mat1))
...@@ -27,8 +37,9 @@ namespace Dec ...@@ -27,8 +37,9 @@ namespace Dec
, indexmap_(std::forward<Map>(map)) , indexmap_(std::forward<Map>(map))
{} {}
/// Interpolate a vector from coarse to fine grid
template <class Vector1, class Vector2> template <class Vector1, class Vector2>
void interpol(Vector1 const& in, Vector2& out, bool add = false) const void prolongate(Vector1 const& in, Vector2& out, bool add = false) const
{ {
assert( std::size_t(in.size()) == mat1_.num_cols() ); assert( std::size_t(in.size()) == mat1_.num_cols() );
out.resize(mat1_.num_rows()); out.resize(mat1_.num_rows());
...@@ -40,41 +51,42 @@ namespace Dec ...@@ -40,41 +51,42 @@ namespace Dec
} }
} }
/// Interpolate a vector from fine to coarse grid
template <class Vector1, class Vector2> template <class Vector1, class Vector2>
void restrict(Vector1 const& in, Vector2& out, bool add = false) const void restrict(Vector1 const& in, Vector2& out, bool add = false) const
{ {
restrict_(in, out, Strategy{}, add); restrictImpl(in, out, Strategy{}, add);
}
// simple injection
template <class Vector1, class Vector2>
void injection(Vector1 const& in, Vector2& out, bool add = false) const
{
assert( std::size_t(in.size()) == mat1_.num_rows() );
out.resize(mat1_.num_cols());
auto stored = in.eval();
if (add) {
for (std::size_t i = 0; i < indexmap_.size(); ++i)
out[i] += stored[indexmap_[i]];
} else {
for (std::size_t i = 0; i < indexmap_.size(); ++i)
out[i] = stored[indexmap_[i]];
}
}
template <class MatrixIn, class MatrixOut>
void coarse(MatrixIn const& in, MatrixOut& out) const
{
// out = (1.0/(1<<dim)) * ((mat1_.transpose() * in).pruned() * mat1_).pruned();
out = mat2_ * in * mat1_;
} }
// // simple injection
public: // implementation of restriction strategies // template <class Vector1, class Vector2>
// void injection(Vector1 const& in, Vector2& out, bool add = false) const
// {
// assert( std::size_t(in.size()) == mat1_.num_rows() );
// out.resize(mat1_.num_cols());
//
// auto stored = in.eval();
// if (add) {
// for (std::size_t i = 0; i < indexmap_.size(); ++i)
// out[i] += stored[indexmap_[i]];
// } else {
// for (std::size_t i = 0; i < indexmap_.size(); ++i)
// out[i] = stored[indexmap_[i]];
// }
// }
// template <class MatrixIn, class MatrixOut>
// void coarse(MatrixIn const& in, MatrixOut& out) const
// {
// // out = (1.0/(1<<dim)) * ((mat1_.transpose() * in).pruned() * mat1_).pruned();
// out = mat2_ * in * mat1_;
// }
private: // implementation of restriction strategies
template <class Vector1, class Vector2> template <class Vector1, class Vector2>
void restrict_(Vector1 const& in, Vector2& out, tag::full_weighting, bool add = false) const void restrictImpl(Vector1 const& in, Vector2& out, tag::full_weighting, bool add = false) const
{ {
assert( std::size_t(in.size()) == mat2_.num_cols() ); assert( std::size_t(in.size()) == mat2_.num_cols() );
out.resize(mat2_.num_rows()); out.resize(mat2_.num_rows());
...@@ -87,7 +99,7 @@ namespace Dec ...@@ -87,7 +99,7 @@ namespace Dec
} }
template <class Vector1, class Vector2> template <class Vector1, class Vector2>
void restrict_(Vector1 const& in, Vector2& out, tag::adjoint, bool add = false) const void restrictImpl(Vector1 const& in, Vector2& out, tag::adjoint, bool add = false) const
{ {
assert( std::size_t(in.size()) == mat2_.num_cols() ); assert( std::size_t(in.size()) == mat2_.num_cols() );
out.resize(mat2_.num_rows()); out.resize(mat2_.num_rows());
...@@ -101,13 +113,23 @@ namespace Dec ...@@ -101,13 +113,23 @@ namespace Dec
} }
template <class Vector1, class Vector2> template <class Vector1, class Vector2>
void restrict_(Vector1 const& in, Vector2& out, tag::injection, bool add = false) const void restrictImpl(Vector1 const& in, Vector2& out, tag::injection, bool add = false) const
{ {
this->injection(in, out, add); assert( std::size_t(in.size()) == mat1_.num_rows() );
out.resize(mat1_.num_cols());
auto stored = in.eval();
if (add) {
for (std::size_t i = 0; i < indexmap_.size(); ++i)
out[i] += stored[indexmap_[i]];
} else {
for (std::size_t i = 0; i < indexmap_.size(); ++i)
out[i] = stored[indexmap_[i]];
}
} }
public: private:
DOFMatrix<float_type> mat1_; DOFMatrix<float_type> mat1_;
DOFMatrix<float_type> mat2_; DOFMatrix<float_type> mat2_;
...@@ -115,6 +137,7 @@ namespace Dec ...@@ -115,6 +137,7 @@ namespace Dec
}; };
#ifndef DOXYGEN
// dummy implementation for grids that do not have an id-set // dummy implementation for grids that do not have an id-set
template <class GridBase, class = void> template <class GridBase, class = void>
class GradTransferFactory class GradTransferFactory
...@@ -133,7 +156,10 @@ namespace Dec ...@@ -133,7 +156,10 @@ namespace Dec
return {mat1, mat2, indexmap}; return {mat1, mat2, indexmap};
} }
}; };
#endif
/// Factory for the interpolation matrices between two grid levels. \relates GridTransfer
template <class GridBase> template <class GridBase>
class GradTransferFactory<GridBase, Void_t<typename GridBase::LocalIdSet>> class GradTransferFactory<GridBase, Void_t<typename GridBase::LocalIdSet>>
{ {
...@@ -223,7 +249,6 @@ namespace Dec ...@@ -223,7 +249,6 @@ namespace Dec
} }
} }
using Matrix = DOFMatrix<float_type>::Matrix; using Matrix = DOFMatrix<float_type>::Matrix;
for (int k = 0; k < mat2.outerSize(); ++k) { for (int k = 0; k < mat2.outerSize(); ++k) {
for (typename Matrix::InnerIterator it(mat2, k); it; ++it) for (typename Matrix::InnerIterator it(mat2, k); it; ++it)
......
...@@ -9,10 +9,37 @@ namespace Dec ...@@ -9,10 +9,37 @@ namespace Dec
{ {
/** /**
* \defgroup operators Operators * \defgroup operators Operators
* \brief Differential operators * \brief Predefined terms representing (differential) operators, to assemble a linear system.
*
* @{ * @{
**/ **/
/*
* Operators are used to simplify the assembling of linear systems. Therefore, standard
* differential operators are provided that can be chained to describe a partial differential
* equation in simple terms.
*
* Basic usage:
* ```
* DOFMatrix<double> mat;
* OperatorXYZ<GridView> op(gridview);
*
* op.build(mat); // assemble operator into matrix.
* ```
*
* The method `build()` creates a new inserter to insert coefficients of one specific
* operator. When combining several operators, use the `assemble()` method. This expects
* an inserter instead of a matrix.
* ```
* op.init(mat);
* { // local scope
* auto ins = mat.inserter();
* op.assemble(ins);
* }
* op.finish(mat);
**/
class OperatorAccessor class OperatorAccessor
{ {
template <class Model, class GridView, int N, int M> friend class Operator; template <class Model, class GridView, int N, int M> friend class Operator;
......
...@@ -13,7 +13,7 @@ ...@@ -13,7 +13,7 @@
#include "operators/FEMLaplace.hpp" #include "operators/FEMLaplace.hpp"
#include "operators/FEMCenterWeight.hpp" #include "operators/FEMCenterWeight.hpp"
#include "linear_algebra/krylov/cg.hpp" #include "linear_algebra/krylov/cg.hpp"
#include "linear_algebra/krylov/diagonal_pc.hpp" #include "linear_algebra/preconditioner/Diagonal.hpp"
#include "ranges/Enumerate.hpp" #include "ranges/Enumerate.hpp"
#include "utility/Timer.hpp" #include "utility/Timer.hpp"
...@@ -52,7 +52,8 @@ namespace Dec ...@@ -52,7 +52,8 @@ namespace Dec
b[0] = 1.0; b[0] = 1.0;
DOFVector<double> weights(grid.leafGridView(), 0, 1.0); DOFVector<double> weights(grid.leafGridView(), 0, 1.0);
auto iter = solver::pcg(A, weights, b, precon::Diagonal<DOFMatrix<double>>{A}); auto precon = precon::Diagonal{};
auto iter = solver::pcg(A, weights, b, precon.compute(A));
msg("nr of iterations: ",iter); msg("nr of iterations: ",iter);
calc_center(grid, weights); calc_center(grid, weights);
......
...@@ -106,12 +106,14 @@ namespace Dec ...@@ -106,12 +106,14 @@ namespace Dec
FixVec<S, M> lower_tri_solve(MatExpr<E0, T0, N, M> const& mat, V const& b); FixVec<S, M> lower_tri_solve(MatExpr<E0, T0, N, M> const& mat, V const& b);
/// \brief Solve a overdetermined system, using a qr factorization approach. /// \brief Solve an overdetermined system, using a qr factorization approach.
/** /**
* Requirement: * Requirement:
* - `T0` models concept `concepts::Arithmetic` * - `T0` models concept `concepts::Arithmetic`
* - `V` models concept `concepts::Vector` with size `N` * - `V` models concept `concepts::Vector` with size `N`
* - `N >= M` * - `N >= M`
*
* Result:
* - `S` := `std::common_type_t<T0, Value_t<V>>` * - `S` := `std::common_type_t<T0, Value_t<V>>`
**/ **/
template <class E0, class T0, class V, std::size_t N, std::size_t M> template <class E0, class T0, class V, std::size_t N, std::size_t M>
......
...@@ -439,7 +439,7 @@ namespace Dec ...@@ -439,7 +439,7 @@ namespace Dec
} }
/** \relates FixMat /** \relates FixMat
* \brief Returns the hermidean of the matrix * \brief Returns the hermitian of the matrix
**/ **/
template <class E, class T, std::size_t N, std::size_t M> template <class E, class T, std::size_t N, std::size_t M>
expression::Hermitian<E,T,N,M> expression::Hermitian<E,T,N,M>
...@@ -525,9 +525,9 @@ namespace Dec ...@@ -525,9 +525,9 @@ namespace Dec
error_exit("Determinant not yet implemented for this matrix!"); error_exit("Determinant not yet implemented for this matrix!");
} }
/** \relates FixMat #ifndef DOXYGEN
* \brief Inverse of a 1x1 matrix
**/ // Inverse of a 1x1 matrix
template <class E, class T> template <class E, class T>
FixMat<T,1,1> inv(MatExpr<E, T, 1, 1> const& mat) FixMat<T,1,1> inv(MatExpr<E, T, 1, 1> const& mat)
{ {
...@@ -535,9 +535,7 @@ namespace Dec ...@@ -535,9 +535,7 @@ namespace Dec
return result; return result;
} }
/** \relates FixMat // Inverse of a 2x2 matrix
* \brief Inverse of a 2x2 matrix
**/
template <class E, class T> template <class E, class T>
FixMat<T,2,2> inv(MatExpr<E, T, 2, 2> const& mat) FixMat<T,2,2> inv(MatExpr<E, T, 2, 2> const& mat)
{ {
...@@ -545,9 +543,7 @@ namespace Dec ...@@ -545,9 +543,7 @@ namespace Dec
return result /= det(mat); return result /= det(mat);
} }
/** \relates FixMat // Inverse of a 3x3 matrix
* \brief Inverse of a 3x3 matrix
**/
template <class E, class T> template <class E, class T>
FixMat<T,3,3> inv(MatExpr<E, T, 3, 3> const& mat) FixMat<T,3,3> inv(MatExpr<E, T, 3, 3> const& mat)
{ {
...@@ -559,6 +555,8 @@ namespace Dec ...@@ -559,6 +555,8 @@ namespace Dec
return result /= det(A); return result /= det(A);
} }
#endif
namespace expression namespace expression
{ {
...@@ -616,7 +614,7 @@ namespace Dec ...@@ -616,7 +614,7 @@ namespace Dec
} }
/** \relates FixMat /** \relates FixMat
* \brief Inverse of a square matrix * \brief Inverse of a square NxN matrix
**/ **/
template <class E, class T, std::size_t N> template <class E, class T, std::size_t N>
FixMat<T,N,N> inv(MatExpr<E, T, N, N> const& mat) FixMat<T,N,N> inv(MatExpr<E, T, N, N> const& mat)
......
...@@ -7,7 +7,19 @@ ...@@ -7,7 +7,19 @@
namespace Dec namespace Dec
{ {
/// Block-Jacobi method, can be used as solver and preconditioner /**
* \addtogroup linear_solver
* @{
**/
/// \brief Block-Jacobi method, can be used as solver and preconditioner
/**
* The (variadic) template parameters `Solver...` must model `LinearSolver`, i.e.
* must provide the methods
* - `compute()` (initialization of the solver)
* - `solve()` (solve with initial solution u=0)
* - `solveWithGuess()` (solve with provided initial solution)
**/
template <class... Solver> template <class... Solver>
struct BlockJacobi struct BlockJacobi
{ {
...@@ -21,7 +33,7 @@ namespace Dec ...@@ -21,7 +33,7 @@ namespace Dec
: solver_(&solver...) : solver_(&solver...)
{} {}
/// Calls `compute` on all solvers, with the diagonal matrix A_ii and stores /// Calls `compute` on all solvers, with the diagonal matrix \f$ A_{ii}\f$ and stores
/// pointer to the matrix. /// pointer to the matrix.
BlockJacobi& compute(Matrix const& mat) BlockJacobi& compute(Matrix const& mat)
{ {
...@@ -34,7 +46,7 @@ namespace Dec ...@@ -34,7 +46,7 @@ namespace Dec
return *this; return *this;
} }
/// Solve A_ii u_i = b_i /// Solve \f$ A_{ii} u_i = b_i\f$
void solve(Vector const& b, Vector& u) const void solve(Vector const& b, Vector& u) const
{ {
forEach(range_<0,N>, [this,&b,&u](auto const _i) forEach(range_<0,N>, [this,&b,&u](auto const _i)
...@@ -43,7 +55,7 @@ namespace Dec ...@@ -43,7 +55,7 @@ namespace Dec
}); });
} }
/// Solve A_ii d = b_i - (A*u)_i, u += d /// Solve \f$ A_{ii} d = b_i - (A\cdot u)_i,\; u += d\f$
void solveWithGuess(Vector const& b, Vector& u) const void solveWithGuess(Vector const& b, Vector& u) const
{ {
assert( mat_ != nullptr && "Call BlockJacobi::compute() first" ); assert( mat_ != nullptr && "Call BlockJacobi::compute() first" );
...@@ -63,11 +75,13 @@ namespace Dec ...@@ -63,11 +75,13 @@ namespace Dec
}; };
/// Generator for Block-Jacobi solver /// Generator for Block-Jacobi solver, \relates BlockJacobi
template <class... Solver> template <class... Solver>
BlockJacobi<Solver...> makeBlockJacobi(Solver&... solver) BlockJacobi<Solver...> makeBlockJacobi(Solver&... solver)
{ {
return {solver...}; return {solver...};
} }
/** @} */
} // end namespace Dec } // end namespace Dec