Skip to content

GitLab

Commit 951bc381 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'master' into feature/multilevel-grid

parents d80b227e d8b987ef
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 951bc381
......@@ -11,9 +11,10 @@
# ignore build directories
build*/
install*/
output/
error.txt
*.out
*.out.*
*.tmp
*.*-swp
\ No newline at end of file
*.*-swp
.gitlab-ci.yml
View file @ 951bc381
......@@ -8,9 +8,7 @@ stages:
variables:
BUILD_DIRECTORY: "./build-cmake/"
COMPILER: UNDEFINED
LINEAR_ALGEBRA_BACKEND: UNDEFINED
CMAKE_BUILD_TYPE: UNDEFINED
COMPILER_FLAGS: "-Wall -Wextra -pedantic"
# ------------------------------------------------------------------------------
# templates
......@@ -25,7 +23,7 @@ variables:
expire_in: 6 h
script:
- dunecontrol --opts=gitlab-ci.$COMPILER.opts --current all
- cmake -DLINEAR_ALGEBRA_BACKEND=$LINEAR_ALGEBRA_BACKEND -DCMAKE_BUILD_TYPE=$CMAKE_BUILD_TYPE -DCMAKE_CXX_FLAGS=$COMPILER_FLAGS $BUILD_DIRECTORY
- cmake -DCMAKE_BUILD_TYPE=$CMAKE_BUILD_TYPE $BUILD_DIRECTORY
# test
.test_template: &test-definition
......@@ -45,88 +43,28 @@ variables:
setup:gcc-Default-Debug:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: Default
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:gcc-Default-Release:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: Default
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
setup:gcc-Eigen-Debug:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: Eigen
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:gcc-Eigen-Release:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: Eigen
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
setup:gcc-PETSc-Debug:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: PETSc
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:gcc-PETSc-Release:
variables:
COMPILER: gcc
LINEAR_ALGEBRA_BACKEND: PETSc
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
setup:clang-Default-Debug:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: Default
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:clang-Default-Release:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: Default
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
setup:clang-Eigen-Debug:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: Eigen
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:clang-Eigen-Release:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: Eigen
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
setup:clang-PETSc-Debug:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: PETSc
CMAKE_BUILD_TYPE: Debug
<<: *setup-definition
setup:clang-PETSc-Release:
variables:
COMPILER: clang
LINEAR_ALGEBRA_BACKEND: PETSc
CMAKE_BUILD_TYPE: Release
<<: *setup-definition
# ------------------------------------------------------------------------------
# test
......@@ -140,26 +78,6 @@ test:gcc-Default-Release:
- setup:gcc-Default-Release
<<: *test-definition
test:gcc-Eigen-Debug:
dependencies:
- setup:gcc-Eigen-Debug
<<: *test-definition
test:gcc-Eigen-Release:
dependencies:
- setup:gcc-Eigen-Release
<<: *test-definition
test:gcc-PETSc-Debug:
dependencies:
- setup:gcc-PETSc-Debug
<<: *test-definition
test:gcc-PETSc-Release:
dependencies:
- setup:gcc-PETSc-Release
<<: *test-definition
test:clang-Default-Debug:
dependencies:
- setup:clang-Default-Debug
......@@ -170,26 +88,6 @@ test:clang-Default-Release:
- setup:clang-Default-Release
<<: *test-definition
test:clang-Eigen-Debug:
dependencies:
- setup:clang-Eigen-Debug
<<: *test-definition
test:clang-Eigen-Release:
dependencies:
- setup:clang-Eigen-Release
<<: *test-definition
test:clang-PETSc-Debug:
dependencies:
- setup:clang-PETSc-Debug
<<: *test-definition
test:clang-PETSc-Release:
dependencies:
- setup:clang-PETSc-Release
<<: *test-definition
# ------------------------------------------------------------------------------
# demo
......@@ -204,26 +102,6 @@ demo:gcc-Default-Release:
- setup:gcc-Default-Release
<<: *demo-definition
demo:gcc-Eigen-Debug:
dependencies:
- setup:gcc-Eigen-Debug
<<: *demo-definition
demo:gcc-Eigen-Release:
dependencies:
- setup:gcc-Eigen-Release
<<: *demo-definition
demo:gcc-PETSc-Debug:
dependencies:
- setup:gcc-PETSc-Debug
<<: *demo-definition
demo:gcc-PETSc-Release:
dependencies:
- setup:gcc-PETSc-Release
<<: *demo-definition
demo:clang-Default-Debug:
dependencies:
- setup:clang-Default-Debug
......@@ -234,23 +112,3 @@ demo:clang-Default-Release:
- setup:clang-Default-Release
<<: *demo-definition
demo:clang-Eigen-Debug:
dependencies:
- setup:clang-Eigen-Debug
<<: *demo-definition
demo:clang-Eigen-Release:
dependencies:
- setup:clang-Eigen-Release
<<: *demo-definition
demo:clang-PETSc-Debug:
dependencies:
- setup:clang-PETSc-Debug
<<: *demo-definition
demo:clang-PETSc-Release:
dependencies:
- setup:clang-PETSc-Release
<<: *demo-definition
README.md
View file @ 951bc381
......@@ -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.
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:
```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();
// 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
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
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;
}
......
consistency.py deleted 100644 → 0
View file @ d80b227e
from sympy import *
from sympy.tensor.array import Array, tensorcontraction, tensorproduct, derive_by_array
def factorial(n):
return 1 if n <= 1 else n*factorial(n-1)
def get(tensor,n):
return tensor[0] if n <= 1 else get(tensor[0],n-1)
def taylor_expansion(f, vec, n=5):
h = symbols('h')
one = sympify('1')
result = f;
df = f;
for i in range(1,n):
df = derive_by_array(df, [x,y])
dfh = df;
for j in range(0,i):
dfh = tensorcontraction(tensorproduct(dfh, vec), (i-j-1,i))
result0 = (one/factorial(i) * dfh).tolist()
result += get(result0, i)
return result + Order(h**(n+1))
def taylor_expansion_2d(f, vec):
fx, fy = f.diff(x), f.diff(y)
f1 = Matrix([fx, fy])
fxx, fxy, fyx, fyy = fx.diff(x), fx.diff(y), fy.diff(x), fy.diff(y)
f2 = Matrix([[fxx, fxy], [fyx, fyy]])
f2x = f2.diff(x)
f2y = f2.diff(y)
h = symbols('h')
one = sympify('1')
return f + vec.dot(f1) + one/2 * vec.dot(f2*vec) + one/6 * vec.dot((f2x*vec[0] + f2y*vec[1])*vec) + Order(h**4)
def laplace_2d(f):
return f.diff(x,2) + f.diff(y,2)
init_printing()
x,y = symbols('x y')
f = Function('f')(x,y)
h = symbols('h')
h1 = h * Matrix([ 1, 0])
h2 = h * Matrix([-1, 0])
h3 = h * Matrix([ 0, 1])
h4 = h * Matrix([ 0, -1])
h5 = h/2 * Matrix([ 1, 1])
h6 = h/2 * Matrix([ 1, -1])
h7 = h/2 * Matrix([-1, -1])
h8 = h/2 * Matrix([-1, 1])
f_v1 = taylor_expansion(f,h1)
f_v2 = taylor_expansion(f,h2)
f_v3 = taylor_expansion(f,h3)
f_v4 = taylor_expansion(f,h4)
f_v5 = taylor_expansion(f,h5)
f_v6 = taylor_expansion(f,h6)
f_v7 = taylor_expansion(f,h7)
f_v8 = taylor_expansion(f,h8)
one = sympify('1')
two = sympify('2')
four = sympify('4')
stencil = simplify( two/(h**2) * (four*f - f_v5 - f_v6 - f_v7 - f_v8) + laplace_2d(f) )
print "stencil:"
print stencil
alpha = symbols('alpha')
tau = alpha*h
tau1 = (h - tau)/sqrt(two)
vol = h**2 - tau1**2
w1 = alpha
w2 = one - w1
correction_factor = (one + two*alpha - alpha**2)/(one + alpha)
weighted_stencil = simplify( one/vol * (four*(w1+w2)*f - w1*f_v1 - w1*f_v2 - w1*f_v3 - w1*f_v4 - w2*f_v5 - w2*f_v6 - w2*f_v7 - w2*f_v8) )
print "weighted_stencil:"
print simplify(weighted_stencil)
vol2 = tau1**2
weighted_stencil2 = simplify( w2/vol2 * (four*f - f_v5 - f_v6 - f_v7 - f_v8) + laplace_2d(f) )
print "weighted_stencil2:"
print simplify(weighted_stencil2)
g = exp(sympify(10) * (x*x + y*y))
print "laplace:"
print laplace_2d(g)
dec.opts
View file @ 951bc381
CMAKE_FLAGS=" \
-DCMAKE_INSTALL_PREFIX:PATH=/opt/software/dune/dune-dec \
-DCMAKE_BUILD_TYPE=RelWithDebInfo \
-DEIGEN3_INCLUDE_DIR:PATH=/opt/software/eigen/include/eigen3"
-DEIGEN3_INCLUDE_DIR:PATH=/opt/software/Eigen/git/include/eigen3"
dune/dec/common/Array.hpp
View file @ 951bc381
......@@ -235,16 +235,6 @@ namespace Dec
}
/** \ingroup vector_norms
* \brief The default norm = alias for \ref two_norm
**/
template <class T, std::size_t N>
auto norm(std::array<T, N> const& x)
{
return two_norm(x);
}
// ----------------------------------------------------------------------------
......
dune/dec/common/FieldMatVec.hpp
View file @ 951bc381
......@@ -173,41 +173,6 @@ namespace Dec
// ----------------------------------------------------------------------------
template <class T,
REQUIRES(concepts::Arithmetic<T>) >
T norm(T const& x)
{
return std::abs(x);
}
template <class T, int D>
T norm(Dune::FieldVector<T,D> const& x)
{
return x.two_norm();
}
template <class T, int N, int M>
T norm(Dune::FieldMatrix<T,N,M> const& x)
{
return x.frobenius_norm();
}
template <class T,
REQUIRES(concepts::Arithmetic<T>) >
T norm(Dune::FieldVector<T,1> const& x)
{
return norm(x[0]);
}
template <class T,
REQUIRES(concepts::Arithmetic<T>) >
T norm(Dune::FieldMatrix<T,1,1> const& x)
{
return norm(x[0][0]);
}
// ----------------------------------------------------------------------------
/// The euklidean distance between two vectors = |lhs-rhs|_2
template <class T, int N>
T distance(Dune::FieldVector<T, N> const& lhs, Dune::FieldVector<T, N> const& rhs)
......
dune/dec/fixmatvec/ReductionOperations.hpp
View file @ 951bc381
......@@ -188,26 +188,6 @@ namespace Dec
}
/** \relates FixMat
* \brief The default norm = alias for \ref two_norm
**/
template <class V,
REQUIRES(concepts::Vector<V>) >
auto norm(V const& x)
{
return two_norm(x);
}
/** \relates FixMat
* \brief Default norm is \ref frobenius_norm
**/
template <class E, class T, std::size_t N, std::size_t M>
T norm(MatExpr<E, T, N, M> const& mat)
{
return frobenius_norm(mat);
}
// ----------------------------------------------------------------------------
......
dune/dec/linear_algebra/block/BlockMatrix.hpp
View file @ 951bc381
......@@ -18,6 +18,7 @@ namespace Dec
using Super = std::array<std::array<DOFMatrix<T>, M>, N>;
public:
using value_type = T;
using size_type = std::size_t;
using Matrix = DOFMatrix<T>;
......@@ -28,9 +29,11 @@ namespace Dec
Expects( n == N && m == M );
}
/// Access the (i,j)'th block
DOFMatrix<T>& operator()(size_type i, size_type j) { return (*this)[i][j]; }
DOFMatrix<T> const& operator()(size_type i, size_type j) const { return (*this)[i][j]; }
/// Return number of non-zeros
std::size_t nonZeros() const
{
std::size_t nnz = 0;
......@@ -40,6 +43,17 @@ namespace Dec
return nnz;
}
/// Sets all matrices to zero
void setZero()
{
for (std::size_t r = 0; r < N; ++r)
for (std::size_t c = 0; c < M; ++c)
(*this)(r, c).setZero();
}
/// Return the global sparse matrix
Matrix const& matrix() const { return matrix_; }
/// Create a sparse-matrix from the blocks
Matrix const& compress()
{
......@@ -69,16 +83,6 @@ namespace Dec
return matrix_;
}
Matrix const& matrix() const { return matrix_; }
/// Sets all matrices to zero
void setZero()
{
for (std::size_t r = 0; r < N; ++r)
for (std::size_t c = 0; c < M; ++c)
(*this)(r, c).setZero();
}
/// Return the row/col block sizes, as pair of vectors
auto blockSizes() const
{
......@@ -107,7 +111,7 @@ namespace Dec
private:
Matrix matrix_;
Matrix matrix_; // global sparse matrix
};
......
dune/dec/linear_algebra/block/BlockVector.hpp
View file @ 951bc381
......@@ -83,7 +83,7 @@ namespace Dec
BlockVector(BlockVector const& that, tag::ressource)
: mapper_(that.mapper_)
, contiguous_(mapper_.size())
, vector_(mapper_.size())
{
for (std::size_t i = 0; i < N; ++i)
(*this)[i].resize(mapper_.size(i));
......@@ -161,23 +161,23 @@ namespace Dec
(*this)[i].setZero();
}
DOFVector<T>& vector() { return contiguous_; }
DOFVector<T> const& vector() const { return contiguous_; }
DOFVector<T>& vector() { return vector_; }
DOFVector<T> const& vector() const { return vector_; }
/// Create a contiguous vector from the blocks
DOFVector<T> const& compress()
{
contiguous_.resize(mapper_.size());
vector_.resize(mapper_.size());
for (std::size_t i = 0; i < N; ++i)
contiguous_.segment(mapper_.shift(i),mapper_.size(i)) = (*this)[i];
vector_.segment(mapper_.shift(i),mapper_.size(i)) = (*this)[i];
return contiguous_;
return vector_;
}
private:
BlockMapper mapper_;
DOFVector<T> contiguous_;
DOFVector<T> vector_; // global vector
};
} // end namespace Dec
dune/dec/linear_algebra/block/Operations.hpp
View file @ 951bc381
......@@ -28,18 +28,12 @@ namespace Dec
}
template <class T, std::size_t N>
inline T norm(BlockVector<T,N> const& vec)
inline T two_norm(BlockVector<T,N> const& vec)
{
using std::sqrt;
return sqrt(unary_dot(vec));
}
template <class T, std::size_t N>
inline T two_norm(BlockVector<T,N> const& vec)
{
return norm(vec);
}
template <class T, std::size_t N>
T dot(BlockVector<T,N> const& a, BlockVector<T,N> const& b)
{
......
dune/dec/linear_algebra/default/Operations.hpp
View file @ 951bc381
......@@ -26,18 +26,12 @@ namespace Dec
}
template <class T>
inline T norm(DOFVector<T> const& vec)
inline T two_norm(DOFVector<T> const& vec)
{
using std::sqrt;
return sqrt(unary_dot(vec));
}
template <class T>
inline T two_norm(DOFVector<T> const& vec)
{
return norm(vec);
}
template <class T>
inline T dot(DOFVector<T> const& vec1, DOFVector<T> const& vec2)
{
......
dune/dec/linear_algebra/eigen/Operations.hpp
View file @ 951bc381
......@@ -12,16 +12,10 @@ namespace Dec
return vec;
}
template <class T>
inline T norm(DOFVector<T> const& vec)
{
return vec.norm();
}
template <class T>
inline T two_norm(DOFVector<T> const& vec)
{
return norm(vec);
return vec.norm();
}