...
  View open merge request
Commits (2)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_FUNCTIONRANGE_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_FUNCTIONRANGE_HH
#include <dune/common/fvector.hh>
#include <dune/common/concept.hh>
#include <dune/functions/functionspacebases/concepts.hh>
namespace Dune {
namespace Functions {
// Function range for vector valued functions
template<typename D, typename R, typename DHasIndexAccess = void, typename RHasIndexAccess = void>
struct FunctionRange {};
template<typename D, typename R>
struct FunctionRange<D, R,
typename std::enable_if< models<Concept::HasIndexAccess, D, int>() >::type,
typename std::enable_if< models<Concept::HasIndexAccess, R, int>() >::type>
{
using Domain = D;
using Range = R;
using RangeField = typename Range::value_type;
using JacobianRange = FieldMatrix<RangeField, Range::dimension, Domain::dimension>;
};
template<typename D, typename R>
struct FunctionRange<D, R,
typename std::enable_if< models<Concept::HasIndexAccess, D, int>() >::type,
typename std::enable_if< not models<Concept::HasIndexAccess, R, int>() >::type>
{
using Domain = D;
using Range = R;
using RangeField = Range;
using JacobianRange = FieldVector<RangeField, Domain::dimension>;
};
template<typename D, typename R>
struct FunctionRange<D, R,
typename std::enable_if< not models<Concept::HasIndexAccess, D, int>() >::type,
typename std::enable_if< models<Concept::HasIndexAccess, R, int>() >::type>
{
using Domain = D;
using Range = R;
using RangeField = typename Range::value_type;
using JacobianRange = FieldVector<RangeField, Range::dimension>;
};
template<typename D, typename R>
struct FunctionRange<D, R,
typename std::enable_if< not models<Concept::HasIndexAccess, D, int>() >::type,
typename std::enable_if< not models<Concept::HasIndexAccess, R, int>() >::type>
{
using Domain = D;
using Range = R;
using RangeField = Range;
using JacobianRange = RangeField;
};
} // namespace Dune::Functions
} // namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONRANGE_HH
......@@ -12,6 +12,7 @@
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/functions/common/treedata.hh>
#include <dune/functions/common/functionrange.hh>
namespace Dune {
namespace Functions {
......@@ -81,6 +82,7 @@ public:
using Domain = typename EntitySet::GlobalCoordinate;
using Range = R;
using JacobianRange = typename FunctionRange<Domain,Range>::JacobianRange;
using LocalDomain = typename EntitySet::LocalCoordinate;
using Element = typename EntitySet::Element;
......@@ -95,11 +97,16 @@ public:
template<class Node>
using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
template<class Node>
using LocalBasisJacobian = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
template<class Node>
using NodeData = typename std::vector<LocalBasisRange<Node>>;
template<class Node>
using NodeDataJacobian = typename std::vector<LocalBasisJacobian<Node>>;
using ShapeFunctionValueContainer = TreeData<SubTree, NodeData, true>;
using ShapeFunctionJacobianContainer = TreeData<SubTree, NodeDataJacobian, true>;
struct LocalEvaluateVisitor
: public TypeTree::TreeVisitor
......@@ -169,12 +176,90 @@ public:
ShapeFunctionValueContainer& shapeFunctionValueContainer_;
};
struct LocalJacobianEvaluateVisitor
: public TypeTree::TreeVisitor
, public TypeTree::DynamicTraversal
{
LocalJacobianEvaluateVisitor(const LocalDomain& x, JacobianRange& y, const LocalIndexSet& localIndexSet, const Vector& coefficients, const NodeToRangeEntry& nodeToRangeEntry, ShapeFunctionJacobianContainer& shapeFunctionJacobianContainer):
x_(x),
y_(y),
localIndexSet_(localIndexSet),
coefficients_(coefficients),
nodeToRangeEntry_(nodeToRangeEntry),
shapeFunctionJacobianContainer_(shapeFunctionJacobianContainer)
{}
template<typename Node, typename TreePath>
void leaf(Node& node, TreePath treePath)
{
using LocalBasisJacobianRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
using MultiIndex = typename LocalIndexSet::MultiIndex;
using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
using RangeBlock = typename std::decay<decltype(nodeToRangeEntry_(node, y_))>::type;
auto&& fe = node.finiteElement();
auto&& localBasis = fe.localBasis();
auto&& shapeFunctionJacobians = shapeFunctionJacobianContainer_[node];
localBasis.evaluateJacobian(x_, shapeFunctionJacobians);
// Require transformation to physical element
auto element = node.element();
auto geometry = element.geometry();
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto referenceShapeFunctionJacobian = shapeFunctionJacobians[i];
jacobianInverseTransposed.mv(referenceShapeFunctionJacobian[0], shapeFunctionJacobians[i][0]);
}
// Get range entry associated to this node
auto&& re = nodeToRangeEntry_(node, y_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
// Get coefficient associated to i-th shape function
auto&& c = coefficients_[multiIndex];
// Get value of i-th shape function
auto&& v = shapeFunctionJacobians[i];
// Notice that the range entry re, the coefficient c, and the shape functions
// value v may all be scalar, vector, matrix, or general container valued.
// The matching of their entries is done via the multistage procedure described
// in the class documentation of DiscreteGlobalBasisFunction.
auto dimC = FlatVectorBackend<CoefficientBlock>::size(c);
auto dimV = FlatVectorBackend<LocalBasisJacobianRange>::size(v);
assert(dimC*dimV == FlatVectorBackend<RangeBlock>::size(re));
for(size_type j=0; j<dimC; ++j)
{
auto&& c_j = FlatVectorBackend<CoefficientBlock>::getEntry(c, j);
for(size_type k=0; k<dimV; ++k)
{
auto&& v_k = FlatVectorBackend<LocalBasisJacobianRange>::getEntry(v, k);
FlatVectorBackend<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
}
}
}
}
const LocalDomain& x_;
JacobianRange& y_;
const LocalIndexSet& localIndexSet_;
const Vector& coefficients_;
const NodeToRangeEntry& nodeToRangeEntry_;
ShapeFunctionJacobianContainer& shapeFunctionJacobianContainer_;
};
public:
using GlobalFunction = DiscreteGlobalBasisFunction;
using Domain = LocalDomain;
using Range = GlobalFunction::Range;
using JacobianRange = GlobalFunction::JacobianRange;
using Element = GlobalFunction::Element;
LocalFunction(const DiscreteGlobalBasisFunction& globalFunction)
......@@ -186,6 +271,7 @@ public:
// and queried for tree indices even in unbound state.
subTree_ = &TypeTree::child(localBasisView_.tree(), globalFunction_->treePath());
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
// localDoFs_.reserve(localBasisView_.maxSize());
}
......@@ -198,6 +284,7 @@ public:
// and queried for tree indices even in unbound state.
subTree_ = &TypeTree::child(localBasisView_.tree(), globalFunction_->treePath());
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
}
LocalFunction operator=(const LocalFunction& other)
......@@ -210,6 +297,7 @@ public:
// Here we assume that the tree can be accessed, traversed,
// and queried for tree indices even in unbound state.
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
}
/**
......@@ -254,6 +342,21 @@ public:
return y;
}
/**
* \brief Evaluate derivative of LocalFunction at bound element.
*/
JacobianRange derivative(const Domain& x) const
{
JacobianRange y;
for (size_t j=0; j<y.size(); ++j)
y[j] = 0.0;
LocalJacobianEvaluateVisitor localJacobianEvaluateVisitor(x, y, localIndexSet_, globalFunction_->dofs(), globalFunction_->nodeToRangeEntry(), shapeFunctionJacobianContainer_);
TypeTree::applyToTree(*subTree_, localJacobianEvaluateVisitor);
return y;
}
const Element& localContext() const
{
return localBasisView_.element();
......@@ -271,6 +374,7 @@ public:
LocalIndexSet localIndexSet_;
mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
mutable ShapeFunctionJacobianContainer shapeFunctionJacobianContainer_;
// std::vector<typename V::value_type> localDoFs_;
const SubTree* subTree_;
};
......
......@@ -5,3 +5,5 @@ dune_add_test(SOURCES analyticgridviewfunctiontest.cc)
dune_add_test(SOURCES discreteglobalbasisfunctiontest.cc)
dune_add_test(SOURCES gridfunctiontest.cc)
dune_add_test(SOURCES derivativetest.cc)
#include <config.h>
#include <dune/common/parallel/mpihelper.hh>
#include<dune/geometry/type.hh>
#include<dune/geometry/quadraturerules.hh>
#include <dune/istl/bvector.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/grid/yaspgrid.hh>
#include "dune/functions/functionspacebases/interpolate.hh"
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <math.h>
#define VECTORVALUED
using namespace Dune;
// Test case:
// Define vector valued finite element function (piecewise bilinear) on three dimensional grid
// by interpolating a given nonlinear, differentiable function. Compare the derivatives of the
// approximation with the exact derivative. Repeat the same test on finer grids and check
// if the interpolation error behaves as expected.
// Compute L2 norm of the difference between a discrete global bassis function and another gridfunctions.
template <class GridView, class FunctionA, class FunctionB>
double l2difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionA(quadPos);
diff -= functionB(quadPos);
#ifdef VECTORVALUED
auto norm = diff.two_norm2();
#else
auto norm = diff * diff;
#endif
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
// Compute L2 norm of the difference between the gradient of a discrete
// global bassis function and another gridfunctions.
template <class GridView, class FunctionA, class FunctionB>
double h1difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionB(quadPos);
// Compute the derivative of a discrete global basis function.
auto grad = functionA.derivative(quadPos);
for (size_t i=0; i<diff.N(); ++i)
#ifdef VECTORVALUED
for (size_t j=0; j<diff.M(); ++j)
diff[i][j] -= grad[i][j];
auto norm = diff.frobenius_norm2();
#else
diff[i] -= grad[i];
auto norm = diff.two_norm2();
#endif
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define a continuous function and its derivative
// f(x,y,z) = [ sin(2πx) * sin(2πy) * sin(2πz) ]
// [ x^2 * exp(y) - z ];
//
// Df(x,y,z) = [ 2π * cos(2πx) * sin(2πy) * sin(2πz) , 2π * sin(2πx) * cos(2πy) * sin(2πz), 2π * sin(2πx) * sin(2πy) * cos(2πz) ]
// [ 2x * exp(y) , x^2 * exp(y) , -1 ];
//
//////////////////////////////////////////////////////////////////////////////////////////////
auto pi = std::acos(-1.0);
auto function = [&pi] (const auto& x)
{
#ifdef VECTORVALUED
FieldVector<double, 2> result(0.0);
#else
double result(0.0);
#endif
#ifdef VECTORVALUED
result[0] = sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[1] = pow(x[0],2) * exp(x[1]) - x[2];
#else
result = sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
#endif
return result;
};
auto derivative = [&pi] (const auto& x)
{
#ifdef VECTORVALUED
FieldMatrix<double, 2, 3> result(0.0);
result[0][0] = 2.*pi * cos(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][1] = 2.*pi * sin(2.*pi*x[0]) * cos(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][2] = 2.*pi * sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * cos(2.*pi*x[2]);
result[1][0] = 2. * x[0] * exp(x[1]);
result[1][1] = pow(x[0],2) * exp(x[1]);
result[1][2] = - 1.;
#else
FieldVector<double, 3> result(0.0);
result[0] = 2.*pi * cos(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[1] = 2.*pi * sin(2.*pi*x[0]) * cos(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[2] = 2.*pi * sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * cos(2.*pi*x[2]);
#endif
return result;
};
//////////////////////////////////////////////////////////////////////////////////////////////
// Generate a three-dimensional grid
//////////////////////////////////////////////////////////////////////////////////////////////
const int dim = 3;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> ur(1.0);
std::array<int,dim> numElements = {10,10,10};
GridType grid(ur,numElements);
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
for (size_t refineLevel = 0; refineLevel < 3; ++refineLevel)
{
//////////////////////////////////////////////////////////////////////////////////////////////
// Define the continuous function and its derivative as grid functions
//////////////////////////////////////////////////////////////////////////////////////////////
auto gridFunction = localFunction(Functions::makeGridViewFunction(function, gridView));
auto gridDerivative = localFunction(Functions::makeGridViewFunction(derivative, gridView));
//////////////////////////////////////////////////////////////////////////////////////////////
// Choose a finite element space (piecewise bilinear)
//////////////////////////////////////////////////////////////////////////////////////////////
using namespace Functions::BasisBuilder;
#ifdef VECTORVALUED
auto basis = makeBasis(gridView, power<2>(pq<1>(),flatInterleaved()));
#else
auto basis = makeBasis(gridView, pq<1>());
#endif
//////////////////////////////////////////////////////////////////////////////////////////////
// Define array type formats
//////////////////////////////////////////////////////////////////////////////////////////////
typedef BlockVector<FieldVector<double,1> > VectorType;
typedef Dune::Functions::HierarchicVectorWrapper<VectorType, double> HierarchicVectorView;
//////////////////////////////////////////////////////////////////////////////////////////////
// Define a coefficient vector corresponding to a discrete function
//////////////////////////////////////////////////////////////////////////////////////////////
VectorType x;
HierarchicVectorView(x).resize(basis);
x = 0;
//////////////////////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
//////////////////////////////////////////////////////////////////////////////////////////////
#ifdef VECTORVALUED
using FunctionRange = FieldVector<double,2>;
#else
using FunctionRange = double;
#endif
auto discreteFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FunctionRange>(basis, TypeTree::hybridTreePath(), HierarchicVectorView(x));
//////////////////////////////////////////////////////////////////////////////////////////////
// Transform discrete functions to local functions
//////////////////////////////////////////////////////////////////////////////////////////////
auto localDiscreteFunction = localFunction(discreteFunction);
//////////////////////////////////////////////////////////////////////////////////////////////
// Interpolate the function
//////////////////////////////////////////////////////////////////////////////////////////////
interpolate(basis, Dune::TypeTree::hybridTreePath(), HierarchicVectorView(x), function);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define VTK writer and output
//////////////////////////////////////////////////////////////////////////////////////////////
// SubsamplingVTKWriter<GridView> vtkWriter(gridView,1);
//#ifdef VECTORVALUED
// vtkWriter.addVertexData(discreteFunction, VTK::FieldInfo("function", VTK::FieldInfo::Type::vector, 2));
//#else
// vtkWriter.addVertexData(discreteFunction, VTK::FieldInfo("function", VTK::FieldInfo::Type::scalar, 1));
//#endif
// vtkWriter.write("derivativetest");
//////////////////////////////////////////////////////////////////////////////////////////////
// Check l2 difference for function and its derivative
//////////////////////////////////////////////////////////////////////////////////////////////
auto l2diff = l2difference(gridView, localDiscreteFunction, gridFunction);
auto h1diff = h1difference(gridView, localDiscreteFunction, gridDerivative);
std::cout << "Refine level: " << refineLevel << std::endl;
std::cout << "L2 interpolation error: " << l2diff << std::endl;
std::cout << "H1-seminorm interpolation error: " << h1diff << std::endl << std::endl;
grid.globalRefine(1);
}
}
// Error handling
catch (Exception e) {
std::cout << e << std::endl;
}