From 304a0bf3b897d2ef4586d7a0d8133a7261716759 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sun, 22 Mar 2020 12:42:49 +0100
Subject: [PATCH] adapt program to non-sparse matrix (debugging), in mmdg.hh:
remove parameter alpha, [bugfix] fix index error in method interfaceFrame
---
dune/mmdg/dg.hh | 95 ++++++++++++++---------
dune/mmdg/mmdg.hh | 191 ++++++++++++++++++++++++----------------------
src/dg.cc | 8 +-
3 files changed, 164 insertions(+), 130 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index f0f3495..344c397 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -26,8 +26,12 @@ public:
static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
using Scalar = typename GridView::ctype;
- using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
- using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+ // using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
+ // using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+
+ using Matrix = Dune::DynamicMatrix<Scalar>;
+ using Vector = Dune::DynamicVector<Scalar>;
+
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<Scalar>>;
@@ -41,21 +45,34 @@ public:
gridView_(gridView), mapper_(mapper), problem_(problem),
dof_((1 + dim) * gridView.size(0))
{
+
+ A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
+ b = Vector(dof_, 0.0); //load vector
+ d = Vector(dof_, 0.0); //solution vector
+
//initialize stiffness matrix A, load vector b and solution vector d
- A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
- b = Vector(dof_);
- d = Vector(dof_);
+ // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ // b = Vector(dof_);
+ // d = Vector(dof_);
}
void operator() (const Scalar mu)
{
//assemble stiffness matrix A and load vector b
assembleSLE(mu);
- A->compress();
+ //A->compress();
- Dune::InverseOperatorResult result;
- Dune::UMFPack<Matrix> solver(*A);
- solver.apply(d, b, result);
+ (*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;
+
+ // Dune::InverseOperatorResult result;
+ // Dune::UMFPack<Matrix> solver(*A);
+ // solver.apply(d, b, result);
//write solution to a vtk file
writeVTKoutput();
@@ -115,10 +132,14 @@ protected:
const int dof) :
gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
{
+ A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
+ b = Vector(dof_, 0.0); //load vector
+ d = Vector(dof_, 0.0); //solution vector
+
//initialize stiffness matrix A, load vector b and solution vector d
- A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
- b = Vector(dof_);
- d = Vector(dof_);
+ // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ // b = Vector(dof_);
+ // d = Vector(dof_);
}
//assemble stiffness matrix A and load vector b
@@ -190,7 +211,7 @@ protected:
{
for (int j = 0; j <= i; j++)
{
- A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
quadratureFactor * problem_.K(qpGlobal)[j][i];
}
}
@@ -303,7 +324,7 @@ protected:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
- A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
+ (*A)[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
@@ -318,7 +339,7 @@ protected:
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
// = 0.5 * int_intersct K[:,i] * normal ds
- A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
for (int j = 0; j <= i; j++)
@@ -329,7 +350,7 @@ protected:
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
// = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
- A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+= mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
@@ -342,11 +363,11 @@ protected:
//exact evaluation of
// int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
- A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
+ (*A)[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
//stiffness matrix A is symmetric
- A->entry(neighborIdxSLE, elemIdxSLE) +=
- A->entry(elemIdxSLE, neighborIdxSLE);
+ (*A)[neighborIdxSLE][elemIdxSLE] +=
+ (*A)[elemIdxSLE][neighborIdxSLE];
for (int i = 0; i < dim; i++)
{
@@ -356,7 +377,7 @@ protected:
//and
// int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
// = -0.5 * int_intersct K[:,i] * normal ds
- A->entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
+ (*A)[elemIdxSLE + i + 1][neighborIdxSLE] +=
-mu * linearIntegrals[i] + 0.5 * KIntegrals[i];
//we use the relations
@@ -365,14 +386,14 @@ protected:
//and
// int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
// = 0.5 * int_intersct K[:,i] * normal ds
- A->entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
+ (*A)[elemIdxSLE][neighborIdxSLE + i + 1] +=
-mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
//stiffness matrix A is symmetric
- A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
- A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
- A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
- A->entry(elemIdxSLE, neighborIdxSLE + 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++)
{
@@ -385,13 +406,13 @@ protected:
//as well as
// int_intersct avg(K*grad(phi_elem,i)) * jump(phi_neighbor,j) ds
// = -0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
- A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+ (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );
//stiffness matrix A is symmetric
- A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
- A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
+ (*A)[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
+ (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
if (i != j)
{
@@ -404,13 +425,13 @@ protected:
//as well as
// int_intersct avg(K*grad(phi_elem,j))*jump(phi_neighbor,i) ds
// = -0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
- A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+ (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );
//stiffness matrix A is symmetric
- A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
- A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
+ (*A)[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
}
}
}
@@ -452,7 +473,7 @@ protected:
//and
// int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
//with K*grad(phi_elem,i) = K[:,i]
- A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - KIntegrals[i];
for (int j = 0; j <= i; j++)
@@ -463,7 +484,7 @@ protected:
// int_intersct (phi_elem,i * K[:,j]
// + phi_elem,j * K[:,i]) * normal ds
//with K*grad(phi_elem,i) = K[:,i]
- A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
mu * quadraticIntregrals[i][j]
- ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
@@ -474,13 +495,13 @@ protected:
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
- A->entry(elemIdxSLE, elemIdxSLE + i + 1) =
- A->entry(elemIdxSLE + i + 1, elemIdxSLE);
+ (*A)[elemIdxSLE][elemIdxSLE + i + 1] =
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE];
for(int j = 0; j < i; j++)
{
- A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
- A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
+ (*A)[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
}
}
}
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index ccaef4d..c830287 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -34,11 +34,21 @@ public:
const void operator() (const Scalar mu, const Scalar xi)
{
assembleSLE(mu, xi);
- Base::A->compress();
+ //Base::A->compress();
- Dune::InverseOperatorResult result;
- Dune::UMFPack<Matrix> solver(*Base::A);
- solver.apply(Base::d, Base::b, result);
+ for (int i = 0; i < Base::dof_; i++)
+ for (int j = 0; j < i; j++)
+ /*assert*/if(std::abs((*Base::A)[i][j] - (*Base::A)[j][i]) >
+ std::numeric_limits<Scalar>::epsilon())
+ std::cout << i << ", " << j << std::endl;
+
+ std::cout << *Base::A << std::endl;
+
+ (*Base::A).solve(Base::d, Base::b);
+
+ // Dune::InverseOperatorResult result;
+ // Dune::UMFPack<Matrix> solver(*Base::A);
+ // solver.apply(Base::d, Base::b, result);
//write solution to a vtk file
writeVTKoutput();
@@ -96,15 +106,10 @@ private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu, const Scalar xi)
{
- //modell parameters //NOTE: what should be the capture of alpha?
- const auto alpha = [this](const Coordinate& pos) -> Scalar
- {
- return Base::problem_.Kperp(pos) / Base::problem_.aperture(pos);
- };
-
- const auto beta = [&alpha, &xi](const Coordinate& pos) -> Scalar
+ //modell parameter
+ const auto beta = [this, &xi](const Coordinate& pos) -> Scalar
{
- return 4 * alpha(pos) / (2 * xi - 1);
+ return 4 * Base::problem_.Kperp(pos) / (2 * xi - 1);
};
//quadrature rules
@@ -157,7 +162,7 @@ private:
//in the total system of linear equations (SLE) Ad = b,
//the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
//and the indices iElemIdxSLE + i + 1, i = 1,...,dim-1, refer to the
- //basis functions (0, phi_iElem,i)
+ //basis functions (0, phi_iElem,i)interfaceinterface
const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
// === coupling entries ===
@@ -174,14 +179,14 @@ private:
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar betaEvaluation = beta(qpGlobal);
- Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
quadratureFactor * betaEvaluation;
for (int i = 0; i < dim - 1; i++)
{
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
- Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
}
}
@@ -202,13 +207,13 @@ private:
{
for (int j = 0; j <= i; j++)
{
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
* (iFrame[j] * qpGlobal);
if (i != j)
{
- Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
* (iFrame[j] * qpGlobal);
}
@@ -224,7 +229,7 @@ private:
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
//evalute the coupling terms
- // int_iElem alpha * jump(phi_elem1,0) * jump(phi_elem2,i) ds
+ // int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut} and i = 0,...,dim
@@ -236,39 +241,39 @@ private:
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar alphaEvaluation = alpha(qpGlobal);
+ const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 =
- quadratureFactor * (alphaEvaluation + betaEvaluation4);
+ quadratureFactor * (KEvaluation + betaEvaluation4);
const Scalar bulkUpdate2 =
- quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+ quadratureFactor * (-KEvaluation + betaEvaluation4);
- Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
- Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
+ (*Base::A)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
+ (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
- Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
- Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
+ (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
+ (*Base::A)[elemOutIdxSLE][elemInIdxSLE] += bulkUpdate2;
for (int i = 0; i < dim; i++)
{
const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
- Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate3;
- Base::A->entry(elemInIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate3;
- Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate3;
- Base::A->entry(elemOutIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate3;
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate3;
+ (*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3;
+ (*Base::A)[elemOutIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate3;
- Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate4;
- Base::A->entry(elemOutIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate4;
- Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate4;
- Base::A->entry(elemInIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate4;
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4;
+ (*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4;
+ (*Base::A)[elemInIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate4;
}
}
//evalute the coupling terms
- // int_iElem alpha * jump(phi_elem1,i) * jump(phi_elem2,j) ds
+ // int_iElem K_perp * jump(phi_elem1,i) * jump(phi_elem2,j) ds
//and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds
//for elem1, elem2 in {elemIn, elemOut} and i,j = 1,...,dim
@@ -280,27 +285,27 @@ private:
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar alphaEvaluation = alpha(qpGlobal);
+ const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 =
- quadratureFactor * (alphaEvaluation + betaEvaluation4);
+ quadratureFactor * (KEvaluation + betaEvaluation4);
const Scalar bulkUpdate2 =
- quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+ quadratureFactor * (-KEvaluation + betaEvaluation4);
for (int i = 0; i < dim; i++)
{
const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
- Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
bulkUpdate3;
- Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate3;
- Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate4;
- Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
bulkUpdate4;
for (int j = 0; j < i; j++)
@@ -308,22 +313,22 @@ private:
const Scalar bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
- Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
bulkUpdate5;
- Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+ (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
bulkUpdate5;
- Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
bulkUpdate5;
- Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+ (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate5;
- Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
bulkUpdate6;
- Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+ (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
bulkUpdate6;
- Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
bulkUpdate6;
- Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+ (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate6;
}
}
@@ -345,10 +350,10 @@ private:
const Scalar couplingUpdate1 = quadratureFactor * betaEvaluation2;
- Base::A->entry(elemInIdxSLE, iElemIdx) -= couplingUpdate1;
- Base::A->entry(iElemIdx, elemInIdxSLE) -= couplingUpdate1;
- Base::A->entry(elemOutIdxSLE, iElemIdx) -= couplingUpdate1;
- Base::A->entry(iElemIdx, elemOutIdxSLE) -= couplingUpdate1;
+ (*Base::A)[elemInIdxSLE][iElemIdx] -= couplingUpdate1;
+ (*Base::A)[iElemIdx][elemInIdxSLE] -= couplingUpdate1;
+ (*Base::A)[elemOutIdxSLE][iElemIdx] -= couplingUpdate1;
+ (*Base::A)[iElemIdx][elemOutIdxSLE] -= couplingUpdate1;
for (int i = 0; i < dim - 1; i++)
{
@@ -356,23 +361,23 @@ private:
const Scalar couplingUpdate3 =
(iFrame[i] * qpGlobal) * couplingUpdate1;
- Base::A->entry(elemInIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
- Base::A->entry(iElemIdx, elemInIdxSLE + i + 1) -= couplingUpdate2;
- Base::A->entry(elemOutIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
- Base::A->entry(iElemIdx, elemOutIdxSLE + i + 1) -= couplingUpdate2;
+ (*Base::A)[elemInIdxSLE + i + 1][iElemIdx] -= couplingUpdate2;
+ (*Base::A)[iElemIdx][elemInIdxSLE + i + 1] -= couplingUpdate2;
+ (*Base::A)[elemOutIdxSLE + i + 1][iElemIdx] -= couplingUpdate2;
+ (*Base::A)[iElemIdx][elemOutIdxSLE + i + 1] -= couplingUpdate2;
- Base::A->entry(elemInIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
- Base::A->entry(iElemIdx + i + 1, elemInIdxSLE) -= couplingUpdate3;
- Base::A->entry(elemOutIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
- Base::A->entry(iElemIdx + i + 1, elemOutIdxSLE) -= couplingUpdate3;
+ (*Base::A)[elemInIdxSLE][iElemIdx + i + 1] -= couplingUpdate3;
+ (*Base::A)[iElemIdx + i + 1][elemInIdxSLE] -= couplingUpdate3;
+ (*Base::A)[elemOutIdxSLE][iElemIdx + i + 1] -= couplingUpdate3;
+ (*Base::A)[iElemIdx + i + 1][elemOutIdxSLE] -= couplingUpdate3;
}
//i = dim - 1
const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
- Base::A->entry(elemInIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
- Base::A->entry(iElemIdx, elemInIdxSLE + dim + 1) -= couplingUpdate4;
- Base::A->entry(elemOutIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
- Base::A->entry(iElemIdx, elemOutIdxSLE + dim + 1) -= couplingUpdate4;
+ (*Base::A)[elemInIdxSLE + dim + 1][iElemIdx] -= couplingUpdate4;
+ (*Base::A)[iElemIdx][elemInIdxSLE + dim + 1] -= couplingUpdate4;
+ (*Base::A)[elemOutIdxSLE + dim + 1][iElemIdx] -= couplingUpdate4;
+ (*Base::A)[iElemIdx][elemOutIdxSLE + dim + 1] -= couplingUpdate4;
}
// === interface entries ===
@@ -418,12 +423,12 @@ private:
for (int j = 0; j <= i; j++)
{
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
quadratureFactor * ( KparDotTau_i * iFrame[j] );
if (i != j)
{
- Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
quadratureFactor * ( KparDotTau_i * iFrame[j] );
}
}
@@ -440,7 +445,7 @@ private:
//evaluate the interface term
// int_intersct mu^Gamma * jump(phi_iElem,0) * jump(phi_iElem,0) dr
- Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu;
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu;
if (intersct.neighbor()) //inner interface edge
{
@@ -459,13 +464,15 @@ private:
const Scalar interfaceUpdate3 =
interfaceUpdate1 - 0.5 * interfaceUpdate2;
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3;
- Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
+ (*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;
@@ -475,9 +482,9 @@ private:
- 0.5 * (interfaceUpdate2 * (center * iFrame[j])
+ (KparDotTau_j * normal) * centerI);
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate4;
- Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
interfaceUpdate4;
}
}
@@ -502,8 +509,8 @@ private:
// -int_intersct avg(Kparallel * grad(phi_neighbor,i))
// * jump(phi_iElem,j) dr
//for i,j = 0,...,dim-1
- Base::A->entry(iElemIdxSLE, neighborIdxSLE) += -mu;
- Base::A->entry(neighborIdxSLE, iElemIdxSLE) += -mu;
+ (*Base::A)[iElemIdxSLE][neighborIdxSLE] += -mu;
+ (*Base::A)[neighborIdxSLE][iElemIdxSLE] += -mu;
for (int i = 0; i < dim - 1; i++)
{
@@ -516,16 +523,16 @@ private:
const Scalar interfaceUpdate4 =
-interfaceUpdate1 - 0.5 * interfaceUpdate2;
- Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
+ (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
interfaceUpdate3;
- Base::A->entry(neighborIdxSLE, iElemIdxSLE + i + 1) +=
+ (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
- Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE) +=
+ (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
-interfaceUpdate4;
- Base::A->entry(iElemIdxSLE, neighborIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
-interfaceUpdate4;
- Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + i + 1] +=
-centerI * interfaceUpdate1;
for (int j = 0; j < i; j++)
@@ -541,13 +548,13 @@ private:
const Scalar interfaceUpdate8 =
-interfaceUpdate5 + 0.5 * interfaceUpdate6;
- Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
interfaceUpdate7;
- Base::A->entry(neighborIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
interfaceUpdate7;
- Base::A->entry(iElemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
interfaceUpdate8;
- Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate8;
}
}
@@ -579,13 +586,15 @@ private:
Base::b[iElemIdxSLE] += dgGamma * (mu * centerI - interfaceUpdate2);
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3;
- Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
+ (*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;
@@ -595,9 +604,9 @@ private:
- (interfaceUpdate2 * (center * iFrame[j])
+ (KparDotTau_j * normal) * centerI);
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate4;
- Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
interfaceUpdate4;
}
}
@@ -620,7 +629,7 @@ private:
for (int i = 0; i < dim - 1; i++)
{
Coordinate tau(0.0);
- for (int j = 0; j < dim - 1; j++)
+ for (int j = 0; j < dim; j++)
{
if (j == i)
{
diff --git a/src/dg.cc b/src/dg.cc
index 744b48c..d2bba7c 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -39,7 +39,9 @@ int main(int argc, char** argv)
//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 GridFactory = Dune::DGFGridFactory<Grid>;
+ using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
+
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Coordinate = Dune::FieldVector<double, dim>;
@@ -54,7 +56,9 @@ int main(int argc, char** argv)
const double mu = std::stod(pt["mu"]);
//create a grid from .dgf file
- GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
+ //GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
+ GridFactory gridFactory( "grids/mmdg2.msh" );
+
const Grid& grid = *gridFactory.grid();
const GridView& gridView = grid.leafGridView();
--
GitLab