diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 014257a2d65ddcffbe86c1f178539fdc4d200348..76def0c56c7a80f6859bd7fe7d646f450560de85 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -27,6 +27,12 @@ public:
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<Scalar>>;
+ using QRuleElem = Dune::QuadratureRule<Scalar, dim>;
+ using QRulesElem = Dune::QuadratureRules<Scalar, dim>;
+ using QRuleEdge = Dune::QuadratureRule<Scalar, dim-1>;
+ using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
+ static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
+ static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
//constructor
DG (const GridView& gridView, const Mapper& mapper,
@@ -40,10 +46,10 @@ public:
d = Vector(dof);
}
- const void operator() (const Scalar K, const Scalar mu)
+ const void operator() (const Scalar mu)
{
//assemble stiffness matrix A and load vector b
- assembleSLE(K, mu);
+ assembleSLE(mu);
Dune::InverseOperatorResult result;
Dune::UMFPack<Matrix> solver(*A);
@@ -61,12 +67,17 @@ public:
private:
//assemble stiffness matrix A and load vector b
- void assembleSLE (const Scalar K, const Scalar mu)
+ void assembleSLE (const Scalar mu)
{
- //quadrature rule
- const Dune::QuadratureRule<Scalar, dim>& rule =
- Dune::QuadratureRules<Scalar, dim>::rule(
- Dune::GeometryTypes::simplex(dim), problem_.quadratureOrder());
+ //quadrature rules
+ const QRuleElem& rule_q =
+ QRulesElem::rule(simplex, problem_.quadratureOrder_q());
+ const QRuleElem& rule_K_Elem =
+ QRulesElem::rule(simplex, problem_.quadratureOrder_K());
+ const QRuleEdge& rule_K_Edge =
+ QRulesEdge::rule(edge, problem_.quadratureOrder_K());
+ const QRuleEdge& rule_Kplus1_Edge =
+ QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1);
//we use the basis
// phi_elem,0 (x) = indicator(elem);
@@ -76,7 +87,6 @@ private:
{
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
- const double elemVol = geo.volume();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
@@ -85,7 +95,7 @@ private:
const int elemIdxSLE = (dim + 1)*elemIdx;
//fill load vector b using a quadrature rule
- for (const auto& ip : rule)
+ for (const auto& ip : rule_q)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
@@ -96,21 +106,37 @@ private:
const Scalar qEvalution(problem_.q(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
- b[elemIdxSLE] += weight * qEvalution * integrationElem;
+ b[elemIdxSLE] += weight * integrationElem * qEvalution;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] +=
- weight * qpGlobal[i] * qEvalution * integrationElem;
+ weight * integrationElem * qpGlobal[i] * qEvalution;
}
}
- for (int i = 0; i < dim; i++)
+ //evaluate
+ // int_elem K*grad(phi_elem,i)*grad(phi_elem,j) dV
+ // = int_elem K[j][i] dV
+ //using a quadrature rule
+ for (const auto& ip : rule_K_Elem)
{
- //exact evaluation of
- // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
- A->entry(elemIdxSLE + i + 1, elemIdxSLE + i + 1) += K * elemVol;
+ //quadrature point in local and global coordinates
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = geo.global(qpLocal);
+
+ const Scalar weight(ip.weight());
+ const Scalar integrationElem(geo.integrationElement(qpLocal));
+
+ for (int i = 0; i < dim; i++)
+ {
+ for (int j = 0; j <= i; j++)
+ {
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ weight * integrationElem * problem_.K(qpGlobal)[j][i];
+ }
+ }
}
//iterate over all intersection with the boundary of elem
@@ -126,17 +152,27 @@ private:
// Dune::QuadratureRules<double,dim-1>::rule(
// intersct.type(), 2);
- //storage for multiply used integral values,
+ //create storage for multiply used integral values
+ //linearIntegrals[i] = int_intersct x[i] ds
+ Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
+
+ //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds
+ //NOTE: is there a better type for a symmetric matrix?
//note that only the lower diagonal and diagonal entries of
//quadraticIntregrals are used
- Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
- //NOTE: is there a better type for a symmetric matrix?
+ //KIntegrals[i] = int_intersct K[:,i]*normal ds
+ Dune::FieldVector<Scalar, dim> KIntegrals(0.0);
+
+ //KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds
+ Dune::FieldMatrix<Scalar, dim, dim> KIntegralsX(0.0);
+
+ //fill vector linearIntegrals and matrix quadraticIntregrals //NOTE
for (int i = 0; i < dim; i++)
{
//use midpoint rule for exact evaluation of
- // int_intersct x_i ds
+ // int_intersct x[i] ds
linearIntegrals[i] = intersctVol * intersctCenter[i];
for (int j = 0; j <= i; j++)
@@ -164,6 +200,47 @@ private:
}
}
+ //fill vector KIntegrals using a quadrature rule
+ for (const auto& ip : rule_K_Edge)
+ {
+ //quadrature point in local and global coordinates
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = intersctGeo.global(qpLocal);
+
+ const Scalar weight(ip.weight());
+ const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+ const Dune::FieldMatrix<Scalar, dim, dim> KEvalution =
+ problem_.K(qpGlobal);
+
+ for (int i = 0; i < dim; i++)
+ {
+ KIntegrals[i] +=
+ weight * integrationElem * (KEvalution[i] * normal);
+ }
+ }
+
+ //fill vector KIntegralsX using a quadrature rule
+ for (const auto& ip : rule_Kplus1_Edge)
+ {
+ //quadrature point in local and global coordinates
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = intersctGeo.global(qpLocal);
+
+ const Scalar weight(ip.weight());
+ const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+ const Dune::FieldMatrix<Scalar, dim, dim> KEvalution =
+ problem_.K(qpGlobal);
+
+ for (int i = 0; i < dim; i++)
+ {
+ for (int j = 0; j < dim; j++)
+ {
+ KIntegralsX[i][j] += weight * integrationElem *
+ qpGlobal[j] * (KEvalution[i] * normal);
+ }
+ }
+ }
+
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
@@ -177,25 +254,24 @@ private:
for (int i = 0; i < dim; i++)
{ //we use the relations
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
- // = mu * int_intersct x_i ds
+ // = mu * int_intersct x[i] ds
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
- // = 0.5 * K * normal[i] * vol(intersct)
+ // = 0.5 * int_intersct K[:,i] * normal ds
A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
- mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
for (int j = 0; j <= i; j++)
{
//we use the relations
// int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
- // = mu * int_intersct x_i*x_j ds
+ // = mu * int_intersct x[i] * x[j] ds
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
- // = 0.5 * K * normal[i] * int_intersct x_j ds
+ // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+= mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]);
+ - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
}
@@ -205,7 +281,7 @@ private:
}
//exact evaluation of
- // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
+ // int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
//stiffness matrix A is symmetric
@@ -215,26 +291,22 @@ private:
for (int i = 0; i < dim; i++)
{
//we use the relations
- // int_intersct mu*jump(phi_elem,i)
- // *jump(phi_neighbor,0) ds
- // = -mu*int_intersct x_i ds
+ // int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,0) ds
+ // = -mu * int_intersct x[i] ds
//and
- // int_intersct avg(K*grad(phi_neighbor,i))
- // *jump(phi_elem,0) ds
- // = 0.5 * K * normal[i] * vol(intersct)
+ // 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) +=
- -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
+ -mu * linearIntegrals[i] + 0.5 * KIntegrals[i];
//we use the relations
- // int_intersct mu*jump(phi_elem,0)
- // *jump(phi_neighbor,i) ds
- // = -mu*int_intersct x_i ds
+ // int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,i) ds
+ // = -mu * int_intersct x[i] ds
//and
- // int_intersct avg(K*grad(phi_neighbor,i))
- // *jump(phi_elem,0) ds
- // = 0.5 * K * normal[i] * vol(intersct)
+ // 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) +=
- -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ -mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
@@ -245,21 +317,17 @@ private:
for (int j = 0; j <= i; j++)
{
//we use the relations
- // int_intersct mu*jump(phi_elem,i)
- // *jump(phi_neighbor,j) ds
- // = -mu*int_intersct x_i*x_j ds
+ // int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,j) ds
+ // = -mu * int_intersct x[i] * x[j] ds
//and
- // int_intersct avg(K*grad(phi_neighbor,j))
- // *jump(phi_elem,i) ds
- // = 0.5 * K * normal[j] * int_intersct x_i ds
+ // int_intersct avg(K*grad(phi_neighbor,j)) * jump(phi_elem,i) ds
+ // = 0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
//as well as
- // int_intersct avg(K*grad(phi_elem,i))
- // *jump(phi_neighbor,j) ds
- // = -0.5 * K * normal[i] * int_intersct x_j ds
+ // 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) +=
-mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[j] * linearIntegrals[i]
- - normal[i] * linearIntegrals[j]);
+ - 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
@@ -267,21 +335,18 @@ private:
if (i != j)
{
- // int_intersct mu*jump(phi_elem,j)
- // *jump(phi_neighbor,i) ds
- // = -mu*int_intersct x_i*x_j ds
+ //we use the relations
+ // int_intersct mu * jump(phi_elem,j) * jump(phi_neighbor,i) ds
+ // = -mu * int_intersct x[i] * x[j] ds
//and
- // int_intersct avg(K*grad(phi_neighbor,i))
- // *jump(phi_elem,j) ds
- // = 0.5 * K * normal[i] * int_intersct x_j ds
+ // int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,j) ds
+ // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
//as well as
- // int_intersct avg(K*grad(phi_elem,j))
- // *jump(phi_neighbor,i) ds
- // = -0.5 * K * normal[j] * int_intersct x_i ds
+ // 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) +=
-mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- - normal[j] * linearIntegrals[i]);
+ - 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
@@ -294,28 +359,25 @@ private:
{
for (int i = 0; i < dim; i++)
{ //we use the relations
- // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
- // = mu * int_intersct x_i ds
+ // int_intersct mu * jump(phi_elem,0) * jump(phi_elem,i) ds
+ // = mu * int_intersct x[i] ds
//and for boundary facets
- // int_intersct avg(K*grad(phi_elem,i))
- // *jump(phi_elem,0) ds
- // = K * normal[i] * vol(intersct)
+ // 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) +=
- mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
for (int j = 0; j <= i; j++)
{
//we use the relations
- // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
- // = mu * int_intersct x_i*x_j ds
+ // int_intersct mu * jump(phi_elem,i) * jump(phi_elem,j) ds
+ // = mu * int_intersct x[i] * x[j] ds
//and for boundary facets
- // int_intersct avg(K*grad(phi_elem,i))
- // *jump(phi_elem,j) ds
- // = 0.5 * K * normal[i] * int_intersct x_j ds
+ // 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) +=
mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]);
+ - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
}
}
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index f8d22f3f88daa7fbf2c1dc2b222e37a26e335017..942c4338083017f5eebbd2c91adcc00e5749db3b 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -36,13 +36,6 @@ class DGProblem
return Scalar(0.0);
}
- //returns the recommended quadrature order to compute an integral
- //over x * q(x)
- virtual int quadratureOrder () const
- {
- return 1;
- }
-
//permeability tensor at position pos
virtual Matrix K (const Vector& pos) const
{
@@ -56,6 +49,21 @@ class DGProblem
//return identity matrix
return permeability;
}
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * q(x)
+ virtual int quadratureOrder_q () const
+ {
+ return 1;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over an entry K[i][j] of the permeability tensor
+ virtual int quadratureOrder_K () const
+ {
+ return 0;
+ }
+
};
#endif
diff --git a/src/dg.cc b/src/dg.cc
index 540ab9c2c19875f40ffa13a335fdbc7703adcf27..8b160fec3ca62e53607e1a442eb1fd2c8e92d804 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -23,8 +23,6 @@
#include <dune/mmdg/problems/poisson.hh>
-//TODO: specify K in problem? exact solution depends on K!
-
int main(int argc, char** argv)
{
try
@@ -48,9 +46,6 @@ int main(int argc, char** argv)
Dune::ParameterTree pt;
Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
- //assume isotropic medium such that permeability tensor reduces to a scalar
- const double K = std::stod(pt["K"]);
-
//assume constant penalty parameter
const double mu = std::stod(pt["mu"]);
@@ -63,8 +58,6 @@ int main(int argc, char** argv)
//element mapper for the leaf grid view
const Mapper mapper(gridView, Dune::mcmgElementLayout());
-// const DGProblem<Dune::FieldVector<double, dim>, double> dummyProblem;
-
DGProblem<Coordinate>* problem; //problem to be solved
//determine problem type
@@ -81,7 +74,7 @@ int main(int argc, char** argv)
const DG dg(gridView, mapper, *problem);
- dg(K, mu);
+ dg(mu);
//print total run time
const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);