Skip to content
Snippets Groups Projects
Commit a1660a23 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

bugfix: we need the derivative of the geodesic interpolation in the gradient... 

bugfix: we need the derivative of the geodesic interpolation in the gradient computation (chain rule), not the derivative of the velocity interpolation

[[Imported from SVN: r4066]]
parent 509c33bf
No related branches found
No related tags found
Hide whitespace changes
Inline Side-by-side
src/localgeodesicfefunction.hh
+ 72
− 113
View file @ a1660a23
......@@ -53,8 +53,8 @@ private:
\todo This is group-specific and should not really be here
*/
static void interpolationVelocityDerivative(const Rotation<3,ctype>& q0, const Rotation<3,ctype>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad);
static void interpolationDerivative(const Rotation<3,ctype>& q0, const Rotation<3,ctype>& q1, double s,
Dune::array<Quaternion<double>,6>& grad);
/** \brief The type of the reference element */
......@@ -73,7 +73,7 @@ derivativeWRTFirstPoint(const Rotation<3,ctype>& a, const Rotation<3,ctype>& b,
// Get the derivative of [a,b](s) wrt 'a'
/** \brief The method actually compute the derivatives wrt to 'a' and 'b'. This is a waste! */
Dune::array<Quaternion<double>,6> grad;
interpolationVelocityDerivative(a,b,s,1,grad);
interpolationDerivative(a,b,s,grad);
// We are really only interested in the first three entries
......@@ -104,126 +104,84 @@ derivativeWRTFirstPoint(const Rotation<3,ctype>& a, const Rotation<3,ctype>& b,
template <int dim, class ctype, class TargetSpace>
void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
interpolationVelocityDerivative(const Rotation<3,ctype>& q0, const Rotation<3,ctype>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad)
interpolationDerivative(const Rotation<3,ctype>& q0, const Rotation<3,ctype>& q1, double s,
Dune::array<Quaternion<double>,6>& grad)
{
// Clear output array
for (int i=0; i<6; i++)
grad[i] = 0;
// Compute q_0^{-1}
Rotation<3,ctype> q0Inv = q0;
q0Inv.invert();
// Compute v = s \exp^{-1} ( q_0^{-1} q_1)
Dune::FieldVector<ctype,3> v = Rotation<3,ctype>::expInv(q0Inv.mult(q1));
v *= s/intervalLength;
Dune::FieldMatrix<ctype,4,3> dExp_v = Rotation<3,ctype>::Dexp(v);
Dune::array<Dune::FieldMatrix<ctype,3,3>, 4> ddExp;
Rotation<3,ctype>::DDexp(v, ddExp);
Dune::FieldMatrix<ctype,3,4> dExpInv = Rotation<3,ctype>::DexpInv(q0Inv.mult(q1));
Dune::FieldMatrix<ctype,4,4> mat(0);
for (int i=0; i<4; i++)
for (int j=0; j<4; j++)
for (int k=0; k<3; k++)
mat[i][j] += 1/intervalLength * dExp_v[i][k] * dExpInv[k][j];
// /////////////////////////////////////////////////
// The derivatives with respect to w^0
// /////////////////////////////////////////////////
for (int i=0; i<3; i++) {
// \partial exp \partial w^1_j at 0
Quaternion<ctype> dw;
for (int j=0; j<4; j++)
dw[j] = 0.5*(i==j); // dExp_v_0[j][i];
// \xi = \exp^{-1} q_0^{-1} q_1
Dune::FieldVector<ctype,3> xi = Rotation<3,ctype>::expInv(q0Inv.mult(q1));
Quaternion<ctype> addend0;
addend0 = 0;
dExp_v.umv(xi,addend0);
addend0 = dw.mult(addend0);
addend0 /= intervalLength;
// \parder{\xi}{w^1_j} = ...
Quaternion<ctype> dwConj = dw;
dwConj.conjugate();
//dwConj[3] -= 2 * dExp_v_0[3][i]; the last row of dExp_v_0 is zero
dwConj = dwConj.mult(q0Inv.mult(q1));
Dune::FieldVector<ctype,3> dxi(0);
Rotation<3,ctype>::DexpInv(q0Inv.mult(q1)).umv(dwConj, dxi);
Quaternion<ctype> vHv;
for (int j=0; j<4; j++) {
vHv[j] = 0;
// vHv[j] = dxi * DDexp * xi
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
vHv[j] += ddExp[j][k][l]*dxi[k]*xi[l];
// The derivatives with respect to w^1
// Compute q_1^{-1}q_0
Rotation<3,ctype> q1InvQ0 = q1;
q1InvQ0.invert();
q1InvQ0 = q1InvQ0.mult(q0);
{
// Compute v = (1-s) \exp^{-1} ( q_1^{-1} q_0)
Dune::FieldVector<ctype,3> v = Rotation<3,ctype>::expInv(q1InvQ0);
v *= (1-s);
Dune::FieldMatrix<ctype,4,3> dExp_v = Rotation<3,ctype>::Dexp(v);
Dune::FieldMatrix<ctype,3,4> dExpInv = Rotation<3,ctype>::DexpInv(q1InvQ0);
Dune::FieldMatrix<ctype,4,4> mat(0);
for (int i=0; i<4; i++)
for (int j=0; j<4; j++)
for (int k=0; k<3; k++)
mat[i][j] += (1-s) * dExp_v[i][k] * dExpInv[k][j];
for (int i=0; i<3; i++) {
Quaternion<ctype> dw;
for (int j=0; j<4; j++)
dw[j] = 0.5 * (i==j); // dExp[j][i] at v=0
dw = q1InvQ0.mult(dw);
mat.umv(dw,grad[i]);
grad[i] = q1.mult(grad[i]);
}
vHv *= s/intervalLength/intervalLength;
// Third addend
mat.umv(dwConj,grad[i]);
// add up
grad[i] += addend0;
grad[i] += vHv;
grad[i] = q0.mult(grad[i]);
}
// /////////////////////////////////////////////////
// The derivatives with respect to w^1
// /////////////////////////////////////////////////
for (int i=3; i<6; i++) {
// \partial exp \partial w^1_j at 0
Quaternion<ctype> dw;
for (int j=0; j<4; j++)
dw[j] = 0.5 * ((i-3)==j); // dw[j] = dExp_v_0[j][i-3];
// \xi = \exp^{-1} q_0^{-1} q_1
Dune::FieldVector<ctype,3> xi = Rotation<3,ctype>::expInv(q0Inv.mult(q1));
// \parder{\xi}{w^1_j} = ...
Dune::FieldVector<ctype,3> dxi(0);
dExpInv.umv(q0Inv.mult(q1.mult(dw)), dxi);
Quaternion<ctype> vHv;
for (int j=0; j<4; j++) {
// vHv[j] = dxi * DDexp * xi
vHv[j] = 0;
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
vHv[j] += ddExp[j][k][l]*dxi[k]*xi[l];
}
vHv *= s/intervalLength/intervalLength;
// ///////////////////////////////////
// second addend
// ///////////////////////////////////
// Compute q_0^{-1}
Rotation<3,ctype> q0InvQ1 = q0;
q0InvQ1.invert();
q0InvQ1 = q0InvQ1.mult(q1);
{
// Compute v = s \exp^{-1} ( q_0^{-1} q_1)
Dune::FieldVector<ctype,3> v = Rotation<3,ctype>::expInv(q0InvQ1);
v *= s;
dw = q0Inv.mult(q1.mult(dw));
Dune::FieldMatrix<ctype,4,3> dExp_v = Rotation<3,ctype>::Dexp(v);
mat.umv(dw,grad[i]);
grad[i] += vHv;
grad[i] = q0.mult(grad[i]);
Dune::FieldMatrix<ctype,3,4> dExpInv = Rotation<3,ctype>::DexpInv(q0InvQ1);
Dune::FieldMatrix<ctype,4,4> mat(0);
for (int i=0; i<4; i++)
for (int j=0; j<4; j++)
for (int k=0; k<3; k++)
mat[i][j] += s * dExp_v[i][k] * dExpInv[k][j];
for (int i=3; i<6; i++) {
Quaternion<ctype> dw;
for (int j=0; j<4; j++)
dw[j] = 0.5 * ((i-3)==j); // dExp[j][i-3] at v=0
dw = q0InvQ1.mult(dw);
mat.umv(dw,grad[i]);
grad[i] = q0.mult(grad[i]);
}
}
}
......@@ -288,6 +246,7 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local)
return result;
}
template <int dim, class ctype, class TargetSpace>
Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::size, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local)
......
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