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

add interpolation for output of exact solution

parent 13f48a97
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 74 additions and 59 deletions
dune/mmdg/dg.hh
+ 74
− 59
View file @ a77ce8dd
......@@ -145,6 +145,8 @@ public:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t" << intersctCenter<< ":\n"
<< mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
if (intersct.neighbor()) //intersct has neighboring element
{
......@@ -161,6 +163,10 @@ public:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
<< mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
for (int j = 0; j <= i; j++)
{
......@@ -174,6 +180,11 @@ public:
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
<< mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
......@@ -201,6 +212,10 @@ public:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][neighborIdxSLE] +=
-mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
//we use the relations
// int_intersct mu*jump(phi_elem,0)
......@@ -212,6 +227,10 @@ public:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE][neighborIdxSLE + i + 1] +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
-0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE][elemIdxSLE + i + 1] +=
......@@ -237,6 +256,11 @@ public:
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]) << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
......@@ -259,6 +283,11 @@ public:
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
......@@ -279,6 +308,11 @@ public:
// = K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
<< intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
for (int j = 0; j <= i; j++)
{
......@@ -293,6 +327,11 @@ public:
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
<< intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
}
......@@ -316,17 +355,42 @@ public:
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
std::cout << b << std::endl;
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;
//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;
};
*/
//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
Dune::DynamicVector<Scalar>
pressure((dim + 1)*gridView_.size(0), 0.0);
Dune::DynamicVector<Scalar> pressure((dim + 1)*gridView_.size(0), 0.0);
Dune::DynamicVector<double>
exactPressure((dim + 1)*gridView_.size(0), 0.0);
for (const auto& elem : elements(gridView_))
{
......@@ -335,6 +399,8 @@ public:
for (int k = 0; k < geo.corners(); k++)
{
exactPressure[elemIdxSLE + k] = exactSolution(geo.corner(k));
//contribution of the basis function
// phi_elem,0 (x) = indicator(elem);
//at the kth corner of elem
......@@ -350,71 +416,20 @@ public:
}
}
Dune::DynamicVector<Scalar>
pressureTEST(gridView_.size(0), 0.0);
for (const auto& elem : elements(gridView_))
{
const int elemIdx = mapper_.index(elem);
const int elemIdxSLE = (dim + 1)*elemIdx;
const auto& center = elem.geometry().center();
pressureTEST[elemIdx] = d[elemIdxSLE];
for (int i = 0; i < dim; i++)
{
pressureTEST[elemIdx] += d[elemIdxSLE + i + 1] * center[i];
}
}
//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);
Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
// std::make_shared<Dune::VTKWriter<GridView>>(gridView_);
// Dune::VTKSequenceWriter<GridView>
// vtkSqWriter(vtkWriter, "pressure", "data", "data");
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<double>>;
VTKFunction* vtkOperator =
new P1Function(gridView_, pressure, "pressure");
//NOTE: why does this only work with a sequence writer?
vtkWriter.addVertexData(
std::shared_ptr<const VTKFunction>(vtkOperator) );
vtkWriter.addCellData(exactPressure, "exactPressure");
vtkWriter.addCellData(pressureTEST, "pressureTEST");
vtkWriter.write("linearpressure");
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, pressure, "pressure")));
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, exactPressure, "exactPressure")));
vtkWriter.pwrite("pressureData", "", "data");
}
const GridView& gridView_;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment