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

correct a wrong sign in the assemblation of the stiffness matrix

parent 11961085
No related branches found
No related tags found
No related merge requests found
#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH
#include <cmath>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
......@@ -44,6 +46,7 @@ public:
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
const double elemVol = geo.volume();
const auto& center = geo.center();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
......@@ -75,17 +78,17 @@ public:
*/
//NOTE: makeshift solution for source term q = -1
b[elemIdxSLE] = -elemVol;
b[elemIdxSLE] += -elemVol;
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] = -elemVol * geo.center()[i];
b[elemIdxSLE + i + 1] += -elemVol * 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[elemIdxSLE + i + 1][elemIdxSLE + i + 1] = K * elemVol;
A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
}
//iterate over all intersection with the boundary of elem
......@@ -141,7 +144,7 @@ public:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A[elemIdxSLE][elemIdxSLE] += mu * intersctGeo.volume();
A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
......@@ -154,8 +157,7 @@ public:
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
// = mu * int_intersct x_i ds
//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)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
......@@ -166,8 +168,7 @@ public:
// int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
// = mu * int_intersct x_i*x_j ds
//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
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+= mu * quadraticIntregrals[i][j]
......@@ -199,7 +200,7 @@ public:
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
A[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
// int_intersct mu*jump(phi_elem,0)
......@@ -310,8 +311,6 @@ public:
}
}
std::cout << A << std::endl;
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
......@@ -349,6 +348,34 @@ public:
}
}
//NOTE: dim = 1
auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
//NOTE: dim = 2
/* auto exactSolution = [](const auto& pos)
{
double solution = 0.0;
for (int m = 1; m <= 21; m = m + 2)
{
for (int n = 1; n <= 21; n = n + 2)
{
solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1])
/ (m * n * (m*m + n*n));
}
}
solution *= -16 / pow(M_PI, 4);
return solution;
};
*/
Dune::DynamicVector<double> exactPressure(gridView_.size(0));
for (const auto& elem : elements(gridView_))
{
exactPressure[mapper_.index(elem)] =
exactSolution(elem.geometry().center());
}
//create output directory if necessary
mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
......@@ -366,6 +393,7 @@ public:
//NOTE: why does this only work with a sequence writer?
vtkSqWriter.addVertexData(
std::shared_ptr<const VTKFunction>(vtkOperator) );
vtkSqWriter.addCellData(exactPressure, "exactPressure");
vtkSqWriter.write(0.0);
}
......
......@@ -3,7 +3,7 @@ DGF
Interval
0
1
20
100
#
Simplex
......
K = 1
mu = 3
mu = 1000
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment