diff --git a/test/interillustration.cc b/test/interillustration.cc
index 16af3cce36208c236a0d0090e192eb720cdb634c..17deabe4a7fc8c781f96aa24a6d49776e654d83c 100644
--- a/test/interillustration.cc
+++ b/test/interillustration.cc
@@ -5,12 +5,17 @@
#include <dune/grid/io/file/vtk.hh>
#include <dune/localfunctions/lagrange/p1.hh>
+#include <dune/localfunctions/lagrange/p2.hh>
+
+#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
+#include <dune/functions/gridfunctions/discretescalarglobalbasisfunction.hh>
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/localtangentfefunction.hh>
+
using namespace Dune;
/** \brief Encapsulates the grid deformation for the GeometryGrid class */
@@ -55,9 +60,13 @@ private:
};
-
+/** \brief
+ *
+ * \param partialDerivative The dof with respect to which we are computing the variation
+ */
template <class TargetSpace, class LocalFEFunctionType>
-void interpolate(std::vector<TargetSpace> coefficients,
+void interpolate(const LocalFEFunctionType& localGeodesicFEFunction,
+ int partialDerivative,
const std::string& nameSuffix)
{
static const int dim = 2;
@@ -76,20 +85,11 @@ void interpolate(std::vector<TargetSpace> coefficients,
shared_ptr<GridType> grid = shared_ptr<GridType>(factory.createGrid());
- grid->globalRefine(6);
+ grid->globalRefine(3);
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- //////////////////////////////////////////////////////////////
- // Set up an interpolation function on the sphere
- //////////////////////////////////////////////////////////////
-
- typedef P1LocalFiniteElement<double,double,dim> LocalFiniteElement;
- LocalFiniteElement localFiniteElement;
-
- LocalFEFunctionType localGeodesicFEFunction(localFiniteElement,coefficients);
-
///////////////////////////////////////////////////////////////
// Sample the interpolating function on the grid
///////////////////////////////////////////////////////////////
@@ -99,6 +99,20 @@ void interpolate(std::vector<TargetSpace> coefficients,
for (auto v=gridView.begin<dim>(); v!=gridView.end<dim>(); ++v)
samples[gridView.indexSet().index(*v)] = localGeodesicFEFunction.evaluate(v->geometry().corner(0)).globalCoordinates();
+ // Sample a variation vector field
+ std::vector<typename TargetSpace::EmbeddedTangentVector> variation0(gridView.size(dim));
+ std::vector<typename TargetSpace::EmbeddedTangentVector> variation1(gridView.size(dim));
+ for (const auto& v : vertices(gridView))
+ {
+ FieldMatrix<double,3,3> derivative;
+ localGeodesicFEFunction.evaluateDerivativeOfValueWRTCoefficient(v.geometry().corner(0),
+ partialDerivative, // derive wrt value at first triangle corner
+ derivative);
+
+ variation0[gridView.indexSet().index(v)] = derivative[1];
+ variation1[gridView.indexSet().index(v)] = derivative[2];
+ }
+
// sample a checkerboard pattern for nicer pictures
uint pattern = 8;
std::vector<double> colors(gridView.size(0));
@@ -120,52 +134,81 @@ void interpolate(std::vector<TargetSpace> coefficients,
// Write a grid with the interpolated positions
///////////////////////////////////////////////////////////////
-
typedef Dune::GeometryGrid<GridType,DeformationFunction<typename GridType::LeafGridView> > DeformedGridType;
DeformationFunction<GridView> deformationFunction(gridView, samples);
// stupid, can't instantiate deformedGrid with a const grid
DeformedGridType deformedGrid(const_cast<GridType&>(*grid), deformationFunction);
+
+ typedef Functions::PQ1NodalBasis<typename DeformedGridType::LeafGridView > FEBasis;
+ FEBasis feBasis(deformedGrid.leafGridView());
+
+ Functions::DiscreteScalarGlobalBasisFunction<decltype(feBasis),decltype(variation0)> variation0Function(feBasis,variation0);
+ Functions::DiscreteScalarGlobalBasisFunction<decltype(feBasis),decltype(variation1)> variation1Function(feBasis,variation1);
+ auto localVariation0Function = localFunction(variation0Function);
+ auto localVariation1Function = localFunction(variation1Function);
+
Dune::VTKWriter<typename DeformedGridType::LeafGridView> vtkWriter(deformedGrid.leafGridView());
vtkWriter.addCellData(colors, "colors");
- vtkWriter.write("sphere-patch-" + nameSuffix);
+ vtkWriter.addVertexData(localVariation0Function, VTK::FieldInfo("variation 0", VTK::FieldInfo::Type::scalar, variation0[0].size()));
+ vtkWriter.addVertexData(localVariation1Function, VTK::FieldInfo("variation 1", VTK::FieldInfo::Type::scalar, variation1[0].size()));
+
+ vtkWriter.write("sphere-patch-" + nameSuffix);
}
int main(int argc, char* argv[]) try
{
static const int dim = 2;
+ typedef UnitVector<double,3> TargetSpace;
//////////////////////////////////////////////////////////////
- // Set up an interpolation function on the sphere
+ // Set up a first-order interpolation function on the sphere
//////////////////////////////////////////////////////////////
- typedef UnitVector<double,3> TargetSpace;
-
std::vector<TargetSpace> coefficients(3);
coefficients[0] = TargetSpace(FieldVector<double,3>({1,0,0}));
coefficients[1] = TargetSpace(FieldVector<double,3>({0,1,0}));
coefficients[2] = TargetSpace(FieldVector<double,3>({0,0,1}));
- typedef P1LocalFiniteElement<double,double,dim> LocalFiniteElement;
- typedef LocalGeodesicFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
- typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalProjectedFEFunctionType;
- typedef GFE::LocalTangentFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalTangentFEFunctionType;
+ typedef LocalGeodesicFEFunction<dim, double, P1LocalFiniteElement<double,double,dim>, TargetSpace> P1LocalGFEFunctionType;
+ //typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalProjectedFEFunctionType;
+ //typedef GFE::LocalTangentFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalTangentFEFunctionType;
+
+ P1LocalFiniteElement<double,double,dim> localFiniteElement;
+ P1LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,coefficients);
+ interpolate<TargetSpace, P1LocalGFEFunctionType>(localGeodesicFEFunction, 0, "riemannian-p1");
+
+ //interpolate<TargetSpace,LocalProjectedFEFunctionType>(coefficients, "projected");
+ //interpolate<TargetSpace,LocalTangentFEFunctionType>(coefficients, "tangent");
+
+ //////////////////////////////////////////////////////////////
+ // Set up a second-order interpolation function on the sphere
+ //////////////////////////////////////////////////////////////
+
+ coefficients.resize(6);
+
+ coefficients[0] = TargetSpace(FieldVector<double,3>({1,0,0}));
+ coefficients[1] = TargetSpace(FieldVector<double,3>({0.5,0.5,0}));
+ coefficients[2] = TargetSpace(FieldVector<double,3>({0,1,0}));
+ coefficients[3] = TargetSpace(FieldVector<double,3>({0.5,0,0.5}));
+ coefficients[4] = TargetSpace(FieldVector<double,3>({0,0.5,0.5}));
+ coefficients[5] = TargetSpace(FieldVector<double,3>({0,0,1}));
- interpolate<TargetSpace,LocalGFEFunctionType>(coefficients, "riemannian");
- interpolate<TargetSpace,LocalProjectedFEFunctionType>(coefficients, "projected");
- interpolate<TargetSpace,LocalTangentFEFunctionType>(coefficients, "tangent");
+ typedef LocalGeodesicFEFunction<dim, double, P2LocalFiniteElement<double,double,dim>, TargetSpace> P2LocalGFEFunctionType;
+ //typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalProjectedFEFunctionType;
+ //typedef GFE::LocalTangentFEFunction<dim, double, LocalFiniteElement, TargetSpace> LocalTangentFEFunctionType;
- std::vector<RealTuple<double,3> > flatCoefficients(3);
+ P2LocalFiniteElement<double,double,dim> p2LocalFiniteElement;
+ P2LocalGFEFunctionType p2LocalGeodesicFEFunction(p2LocalFiniteElement,coefficients);
+ interpolate<TargetSpace, P2LocalGFEFunctionType>(p2LocalGeodesicFEFunction, 0, "riemannian-p2-vertex");
+ interpolate<TargetSpace, P2LocalGFEFunctionType>(p2LocalGeodesicFEFunction, 1, "riemannian-p2-edge");
- flatCoefficients[0] = RealTuple<double,3>(FieldVector<double,3>({0,0,0}));
- flatCoefficients[1] = RealTuple<double,3>(FieldVector<double,3>({1,0,0}));
- flatCoefficients[2] = RealTuple<double,3>(FieldVector<double,3>({0,1,0}));
+ //interpolate<TargetSpace,LocalProjectedFEFunctionType>(coefficients, "projected");
+ //interpolate<TargetSpace,LocalTangentFEFunctionType>(coefficients, "tangent");
- typedef GFE::LocalProjectedFEFunction<dim, double, LocalFiniteElement, RealTuple<double,3> > RealTupleProjectedFEFunctionType;
- interpolate<RealTuple<double,3>, RealTupleProjectedFEFunctionType>(flatCoefficients, "flat");
} catch (Exception e) {
std::cout << e << std::endl;