diff --git a/dune/gfe/cosseratamirameshwriter.hh b/dune/gfe/cosseratamirameshwriter.hh
index 19b6336de1982ad688277f1688a51c04ed8c1e56..458705a1c9cc0b919f2e2c5d3051d362fcb7123b 100644
--- a/dune/gfe/cosseratamirameshwriter.hh
+++ b/dune/gfe/cosseratamirameshwriter.hh
@@ -27,7 +27,7 @@ class CosseratAmiraMeshWriter
public:
DeformationFunction(const HostGridView& gridView,
- const std::vector<RigidBodyMotion<3> >& deformedPosition)
+ const std::vector<RigidBodyMotion<double,3> >& deformedPosition)
: gridView_(gridView),
deformedPosition_(deformedPosition)
{}
@@ -57,14 +57,14 @@ class CosseratAmiraMeshWriter
HostGridView gridView_;
- const std::vector<RigidBodyMotion<3> > deformedPosition_;
+ const std::vector<RigidBodyMotion<double,3> > deformedPosition_;
};
public:
static void write(const GridType& grid,
- const std::vector<RigidBodyMotion<3> >& configuration,
+ const std::vector<RigidBodyMotion<double,3> >& configuration,
const std::string& filePrefix)
{
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index cc72f3b64cb657ea6e3ff1f8063a7a3810dfa598..17cf5f1a88e88fc8d8b340ba9b828a7210e416d9 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -15,11 +15,11 @@
template<class GridView, class LocalFiniteElement, int dim>
class CosseratEnergyLocalStiffness
- : public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,RigidBodyMotion<dim> >
+ : public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,RigidBodyMotion<double,dim> >
{
// grid types
typedef typename GridView::Grid::ctype DT;
- typedef RigidBodyMotion<dim> TargetSpace;
+ typedef RigidBodyMotion<double,dim> TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -74,7 +74,7 @@ class CosseratEnergyLocalStiffness
public: // for testing
/** \brief Compute the derivative of the rotation, but wrt matrix coordinates */
- static void computeDR(const RigidBodyMotion<3>& value,
+ static void computeDR(const RigidBodyMotion<double,3>& value,
const Dune::FieldMatrix<double,7,gridDim>& derivative,
Tensor3<double,3,3,3>& DR)
{
@@ -201,7 +201,7 @@ typename CosseratEnergyLocalStiffness<GridView,LocalFiniteElement,dim>::RT
CosseratEnergyLocalStiffness<GridView,LocalFiniteElement,dim>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const std::vector<RigidBodyMotion<dim> >& localSolution) const
+ const std::vector<RigidBodyMotion<double,dim> >& localSolution) const
{
RT energy = 0;
@@ -226,7 +226,7 @@ energy(const Entity& element,
double weight = quad[pt].weight() * integrationElement;
// The value of the local function
- RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ RigidBodyMotion<double,dim> value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function defined on the reference element
Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
@@ -314,7 +314,7 @@ energy(const Entity& element,
const double integrationElement = it->geometry().integrationElement(quad[pt].position());
// The value of the local function
- RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ RigidBodyMotion<double,dim> value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
Dune::FieldVector<double,3> neumannValue;
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 6a19bf65c141a3f22f17c137164692bb486fa5f1..ae3238c29532ab5e05bead07c9be5c43715df97f 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -560,9 +560,9 @@ evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>
\tparam ctype Type used for coordinates on the reference element
*/
template <int dim, class ctype, class LocalFiniteElement>
-class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,RigidBodyMotion<3> >
+class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,RigidBodyMotion<double,3> >
{
- typedef RigidBodyMotion<3> TargetSpace;
+ typedef RigidBodyMotion<double,3> TargetSpace;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
@@ -582,11 +582,11 @@ public:
for (size_t i=0; i<coefficients.size(); i++)
translationCoefficients_[i] = coefficients[i].r;
- std::vector<Rotation<3,ctype> > orientationCoefficients(coefficients.size());
+ std::vector<Rotation<ctype,3> > orientationCoefficients(coefficients.size());
for (size_t i=0; i<coefficients.size(); i++)
orientationCoefficients[i] = coefficients[i].q;
- orientationFEFunction_ = std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> > > (new LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> >(localFiniteElement,orientationCoefficients));
+ orientationFEFunction_ = std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<double,3> > > (new LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<double,3> >(localFiniteElement,orientationCoefficients));
}
@@ -765,7 +765,7 @@ private:
std::vector<Dune::FieldVector<double,3> > translationCoefficients_;
- std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> > > orientationFEFunction_;
+ std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<double,3> > > orientationFEFunction_;
};
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index 9fa2e96e847942862b0cbd937969cc69940cbcd6..bf6b08f013ff21c50cf149f520477a64e790ae12 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -13,7 +13,7 @@
template <class GridView>
void RodAssembler<GridView,3>::
-assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
+assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
using namespace Dune;
@@ -34,7 +34,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
static const int nDofs = 2;
// Extract local solution
- std::vector<RigidBodyMotion<3> > localSolution(nDofs);
+ std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[this->basis_.index(*it,i)];
@@ -58,7 +58,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
template <class GridView>
void RodAssembler<GridView,3>::
-getStrain(const std::vector<RigidBodyMotion<3> >& sol,
+getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const
{
using namespace Dune;
@@ -84,7 +84,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
int numOfBaseFct = localFiniteElement.localCoefficients().size();
- std::vector<RigidBodyMotion<3> > localSolution(2);
+ std::vector<RigidBodyMotion<double,3> > localSolution(2);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -121,7 +121,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
template <class GridView>
void RodAssembler<GridView,3>::
-getStress(const std::vector<RigidBodyMotion<3> >& sol,
+getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& stress) const
{
// Get the strain
@@ -144,7 +144,7 @@ template <class GridView>
template <class PatchGridView>
Dune::FieldVector<double,6> RodAssembler<GridView,3>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
- const std::vector<RigidBodyMotion<3> >& sol) const
+ const std::vector<RigidBodyMotion<double,3> >& sol) const
{
using namespace Dune;
@@ -171,11 +171,11 @@ getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
double pos = it->geometryInInside().corner(0);
- std::vector<RigidBodyMotion<3> > localSolution(2);
+ std::vector<RigidBodyMotion<double,3> > localSolution(2);
localSolution[0] = sol[indexSet.subIndex(*it->inside(),0,1)];
localSolution[1] = sol[indexSet.subIndex(*it->inside(),1,1)];
- std::vector<RigidBodyMotion<3> > localRefConf(2);
+ std::vector<RigidBodyMotion<double,3> > localRefConf(2);
localRefConf[0] = dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->referenceConfiguration_[indexSet.subIndex(*it->inside(),0,1)];
localRefConf[1] = dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->referenceConfiguration_[indexSet.subIndex(*it->inside(),1,1)];
@@ -222,7 +222,7 @@ getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
template <class GridView>
void RodAssembler<GridView,2>::
-assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
+assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix)
{
Dune::MatrixIndexSet neighborsPerVertex;
@@ -240,7 +240,7 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
const int numOfBaseFct = 2;
// Extract local solution
- std::vector<RigidBodyMotion<2> > localSolution(numOfBaseFct);
+ std::vector<RigidBodyMotion<double,2> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(*it,i)];
@@ -273,7 +273,7 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
template <class GridView>
void RodAssembler<GridView,2>::
getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<2> >& localSolution,
+ const std::vector<RigidBodyMotion<double,2> >& localSolution,
Dune::Matrix<MatrixBlock>& localMat) const
{
/* ndof is the number of vectors of the element */
@@ -420,7 +420,7 @@ getLocalMatrix( EntityType &entity,
template <class GridView>
void RodAssembler<GridView,2>::
-assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
+assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
if (sol.size()!=this->basis_.size())
@@ -439,7 +439,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
const int numOfBaseFct = localFiniteElement.localBasis().size();
- RigidBodyMotion<2> localSolution[numOfBaseFct];
+ RigidBodyMotion<double,2> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(*it,i)];
@@ -519,7 +519,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
template <class GridView>
double RodAssembler<GridView,2>::
-computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
+computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const
{
double energy = 0;
@@ -537,7 +537,7 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
int numOfBaseFct = localFiniteElement.localBasis().size();
- RigidBodyMotion<2> localSolution[numOfBaseFct];
+ RigidBodyMotion<double,2> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(*it,i)];
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index 1dee527af60933dbdf451071cd0a253c748e2276..5d4084e7618bff9f50dd0cf9bb63b6987e114c1d 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -24,7 +24,7 @@ class RodAssembler
/** \brief The FEM operator for an extensible, shearable rod in 3d
*/
template <class GridView>
-class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >
+class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >
{
//typedef typename GridType::template Codim<0>::Entity EntityType;
@@ -46,9 +46,9 @@ public:
//! ???
RodAssembler(const GridView &gridView,
RodLocalStiffness<GridView,double>* localStiffness)
- : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >(gridView,localStiffness)
+ : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >(gridView,localStiffness)
{
- std::vector<RigidBodyMotion<3> > referenceConfiguration(gridView.size(gridDim));
+ std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(gridDim));
typename GridView::template Codim<gridDim>::Iterator it = gridView.template begin<gridDim>();
typename GridView::template Codim<gridDim>::Iterator endIt = gridView.template end<gridDim>();
@@ -60,23 +60,23 @@ public:
referenceConfiguration[idx].r[0] = 0;
referenceConfiguration[idx].r[1] = 0;
referenceConfiguration[idx].r[2] = it->geometry().corner(0)[0];
- referenceConfiguration[idx].q = Rotation<3,double>::identity();
+ referenceConfiguration[idx].q = Rotation<double,3>::identity();
}
dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->setReferenceConfiguration(referenceConfiguration);
}
- std::vector<RigidBodyMotion<3> > getRefConfig()
+ std::vector<RigidBodyMotion<double,3> > getRefConfig()
{ return dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->referenceConfiguration_;
}
- void assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
+ void assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
- void getStrain(const std::vector<RigidBodyMotion<3> >& sol,
+ void getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const;
- void getStress(const std::vector<RigidBodyMotion<3> >& sol,
+ void getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& stress) const;
/** \brief Return resultant force across boundary in canonical coordinates
@@ -84,7 +84,7 @@ public:
\note Linear run-time in the size of the grid */
template <class PatchGridView>
Dune::FieldVector<double,6> getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
- const std::vector<RigidBodyMotion<3> >& sol) const;
+ const std::vector<RigidBodyMotion<double,3> >& sol) const;
}; // end class
@@ -92,7 +92,7 @@ public:
/** \brief The FEM operator for a 2D extensible, shearable rod
*/
template <class GridView>
-class RodAssembler<GridView,2> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<2> >
+class RodAssembler<GridView,2> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >
{
typedef typename GridView::template Codim<0>::Entity EntityType;
@@ -118,7 +118,7 @@ public:
//! ???
RodAssembler(const GridView &gridView)
- : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<2> >(gridView,NULL)
+ : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >(gridView,NULL)
{
B = 1;
A1 = 1;
@@ -135,20 +135,20 @@ public:
/** \brief Assemble the tangent stiffness matrix and the right hand side
*/
- void assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
+ void assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix);
- void assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
+ void assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Compute the energy of a deformation state */
- double computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const;
+ double computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const;
protected:
/** \brief Compute the element tangent stiffness matrix */
void getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<2> >& localSolution,
+ const std::vector<RigidBodyMotion<double,2> >& localSolution,
Dune::Matrix<MatrixBlock>& mat) const;
}; // end class
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index 927c6584c8f12b6c6446ca990e7139aa2aae7f96..a526da4bf800d98bc699e7abcbdfcbed38450153 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -13,9 +13,9 @@
template<class GridView, class RT>
class RodLocalStiffness
- : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<3> >
+ : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<RT,3> >
{
- typedef RigidBodyMotion<3> TargetSpace;
+ typedef RigidBodyMotion<RT,3> TargetSpace;
// grid types
typedef typename GridView::Grid::ctype DT;
@@ -33,7 +33,7 @@ class RodLocalStiffness
public:
/** \brief The stress-free configuration */
- std::vector<RigidBodyMotion<3> > referenceConfiguration_;
+ std::vector<RigidBodyMotion<RT,3> > referenceConfiguration_;
public:
@@ -93,36 +93,36 @@ public:
- void setReferenceConfiguration(const std::vector<RigidBodyMotion<3> >& referenceConfiguration) {
+ void setReferenceConfiguration(const std::vector<RigidBodyMotion<RT,3> >& referenceConfiguration) {
referenceConfiguration_ = referenceConfiguration;
}
/** \brief Local element energy for a P1 element */
virtual RT energy (const Entity& e,
- const Dune::array<RigidBodyMotion<3>, dim+1>& localSolution) const;
+ const Dune::array<RigidBodyMotion<RT,3>, dim+1>& localSolution) const;
virtual RT energy (const Entity& e,
const typename P1NodalBasis<GridView>::LocalFiniteElement& localFiniteElement,
- const std::vector<RigidBodyMotion<3> >& localSolution) const
+ const std::vector<RigidBodyMotion<RT,3> >& localSolution) const
{
assert(localSolution.size()==2);
- Dune::array<RigidBodyMotion<3>, 2> localSolutionArray = {localSolution[0], localSolution[1]};
+ Dune::array<RigidBodyMotion<RT,3>, 2> localSolutionArray = {localSolution[0], localSolution[1]};
return energy(e,localSolutionArray);
}
/** \brief Assemble the element gradient of the energy functional */
void assembleGradient(const Entity& element,
- const std::vector<RigidBodyMotion<3> >& solution,
+ const std::vector<RigidBodyMotion<RT,3> >& solution,
std::vector<Dune::FieldVector<double,6> >& gradient) const;
- Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
+ Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
const Entity& element,
const Dune::FieldVector<double,1>& pos) const;
protected:
void getLocalReferenceConfiguration(const Entity& element,
- std::vector<RigidBodyMotion<3> >& localReferenceConfiguration) const {
+ std::vector<RigidBodyMotion<RT,3> >& localReferenceConfiguration) const {
int numOfBaseFct = element.template count<dim>();
localReferenceConfiguration.resize(numOfBaseFct);
@@ -131,14 +131,14 @@ protected:
localReferenceConfiguration[i] = referenceConfiguration_[gridView_.indexSet().subIndex(element,i,dim)];
}
- static void interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
+ static void interpolationDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
Dune::array<Quaternion<double>,6>& grad);
- static void interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
+ static void interpolationVelocityDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad);
template <class T>
- static Dune::FieldVector<T,3> darboux(const Rotation<3,T>& q, const Dune::FieldVector<T,4>& q_s)
+ static Dune::FieldVector<T,3> darboux(const Rotation<T,3>& q, const Dune::FieldVector<T,4>& q_s)
{
Dune::FieldVector<double,3> u; // The Darboux vector
@@ -154,12 +154,12 @@ protected:
template <class GridType, class RT>
RT RodLocalStiffness<GridType, RT>::
energy(const Entity& element,
- const Dune::array<RigidBodyMotion<3>, dim+1>& localSolution
+ const Dune::array<RigidBodyMotion<RT,3>, dim+1>& localSolution
) const
{
RT energy = 0;
- std::vector<RigidBodyMotion<3> > localReferenceConfiguration;
+ std::vector<RigidBodyMotion<RT,3> > localReferenceConfiguration;
getLocalReferenceConfiguration(element, localReferenceConfiguration);
// ///////////////////////////////////////////////////////////////////////////////
@@ -172,7 +172,7 @@ energy(const Entity& element,
= Dune::QuadratureRules<double, 1>::rule(element.type(), shearQuadOrder);
// hack: convert from std::array to std::vector
- std::vector<RigidBodyMotion<3> > localSolutionVector(localSolution.begin(), localSolution.end());
+ std::vector<RigidBodyMotion<RT,3> > localSolutionVector(localSolution.begin(), localSolution.end());
for (size_t pt=0; pt<shearingQuad.size(); pt++) {
@@ -221,7 +221,7 @@ energy(const Entity& element,
template <class GridType, class RT>
void RodLocalStiffness<GridType, RT>::
-interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
+interpolationDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
Dune::array<Quaternion<double>,6>& grad)
{
// Clear output array
@@ -231,18 +231,18 @@ interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, doub
// The derivatives with respect to w^0
// Compute q_1^{-1}q_0
- Rotation<3,RT> q1InvQ0 = q1;
+ Rotation<RT,3> q1InvQ0 = q1;
q1InvQ0.invert();
q1InvQ0 = q1InvQ0.mult(q0);
{
// Compute v = (1-s) \exp^{-1} ( q_1^{-1} q_0)
- Dune::FieldVector<RT,3> v = Rotation<3,RT>::expInv(q1InvQ0);
+ Dune::FieldVector<RT,3> v = Rotation<RT,3>::expInv(q1InvQ0);
v *= (1-s);
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<3,RT>::Dexp(v);
+ Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(v);
- Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<3,RT>::DexpInv(q1InvQ0);
+ Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<RT,3>::DexpInv(q1InvQ0);
Dune::FieldMatrix<RT,4,4> mat(0);
for (int i=0; i<4; i++)
@@ -268,18 +268,18 @@ interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, doub
// The derivatives with respect to w^1
// Compute q_0^{-1}
- Rotation<3,RT> q0InvQ1 = q0;
+ Rotation<RT,3> q0InvQ1 = q0;
q0InvQ1.invert();
q0InvQ1 = q0InvQ1.mult(q1);
{
// Compute v = s \exp^{-1} ( q_0^{-1} q_1)
- Dune::FieldVector<RT,3> v = Rotation<3,RT>::expInv(q0InvQ1);
+ Dune::FieldVector<RT,3> v = Rotation<RT,3>::expInv(q0InvQ1);
v *= s;
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<3,RT>::Dexp(v);
+ Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(v);
- Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<3,RT>::DexpInv(q0InvQ1);
+ Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<RT,3>::DexpInv(q0InvQ1);
Dune::FieldMatrix<RT,4,4> mat(0);
for (int i=0; i<4; i++)
@@ -306,7 +306,7 @@ interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, doub
template <class GridType, class RT>
void RodLocalStiffness<GridType, RT>::
-interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
+interpolationVelocityDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad)
{
// Clear output array
@@ -314,20 +314,20 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
grad[i] = 0;
// Compute q_0^{-1}
- Rotation<3,RT> q0Inv = q0;
+ Rotation<RT,3> q0Inv = q0;
q0Inv.invert();
// Compute v = s \exp^{-1} ( q_0^{-1} q_1)
- Dune::FieldVector<RT,3> v = Rotation<3,RT>::expInv(q0Inv.mult(q1));
+ Dune::FieldVector<RT,3> v = Rotation<RT,3>::expInv(q0Inv.mult(q1));
v *= s/intervalLength;
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<3,RT>::Dexp(v);
+ Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(v);
Dune::array<Dune::FieldMatrix<RT,3,3>, 4> ddExp;
- Rotation<3,RT>::DDexp(v, ddExp);
+ Rotation<RT,3>::DDexp(v, ddExp);
- Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<3,RT>::DexpInv(q0Inv.mult(q1));
+ Dune::FieldMatrix<RT,3,4> dExpInv = Rotation<RT,3>::DexpInv(q0Inv.mult(q1));
Dune::FieldMatrix<RT,4,4> mat(0);
for (int i=0; i<4; i++)
@@ -347,7 +347,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
dw[j] = 0.5*(i==j); // dExp_v_0[j][i];
// \xi = \exp^{-1} q_0^{-1} q_1
- Dune::FieldVector<RT,3> xi = Rotation<3,RT>::expInv(q0Inv.mult(q1));
+ Dune::FieldVector<RT,3> xi = Rotation<RT,3>::expInv(q0Inv.mult(q1));
Quaternion<RT> addend0;
addend0 = 0;
@@ -362,7 +362,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
dwConj = dwConj.mult(q0Inv.mult(q1));
Dune::FieldVector<RT,3> dxi(0);
- Rotation<3,RT>::DexpInv(q0Inv.mult(q1)).umv(dwConj, dxi);
+ Rotation<RT,3>::DexpInv(q0Inv.mult(q1)).umv(dwConj, dxi);
Quaternion<RT> vHv;
for (int j=0; j<4; j++) {
@@ -397,7 +397,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
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<RT,3> xi = Rotation<3,RT>::expInv(q0Inv.mult(q1));
+ Dune::FieldVector<RT,3> xi = Rotation<RT,3>::expInv(q0Inv.mult(q1));
// \parder{\xi}{w^1_j} = ...
Dune::FieldVector<RT,3> dxi(0);
@@ -432,7 +432,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
template <class GridType, class RT>
Dune::FieldVector<double, 6> RodLocalStiffness<GridType, RT>::
-getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
+getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
const Entity& element,
const Dune::FieldVector<double,1>& pos) const
{
@@ -475,10 +475,10 @@ getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
r_s[i] = localSolution[0].r[i]*shapeGrad[0][0] + localSolution[1].r[i]*shapeGrad[1][0];
// Interpolate the rotation at the quadrature point
- Rotation<3,double> q = Rotation<3,double>::interpolate(localSolution[0].q, localSolution[1].q, pos);
+ Rotation<RT,3> q = Rotation<RT,3>::interpolate(localSolution[0].q, localSolution[1].q, pos);
// Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> q_s = Rotation<3,double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
+ Quaternion<double> q_s = Rotation<RT,3>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
pos);
// Transformation from the reference element
q_s *= inv[0][0];
@@ -506,12 +506,12 @@ getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
template <class GridType, class RT>
void RodLocalStiffness<GridType, RT>::
assembleGradient(const Entity& element,
- const std::vector<RigidBodyMotion<3> >& solution,
+ const std::vector<RigidBodyMotion<RT,3> >& solution,
std::vector<Dune::FieldVector<double,6> >& gradient) const
{
using namespace Dune;
- std::vector<RigidBodyMotion<3> > localReferenceConfiguration;
+ std::vector<RigidBodyMotion<RT,3> > localReferenceConfiguration;
getLocalReferenceConfiguration(element, localReferenceConfiguration);
// Extract local solution on this element
@@ -568,7 +568,7 @@ assembleGradient(const Entity& element,
r_s[i] = solution[0].r[i]*shapeGrad[0] + solution[1].r[i]*shapeGrad[1];
// Interpolate current rotation at this quadrature point
- Rotation<3,double> q = Rotation<3,double>::interpolate(solution[0].q, solution[1].q,quadPos[0]);
+ Rotation<RT,3> q = Rotation<RT,3>::interpolate(solution[0].q, solution[1].q,quadPos[0]);
// The current strain
FieldVector<double,blocksize> strain = getStrain(solution, element, quadPos);
@@ -637,10 +637,10 @@ assembleGradient(const Entity& element,
double weight = bendingQuad[pt].weight() * integrationElement;
// Interpolate current rotation at this quadrature point
- Rotation<3,double> q = Rotation<3,double>::interpolate(solution[0].q, solution[1].q,quadPos[0]);
+ Rotation<RT,3> q = Rotation<RT,3>::interpolate(solution[0].q, solution[1].q,quadPos[0]);
// Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> q_s = Rotation<3,double>::interpolateDerivative(solution[0].q, solution[1].q,
+ Quaternion<double> q_s = Rotation<RT,3>::interpolateDerivative(solution[0].q, solution[1].q,
quadPos);
// Transformation from the reference element
q_s *= inv[0][0];