implement frame-invariance test for the energy

[[Imported from SVN: r8055]]

test/cosseratenergytest.cc

@@ -4,6 +4,8 @@
 
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
+#include <dune/fufem/functions/constantfunction.hh>
+
 #include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/cosseratenergystiffness.hh>
 
@@ -102,10 +104,115 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
     }
 }
 
+//////////////////////////////////////////////////////////////////////////////////////
+//   Test invariance of the energy functional under rotations
+//////////////////////////////////////////////////////////////////////////////////////
+
+template <class GridType>
+void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients) {
+
+    PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
+    typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
+    
+    //LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
+
+    ParameterTree materialParameters;
+    materialParameters["thickness"] = "0.1";
+    materialParameters["mu"] = "3.8462e+05";
+    materialParameters["lambda"] = "2.7149e+05";
+    materialParameters["mu_c"] = "3.8462e+05";
+    materialParameters["L_c"] = "0.1";
+    materialParameters["q"] = "2.5";
+
+    ConstantFunction<Dune::FieldVector<double,GridType::dimension>, Dune::FieldVector<double,3> > zeroFunction(Dune::FieldVector<double,3>(0));
+    
+    CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3> assembler(materialParameters,
+                                                                                                 NULL,
+                                                                                                 &zeroFunction);
+    std::vector<TargetSpace> rotatedCoefficients(coefficients.size());
+    
+    std::vector<Rotation<3> > testRotations;
+    ValueFactory<Rotation<3> >::get(testRotations);
 
+    for (size_t i=0; i<testRotations.size(); i++) {
 
+        /////////////////////////////////////////////////////////////////////////
+        //  Multiply the given configuration by the test rotation.
+        //  The energy should remain unchanged.
+        /////////////////////////////////////////////////////////////////////////
+        FieldMatrix<double,3,3> matrix;
+        testRotations[i].matrix(matrix);
+        
+        for (size_t j=0; j<coefficients.size(); j++) {
+            FieldVector<double,3> tmp;
+            matrix.mv(coefficients[j].r, tmp);
+            rotatedCoefficients[j].r = tmp;
+            
+            rotatedCoefficients[j].q = testRotations[i].mult(rotatedCoefficients[j].q);
+        }
+        
+        std::cout << "energy: " << assembler.energy(*grid->template leafbegin<0>(), 
+                                                    feCache.get(grid->template leafbegin<0>()->type()),
+                                                    rotatedCoefficients) << std::endl;
 
-int main(int argc, char** argv)
+    }
+
+}
+
+
+template <int domainDim>
+void testFrameInvariance()
+{
+    // ////////////////////////////////////////////////////////
+    //   Make a test grid consisting of a single simplex
+    // ////////////////////////////////////////////////////////
+
+    typedef UGGrid<domainDim> GridType;
+
+    GridFactory<GridType> factory;
+
+    FieldVector<double,dim> pos(0);
+    factory.insertVertex(pos);
+
+    for (int i=0; i<domainDim+1; i++) {
+        pos = 0;
+        pos[i] = 1;
+        factory.insertVertex(pos);
+    }
+
+    std::vector<unsigned int> v(domainDim+1);
+    for (int i=0; i<domainDim+1; i++)
+        v[i] = i;
+    factory.insertElement(GeometryType(GeometryType::simplex,dim), v);
+
+    const GridType* grid = factory.createGrid();
+    
+    // //////////////////////////////////////////////////////////
+    //  Test whether the energy is invariant under isometries
+    // //////////////////////////////////////////////////////////
+
+    std::vector<TargetSpace> testPoints;
+    ValueFactory<TargetSpace>::get(testPoints);
+
+    // Set up elements of SE(3)
+    std::vector<TargetSpace> coefficients(dim+1);
+
+    MultiIndex index(dim+1, testPoints.size());
+    int numIndices = index.cycle();
+
+    for (int i=0; i<numIndices; i++, ++index) {
+        
+        for (int j=0; j<dim+1; j++)
+            coefficients[j] = testPoints[index[j]];
+
+        testEnergy<GridType>(grid, coefficients);
+        
+    }
+    
+}
+
+
+int main(int argc, char** argv) try
 {
     const int domainDim = 2;
     std::cout << " --- Testing Rotation<3>, domain dimension: " << domainDim << " ---" << std::endl;
@@ -131,4 +238,14 @@ int main(int argc, char** argv)
                 
     }
     
-}
+    //////////////////////////////////////////////////////////////////////////////////////
+    //   Test invariance of the energy functional under rotations
+    //////////////////////////////////////////////////////////////////////////////////////
+    
+    testFrameInvariance<2>();
+    
+} catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }