Introduce new class CosseratEnergyLocalStiffness.

It will implement the energy of the elastic Cosserat modell
for continua and shells proposed by Patrizio Neff.  So far
it is only a copy of the harmonic energy, though.

[[Imported from SVN: r7417]]

cosserat-continuum.cc

@@ -2,7 +2,6 @@
 
 #include <fenv.h>
 
-//#define LAPLACE_DEBUG
 //#define HARMONIC_ENERGY_FD_GRADIENT
 
 #define RIGIDBODYMOTION3
@@ -25,7 +24,7 @@
 #include <dune/gfe/rotation.hh>
 #include <dune/gfe/unitvector.hh>
 #include <dune/gfe/realtuple.hh>
-#include <dune/gfe/harmonicenergystiffness.hh>
+#include <dune/gfe/cosseratenergystiffness.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
 
@@ -182,7 +181,7 @@ int main (int argc, char *argv[]) try
     //   Create an assembler for the Harmonic Energy Functional
     // ////////////////////////////////////////////////////////////
 
-    HarmonicEnergyLocalStiffness<GridType::LeafGridView,TargetSpace> harmonicEnergyLocalStiffness;
+    CosseratEnergyLocalStiffness<GridType::LeafGridView,TargetSpace> harmonicEnergyLocalStiffness;
 
     GeodesicFEAssembler<GridType::LeafGridView,TargetSpace> assembler(grid.leafView(),
                                                                       &harmonicEnergyLocalStiffness);

dune/gfe/cosseratenergystiffness.hh

+#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
+#define COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/quadraturerules.hh>
+
+#include "localgeodesicfestiffness.hh"
+#include "localgeodesicfefunction.hh"
+
+
+template<class GridView, class TargetSpace>
+class CosseratEnergyLocalStiffness 
+    : public LocalGeodesicFEStiffness<GridView,TargetSpace>
+{
+    // grid types
+    typedef typename GridView::Grid::ctype DT;
+    typedef typename TargetSpace::ctype RT;
+    typedef typename GridView::template Codim<0>::Entity Entity;
+    
+    // some other sizes
+    enum {gridDim=GridView::dimension};
+
+public:
+    
+    //! Dimension of a tangent space
+    enum { blocksize = TargetSpace::TangentVector::size };
+
+    /** \brief Assemble the energy for a single element */
+    RT energy (const Entity& e,
+               const std::vector<TargetSpace>& localSolution) const;
+
+};
+
+template <class GridView, class TargetSpace>
+typename CosseratEnergyLocalStiffness<GridView, TargetSpace>::RT CosseratEnergyLocalStiffness<GridView, TargetSpace>::
+energy(const Entity& element,
+       const std::vector<TargetSpace>& localSolution) const
+{
+    RT energy = 0;
+
+    assert(element.type().isSimplex());
+    
+    LocalGeodesicFEFunction<gridDim, double, TargetSpace> localGeodesicFEFunction(localSolution);
+
+    int quadOrder = 1;//gridDim;
+
+    const Dune::QuadratureRule<double, gridDim>& quad 
+        = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
+    
+    for (size_t pt=0; pt<quad.size(); pt++) {
+        
+        // Local position of the quadrature point
+        const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
+        
+        const double integrationElement = element.geometry().integrationElement(quadPos);
+
+        const Dune::FieldMatrix<double,gridDim,gridDim>& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+        
+        double weight = quad[pt].weight() * integrationElement;
+
+        // The derivative of the local function defined on the reference element
+        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+
+        // The derivative of the function defined on the actual element
+        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> derivative(0);
+
+        for (size_t comp=0; comp<referenceDerivative.N(); comp++)
+            jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
+
+        // Add the local energy density
+        // The Frobenius norm is the correct norm here if the metric of TargetSpace is the identity.
+        // (And if the metric of the domain space is the identity, which it always is here.)
+        energy += weight * derivative.frobenius_norm2();
+
+    }
+
+    return 0.5 * energy;
+}
+
+#endif
+