Visualize test function vector fields for second order test functions

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;