globally refine a grid and transfer a given function to the next level. Works...

globally refine a grid and transfer a given function to the next level.  Works for all TargetSpaces and replace rodRefine, which did the same only for RigidBodyMotions

[[Imported from SVN: r5771]]

src/Makefile.am

@@ -6,7 +6,7 @@ AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/..
 
 srcincludedir = $(includedir)/dune/common
 srcinclude_HEADERS = averagedistanceassembler.hh localgeodesicfefunction.hh \
-                     quaternion.hh rodrefine.hh averageinterface.hh localstiffness.hh \
+                     quaternion.hh geodesicfefunctionadaptor.hh averageinterface.hh localstiffness.hh \
                      localgeodesicfestiffness.hh riemanniantrsolver.hh rodwriter.hh \
                      geodesicdifference.hh makestraightrod.hh rigidbodymotion.hh \
                      rotation.hh geodesicfeassembler.hh maxnormtrustregion.hh \

src/rodrefine.hh → src/geodesicfefunctionadaptor.hh

-#ifndef ROD_REFINE_HH
-#define ROD_REFINE_HH
+#ifndef GEODESIC_FE_FUNCTION_ADAPTOR_HH
+#define GEODESIC_FE_FUNCTION_ADAPTOR_HH
 
 #include <vector>
 #include <map>
 
-#include "rigidbodymotion.hh"
+#include "localgeodesicfefunction.hh"
 
-template <class GridType>
-void globalRodRefine(GridType& grid, std::vector<RigidBodyMotion<3> >& x)
+template <class GridType, class TargetSpace>
+void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
 {
-    // Should be a static assertion
-    assert(GridType::dimension == 1);
+    const int dim = GridType::dimension;
 
     typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
-    typedef typename GridType::template Codim<1>::LeafIterator VertexIterator;
+    typedef typename GridType::template Codim<dim>::LeafIterator VertexIterator;
 
     // /////////////////////////////////////////////////////
     //   Save leaf p1 data in a map
     // /////////////////////////////////////////////////////
 
-    const typename GridType::Traits::LocalIdSet&    idSet    = grid.localIdSet();
+    const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
     const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
 
-    std::map<typename GridType::Traits::LocalIdSet::IdType, RigidBodyMotion<3> > dofMap;
+    std::map<typename GridType::Traits::LocalIdSet::IdType, TargetSpace> dofMap;
 
-    VertexIterator vIt    = grid.template leafbegin<1>();
-    VertexIterator vEndIt = grid.template leafend<1>();
+    VertexIterator vIt    = grid.template leafbegin<dim>();
+    VertexIterator vEndIt = grid.template leafend<dim>();
 
     for (; vIt!=vEndIt; ++vIt)
         dofMap.insert(std::make_pair(idSet.id(*vIt), x[indexSet.index(*vIt)]));
@@ -43,28 +42,36 @@ void globalRodRefine(GridType& grid, std::vector<RigidBodyMotion<3> >& x)
     //   Restore and interpolate the data
     // /////////////////////////////////////////////////////
 
-    x.resize(grid.size(1));
+    x.resize(grid.size(dim));
 
     ElementIterator eIt    = grid.template leafbegin<0>();
     ElementIterator eEndIt = grid.template leafend<0>();
 
     for (; eIt!=eEndIt; ++eIt) {
 
-        for (int i=0; i<2; i++) {
+        // Set up a local gfe function on the father element
+        std::vector<TargetSpace> coefficients(eIt->father()->template count<dim>());
 
-            if (dofMap.find(idSet.subId(*eIt,i,1)) != dofMap.end()) {
+        for (int i=0; i<eIt->father()->template count<dim>(); i++)
+            coefficients[i] = dofMap.find(idSet.subId(*eIt->father(),i,dim))->second;
 
-                x[indexSet.subIndex(*eIt,i,1)] = dofMap[idSet.subId(*eIt,i,1)];
+        LocalGeodesicFEFunction<dim,double,TargetSpace> fatherFunction(coefficients);
+
+        // The embedding of this element into the father geometry
+        const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
+
+        for (int i=0; i<eIt->template count<dim>(); i++) {
+
+            if (dofMap.find(idSet.subId(*eIt,i,dim)) != dofMap.end()) {
+
+                // If the vertex exists on the coarser level we take the value from there.
+                // That should be faster and more accurate than interpolating
+                x[indexSet.subIndex(*eIt,i,dim)] = dofMap[idSet.subId(*eIt,i,dim)];
 
             } else {
-                // Interpolate
-                RigidBodyMotion<3> p0 = dofMap[idSet.subId(*eIt->father(),0,1)];
-                RigidBodyMotion<3> p1 = dofMap[idSet.subId(*eIt->father(),1,1)];
 
-                x[indexSet.subIndex(*eIt,i,1)].r = (p0.r + p1.r);
-                x[indexSet.subIndex(*eIt,i,1)].r *= 0.5;
-                x[indexSet.subIndex(*eIt,i,1)].q
-                    = Rotation<3,double>::interpolate(p0.q, p1.q, 0.5);
+                // Interpolate
+                x[indexSet.subIndex(*eIt,i,dim)] = fatherFunction.evaluate(geometryInFather.corner(i));
 
             }