From 7be1ade759c3f3b9b059430329d79cee207eb778 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 19 Oct 2010 14:39:48 +0000
Subject: [PATCH] =?UTF-8?q?Der=20Code=20zum=20analytischen=20Berechnen=20d?=
=?UTF-8?q?es=20Gradienten=20des=20Energiefunktionals=20=C3=BCbersetzt=20j?=
=?UTF-8?q?etzt=20und=20st=C3=BCrzt=20nicht=20ab.=20=20Es=20kommt=20allerd?=
=?UTF-8?q?ings=20noch=20Unfug=20raus.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
[[Imported from SVN: r6451]]
---
dune/gfe/harmonicenergystiffness.hh | 19 ++++---
dune/gfe/localgeodesicfefunction.hh | 81 ++++++++++++++++-------------
2 files changed, 57 insertions(+), 43 deletions(-)
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 0113dafc..d111784d 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -1,6 +1,8 @@
#ifndef HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
#define HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
+//#define HARMONIC_ENERGY_FD_GRADIENT
+
#include <dune/common/fmatrix.hh>
#include <dune/grid/common/quadraturerules.hh>
@@ -27,12 +29,13 @@ public:
/** \brief Assemble the energy for a single element */
RT energy (const Entity& e,
const std::vector<TargetSpace>& localSolution) const;
-
+
+#ifndef HARMONIC_ENERGY_FD_GRADIENT
/** \brief Assemble the gradient of the energy functional on one element */
virtual void assembleEmbeddedGradient(const Entity& element,
const std::vector<TargetSpace>& solution,
std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
-
+#endif
};
template <class GridView, class TargetSpace>
@@ -97,7 +100,7 @@ energy(const Entity& element,
return 0.5 * energy;
}
-
+#ifndef HARMONIC_ENERGY_FD_GRADIENT
template <class GridView, class TargetSpace>
void HarmonicEnergyLocalStiffness<GridView, TargetSpace>::
assembleEmbeddedGradient(const Entity& element,
@@ -141,7 +144,7 @@ assembleEmbeddedGradient(const Entity& element,
// loop over all the element's degrees of freedom and compute the gradient wrt it
for (size_t i=0; i<localSolution.size(); i++) {
- Dune::array<Dune::FieldMatrix<double,gridDim,TargetSpace::EmbeddedTangentVector::size>, TargetSpace::EmbeddedTangentVector::size> derivativeDerivative;
+ Tensor3<double, TargetSpace::EmbeddedTangentVector::size, gridDim,TargetSpace::EmbeddedTangentVector::size> derivativeDerivative;
localGeodesicFEFunction.evaluateDerivativeOfGradientWRTCoefficient(quadPos, i, derivativeDerivative);
for (int j=0; j<derivative.rows; j++) {
@@ -154,11 +157,11 @@ assembleEmbeddedGradient(const Entity& element,
}
-
-
}
- }
-}
+ }
+
+}
+#endif
#endif
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 1f5b4943..c2fe8147 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -44,6 +44,17 @@ Dune::FieldMatrix<K,m,n> operator- ( const Dune::FieldMatrix<K, m, n> &A, const
return ret;
}
+#if 0
+template< class K, int m, int n>
+void transpose(Dune::FieldMatrix<K, m, n> &A)
+{
+ for( size_type i = 0; i < m; ++i )
+ for( size_type j = 0; j < i; ++j )
+ std::swap(A[i][j], A[j][i]);
+}
+#endif
+
+
/** \brief A function defined by simplicial geodesic interpolation
from the reference element to a Riemannian manifold.
@@ -56,7 +67,7 @@ class LocalGeodesicFEFunction
{
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
- static const int targetDim = EmbeddedTangentVector::size;
+ static const int embeddedDim = EmbeddedTangentVector::size;
public:
@@ -77,7 +88,7 @@ public:
/** \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,
- Dune::array<Dune::FieldMatrix<double,dim,TargetSpace::EmbeddedTangentVector::size>, TargetSpace::EmbeddedTangentVector::size>& result) const;
+ Tensor3<double, TargetSpace::EmbeddedTangentVector::size,dim,TargetSpace::EmbeddedTangentVector::size>& result) const;
private:
@@ -103,22 +114,22 @@ private:
return result;
}
- static Dune::FieldMatrix<double,targetDim,targetDim> pseudoInverse(const Dune::FieldMatrix<double,targetDim,targetDim>& A)
+ static Dune::FieldMatrix<double,embeddedDim,embeddedDim> pseudoInverse(const Dune::FieldMatrix<double,embeddedDim,embeddedDim>& A)
{
- Dune::FieldMatrix<double,targetDim,targetDim> U = A;
- Dune::FieldVector<double,targetDim> W;
- Dune::FieldMatrix<double,targetDim,targetDim> V;
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim> U = A;
+ Dune::FieldVector<double,embeddedDim> W;
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim> V;
svdcmp(U,W,V);
// pseudoInv = V W^{-1} U^T
- Dune::FieldMatrix<double,targetDim,targetDim> UT;
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim> UT;
- for (int i=0; i<targetDim; i++)
- for (int j=0; j<targetDim; j++)
+ for (int i=0; i<embeddedDim; i++)
+ for (int j=0; j<embeddedDim; j++)
UT[i][j] = U[j][i];
- for (int i=0; i<targetDim; i++) {
+ for (int i=0; i<embeddedDim; i++) {
if (std::abs(W[i]) > 1e-12) // Diagonal may be zero, that's why we're using the pseudo inverse
UT[i] /= W[i];
else
@@ -162,16 +173,14 @@ template <int dim, class ctype, class TargetSpace>
Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::size, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
- const int targetDim = EmbeddedTangentVector::size;
-
- Dune::FieldMatrix<ctype, targetDim, dim> result;
+ Dune::FieldMatrix<ctype, embeddedDim, dim> result;
#if 0 // this is probably faster than the general implementation, but we leave it out for testing purposes
if (dim==1) {
EmbeddedTangentVector tmp = TargetSpace::interpolateDerivative(coefficients_[0], coefficients_[1], local[0]);
- for (int i=0; i<targetDim; i++)
+ for (int i=0; i<embeddedDim; i++)
result[i][0] = tmp[i];
}
@@ -189,23 +198,23 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
Dune::FieldMatrix<ctype,dim+1,dim> B = referenceToBarycentricLinearPart();
// compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
- Dune::FieldMatrix<ctype,targetDim,dim+1> dFdw;
+ Dune::FieldMatrix<ctype,embeddedDim,dim+1> dFdw;
for (int i=0; i<dim+1; i++) {
- Dune::FieldVector<ctype,targetDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
- for (int j=0; j<targetDim; j++)
+ Dune::FieldVector<ctype,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
+ for (int j=0; j<embeddedDim; j++)
dFdw[j][i] = tmp[j];
}
dFdw *= -1;
// multiply the two previous matrices: the result is the right hand side
- Dune::FieldMatrix<ctype,targetDim,dim> RHS = dFdw * B;
+ Dune::FieldMatrix<ctype,embeddedDim,dim> RHS = dFdw * B;
// the actual system matrix
std::vector<ctype> w = barycentricCoordinates(local);
AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
- Dune::FieldMatrix<ctype,targetDim,targetDim> dFdq(0);
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleHessian(q,dFdq);
// ////////////////////////////////////
@@ -215,18 +224,18 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
// 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,targetDim,targetDim> dFdqPseudoInv = pseudoInverse(dFdq);
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq);
for (int i=0; i<dim; i++) {
- Dune::FieldVector<ctype,targetDim> rhs, x;
- for (int j=0; j<targetDim; j++)
+ Dune::FieldVector<ctype,embeddedDim> rhs, x;
+ for (int j=0; j<embeddedDim; j++)
rhs[j] = RHS[j][i];
//dFdq.solve(x, rhs);
dFdqPseudoInv.mv(rhs,x);
- for (int j=0; j<targetDim; j++)
+ for (int j=0; j<embeddedDim; j++)
result[j][i] = x[j];
}
@@ -265,9 +274,9 @@ template <int dim, class ctype, class TargetSpace>
void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
- Dune::array<Dune::FieldMatrix<double,dim,TargetSpace::EmbeddedTangentVector::size>, TargetSpace::EmbeddedTangentVector::size>& result) const
+ Tensor3<double, TargetSpace::EmbeddedTangentVector::size,dim,TargetSpace::EmbeddedTangentVector::size>& result) const
{
- const int targetDim = EmbeddedTangentVector::size;
+ const int embeddedDim = EmbeddedTangentVector::size;
// the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
@@ -276,10 +285,10 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
Dune::FieldMatrix<ctype,dim+1,dim> B = referenceToBarycentricLinearPart();
// compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w
- Dune::FieldMatrix<ctype,targetDim,dim+1> dFdw;
+ Dune::FieldMatrix<ctype,embeddedDim,dim+1> dFdw;
for (int i=0; i<dim+1; i++) {
- Dune::FieldVector<ctype,targetDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
- for (int j=0; j<targetDim; j++)
+ Dune::FieldVector<ctype,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
+ for (int j=0; j<embeddedDim; j++)
dFdw[j][i] = tmp[j];
}
@@ -287,29 +296,31 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
std::vector<ctype> w = barycentricCoordinates(local);
AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
- Dune::FieldMatrix<ctype,targetDim,targetDim> dFdq(0);
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleHessian(q,dFdq);
- Tensor3<double,dim+1,targetDim,dim+1> dcDwF;
- for (size_t i=0; i<dcDwF.size(); i++)
- dcDwF[i] = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[i], q);
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> mixedDerivative = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
+ Tensor3<double,embeddedDim,embeddedDim,dim+1> dpDwF(0);
+ for (int i=0; i<embeddedDim; i++)
+ for (int j=0; j<embeddedDim; j++)
+ dpDwF[i][j][coefficient] = mixedDerivative[i][j];
// 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,targetDim,targetDim> dFdqPseudoInv = pseudoInverse(dFdq);
+ Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq);
// Put it all together
for (size_t i=0; i<result.size(); i++) {
//
- Tensor3<double,targetDim,targetDim,targetDim> dcDqF
+ Tensor3<double,embeddedDim,embeddedDim,embeddedDim> dcDqF
= TargetSpace::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(coefficients_[i], q);
- result[i] = dFdqPseudoInv * ( dcDqF[i] * dFdqPseudoInv * dFdw - dcDwF[i]) * B;
+ result[i] = dFdqPseudoInv * ( dcDqF[i] * dFdqPseudoInv * dFdw - dpDwF[i]) * B;
}
--
GitLab