diff --git a/src/samplefunctional.hh b/src/samplefunctional.hh index a32e96645407c98d8ac43d9387cccef888b568c8..23d9fd859e74293771c3c56aface31696949f473 100644 --- a/src/samplefunctional.hh +++ b/src/samplefunctional.hh @@ -39,7 +39,8 @@ template <int dim> class SampleFunctional { // negative_projection() divides by it if (x.two_norm2() == 0.0) { // If there is a direction of descent, this is it - SmallVector const d = smoothGradient(x); + SmallVector d; + smoothGradient(x, d); Interval<double> D; phi.directionalSubDiff(x, d, D); @@ -56,8 +57,10 @@ template <int dim> class SampleFunctional { return; } - SmallVector const pg = upperGradient(x); - SmallVector const mg = lowerGradient(x); + SmallVector pg; + upperGradient(x, pg); + SmallVector mg; + lowerGradient(x, mg); double const pgx = pg * x; double const mgx = mg * x; if (pgx >= 0 && mgx >= 0) { @@ -73,7 +76,9 @@ template <int dim> class SampleFunctional { } else { // Includes the case that pg points in direction x and mg // points in direction -x. The projection will then be zero. - ret = negative_projection(smoothGradient(x), x); + SmallVector d; + smoothGradient(x, d); + ret = negative_projection(d, x); dverb << "## Directional derivative (as per scalar product w/ " "semigradient): " << -(ret * ret) << " (coordinates of the restriction)" << std::endl; @@ -87,25 +92,23 @@ template <int dim> class SampleFunctional { private: // Gradient of the smooth part - SmallVector smoothGradient(SmallVector const x) const { - SmallVector y; + void smoothGradient(SmallVector const &x, SmallVector &y) const { A.mv(x, y); // y = Av y -= b; // y = Av - b - return y; } - SmallVector upperGradient(SmallVector const x) const { - SmallVector y; + void upperGradient(SmallVector const &x, SmallVector &y) const { phi.upperGradient(x, y); - y += smoothGradient(x); - return y; + SmallVector z; + smoothGradient(x, z); + y += z; } - SmallVector lowerGradient(SmallVector const x) const { - SmallVector y; + void lowerGradient(SmallVector const &x, SmallVector &y) const { phi.lowerGradient(x, y); - y += smoothGradient(x); - return y; + SmallVector z; + smoothGradient(x, z); + y += z; } // No normalising is done!