From ce928b3310dcb1ec2c0d7f142789809ef9c8ca9d Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 1 Jan 2012 16:37:29 +0000
Subject: [PATCH] Helper method which computes the energy norm (squared) of the
difference of the fe functions on two different refinements of the same base
grid.
[[Imported from SVN: r8331]]
---
dune/gfe/gfedifferencenormsquared.hh | 288 +++++++++++++++++++++++++++
1 file changed, 288 insertions(+)
create mode 100644 dune/gfe/gfedifferencenormsquared.hh
diff --git a/dune/gfe/gfedifferencenormsquared.hh b/dune/gfe/gfedifferencenormsquared.hh
new file mode 100644
index 00000000..9b452e07
--- /dev/null
+++ b/dune/gfe/gfedifferencenormsquared.hh
@@ -0,0 +1,288 @@
+#ifndef GFE_DIFFERENCE_NORM_SQUARED_HH
+#define GFE_DIFFERENCE_NORM_SQUARED_HH
+
+/** \file
+ \brief Helper method which computes the energy norm (squared) of the difference of
+ the fe functions on two different refinements of the same base grid.
+*/
+
+#include <dune/geometry/type.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/fufem/assemblers/localoperatorassembler.hh>
+
+#include <dune/gfe/globalgeodesicfefunction.hh>
+
+template <class Basis, class TargetSpace>
+class GFEDifferenceNormSquared {
+
+ typedef typename Basis::GridView GridView;
+ typedef typename Basis::GridView::Grid GridType;
+
+ // Grid dimension
+ const static int dim = GridType::dimension;
+
+ /** \brief Check whether given local coordinates are contained in the reference element
+ \todo This method exists in the Dune grid interface! But we need the eps.
+ */
+ static bool checkInside(const Dune::GeometryType& type,
+ const Dune::FieldVector<double, dim> &loc,
+ double eps)
+ {
+ switch (type.dim()) {
+
+ case 0: // vertex
+ return false;
+ case 1: // line
+ return -eps <= loc[0] && loc[0] <= 1+eps;
+ case 2:
+
+ if (type.isSimplex()) {
+ return -eps <= loc[0] && -eps <= loc[1] && (loc[0]+loc[1])<=1+eps;
+ } else if (type.isCube()) {
+ return -eps <= loc[0] && loc[0] <= 1+eps
+ && -eps <= loc[1] && loc[1] <= 1+eps;
+ } else {
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+ }
+
+ case 3:
+ if (type.isSimplex()) {
+ return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+ && (loc[0]+loc[1]+loc[2]) <= 1+eps;
+ } else if (type.isPyramid()) {
+ return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+ && (loc[0]+loc[2]) <= 1+eps
+ && (loc[1]+loc[2]) <= 1+eps;
+ } else if (type.isPrism()) {
+ return -eps <= loc[0] && -eps <= loc[1]
+ && (loc[0]+loc[1])<= 1+eps
+ && -eps <= loc[2] && loc[2] <= 1+eps;
+ } else if (type.isCube()) {
+ return -eps <= loc[0] && loc[0] <= 1+eps
+ && -eps <= loc[1] && loc[1] <= 1+eps
+ && -eps <= loc[2] && loc[2] <= 1+eps;
+ }else {
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+ }
+ default:
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+ }
+
+ }
+
+
+ /** \brief Coordinate function in one variable, constant in the others
+
+ This is used to extract the positions of the Lagrange nodes.
+ */
+ struct CoordinateFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,1> >
+ {
+ CoordinateFunction(int d)
+ : d_(d)
+ {}
+
+ void evaluate(const Dune::FieldVector<double, dim>& x, Dune::FieldVector<double,1>& out) const {
+ out[0] = x[d_];
+ }
+
+ //
+ int d_;
+ };
+
+public:
+
+ static void interpolate(const GridType& sourceGrid,
+ const std::vector<TargetSpace>& source,
+ const GridType& targetGrid,
+ std::vector<TargetSpace>& target)
+ {
+ // Create a leaf function, which we need to be able to call 'evalall()'
+ Basis sourceBasis(sourceGrid.leafView());
+ GlobalGeodesicFEFunction<Basis,TargetSpace> sourceFunction(sourceBasis, source);
+
+ Basis targetBasis(targetGrid.leafView());
+
+ // ///////////////////////////////////////////////////////////////////////////////////////////
+ // Prolong the adaptive solution onto the uniform grid in order to make it comparable
+ // ///////////////////////////////////////////////////////////////////////////////////////////
+
+ target.resize(targetBasis.size());
+
+ // handle each dof only once
+ Dune::BitSetVector<1> handled(targetBasis.size(), false);
+
+ typename GridView::template Codim<0>::Iterator eIt = targetGrid.template leafbegin<0>();
+ typename GridView::template Codim<0>::Iterator eEndIt = targetGrid.template leafend<0>();
+
+ // Loop over all dofs by looping over all elements
+ for (; eIt!=eEndIt; ++eIt) {
+
+ const typename Basis::LocalFiniteElement& targetLFE = targetBasis.getLocalFiniteElement(*eIt);
+
+ // Generate position of the Lagrange nodes
+ std::vector<Dune::FieldVector<double,dim> > lagrangeNodes(targetLFE.localBasis().size());
+
+ for (int i=0; i<dim; i++) {
+ CoordinateFunction lFunction(i);
+ std::vector<Dune::FieldVector<double,1> > coordinates;
+ targetLFE.localInterpolation().interpolate(lFunction, coordinates);
+
+ for (size_t j=0; j<coordinates.size(); j++)
+ lagrangeNodes[j][i] = coordinates[j];
+
+ }
+
+ for (size_t i=0; i<targetLFE.localCoefficients().size(); i++) {
+
+ typename GridType::template Codim<0>::EntityPointer element(eIt);
+
+ Dune::FieldVector<double,dim> pos = lagrangeNodes[i];
+
+ if (handled[targetBasis.index(*eIt,i)][0])
+ continue;
+
+ handled[targetBasis.index(*eIt,i)] = true;
+
+ assert(checkInside(element->type(), pos, 1e-7));
+
+ // ////////////////////////////////////////////////////////////////////
+ // Get an element on the coarsest grid which contains the vertex and
+ // its local coordinates there
+ // ////////////////////////////////////////////////////////////////////
+ while (element->level() != 0){
+
+ pos = element->geometryInFather().global(pos);
+ element = element->father();
+
+ assert(checkInside(element->type(), pos, 1e-7));
+ }
+
+ // ////////////////////////////////////////////////////////////////////
+ // Find the corresponding coarse grid element on the adaptive grid.
+ // This is a linear algorithm, but we expect the coarse grid to be small.
+ // ////////////////////////////////////////////////////////////////////
+ Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> uniformP0Mapper (targetGrid, 0);
+ Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,Dune::MCMGElementLayout> adaptiveP0Mapper(sourceGrid, 0);
+ int coarseIndex = uniformP0Mapper.map(*element);
+
+ typename GridType::template Codim<0>::LevelIterator adaptEIt = sourceGrid.template lbegin<0>(0);
+ typename GridType::template Codim<0>::LevelIterator adaptEEndIt = sourceGrid.template lend<0>(0);
+
+ for (; adaptEIt!=adaptEEndIt; ++adaptEIt)
+ if (adaptiveP0Mapper.map(*adaptEIt) == coarseIndex)
+ break;
+
+ element = adaptEIt;
+
+ // ////////////////////////////////////////////////////////////////////////
+ // Find a corresponding point (not necessarily vertex) on the leaf level
+ // of the adaptive grid.
+ // ////////////////////////////////////////////////////////////////////////
+
+ while (!element->isLeaf()) {
+
+
+ typename GridType::template Codim<0>::Entity::HierarchicIterator hIt = element->hbegin(element->level()+1);
+ typename GridType::template Codim<0>::Entity::HierarchicIterator hEndIt = element->hend(element->level()+1);
+
+ Dune::FieldVector<double,dim> childPos;
+ assert(checkInside(element->type(), pos, 1e-7));
+
+ for (; hIt!=hEndIt; ++hIt) {
+
+ childPos = hIt->geometryInFather().local(pos);
+ if (checkInside(hIt->type(), childPos, 1e-7))
+ break;
+
+ }
+
+ assert(hIt!=hEndIt);
+ element = hIt;
+ pos = childPos;
+
+ }
+
+ // ////////////////////////////////////////////////////////////////////////
+ // Sample adaptive function
+ // ////////////////////////////////////////////////////////////////////////
+ sourceFunction.evaluateLocal(*element, pos,
+ target[targetBasis.index(*eIt,i)]);
+
+ }
+
+ }
+
+ }
+
+
+ template <int blocksize>
+ static double computeNormSquared(const GridType& grid,
+ const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& x,
+ const LocalOperatorAssembler<GridType,
+ typename Basis::LocalFiniteElement,
+ typename Basis::LocalFiniteElement,
+ Dune::FieldMatrix<double,1,1> >* localStiffness)
+ {
+ // ///////////////////////////////////////////////////////////////////////////////////////////
+ // Compute the energy norm
+ // ///////////////////////////////////////////////////////////////////////////////////////////
+
+ double energyNormSquared = 0;
+ Basis basis(grid.leafView());
+ Dune::Matrix<Dune::FieldMatrix<double,1,1> > localMatrix;
+
+ typename GridType::template Codim<0>::LeafIterator eIt = grid.template leafbegin<0>();
+ typename GridType::template Codim<0>::LeafIterator eEndIt = grid.template leafend<0>();
+
+ for (; eIt!=eEndIt; ++eIt) {
+
+ size_t nLocalDofs = basis.getLocalFiniteElement(*eIt).localCoefficients().size();
+
+ localMatrix.setSize(nLocalDofs,nLocalDofs);
+
+ localStiffness->assemble(*eIt, localMatrix,
+ basis.getLocalFiniteElement(*eIt),
+ basis.getLocalFiniteElement(*eIt));
+
+ for (size_t i=0; i<nLocalDofs; i++)
+ for (size_t j=0; j<nLocalDofs; j++)
+ energyNormSquared += localMatrix[i][j][0][0] * (x[basis.index(*eIt,i)] * x[basis.index(*eIt,j)]);
+
+ }
+
+ return energyNormSquared;
+ }
+
+ static double compute(const GridType& uniformGrid,
+ const std::vector<TargetSpace>& uniformSolution,
+ const GridType& adaptiveGrid,
+ const std::vector<TargetSpace>& adaptiveSolution,
+ const LocalOperatorAssembler<GridType,
+ typename Basis::LocalFiniteElement,
+ typename Basis::LocalFiniteElement,
+ Dune::FieldMatrix<double,1,1> >* localStiffness)
+ {
+ std::vector<TargetSpace> uniformAdaptiveSolution;
+
+ interpolate(adaptiveGrid, adaptiveSolution, uniformGrid, uniformAdaptiveSolution);
+
+ // Compute difference in the embedding space
+ Dune::BlockVector<typename TargetSpace::CoordinateType> difference(uniformSolution.size());
+
+ for (size_t i=0; i<difference.size(); i++)
+ difference[i] = uniformAdaptiveSolution[i].globalCoordinates() - uniformSolution[i].globalCoordinates();
+
+ //std::cout << "new difference:\n" << difference << std::endl;
+// std::cout << "new projected:\n" << std::endl;
+// for (size_t i=0; i<difference.size(); i++)
+// std::cout << uniformAdaptiveSolution[i].globalCoordinates() << std::endl;
+
+ return computeNormSquared(uniformGrid, difference, localStiffness);
+ }
+};
+
+#endif
--
GitLab