Commit d1a5767b authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

Extend and improve IntegralFunctionalDirectionalRestriction

* Export origin and direction
* Cache and use evaluation matrix times origin
* Implement `operator()`
parent b7e0f05b
......@@ -384,13 +384,31 @@ public:
// Premultiply evaluationMatrix and direction vector.
// This is used later in the subDifferential method.
f_->evaluationMatrix().mv(*direction_, evaluationMatrixTimesDirection_);
f_->evaluationMatrix().mv(*origin_, evaluationMatrixTimesOrigin_);
}
Range operator()(const Scalar& v) const
/** \brief Compute the value of the functional a given point
*
* \param z Where to evaluate the subdifferential
* (The coordinate is a single number because we are in a one-dimensional space).
*/
Range operator()(const Scalar& z) const
{
DUNE_THROW(NotImplemented, "DirectionalRestriction::operator()");
Range result(0.0);
return result;
// The global position of where the directional derivative
// is to be evaluated is v = origin_ + z*direction. However,
// in the following we only need B*(origin + z*direction) = Bv
auto BXv = evaluationMatrixTimesDirection_;
BXv *= z;
BXv += evaluationMatrixTimesOrigin_;
auto&& density = f_->density();
auto&& weights = f_->quadratureWeights();
Range value = 0;
for (std::size_t i=0; i<f_->evaluationMatrix().N(); ++i)
value += density(BXv[i])*weights[i];
return value;;
}
/** \brief The set of points where the functional takes finite values
......@@ -415,26 +433,21 @@ public:
*/
Solvers::Interval<double> subDifferential(Scalar z) const
{
// The global position of where the directional derivative is to be evaluated
auto evaluationPoint = *origin_;
evaluationPoint.axpy(z,*direction_);
GlobalVector gradient;
Solvers::resizeInitializeZero(gradient,(*origin_));
// The global position of where the directional derivative
// is to be evaluated is v = origin_ + z*direction. However,
// in the following we only need B*(origin + z*direction) = Bv
auto BXv = evaluationMatrixTimesDirection_;
BXv *= z;
BXv += evaluationMatrixTimesOrigin_;
typename Functional::Density::Domain BiXv;
auto densityDerivative = derivative(f_->density());
auto&& weights = f_->quadratureWeights();
double directionalDerivative = 0;
for (std::size_t i=0; i<f_->evaluationMatrix().N(); ++i)
{
BiXv = 0.0;
Functional::addRowXVector(f_->evaluationMatrix()[i], evaluationPoint, BiXv);
auto densityDerivative = derivative(f_->density());
auto Dgamma = densityDerivative(BiXv);
Dgamma *= f_->quadratureWeights()[i];
auto Dgamma = densityDerivative(BXv[i]);
Dgamma *= weights[i];
directionalDerivative += Dgamma * evaluationMatrixTimesDirection_[i];
}
......@@ -451,15 +464,26 @@ public:
s = subDifferential(z);
}
const GlobalVector& origin() const
{
return *origin_;
}
const GlobalVector& direction() const
{
return *direction_;
}
protected:
const Functional* f_;
const GlobalVector* origin_;
const GlobalVector* direction_;
// The evaluation matrix multiplied by the direction vector
// The evaluation matrix multiplied by the origin and direction vector
// This product is used in the subDifferential method,
// and it pays to precompute it.
BlockVector<typename Functional::Density::Domain> evaluationMatrixTimesDirection_;
BlockVector<typename Functional::Density::Domain> evaluationMatrixTimesOrigin_;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment