Skip to content
Snippets Groups Projects
Commit 295dac6f authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

use BCRSMatrix and UMFPack solver

parent a0c1abd6
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune.module
+ 1
− 1
View file @ 295dac6f
......@@ -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
dune/mmdg/dg.hh
+ 44
− 34
View file @ 295dac6f
......@@ -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();
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment