From 8a7cc91875c4c1df64c3810fb83355c073756dc7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Wed, 1 Apr 2020 20:37:21 +0200
Subject: [PATCH] back to sparse
---
dune/mmdg/dg.hh | 93 +++++++++--------------
dune/mmdg/mmdg.hh | 170 +++++++++++++++++++++---------------------
src/mmdgAnalysis.py | 6 +-
src/parameterMMDG.ini | 4 +-
4 files changed, 126 insertions(+), 147 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index b3372b7..adfb45e 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -28,11 +28,8 @@ public:
using Scalar = typename GridView::ctype;
using Coordinate = typename Problem::CoordinateType;
using CoordinateMatrix = typename Problem::Matrix;
- // 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 Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
+ using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<Scalar>>;
@@ -46,34 +43,21 @@ 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 mu0)
{
//assemble stiffness matrix A and load vector b
assembleSLE(mu0);
- //A->compress();
+ A->compress();
- (*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);
+ Dune::InverseOperatorResult result;
+ Dune::UMFPack<Matrix> solver(*A);
+ solver.apply(d, b, result);
//write solution to a vtk file
writeVTKoutput();
@@ -130,14 +114,10 @@ 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
@@ -199,7 +179,7 @@ protected:
{
for (int j = 0; j <= i; j++)
{
- (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
weight * problem_.K(qpGlobal)[j][i];
}
}
@@ -301,7 +281,7 @@ protected:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
- (*A)[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
+ A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
@@ -316,7 +296,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)[elemIdxSLE + i + 1][elemIdxSLE] +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
for (int j = 0; j <= i; j++)
@@ -327,7 +307,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)[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+= mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
@@ -340,10 +320,11 @@ protected:
//exact evaluation of
// int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
- (*A)[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
+ A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
//stiffness matrix A is symmetric
- (*A)[neighborIdxSLE][elemIdxSLE] += (*A)[elemIdxSLE][neighborIdxSLE];
+ A->entry(neighborIdxSLE, elemIdxSLE) +=
+ A->entry(elemIdxSLE, neighborIdxSLE);
for (int i = 0; i < dim; i++)
{
@@ -353,7 +334,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)[elemIdxSLE + i + 1][neighborIdxSLE] +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
-mu * linearIntegrals[i] + 0.5 * KIntegrals[i];
//we use the relations
@@ -362,14 +343,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)[elemIdxSLE][neighborIdxSLE + i + 1] +=
+ A->entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
-mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
//stiffness matrix A is symmetric
- (*A)[neighborIdxSLE][elemIdxSLE + i + 1] +=
- (*A)[elemIdxSLE + i + 1][neighborIdxSLE];
- (*A)[neighborIdxSLE + i + 1][elemIdxSLE] +=
- (*A)[elemIdxSLE][neighborIdxSLE + i + 1];
+ A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
+ A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
+ A->entry(elemIdxSLE, neighborIdxSLE + i + 1);
for (int j = 0; j <= i; j++)
{
@@ -382,13 +363,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)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );
//stiffness matrix A is symmetric
- (*A)[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
- (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
+ A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
if (i != j)
{
@@ -401,13 +382,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)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
+ A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );
//stiffness matrix A is symmetric
- (*A)[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
- (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
+ A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
}
}
}
@@ -445,7 +426,7 @@ protected:
//and
// int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
//with K*grad(phi_elem,i) = K[:,i]
- (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - KIntegrals[i];
for (int j = 0; j <= i; j++)
@@ -456,7 +437,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)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
mu * quadraticIntregrals[i][j]
- ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
@@ -467,13 +448,13 @@ protected:
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
- (*A)[elemIdxSLE][elemIdxSLE + i + 1] =
- (*A)[elemIdxSLE + i + 1][elemIdxSLE];
+ A->entry(elemIdxSLE, elemIdxSLE + i + 1) =
+ A->entry(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];
+ A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
}
}
}
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index f58bbf0..73890ca 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -38,13 +38,11 @@ public:
const void operator() (const Scalar mu0, const Scalar xi)
{
assembleSLE(mu0, xi);
- //Base::A->compress();
+ Base::A->compress();
- (*Base::A).solve(Base::d, Base::b);
-
- // Dune::InverseOperatorResult result;
- // Dune::UMFPack<Matrix> solver(*Base::A);
- // solver.apply(Base::d, Base::b, result);
+ Dune::InverseOperatorResult result;
+ Dune::UMFPack<Matrix> solver(*Base::A);
+ solver.apply(Base::d, Base::b, result);
//write solution to a vtk file
writeVTKoutput();
@@ -203,25 +201,25 @@ private:
//quadrature for
// int_iElem beta * phi_iElem,0 * phi_iElem,0 ds
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] += weight * betaEvaluation;
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * betaEvaluation;
//quadrature for
// int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,0) ds
//for elem1, elem2 in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
- (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
- (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
- (*Base::A)[elemOutIdxSLE][elemInIdxSLE] += bulkUpdate2;
+ Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
+ Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
+ Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
+ Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
//quadrature for
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,j) ds
//for elem in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE][iElemIdxSLE] -= couplingUpdate1;
- (*Base::A)[iElemIdxSLE][elemInIdxSLE] -= couplingUpdate1;
- (*Base::A)[elemOutIdxSLE][iElemIdxSLE] -= couplingUpdate1;
- (*Base::A)[iElemIdxSLE][elemOutIdxSLE] -= couplingUpdate1;
+ Base::A->entry(elemInIdxSLE, iElemIdxSLE) -= couplingUpdate1;
+ Base::A->entry(iElemIdxSLE, elemInIdxSLE) -= couplingUpdate1;
+ Base::A->entry(elemOutIdxSLE, iElemIdxSLE) -= couplingUpdate1;
+ Base::A->entry(iElemIdxSLE, elemOutIdxSLE) -= couplingUpdate1;
for (int i = 0; i < dim; i++)
{
@@ -238,16 +236,16 @@ private:
//quadrature for
// int_iElem beta * phi_iElem,0 * phi_iElem,i ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += interfaceUpdate;
- (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += interfaceUpdate;
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate;
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,0) ds
//for elem in {elemIn, elemOut} using a quadrature rule
- (*Base::A)[elemInIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
- (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE] -= couplingUpdate3;
- (*Base::A)[elemOutIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
- (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE] -= couplingUpdate3;
+ Base::A->entry(elemInIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
+ Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE) -= couplingUpdate3;
+ Base::A->entry(elemOutIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
+ Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE) -= couplingUpdate3;
}
//quadrature for
@@ -255,22 +253,22 @@ private:
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut}
- (*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)[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;
+ 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->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;
//quadrature for
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
- (*Base::A)[iElemIdxSLE][elemInIdxSLE + i + 1] -= couplingUpdate2;
- (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
- (*Base::A)[iElemIdxSLE][elemOutIdxSLE + i + 1] -= couplingUpdate2;
+ Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
+ Base::A->entry(iElemIdxSLE, elemInIdxSLE + i + 1) -= couplingUpdate2;
+ Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
+ Base::A->entry(iElemIdxSLE, elemOutIdxSLE + i + 1) -= couplingUpdate2;
}
};
@@ -311,19 +309,19 @@ private:
{
//quadrature for
// int_iElem beta * (phi_iElem,i)^2 ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate1 * qpI;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + i + 1] -=
+ Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE + i + 1) -=
couplingUpdate2;
- (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE + i + 1] -=
+ Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE + i + 1) -=
couplingUpdate2;
- (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + i + 1] -=
+ Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE + i + 1) -=
couplingUpdate2;
- (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE + i + 1] -=
+ Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE + i + 1) -=
couplingUpdate2;
}
@@ -332,13 +330,13 @@ private:
//and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
+ Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
bulkUpdate3;
- (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
+ Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate3;
- (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
+ Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate4;
- (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
+ Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
bulkUpdate4;
for (int j = 0; j < i; j++)
@@ -359,21 +357,21 @@ private:
//quadrature for
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate2;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate2;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
//for elem in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + j + 1][iElemIdxSLE + i + 1] -=
+ Base::A->entry(elemInIdxSLE + j + 1, iElemIdxSLE + i + 1) -=
couplingUpdate4;
- (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE + j + 1] -=
+ Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE + j + 1) -=
couplingUpdate4;
- (*Base::A)[elemOutIdxSLE + j + 1][iElemIdxSLE + i + 1] -=
+ Base::A->entry(elemOutIdxSLE + j + 1, iElemIdxSLE + i + 1) -=
couplingUpdate4;
- (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE + j + 1] -=
+ Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE + j + 1) -=
couplingUpdate4;
}
@@ -382,33 +380,33 @@ private:
//and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds
//for elem1, elem2 in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
+ Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
bulkUpdate5;
- (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
+ Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
bulkUpdate5;
- (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
+ Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
bulkUpdate5;
- (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
+ Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate5;
- (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
+ Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
bulkUpdate6;
- (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
+ Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
bulkUpdate6;
- (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
+ Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
bulkUpdate6;
- (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
+ Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate6;
//quadrature for
// -int_iElem beta * phi_iElem,j * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
- (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE + j + 1) -=
couplingUpdate3;
- (*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -=
+ Base::A->entry(iElemIdxSLE + j + 1, elemInIdxSLE + i + 1) -=
couplingUpdate3;
- (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE + j + 1) -=
couplingUpdate3;
- (*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -=
+ Base::A->entry(iElemIdxSLE + j + 1, elemOutIdxSLE + i + 1) -=
couplingUpdate3;
}
}
@@ -457,7 +455,7 @@ private:
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,0)) * grad(d * phi,iElem,0) ds
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] += weight * interfaceUpdate1;
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * interfaceUpdate1;
for (int i = 0; i < dim - 1 ; i++)
{
@@ -472,13 +470,13 @@ private:
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,0) ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += interfaceUpdate3;
- (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += interfaceUpdate3;
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate3;
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate3;
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,i) ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i])
+ qpI * (qpI * interfaceUpdate1 + 2.0 * interfaceUpdate2) );
@@ -493,9 +491,9 @@ private:
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,j) ds
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
}
}
@@ -548,7 +546,7 @@ private:
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
@@ -568,9 +566,9 @@ private:
//and
// -int_intersct avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_iElem,0) dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate4;
- (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
//quadrature for
@@ -578,7 +576,7 @@ private:
//and
// -int_intersct avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_iElem,i) dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * dEvaluation * qpI * (interfaceUpdate2
- interfaceUpdate3 - qpI * interfaceUpdate1);
@@ -596,9 +594,9 @@ private:
//and
// -int_intersct avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_iElem,j) dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate5;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
@@ -620,9 +618,9 @@ private:
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_neighbor,0) dr
- (*Base::A)[iElemIdxSLE][neighborIdxSLE] -=
+ Base::A->entry(iElemIdxSLE, neighborIdxSLE) -=
weight * mu * dEvaluation2;
- (*Base::A)[neighborIdxSLE][iElemIdxSLE] -=
+ Base::A->entry(neighborIdxSLE, iElemIdxSLE) -=
weight * mu * dEvaluation2;
for (int i = 0; i < dim - 1; i++)
@@ -647,9 +645,9 @@ private:
//and
// -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_neighbor,0) dr
- (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
+ Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
interfaceUpdate3a;
- (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
+ Base::A->entry(neighborIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate3a;
//quadrature for
@@ -658,9 +656,9 @@ private:
//and.
// -int_intersct d * avg(Kparallel * grad(d * phi_neighbor,i))
// * jump(phi_iElem,0) dr
- (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
+ Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate3b;
- (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE, neighborIdxSLE + i + 1) +=
interfaceUpdate3b;
for (int j = 0; j < dim - 1; j++)
@@ -683,9 +681,9 @@ private:
//resp.
// -int_intersct d * avg(Kparallel * grad(d * phi_neighbor,i))
// * jump(phi_iElem,j) dr
- (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
interfaceUpdate5;
- (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(neighborIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
@@ -721,7 +719,7 @@ private:
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
@@ -741,9 +739,9 @@ private:
//and
// -int_intersct d * (Kparallel * grad(d * phi_iElem,i))
// * phi_iElem,0 * normal dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate4;
- (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
//quadrature for
@@ -751,7 +749,7 @@ private:
//and
// -int_intersct d * (Kparallel * grad(d * phi_iElem,i))
// * phi_iElem,i * normal dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * qpI * dEvaluation *
(interfaceUpdate2 - 2.0 * interfaceUpdate3
- qpI * interfaceUpdate1);
@@ -771,9 +769,9 @@ private:
//and
// -int_intersct d * (Kparallel * grad(d * phi_iElem,i))
// * phi_iElem,j * normal dr
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate5;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
diff --git a/src/mmdgAnalysis.py b/src/mmdgAnalysis.py
index e8312ed..c5b972e 100644
--- a/src/mmdgAnalysis.py
+++ b/src/mmdgAnalysis.py
@@ -61,9 +61,9 @@ if not exists("plots"):
numberOfMeans = 1
# names of the grid files sorted by refinement
-gridfiles = ["mmdg2_1", "mmdg2_5e-1", "mmdg2_2e-1", "mmdg2_8e-2", "mmdg2_5e-2"]
-# gridfiles = ["mmdg3_1", "mmdg3_5e-1", "mmdg3_2e-1", "mmdg3_8e-2", "mmdg3_5e-2"]
-# gridfiles = ["mmdg5_1", "mmdg5_5e-1", "mmdg5_2e-1", "mmdg5_8e-2", "mmdg5_5e-2"]
+gridfiles = ["mmdg2_1", "mmdg2_5e-1", "mmdg2_2e-1", "mmdg2_8e-2", "mmdg2_5e-2", "mmdg2_3e-2", "mmdg2_2e-2", "mmdg2_1e-2", "mmdg2_9e-3", "mmdg2_8e-3"]
+# gridfiles = ["mmdg3_1", "mmdg3_5e-1", "mmdg3_2e-1", "mmdg3_8e-2", "mmdg3_5e-2", "mmdg3_3e-2", "mmdg3_2e-2", "mmdg3_1e-2", "mmdg3_9e-3", "mmdg3_8e-3"]
+# gridfiles = ["mmdg5_1", "mmdg5_5e-1", "mmdg5_2e-1", "mmdg5_8e-2", "mmdg5_5e-2", "mmdg5_3e-2", "mmdg5_2e-2", "mmdg5_1e-2", "mmdg5_9e-3", "mmdg5_8e-3"]
dim = 2
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
index 69a5f03..03de01e 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -1,5 +1,5 @@
mu0 = 1000
-xi = 0.75
-d = 1
+#xi = 0.75
+d = 0.01
problem = mmdg2
gridfile = mmdg2_5e-2
--
GitLab