Downsample higher-order Lagrangian functions to first-order ones.

That way, they can at least be looked at.  You won't see the intra-element
details, but frequently they are not even all this important.  Proper
subsampling of higher-order functions is more difficult and has to wait
for another day.

[[Imported from SVN: r9665]]

dune/gfe/cosseratvtkwriter.hh

@@ -63,6 +63,25 @@ class CosseratVTKWriter
 
     };
 
+    template <typename Basis1, typename Basis2>
+    static void downsample(const Basis1& basis1, const std::vector<RigidBodyMotion<double,3> >& v1,
+                           const Basis2& basis2,       std::vector<RigidBodyMotion<double,3> >& v2)
+    {
+      // Embed v1 into R^7
+      std::vector<Dune::FieldVector<double,7> > v1Embedded(v1.size());
+      for (size_t i=0; i<v1.size(); i++)
+        v1Embedded[i] = v1[i].globalCoordinates();
+
+      // Interpolate
+      BasisGridFunction<Basis1, std::vector<Dune::FieldVector<double,7> > > function(basis1, v1Embedded);
+      std::vector<Dune::FieldVector<double,7> > v2Embedded;
+      Functions::interpolate(basis2, v2Embedded, function);
+
+      // Copy back from R^7 into RigidBodyMotions
+      v2.resize(v2Embedded.size());
+      for (size_t i=0; i<v2.size(); i++)
+        v2[i] = RigidBodyMotion<double,3>(v2Embedded[i]);
+    }
     
 public:
 
@@ -120,15 +139,30 @@ public:
 
         const GridType& grid = basis.getGridView().grid();
 
+        //////////////////////////////////////////////////////////////////////////////////
+        //  Downsample the function onto a P1-space.  That's all we can visualize today.
+        //////////////////////////////////////////////////////////////////////////////////
+
+        typedef P1NodalBasis<GridView,double> P1Basis;
+        P1Basis p1Basis(basis.getGridView());
+
+        std::vector<RigidBodyMotion<double,3> > downsampledConfig;
+
+        downsample(basis, configuration, p1Basis, downsampledConfig);
+
+        //////////////////////////////////////////////////////////////////////////////////
+        //  Deform the grid according to the position information in 'configuration'
+        //////////////////////////////////////////////////////////////////////////////////
+
         typedef Dune::GeometryGrid<GridType,DeformationFunction<GridView> > DeformedGridType;
 
-        DeformationFunction<typename GridType::LeafGridView> deformationFunction(grid.leafGridView(), configuration);
+        DeformationFunction<typename GridType::LeafGridView> deformationFunction(grid.leafGridView(), downsampledConfig);
 
         // stupid, can't instantiate deformedGrid with a const grid
         DeformedGridType deformedGrid(const_cast<GridType&>(grid), deformationFunction);
 
-        typedef P1NodalBasis<typename DeformedGridType::LeafGridView,double> P1Basis;
-        P1Basis p1Basis(deformedGrid.leafGridView());
+        typedef P1NodalBasis<typename DeformedGridType::LeafGridView,double> DeformedP1Basis;
+        DeformedP1Basis deformedP1Basis(deformedGrid.leafGridView());
 
         Dune::VTKWriter<typename DeformedGridType::LeafGridView> vtkWriter(deformedGrid.leafGridView());
 
@@ -139,16 +173,16 @@ public:
 
         for (int i=0; i<3; i++) {
 
-            directors[i].resize(configuration.size());
-            for (size_t j=0; j<configuration.size(); j++)
-                directors[i][j] = configuration[j].q.director(i);
+            directors[i].resize(downsampledConfig.size());
+            for (size_t j=0; j<downsampledConfig.size(); j++)
+                directors[i][j] = downsampledConfig[j].q.director(i);
 
             std::stringstream iAsAscii;
             iAsAscii << i;
 
-            Dune::shared_ptr<VTKBasisGridFunction<P1Basis,CoefficientType> > vtkDirector
-               = Dune::shared_ptr<VTKBasisGridFunction<P1Basis,CoefficientType> >
-                   (new VTKBasisGridFunction<P1Basis,CoefficientType>(p1Basis, directors[i], "director"+iAsAscii.str()));
+            Dune::shared_ptr<VTKBasisGridFunction<DeformedP1Basis,CoefficientType> > vtkDirector
+               = Dune::shared_ptr<VTKBasisGridFunction<DeformedP1Basis,CoefficientType> >
+                   (new VTKBasisGridFunction<DeformedP1Basis,CoefficientType>(deformedP1Basis, directors[i], "director"+iAsAscii.str()));
             vtkWriter.addVertexData(vtkDirector);
         }