Don't ever store the derivative -- only compute its norm

The main reason for this is that it is apparently needed for situations
where the grid dimension is not the world dimension.  It is a bit
shorter two.  I was hoping that it would also be faster, because we
eliminate an intermediate matrix object.  However, preliminary tests
have not shown any differences.

dune/gfe/harmonicenergystiffness.hh

@@ -58,16 +58,16 @@ energy(const typename Basis::LocalView& localView,
         // The derivative of the local function defined on the reference element
         auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos);
 
-        // The derivative of the function defined on the actual element
-        typename LocalInterpolationRule::DerivativeType 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();
+        // Compute the Frobenius norm squared of the derivative.  This is the correct term
+        // if both domain and target space use the metric inherited from an embedding.
+        for (size_t i=0; i<jacobianInverseTransposed.N(); i++)
+          for (int j=0; j<TargetSpace::embeddedDim; j++)
+          {
+            RT entry = 0;
+            for (size_t k=0; k<jacobianInverseTransposed.M(); k++)
+              entry += jacobianInverseTransposed[i][k] * referenceDerivative[j][k];
+            energy += weight * entry * entry;
+          }
 
     }