diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index d77c032191a3ad0c5d12fe0fa847071bb9b12f78..03387bd35f0d6335eb673d0a6304619676be3625 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -27,6 +27,7 @@ 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>>;
@@ -79,7 +80,7 @@ public:
}
//compute L2 error using a quadrature rule
- Scalar computeL2error(const int quadratureOrder = 5) const
+ Scalar computeL2error(const int quadratureOrder = 5) //const
{
if (!problem_.hasExactSolution())
{
@@ -96,26 +97,23 @@ public:
const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);
//determine L2-error on element elem using a quadrature rule
- for (const auto& ip : qrule)
- {
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = geo.global(qpLocal);
+ const auto& qlambda =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
+ {
+ Scalar numericalSolutionAtQP = d[elemIdxSLE];
- const Scalar quadratureFactor =
- ip.weight() * geo.integrationElement(qpLocal);
+ for (int i = 0; i < dim; i++)
+ {
+ numericalSolutionAtQP += d[elemIdxSLE + i + 1] * qpGlobal[i];
+ }
- Scalar numericalSolutionAtQP = d[elemIdxSLE];
+ const Scalar errorEvaluation =
+ numericalSolutionAtQP - problem_.exactSolution(qpGlobal);
- for (int i = 0; i < dim; i++)
- {
- numericalSolutionAtQP += d[elemIdxSLE + i + 1] * qpGlobal[i];
- }
-
- const Scalar errorEvaluation =
- numericalSolutionAtQP - problem_.exactSolution(qpGlobal);
+ error += weight * errorEvaluation * errorEvaluation;
+ };
- error += quadratureFactor * errorEvaluation * errorEvaluation;
- }
+ applyQuadratureRule(geo, qrule, qlambda);
}
return std::sqrt(error);
@@ -174,49 +172,42 @@ protected:
//basis function phi_elem,i
const int elemIdxSLE = (dim + 1)*elemIdx;
- //fill load vector b using a quadrature rule
- for (const auto& ip : rule_q)
- {
- //quadrature point in local and global coordinates
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = geo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * geo.integrationElement(qpLocal);
- const Scalar qEvaluation(problem_.q(qpGlobal));
+ //lambda to fill load vector b using a quadrature rule
+ const auto lambda_q =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
+ {
+ const Scalar qEvaluation(problem_.q(qpGlobal));
- //quadrature for int_elem q*phi_elem,0 dV
- b[elemIdxSLE] += quadratureFactor * qEvaluation;
+ //quadrature for int_elem q*phi_elem,0 dV
+ b[elemIdxSLE] += weight * qEvaluation;
- //quadrature for int_elem q*phi_elem,i dV
- for (int i = 0; i < dim; i++)
- {
- b[elemIdxSLE + i + 1] += quadratureFactor * qpGlobal[i] * qEvaluation;
- }
- }
+ //quadrature for int_elem q*phi_elem,i dV
+ for (int i = 0; i < dim; i++)
+ {
+ b[elemIdxSLE + i + 1] += weight * qpGlobal[i] * qEvaluation;
+ }
+ };
- //evaluate
+ //lambda to 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)
- {
- //quadrature point in local and global coordinates
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = geo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * geo.integrationElement(qpLocal);
-
- for (int i = 0; i < dim; i++)
+ const auto lambda_K_Elem =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- for (int j = 0; j <= i; j++)
+ for (int i = 0; i < dim; i++)
{
- (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
- quadratureFactor * problem_.K(qpGlobal)[j][i];
+ for (int j = 0; j <= i; j++)
+ {
+ (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ weight * problem_.K(qpGlobal)[j][i];
+ }
}
- }
- }
+ };
+
+ //apply quadrature rules on elem
+ applyQuadratureRule(geo, rule_q, lambda_q);
+ applyQuadratureRule(geo, rule_K_Elem, lambda_K_Elem);
//iterate over all intersection with the boundary of elem
for (const auto& intersct : intersections(gridView_, elem))
@@ -239,90 +230,72 @@ protected:
//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::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
-
+ //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds //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);
+ Coordinate linearIntegrals(0.0);
+ CoordinateMatrix quadraticIntregrals(0.0);
+ Coordinate KIntegrals(0.0);
+ CoordinateMatrix KIntegralsX(0.0);
- //fill vector linearIntegrals and matrix quadraticIntregrals
- for (const auto& ip : rule_2ndOrder)
+ //use midpoint rule to fill vector linearIntegrals (exact evaluation)
+ for (int i = 0; i < dim; i++)
{
- //quadrature point in local and global coordinates
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = intersctGeo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * intersctGeo.integrationElement(qpLocal);
+ linearIntegrals[i] = intersctVol * intersctCenter[i];
+ }
- for (int i = 0; i < dim; i++)
+ //lambda to fill matrix quadraticIntregrals using a quadrature rule
+ const auto lambda_2ndOrder =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- //use midpoint rule for exact evaluation of
- // int_intersct x[i] ds //TODO
- linearIntegrals[i] = intersctVol * intersctCenter[i];
-
- for (int j = 0; j <= i; j++)
+ for (int i = 0; i < dim; i++)
{
- quadraticIntregrals[i][j] +=
- quadratureFactor * qpGlobal[i] * qpGlobal[j];
+ for (int j = 0; j <= i; j++)
+ {
+ quadraticIntregrals[i][j] += weight * qpGlobal[i] * qpGlobal[j];
}
}
- }
+ };
- //quadraticIntregrals is symmetric
- for (int i = 0; i < dim; i++)
+ //lambda to fill vector KIntegrals using a quadrature rule
+ const auto lambda_K_Edge =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- for (int j = 0; j < i; j++)
+ const CoordinateMatrix KEvaluation = problem_.K(qpGlobal);
+
+ for (int i = 0; i < dim; i++)
{
- quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
+ KIntegrals[i] += weight * (KEvaluation[i] * normal);
}
- }
+ };
- //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 quadratureFactor =
- ip.weight() * intersctGeo.integrationElement(qpLocal);
- const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
- problem_.K(qpGlobal);
-
- for (int i = 0; i < dim; i++)
+ //lambda to fill vector KIntegralsX using a quadrature rule
+ const auto lambda_Kplus1_Edge =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- KIntegrals[i] += quadratureFactor * (KEvaluation[i] * normal);
- }
- }
+ const CoordinateMatrix KEvaluation = problem_.K(qpGlobal);
- //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);
+ for (int i = 0; i < dim; i++)
+ {
+ for (int j = 0; j < dim; j++)
+ {
+ KIntegralsX[i][j] += weight * qpGlobal[j]
+ * (KEvaluation[i] * normal);
+ }
+ }
+ };
- const Scalar quadratureFactor =
- ip.weight() * intersctGeo.integrationElement(qpLocal);
- const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
- problem_.K(qpGlobal);
+ //apply quadrature rules in intersct
+ applyQuadratureRule(intersctGeo, rule_2ndOrder, lambda_2ndOrder);
+ applyQuadratureRule(intersctGeo, rule_K_Edge, lambda_K_Edge);
+ applyQuadratureRule(intersctGeo, rule_Kplus1_Edge, lambda_Kplus1_Edge);
- for (int i = 0; i < dim; i++)
+ //quadraticIntregrals is symmetric
+ for (int i = 0; i < dim; i++)
+ {
+ for (int j = 0; j < i; j++)
{
- for (int j = 0; j < dim; j++)
- {
- KIntegralsX[i][j] += quadratureFactor *
- qpGlobal[j] * (KEvaluation[i] * normal);
- }
+ quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
}
}
@@ -370,8 +343,7 @@ protected:
(*A)[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
//stiffness matrix A is symmetric
- (*A)[neighborIdxSLE][elemIdxSLE] +=
- (*A)[elemIdxSLE][neighborIdxSLE];
+ (*A)[neighborIdxSLE][elemIdxSLE] += (*A)[elemIdxSLE][neighborIdxSLE];
for (int i = 0; i < dim; i++)
{
@@ -442,33 +414,29 @@ protected:
}
else //boundary facet
{
- //fill load vector b using a quadrature rule (boundary terms)
- for (const auto& ip : rule_boundary)
- {
- //quadrature point in local and global coordinates
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = intersctGeo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * intersctGeo.integrationElement(qpLocal);
-
- const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
- const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
- problem_.K(qpGlobal);
+ //lambda to fill load vector b with boundary entries using a
+ //quadrature rule
+ const auto lambda_boundary =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
+ {
+ const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
+ const CoordinateMatrix KEvaluation = problem_.K(qpGlobal);
- //quadrature for int_intersct mu * g * phi_elem,0 ds
- b[elemIdxSLE] += quadratureFactor * mu * boundaryEvaluation;
+ //quadrature for int_intersct mu * g * phi_elem,0 ds
+ b[elemIdxSLE] += weight * mu * boundaryEvaluation;
- //quadrature for
- // int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
- //with K*grad(phi_elem,i) = K[:,i]
- for (int i = 0; i < dim; i++)
- {
- b[elemIdxSLE + i + 1] +=
- quadratureFactor * (mu * qpGlobal[i] -
+ //quadrature for
+ // int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
+ //with K*grad(phi_elem,i) = K[:,i]
+ for (int i = 0; i < dim; i++)
+ {
+ b[elemIdxSLE + i + 1] += weight * (mu * qpGlobal[i] -
(KEvaluation[i] * normal)) * boundaryEvaluation;
- }
- }
+ }
+ };
+
+ //apply quadrature on intersct
+ applyQuadratureRule(intersctGeo, rule_boundary, lambda_boundary);
for (int i = 0; i < dim; i++)
{
@@ -587,6 +555,26 @@ protected:
return eigenvalues.infinity_norm();
}
+ //apply the quadrature rule "rule" on the given geometry as specified in the
+ //lambda function "lambda", that takes the quadrature point as global
+ //coordinate and the quadrature weight (resp. weight * integrationElement) as
+ //arguments
+ void applyQuadratureRule (const auto& geometry, const auto& rule,
+ const auto& lambda)
+ {
+ for (const auto& ip : rule)
+ {
+ //quadrature point in local and global coordinates
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = geometry.global(qpLocal);
+
+ const Scalar quadratureFactor =
+ ip.weight() * geometry.integrationElement(qpLocal);
+
+ lambda(qpGlobal, quadratureFactor);
+ }
+ }
+
std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector