complete implementation, not tested beyond compilation yet, though

[[Imported from SVN: r4059]]

src/harmonicenergystiffness.hh

 #ifndef HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
 #define HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-#include <dune/disc/operators/localstiffness.hh>
+#include <dune/grid/common/quadraturerules.hh>
+
+#include "localgeodesicfestiffness.hh"
+#include "localgeodesicfefunction.hh"
 
 template<class GridView, class TargetSpace>
 class HarmonicEnergyLocalStiffness 
-    : public Dune::LocalGeodesicFEStiffness<GridView,TargetSpace>
+    : 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
@@ -23,14 +24,18 @@ public:
     //! Dimension of a tangent space
     enum { blocksize = TargetSpace::TangentVector::size };
 
+#if 0
     // types for matrics, vectors and boundary conditions
     typedef Dune::FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
     typedef Dune::FieldVector<RT,m> VBlockType;   // one entry in the global vectors
     typedef Dune::array<Dune::BoundaryConditions::Flags,m> BCBlockType;     // componentwise boundary conditions
+#endif
 
+#if 0
     //! Default Constructor
     HarmonicEnergyLocalStiffness ()
     {}
+#endif
 
     /** \brief Assemble the energy for a single element */
     RT energy (const Entity& e,
@@ -38,19 +43,19 @@ public:
 
 };
 
-template <class GridType, class RT>
-RT RodLocalStiffness<GridType, RT>::
+template <class GridView, class TargetSpace>
+typename HarmonicEnergyLocalStiffness<GridView, TargetSpace>::RT HarmonicEnergyLocalStiffness<GridView, TargetSpace>::
 energy(const Entity& element,
-       const Dune::array<RigidBodyMotion<3>,2>& localSolution) const
+       const std::vector<TargetSpace>& localSolution) const
 {
     RT energy = 0;
     
-    assert(element.type().isSimplex());
+    LocalGeodesicFEFunction<gridDim, double, TargetSpace> localGeodesicFEFunction(element.type(), localSolution);
 
     int quadOrder = gridDim;
 
-    const Dune::QuadratureRule<double, 1>& quad 
-        = Dune::QuadratureRules<double, 1>::rule(element.type(), quadOrder);
+    const Dune::QuadratureRule<double, gridDim>& quad 
+        = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
     
     for (size_t pt=0; pt<quad.size(); pt++) {
         
@@ -58,17 +63,24 @@ energy(const Entity& element,
         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;
 
-        for (int comp=0; comp<4; comp++) {
+        // 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 (int dir=0; dir<gridDim; dir++) {
+        for (int comp=0; comp<4; comp++)
+            jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
 
-                double derivative = localGeodesicFunction.evaluateDerivative(comp, dir);
-                energy += weight * derivative * derivative;
+        for (int comp=0; comp<4; comp++) {
 
-            }
+            for (int dir=0; dir<gridDim; dir++)
+                energy += weight * derivative[comp][dir] * derivative[comp][dir];
 
         }