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

[bugfix] stiffness matrix is now symmetric

parent 2d6f739c
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/mmdg/dg.hh
+ 85
− 80
View file @ 11961085
#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
template<class GridView, class Mapper, class Problem>
......@@ -38,6 +45,12 @@ public:
const auto& geo = elem.geometry();
const double elemVol = geo.volume();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
//and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
//basis function phi_elem,i
const int elemIdxSLE = (dim + 1)*elemIdx;
/* //TODO: can be done outside of the loop?
const Dune::QuadratureRule<double,dim>& rule =
Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
......@@ -51,29 +64,28 @@ public:
const auto weight = ip.weight();
//quadrature for int_elem q*phi_elem,0 dV
b[(1 + dim)*elemIdx] += weight * q(qp);
b[elemIdxSLE] += weight * q(qp);
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
}
}
*/
//NOTE: makeshift solution for source term q = -1
b[(1 + dim)*elemIdx] = -elemVol;
b[elemIdxSLE] = -elemVol;
for (int i = 0; i < dim; i++)
{
b[(1 + dim)*elemIdx + i + 1] = -elemVol * geo.center()[i];
b[elemIdxSLE + i + 1] = -elemVol * geo.center()[i];
}
for (int i = 0; i < dim; i++)
{
//exact evaluation of
// int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + i + 1] =
K * elemVol;
A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] = K * elemVol;
}
//iterate over all intersection with the boundary of elem
......@@ -122,18 +134,20 @@ public:
ip.weight() * qp[i] * qp[j] * intersctVol;
}
*/
//NOTE: probably unnecessary
quadraticIntregrals[i][j] = quadraticIntregrals[j][i];
}
}
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] +=
mu * intersctGeo.volume();
A[elemIdxSLE][elemIdxSLE] += mu * intersctGeo.volume();
if (intersct.neighbor()) //intersct has neighboring element
{
//index of the neighboring element
const int neighborIdx = mapper_.index(intersct.outside());
const int neighborIdxSLE = (dim + 1)*neighborIdx;
for (int i = 0; i < dim; i++)
{ //we use the relations
......@@ -143,13 +157,8 @@ public:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
mu * linearIntegrals[i]
- 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric
A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
{
......@@ -160,15 +169,10 @@ public:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1]
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric
A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
+= A[(1 + dim)*elemIdx + i + 1]
[(1 + dim)*elemIdx + j + 1];
}
}
......@@ -179,12 +183,10 @@ public:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] +=
-mu * intersctVol;
A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
//stiffness matrix A is symmetric
A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] +=
A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx];
A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
for (int i = 0; i < dim; i++)
{
......@@ -196,9 +198,8 @@ public:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] +=
-mu * linearIntegrals[i]
- 0.5 * K * normal[i] * intersctVol;
A[elemIdxSLE + i + 1][neighborIdxSLE] +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
//we use the relations
// int_intersct mu*jump(phi_elem,0)
......@@ -208,15 +209,14 @@ public:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] +=
-mu * linearIntegrals[i]
- 0.5 * K * normal[i] * intersctVol;
A[elemIdxSLE][neighborIdxSLE + i + 1] +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric
A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] +=
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx];
A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] +=
A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1];
A[neighborIdxSLE][elemIdxSLE + i + 1] +=
A[elemIdxSLE + i + 1][neighborIdxSLE];
A[neighborIdxSLE + i + 1][elemIdxSLE] +=
A[elemIdxSLE][neighborIdxSLE + i + 1];
for (int j = 0; j <= i; j++)
{
......@@ -232,38 +232,37 @@ public:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_neighbor,j) ds
// = -0.5 * K * normal[i] * int_intersct x_j ds
A[(1 + dim)*elemIdx + i + 1]
[(1 + dim)*neighborIdx + j + 1] +=
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
// int_intersct mu*jump(phi_elem,j)
// *jump(phi_neighbor,i) ds
// = -mu*int_intersct x_i*x_j ds
//and
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
//as well as
// int_intersct avg(K*grad(phi_elem,j))
// *jump(phi_neighbor,i) ds
// = -0.5 * K * normal[j] * int_intersct x_i ds
A[(1 + dim)*elemIdx + j + 1]
[(1 + dim)*neighborIdx + i + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric
A[(1 + dim)*neighborIdx + j + 1]
[(1 + dim)*elemIdx + i + 1] +=
A[(1 + dim)*elemIdx + i + 1]
[(1 + dim)*neighborIdx + j + 1];
A[(1 + dim)*neighborIdx + i + 1]
[(1 + dim)*elemIdx + j + 1] +=
A[(1 + dim)*elemIdx + j + 1]
[(1 + dim)*neighborIdx + i + 1];
A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
if (i != j)
{
// int_intersct mu*jump(phi_elem,j)
// *jump(phi_neighbor,i) ds
// = -mu*int_intersct x_i*x_j ds
//and
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
//as well as
// 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] +=
-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];
}
}
}
}
......@@ -277,13 +276,8 @@ public:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds
// = K * normal[i] * vol(intersct)
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
mu * linearIntegrals[i]
- 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric
A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
{
......@@ -294,20 +288,26 @@ public:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
A[(1 + dim)*elemIdx + i + 1]
[(1 + dim)*elemIdx + j + 1] +=
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric
A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
+= A[(1 + dim)*elemIdx + i + 1]
[(1 + dim)*elemIdx + j + 1];
}
}
}
}
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
A[elemIdxSLE][elemIdxSLE + i + 1] = A[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];
}
}
}
std::cout << A << std::endl;
......@@ -315,6 +315,12 @@ public:
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
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;
//storage for pressure data, for each element we store the
//pressure at the corners of the element, the VTKFunction will
//create the output using a linear interpolation for each element
......@@ -323,7 +329,7 @@ public:
for (const auto& elem : elements(gridView_))
{
const int elemIdx = mapper_.index(elem);
const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
const auto& geo = elem.geometry();
for (int k = 0; k < geo.corners(); k++)
......@@ -331,15 +337,14 @@ public:
//contribution of the basis function
// phi_elem,0 (x) = indicator(elem);
//at the kth corner of elem
pressure[(dim + 1)*elemIdx + k] = d[(1 + dim)*elemIdx];
pressure[elemIdxSLE + k] = d[elemIdxSLE];
for (int i = 0; i < dim; i++)
{
//contribution of the basis function
// phi_elem,i (x) = x[i]*indicator(elem);
//at the kth corner of elem
pressure[(dim + 1)*elemIdx + k] +=
d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i];
pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
}
}
}
......
src/dg.cc
+ 1
− 6
View file @ 11961085
......@@ -10,17 +10,12 @@
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/mmesh/mmesh.hh>
#include <dune/mmdg/dg.hh>
......@@ -38,12 +33,12 @@ int main(int argc, char** argv)
static constexpr int dim = GRIDDIM;
//use YaspGrid if dim = 1 and MMesh otherwise
using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
Dune::MovingMesh<dim>>::type;
using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Coordinate = Dune::FieldVector<double, dim>;
//Message Passing Interface (MPI) for parallel computation
Dune::MPIHelper::instance(argc, argv);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment