diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 4390d219c296afdb1153369168c20d258de52138..f352e238227638252af237c20fcb6b7ab31c65e3 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -63,6 +63,7 @@ public: writeVTKoutput(); } + //compute L2 error using a quadrature rule Scalar computeL2error(const int quadratureOrder = 5) const { if (!problem_.hasExactSolution()) diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index 69a48fe128f643eb37bd60d7ddf86d33f97e46f9..62d2d0b0c325eeec680f291656c60b170ff1c9e5 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -44,6 +44,51 @@ public: writeVTKoutput(); } + //compute L2 error using quadrature rules + // norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) ) + Scalar computeL2error(const int quadratureOrder = 5) const + { + const Scalar bulkError = Base::computeL2error(quadratureOrder); + Scalar interfaceError(0.0); + + for (const auto& iElem : elements(iGridView_)) + { + const auto& iGeo = iElem.geometry(); + const int iElemIdxSLE = + (dim + 1) * Base::gridView_.size(0) + dim * iMapper_.index(iElem); + const QRuleInterface& qrule = + QRulesInterface::rule(iSimplex, quadratureOrder); + const InterfaceBasis iFrame = interfaceFrame( + Base::grid_.asIntersection(iElem).centerUnitOuterNormal()); + + //determine L2-error on the interface element iElem using a + //quadrature rule + for (const auto& ip : qrule) + { + const auto& qpLocal = ip.position(); + const auto& qpGlobal = iGeo.global(qpLocal); + + const Scalar quadratureFactor = + ip.weight() * iGeo.integrationElement(qpLocal); + + Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE]; + + for (int i = 0; i < dim - 1; i++) + { + numericalSolutionAtQP += + Base::d[iElemIdxSLE + i + 1] * (qpGlobal * iFrame[i]); + } + + const Scalar errorEvaluation = numericalSolutionAtQP + - Base::problem_.exactInterfaceSolution(qpGlobal); + + interfaceError += quadratureFactor * errorEvaluation * errorEvaluation; + } + } + + return sqrt(bulkError * bulkError + interfaceError); + } + const InterfaceGridView& iGridView_; const InterfaceMapper& iMapper_; diff --git a/src/mmdg.cc b/src/mmdg.cc index 4d567e270cc203f4a28909349221f5cdf073e78a..e2d4348fb89b91e95047dc57c48bdc8dd7cbf06a 100644 --- a/src/mmdg.cc +++ b/src/mmdg.cc @@ -83,7 +83,7 @@ int main(int argc, char** argv) //print total run time const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); std::cout << "Finished with a total run time of " << timeTotal << - " Seconds.\n"; + " Seconds.\nL2-error: " << mmdg.computeL2error() << ".\n"; return 0; }