diff --git a/src/localgeodesicfefunction.hh b/src/localgeodesicfefunction.hh
index db9112ad4cdca3a4bf9c6ea772030541d4820042..e1e50aae9c5acd12bfc9c9f48414de8db10554b5 100644
--- a/src/localgeodesicfefunction.hh
+++ b/src/localgeodesicfefunction.hh
@@ -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)