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

[bugfix] fix index error in interface vtk output

parent b11006cd
No related branches found
No related tags found
No related merge requests found
......@@ -46,6 +46,8 @@ public:
(*Base::A).solve(Base::d, Base::b);
std::cout << "\n\n" << Base::b << "\n\n" << Base::d << "\n\n";
// Dune::InverseOperatorResult result;
// Dune::UMFPack<Matrix> solver(*Base::A);
// solver.apply(Base::d, Base::b, result);
......@@ -471,8 +473,6 @@ private:
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
centerI * (interfaceUpdate1 - interfaceUpdate2);
std::cout << iElemIdxSLE + i + 1 << "\t" << interfaceUpdate3 << "\t" << centerI * (interfaceUpdate1 - interfaceUpdate2);
for (int j = 0; j < i; i++)
{ //NOTE: only relevant for n=3, requires quadrature rule for n=3
Coordinate KparDotTau_j;
......@@ -593,8 +593,6 @@ private:
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2);
std::cout << iElemIdxSLE + i + 1 << "\t" << iFrame[i] << "\t" << center * iFrame[i] << "\t" << interfaceUpdate3 << "\t" << centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2) << "\n\n";
for (int j = 0; j < i; i++)
{ //NOTE: only relevant for n=3, requires quadrature rule for n=3
Coordinate KparDotTau_j;
......@@ -670,7 +668,8 @@ private:
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdxSLE = bulkDOF + dim * iMapper_.index(iElem);
const int iElemIdxOutput = dim * iMapper_.index(iElem);
const int iElemIdxSLE = bulkDOF + iElemIdxOutput;
const auto& iGeo = iElem.geometry();
//get basis of the interface coordinate system for evaluation of
......@@ -682,21 +681,21 @@ private:
{
if (Base::problem_.hasExactSolution())
{
exactIPressure[iElemIdxSLE + k] =
exactIPressure[iElemIdxOutput + k] =
Base::problem_.exactInterfaceSolution(iGeo.corner(k));
}
//contribution of the basis function
// phi_iElem,0 (x) = indicator(iElem);
//at the kth corner of iElem
iPressure[iElemIdxSLE + k] = Base::d[iElemIdxSLE];
iPressure[iElemIdxOutput + k] = Base::d[iElemIdxSLE];
for (int i = 0; i < dim - 1; i++)
{
//contribution of the basis function
// phi_iElem,i (x) = (x * tau[i]) * indicator(iElem);
//at the kth corner of iElem
iPressure[iElemIdxSLE + k] +=
iPressure[iElemIdxOutput + k] +=
Base::d[iElemIdxSLE + i + 1] * (iGeo.corner(k) * iFrame[i]);
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment