Skip to content
Snippets Groups Projects
Commit fbef589c authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

add method to mmdg to compute L2 error

parent dbffdbb9
No related branches found
No related tags found
No related merge requests found
...@@ -63,6 +63,7 @@ public: ...@@ -63,6 +63,7 @@ public:
writeVTKoutput(); writeVTKoutput();
} }
//compute L2 error using a quadrature rule
Scalar computeL2error(const int quadratureOrder = 5) const Scalar computeL2error(const int quadratureOrder = 5) const
{ {
if (!problem_.hasExactSolution()) if (!problem_.hasExactSolution())
......
...@@ -44,6 +44,51 @@ public: ...@@ -44,6 +44,51 @@ public:
writeVTKoutput(); 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 InterfaceGridView& iGridView_;
const InterfaceMapper& iMapper_; const InterfaceMapper& iMapper_;
......
...@@ -83,7 +83,7 @@ int main(int argc, char** argv) ...@@ -83,7 +83,7 @@ int main(int argc, char** argv)
//print total run time //print total run time
const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
std::cout << "Finished with a total run time of " << timeTotal << std::cout << "Finished with a total run time of " << timeTotal <<
" Seconds.\n"; " Seconds.\nL2-error: " << mmdg.computeL2error() << ".\n";
return 0; return 0;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment