Commit 1a8a1c71 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Extend and improve QuadraticDirectionalRestriction

* Export origin and direction
* Implement `operator()`. This requires computation of the constant
  part which is done lazily, because it's not needed in many cases
  to find minimizers.
parent 8e497364
......@@ -152,7 +152,11 @@ public:
using Range = R;
QuadraticDirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const GlobalVector& origin, const GlobalVector& direction)
QuadraticDirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const GlobalVector& origin, const GlobalVector& direction) :
globalMatrix_(&matrix),
globalLinearTerm_(&linearTerm),
origin_(&origin),
direction_(&direction)
{
GlobalVector temp = linearTerm;
matrix.mv(direction, temp);
......@@ -161,9 +165,15 @@ public:
linearPart_ = linearTerm*direction - temp*origin;
}
Range operator()(const Vector& v) const
Range operator()(const Vector& z) const
{
DUNE_THROW(Dune::NotImplemented, "Evaluation of QuadraticDirectionalRestriction not implemented");
if (not constantPart_.has_value())
{
GlobalVector temp = *origin_;
globalMatrix_->mv(*origin_, temp);
constantPart_ = .5 * (temp * (*origin_));
}
return 0.5 * quadraticPart_*z*z - linearPart_*z + constantPart_.value();
}
const Matrix& quadraticPart() const
......@@ -176,9 +186,24 @@ public:
return linearPart_;
}
const GlobalVector& origin() const
{
return *origin_;
}
const GlobalVector& direction() const
{
return *direction_;
}
protected:
Matrix quadraticPart_;
Vector linearPart_;
const GlobalMatrix* globalMatrix_;
const GlobalVector* globalLinearTerm_;
const GlobalVector* origin_;
const GlobalVector* direction_;
mutable std::optional<R> constantPart_;
};
......
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