diff --git a/dune.module b/dune.module
index 58ea718f31ce5cf7c3a772768e36fafcfc563de1..57f4d63a71c28d3e24d78b28e92466da3d288893 100644
--- a/dune.module
+++ b/dune.module
@@ -7,4 +7,4 @@ Module: dune-mmdg
Version: 0.1
Maintainer: maximilian.hoerl@mathematik.uni-stuttgart.de
#depending on
-Depends: dune-common dune-geometry dune-mmesh
+Depends: dune-common dune-geometry dune-mmesh dune-istl
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 67227159f95b204edecff0331c6507ee45f3715f..d115dfa2a5af452db9cdd78380620c9589efa1a9 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -10,6 +10,10 @@
#include <dune/geometry/quadraturerules.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/solver.hh>
+#include <dune/istl/umfpack.hh>
+
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
template<class GridView, class Mapper, class Problem>
@@ -18,8 +22,8 @@ class DG
public:
using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
- using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS
- using Vector = Dune::DynamicVector<Scalar>;
+ using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
+ using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
//constructor
DG (const GridView& gridView, const Mapper& mapper,
@@ -27,9 +31,10 @@ public:
gridView_(gridView), mapper_(mapper), problem_(problem),
dof((1 + dim) * gridView.size(0))
{
- A = Matrix(dof, dof, 0.0); //initialize stiffness matrix
- b = Vector(dof, 0.0); //initialize load vector
- d = Vector(dof, 0.0); //initialize solution vector
+ //initialize stiffness matrix A, load vector b and solution vector d
+ A = Matrix(dof, dof, 4, 0.1, Matrix::implicit);
+ b = Vector(dof);
+ d = Vector(dof);
}
const void operator() (const Scalar K, const Scalar mu)
@@ -37,8 +42,9 @@ public:
//assemble stiffness matrix A and load vector b
assembleSLE(K, mu);
- //NOTE: what would be an appropiate solver here?
- A.solve(d, b);
+ Dune::InverseOperatorResult result;
+ Dune::UMFPack<Matrix> solver(A);
+ solver.apply(d, b, result);
//storage for pressure data, for each element we store the
//pressure at the corners of the element, the VTKFunction will
@@ -150,7 +156,7 @@ private:
{
//exact evaluation of
// int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
- A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE + i + 1) += K * elemVol;
}
//iterate over all intersection with the boundary of elem
@@ -206,7 +212,7 @@ private:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
- A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
+ A.entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
@@ -221,7 +227,7 @@ private:
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A[elemIdxSLE + i + 1][elemIdxSLE] +=
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
@@ -232,7 +238,7 @@ private:
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
- A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
@@ -246,10 +252,11 @@ private:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
- A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
+ A.entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
//stiffness matrix A is symmetric
- A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
+ A.entry(neighborIdxSLE, elemIdxSLE) +=
+ A.entry(elemIdxSLE, neighborIdxSLE);
for (int i = 0; i < dim; i++)
{
@@ -261,7 +268,7 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A[elemIdxSLE + i + 1][neighborIdxSLE] +=
+ A.entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
-mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
//we use the relations
@@ -272,14 +279,14 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A[elemIdxSLE][neighborIdxSLE + i + 1] +=
+ A.entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric
- A[neighborIdxSLE][elemIdxSLE + i + 1] +=
- A[elemIdxSLE + i + 1][neighborIdxSLE];
- A[neighborIdxSLE + i + 1][elemIdxSLE] +=
- A[elemIdxSLE][neighborIdxSLE + i + 1];
+ A.entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
+ A.entry(elemIdxSLE + i + 1, neighborIdxSLE);
+ A.entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
+ A.entry(elemIdxSLE, neighborIdxSLE + i + 1);
for (int j = 0; j <= i; j++)
{
@@ -295,14 +302,14 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_neighbor,j) ds
// = -0.5 * K * normal[i] * int_intersct x_j ds
- A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
+ A.entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
//stiffness matrix A is symmetric
- A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
- A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
+ A.entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
+ A.entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
if (i != j)
{
@@ -317,14 +324,14 @@ private:
// int_intersct avg(K*grad(phi_elem,j))
// *jump(phi_neighbor,i) ds
// = -0.5 * K * normal[j] * int_intersct x_i ds
- A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
+ A.entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric
- A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
- A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
+ A.entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ A.entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
}
}
}
@@ -339,7 +346,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds
// = K * normal[i] * vol(intersct)
- A[elemIdxSLE + i + 1][elemIdxSLE] +=
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
@@ -351,7 +358,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
- A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
@@ -363,22 +370,25 @@ private:
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
- A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
+ A.entry(elemIdxSLE, elemIdxSLE + i + 1) =
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE);
for(int j = 0; j < i; j++)
{
- A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
- A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
+ A.entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
+ A.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
}
}
}
//NOTE: check if A is symmetric
- for (int i = 0; i < dof; i++)
- for (int j = 0; j < i; j++)
- /*assert*/if(std::abs(A[i][j] - A[j][i]) >
+ /* for (int i = 0; i < dof; i++)
+ for (int j = 0; j < i; j++) */
+ /*assert*//*if(std::abs(A[i][j] - A[j][i]) >
std::numeric_limits<Scalar>::epsilon())
- std::cout << i << ", " << j << std::endl;
+ std::cout << i << ", " << j << std::endl; */
+
+ A.compress();
}