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

Merge branch 'feature/sparse'

parents 105ccc22 e1a3041b
No related branches found
No related tags found
No related merge requests found
...@@ -7,4 +7,4 @@ Module: dune-mmdg ...@@ -7,4 +7,4 @@ Module: dune-mmdg
Version: 0.1 Version: 0.1
Maintainer: maximilian.hoerl@mathematik.uni-stuttgart.de Maintainer: maximilian.hoerl@mathematik.uni-stuttgart.de
#depending on #depending on
Depends: dune-common dune-geometry dune-mmesh Depends: dune-common dune-geometry dune-mmesh dune-istl
...@@ -10,6 +10,10 @@ ...@@ -10,6 +10,10 @@
#include <dune/geometry/quadraturerules.hh> #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> #include <dune/mmdg/nonconformingp1vtkfunction.hh>
template<class GridView, class Mapper, class Problem> template<class GridView, class Mapper, class Problem>
...@@ -18,8 +22,8 @@ class DG ...@@ -18,8 +22,8 @@ class DG
public: public:
using Scalar = typename GridView::ctype; using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension; static constexpr int dim = GridView::dimension;
using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
using Vector = Dune::DynamicVector<Scalar>; using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
using VTKFunction = Dune::VTKFunction<GridView>; using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView, using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<double>>; Dune::DynamicVector<double>>;
...@@ -30,9 +34,10 @@ public: ...@@ -30,9 +34,10 @@ public:
gridView_(gridView), mapper_(mapper), problem_(problem), gridView_(gridView), mapper_(mapper), problem_(problem),
dof((1 + dim) * gridView.size(0)) dof((1 + dim) * gridView.size(0))
{ {
A = Matrix(dof, dof, 0.0); //initialize stiffness matrix //initialize stiffness matrix A, load vector b and solution vector d
b = Vector(dof, 0.0); //initialize load vector A = std::make_shared<Matrix>(dof, dof, 10, 1, Matrix::implicit); //TODO
d = Vector(dof, 0.0); //initialize solution vector b = Vector(dof);
d = Vector(dof);
} }
const void operator() (const Scalar K, const Scalar mu) const void operator() (const Scalar K, const Scalar mu)
...@@ -40,8 +45,9 @@ public: ...@@ -40,8 +45,9 @@ public:
//assemble stiffness matrix A and load vector b //assemble stiffness matrix A and load vector b
assembleSLE(K, mu); assembleSLE(K, mu);
//NOTE: what would be an appropiate solver here? Dune::InverseOperatorResult result;
A.solve(d, b); Dune::UMFPack<Matrix> solver(*A);
solver.apply(d, b, result);
//write solution to a vtk file //write solution to a vtk file
writeVTKoutput(); writeVTKoutput();
...@@ -108,7 +114,7 @@ private: ...@@ -108,7 +114,7 @@ private:
{ {
//exact evaluation of //exact evaluation of
// int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV // 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 //iterate over all intersection with the boundary of elem
...@@ -164,7 +170,7 @@ private: ...@@ -164,7 +170,7 @@ private:
//exact evaluation of //exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds // 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 if (intersct.neighbor()) //intersct has neighboring element
{ {
...@@ -179,7 +185,7 @@ private: ...@@ -179,7 +185,7 @@ private:
//and //and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct) // = 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; mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++) for (int j = 0; j <= i; j++)
...@@ -190,7 +196,7 @@ private: ...@@ -190,7 +196,7 @@ private:
//and //and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_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] += mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j] - 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]); + normal[j] * linearIntegrals[i]);
...@@ -204,10 +210,11 @@ private: ...@@ -204,10 +210,11 @@ private:
//exact evaluation of //exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds // 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 //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++) for (int i = 0; i < dim; i++)
{ {
...@@ -219,7 +226,7 @@ private: ...@@ -219,7 +226,7 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i)) // int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds // *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct) // = 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; -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
//we use the relations //we use the relations
...@@ -230,14 +237,14 @@ private: ...@@ -230,14 +237,14 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i)) // int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds // *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct) // = 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; -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric //stiffness matrix A is symmetric
A[neighborIdxSLE][elemIdxSLE + i + 1] += A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
A[elemIdxSLE + i + 1][neighborIdxSLE]; A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
A[neighborIdxSLE + i + 1][elemIdxSLE] += A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
A[elemIdxSLE][neighborIdxSLE + i + 1]; A->entry(elemIdxSLE, neighborIdxSLE + i + 1);
for (int j = 0; j <= i; j++) for (int j = 0; j <= i; j++)
{ {
...@@ -253,14 +260,14 @@ private: ...@@ -253,14 +260,14 @@ private:
// int_intersct avg(K*grad(phi_elem,i)) // int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_neighbor,j) ds // *jump(phi_neighbor,j) ds
// = -0.5 * K * normal[i] * int_intersct x_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] -mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i] - 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]); - normal[i] * linearIntegrals[j]);
//stiffness matrix A is symmetric //stiffness matrix A is symmetric
A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] += A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1]; A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
if (i != j) if (i != j)
{ {
...@@ -275,14 +282,14 @@ private: ...@@ -275,14 +282,14 @@ private:
// int_intersct avg(K*grad(phi_elem,j)) // int_intersct avg(K*grad(phi_elem,j))
// *jump(phi_neighbor,i) ds // *jump(phi_neighbor,i) ds
// = -0.5 * K * normal[j] * int_intersct x_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] -mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j] - 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]); - normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric //stiffness matrix A is symmetric
A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] += A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1]; A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
} }
} }
} }
...@@ -297,7 +304,7 @@ private: ...@@ -297,7 +304,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i)) // int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds // *jump(phi_elem,0) ds
// = K * normal[i] * vol(intersct) // = 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; mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++) for (int j = 0; j <= i; j++)
...@@ -309,7 +316,7 @@ private: ...@@ -309,7 +316,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i)) // int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds // *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_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] mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j] - 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]); + normal[j] * linearIntegrals[i]);
...@@ -321,22 +328,18 @@ private: ...@@ -321,22 +328,18 @@ private:
//stiffness matrix A is symmetric //stiffness matrix A is symmetric
for (int i = 0; i < dim; i++) 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++) for(int j = 0; j < i; j++)
{ {
A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] = A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]; A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
} }
} }
} }
//NOTE: check if A is symmetric A->compress();
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;
} }
//writes the solution to a vtk file //writes the solution to a vtk file
...@@ -392,7 +395,7 @@ private: ...@@ -392,7 +395,7 @@ private:
} }
Matrix A; //stiffness matrix std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector Vector b; //load vector
Vector d; //solution vector Vector d; //solution vector
}; };
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment