From ce4902bdc126d2df0d26d98818768eebdfa58a14 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 18 Apr 2011 10:36:06 +0000
Subject: [PATCH] Implement the derivative of a gfe function value wrt the
coefficients.
Also: fix more bugs in the implementation of the derivative of gfe
function gradients wrt the coefficients. This appears to work now,
but inexact evaluations remain an issue.
[[Imported from SVN: r7180]]
---
dune/gfe/localgeodesicfefunction.hh | 117 ++++++++++++++++++++++++++--
1 file changed, 109 insertions(+), 8 deletions(-)
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index a33dd594..e624d897 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -87,6 +87,16 @@ public:
/** \brief For debugging: Evaluate the derivative of the function using a finite-difference approximation*/
Dune::FieldMatrix<ctype, EmbeddedTangentVector::size, dim> evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const;
+ /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+ void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
+ int coefficient,
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim>& derivative) const;
+
+ /** \brief Evaluate the derivative of the function value with respect to a coefficient */
+ void evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
+ int coefficient,
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim>& derivative) const;
+
/** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
void evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
@@ -158,6 +168,15 @@ private:
return dFdw;
}
+ Tensor3<ctype,embeddedDim,embeddedDim,embeddedDim> computeDqDqF(const std::vector<ctype>& w, const TargetSpace& q) const
+ {
+ Tensor3<ctype,embeddedDim,embeddedDim,embeddedDim> result;
+ result = 0;
+ for (int i=0; i<dim+1; i++)
+ result.axpy(w[i], TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i],q));
+ return result;
+ }
+
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
@@ -283,6 +302,66 @@ evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
return result;
}
+template <int dim, class ctype, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
+ int coefficient,
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim>& result) const
+{
+ // the function value at the point where we are evaluating the derivative
+ TargetSpace q = evaluate(local);
+
+ // dFdq
+ std::vector<ctype> w = barycentricCoordinates(local);
+ AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
+ assembler.assembleEmbeddedHessian(q,dFdq);
+
+ // dFDq is not invertible, if the target space is embedded into a higher-dimensional
+ // Euclidean space. Therefore we use its pseudo inverse. I don't think that is the
+ // best way, though.
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq);
+
+ //
+ result = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
+ result *= -w[coefficient];
+
+ result = result * dFdqPseudoInv;
+}
+
+template <int dim, class ctype, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
+ int coefficient,
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim>& result) const
+{
+ double eps = 1e-6;
+ for (int j=0; j<TargetSpace::EmbeddedTangentVector::size; j++) {
+
+ std::vector<TargetSpace> cornersPlus = coefficients_;
+ std::vector<TargetSpace> cornersMinus = coefficients_;
+ typename TargetSpace::CoordinateType aPlus = coefficients_[coefficient].globalCoordinates();
+ typename TargetSpace::CoordinateType aMinus = coefficients_[coefficient].globalCoordinates();
+ aPlus[j] += eps;
+ aMinus[j] -= eps;
+ cornersPlus[coefficient] = TargetSpace(aPlus);
+ cornersMinus[coefficient] = TargetSpace(aMinus);
+ LocalGeodesicFEFunction<dim,double,TargetSpace> fPlus(cornersPlus);
+ LocalGeodesicFEFunction<dim,double,TargetSpace> fMinus(cornersMinus);
+
+ TargetSpace hPlus = fPlus.evaluate(local);
+ TargetSpace hMinus = fMinus.evaluate(local);
+
+ result[j] = hPlus.globalCoordinates();
+ result[j] -= hMinus.globalCoordinates();
+ result[j] /= 2*eps;
+
+ }
+
+}
+
+
template <int dim, class ctype, class TargetSpace>
void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
@@ -325,20 +404,42 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
dvDqF = w[coefficient] * dvDqF;
// Put it all together
-#if 0
- for (size_t i=0; i<result.size(); i++)
- result[i] = dFdqPseudoInv * ( dvDqF[i] * dFdqPseudoInv * dFdw - dvDwF[i]) * B;
-#else
- Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::size, dim> derivative = evaluateDerivative(local);
- Tensor3<double, TargetSpace::EmbeddedTangentVector::size,TargetSpace::EmbeddedTangentVector::size,dim> foo;
- foo = -1 * dvDwF * B - dvDqF*derivative;
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim> dvq;
+ evaluateDerivativeOfValueWRTCoefficient(local,coefficient,dvq);
+
+ Dune::FieldMatrix<ctype, embeddedDim, dim> derivative = evaluateDerivative(local);
+
+ Tensor3<double, embeddedDim,embeddedDim,embeddedDim> dqdqF;
+ dqdqF = computeDqDqF(w,q);
+
+ Tensor3<double, embeddedDim,embeddedDim,dim+1> dqdwF;
+
+ for (int k=0; k<dim+1; k++) {
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
+ for (int i=0; i<embeddedDim; i++)
+ for (int j=0; j<embeddedDim; j++)
+ dqdwF[i][j][k] = hesse[i][j];
+ }
+
+ Tensor3<double, embeddedDim,embeddedDim,dim+1> dqdwF_times_dvq;
+ for (int i=0; i<embeddedDim; i++)
+ for (int j=0; j<embeddedDim; j++)
+ for (int k=0; k<dim+1; k++) {
+ dqdwF_times_dvq[i][j][k] = 0;
+ for (int l=0; l<embeddedDim; l++)
+ dqdwF_times_dvq[i][j][k] += dqdwF[l][j][k] * dvq[l][i];
+ }
+
+ Tensor3<double, embeddedDim,embeddedDim,dim> foo;
+ foo = -1 * dvDqF*derivative - (dvq*dqdqF)*derivative - dvDwF * B - dqdwF_times_dvq*B;
+
result = 0;
for (int i=0; i<embeddedDim; i++)
for (int j=0; j<embeddedDim; j++)
for (int k=0; k<dim; k++)
for (int l=0; l<embeddedDim; l++)
result[i][j][k] += dFdqPseudoInv[j][l] * foo[i][l][k];
-#endif
+
}
--
GitLab