From 15aac20226b7208b9281ef7a01dbdf3d7a42ba81 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 11 Feb 2016 10:52:39 +0100
Subject: [PATCH] Refactor code to make it more flexible
---
src/compute-disc-error.cc | 231 ++++++++++++++++++++++++--------------
1 file changed, 144 insertions(+), 87 deletions(-)
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index 47b9ab09..4ccd6b9f 100644
--- a/src/compute-disc-error.cc
+++ b/src/compute-disc-error.cc
@@ -39,9 +39,9 @@ const int blocksize = TargetSpace::TangentVector::dimension;
using namespace Dune;
template <class GridView, int order>
-void measureEOC(const GridView gridView,
- const GridView referenceGridView,
- const ParameterTree& parameterSet)
+void measureDiscreteEOC(const GridView gridView,
+ const GridView referenceGridView,
+ const ParameterTree& parameterSet)
{
typedef std::vector<TargetSpace> SolutionType;
@@ -51,6 +51,11 @@ void measureEOC(const GridView gridView,
typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
+ FEBasis referenceFEBasis(referenceGridView);
+
+ using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
+ FufemFEBasis fufemReferenceFEBasis(referenceFEBasis);
+ FufemFEBasis fufemFEBasis(feBasis);
//////////////////////////////////////////////////////////////////////////////////
// Read the data whose error is to be measured
@@ -70,108 +75,132 @@ void measureEOC(const GridView gridView,
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(embeddedX[i]);
+ // The numerical solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
+
+ ///////////////////////////////////////////////////////////////////////////
+ // Read the reference configuration
+ ///////////////////////////////////////////////////////////////////////////
+ EmbeddedVectorType embeddedReferenceX(referenceFEBasis.size());
+ inFile.open(parameterSet.get<std::string>("referenceData"), std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("referenceData") << " could not be opened.");
+ GenericVector::readBinary(inFile, embeddedReferenceX);
+
+ SolutionType referenceX(embeddedReferenceX.size());
+ for (size_t i=0; i<referenceX.size(); i++)
+ referenceX[i] = TargetSpace(embeddedReferenceX[i]);
+
+ // The reference solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> referenceSolution(fufemReferenceFEBasis, referenceX);
+
/////////////////////////////////////////////////////////////////
// Measure the discretization error
/////////////////////////////////////////////////////////////////
- if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
+ double l2ErrorSquared = 0;
+ double h1ErrorSquared = 0;
+
+ HierarchicSearch<typename GridView::Grid,typename GridView::IndexSet> hierarchicSearch(gridView.grid(), gridView.indexSet());
+
+ for (const auto& rElement : elements(referenceGridView))
{
- using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
- FufemFEBasis fufemFEBasis(feBasis);
+ const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
+
+ for (const auto& qp : quadRule)
+ {
+ auto integrationElement = rElement.geometry().integrationElement(qp.position());
- // Read reference solution and its derivative into a PythonFunction
- typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
+ auto globalPos = rElement.geometry().global(qp.position());
- Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
- auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
+ auto element = hierarchicSearch.findEntity(globalPos);
+ auto localPos = element.geometry().local(globalPos);
- // The numerical solution, as a grid function
- GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
+ auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
- // QuadratureRule for the integral of the L^2 error
- QuadratureRuleKey quadKey(dim,6);
+ l2ErrorSquared += integrationElement * qp.weight() * diff.two_norm2();
- // Compute the embedded L^2 error
- double l2Error = DiscretizationError<GridView>::computeL2Error(&numericalSolution,
- referenceSolution.get(),
- quadKey);
+ auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
- // Compute the embedded H^1 error
- double h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
- &numericalSolution,
- referenceSolution.get(),
- quadKey);
+ h1ErrorSquared += integrationElement * qp.weight() * derDiff.frobenius_norm2();
- std::cout << "elements: " << gridView.size(0)
- << " "
- << "L^2 error: " << l2Error
- << " ";
- std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
+ }
}
- if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
- {
- FEBasis referenceFEBasis(referenceGridView);
+ std::cout << "levels: " << gridView.grid().maxLevel()+1
+ << " "
+ << "L^2 error: " << std::sqrt(l2ErrorSquared)
+ << " "
+ << "H^1 error: " << std::sqrt(l2ErrorSquared + h1ErrorSquared)
+ << std::endl;
+}
- // Reference configuration
- EmbeddedVectorType embeddedReferenceX(referenceFEBasis.size());
- inFile.open(parameterSet.get<std::string>("referenceData"), std::ios_base::binary);
- if (not inFile)
- DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("referenceData") << " could not be opened.");
- GenericVector::readBinary(inFile, embeddedReferenceX);
+template <class GridView, int order>
+void measureAnalyticalEOC(const GridView gridView,
+ const ParameterTree& parameterSet)
+{
+ typedef std::vector<TargetSpace> SolutionType;
- SolutionType referenceX(embeddedReferenceX.size());
- for (size_t i=0; i<referenceX.size(); i++)
- referenceX[i] = TargetSpace(embeddedReferenceX[i]);
+ //////////////////////////////////////////////////////////////////////////////////
+ // Construct the scalar function space bases corresponding to the GFE space
+ //////////////////////////////////////////////////////////////////////////////////
- using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
- FufemFEBasis fufemReferenceFEBasis(referenceFEBasis);
- FufemFEBasis fufemFEBasis(feBasis);
+ typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
+ FEBasis feBasis(gridView);
- // The numerical solution, as a grid function
- GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> referenceSolution(fufemReferenceFEBasis, referenceX);
+ //////////////////////////////////////////////////////////////////////////////////
+ // Read the data whose error is to be measured
+ //////////////////////////////////////////////////////////////////////////////////
- // The numerical solution, as a grid function
- GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
+ // Input data
+ typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
- double l2ErrorSquared = 0;
- double h1ErrorSquared = 0;
+ EmbeddedVectorType embeddedX(feBasis.size());
+ std::ifstream inFile(parameterSet.get<std::string>("simulationData"), std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("simulationData") << " could not be opened.");
+ GenericVector::readBinary(inFile, embeddedX);
+ inFile.close();
- HierarchicSearch<typename GridView::Grid,typename GridView::IndexSet> hierarchicSearch(gridView.grid(), gridView.indexSet());
+ SolutionType x(embeddedX.size());
+ for (size_t i=0; i<x.size(); i++)
+ x[i] = TargetSpace(embeddedX[i]);
- std::cout << "l2: " << l2ErrorSquared << std::endl;
- for (const auto& rElement : elements(referenceGridView))
- {
- const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
+ /////////////////////////////////////////////////////////////////
+ // Measure the discretization error
+ /////////////////////////////////////////////////////////////////
- for (const auto& qp : quadRule)
- {
- auto integrationElement = rElement.geometry().integrationElement(qp.position());
+ using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
+ FufemFEBasis fufemFEBasis(feBasis);
- auto globalPos = rElement.geometry().global(qp.position());
+ // Read reference solution and its derivative into a PythonFunction
+ typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
- auto element = hierarchicSearch.findEntity(globalPos);
- auto localPos = element.geometry().local(globalPos);
+ Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
+ auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
- auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
+ // The numerical solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
- l2ErrorSquared += integrationElement * qp.weight() * diff.two_norm2();
+ // QuadratureRule for the integral of the L^2 error
+ QuadratureRuleKey quadKey(dim,6);
- auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
+ // Compute the embedded L^2 error
+ double l2Error = DiscretizationError<GridView>::computeL2Error(&numericalSolution,
+ referenceSolution.get(),
+ quadKey);
- h1ErrorSquared += integrationElement * qp.weight() * derDiff.frobenius_norm2();
+ // Compute the embedded H^1 error
+ double h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
+ &numericalSolution,
+ referenceSolution.get(),
+ quadKey);
- }
- }
- std::cout << "l2: " << l2ErrorSquared << std::endl;
-
- std::cout << "levels: " << gridView.grid().maxLevel()+1
- << " "
- << "L^2 error: " << std::sqrt(l2ErrorSquared)
- << " "
- << "H^1 error: " << std::sqrt(l2ErrorSquared + h1ErrorSquared)
- << std::endl;
- }
+ std::cout << "elements: " << gridView.size(0)
+ << " "
+ << "L^2 error: " << l2Error
+ << " ";
+ std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
}
@@ -231,25 +260,53 @@ int main (int argc, char *argv[]) try
referenceGrid->globalRefine(parameterSet.get<int>("numReferenceLevels")-1);
// Do the actual measurement
+
const int order = parameterSet.get<int>("order");
- switch (order)
+
+ if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
+ {
+ switch (order)
+ {
+ case 1:
+ measureDiscreteEOC<GridType::LeafGridView,1>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ case 2:
+ measureDiscreteEOC<GridType::LeafGridView,2>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ case 3:
+ measureDiscreteEOC<GridType::LeafGridView,3>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ default:
+ DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
+ }
+ }
+
+ if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
{
- case 1:
- measureEOC<GridType::LeafGridView,1>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
- break;
+ switch (order)
+ {
+ case 1:
+ measureAnalyticalEOC<GridType::LeafGridView,1>(grid->leafGridView(), parameterSet);
+ break;
- case 2:
- measureEOC<GridType::LeafGridView,2>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
- break;
+ case 2:
+ measureAnalyticalEOC<GridType::LeafGridView,2>(grid->leafGridView(), parameterSet);
+ break;
- case 3:
- measureEOC<GridType::LeafGridView,3>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
- break;
+ case 3:
+ measureAnalyticalEOC<GridType::LeafGridView,3>(grid->leafGridView(), parameterSet);
+ break;
- default:
- DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
+ default:
+ DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
+ }
}
+
+
return 0;
}
catch (Exception e)
--
GitLab