From f901bfba01b3da7f24c21ef40f68aa4c2589be97 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 24 May 2018 11:44:32 +0200
Subject: [PATCH] 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 | 20 ++++++++++----------
1 file changed, 10 insertions(+), 10 deletions(-)
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index fa3cce66..253bd102 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/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;
+ }
}
--
GitLab