diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 212b40b1022b740ad2aa9b13511a0f948e6f855f..197a1d1133413e4d3a62c8c04800a033d114cfcc 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -227,10 +227,11 @@ public:
RT curvatureEnergy(const Tensor3<field_type,3,3,gridDim>& DR) const
{
+ using std::pow;
#ifdef DONT_USE_CURL
- return mu_ * std::pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
+ return mu_ * pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
#else
- return mu_ * std::pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
+ return mu_ * pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
#endif
}
diff --git a/dune/gfe/cosseratrodenergy.hh b/dune/gfe/cosseratrodenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a7c1386d4a40aebe08dfec1eba33ad9b8fed999e
--- /dev/null
+++ b/dune/gfe/cosseratrodenergy.hh
@@ -0,0 +1,458 @@
+#ifndef DUNE_GFE_COSSERATRODENERGY_HH
+#define DUNE_GFE_COSSERATRODENERGY_HH
+
+#include <array>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 8)
+#include <dune/common/transpose.hh>
+#endif
+
+#include <dune/istl/matrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+
+namespace Dune::GFE {
+
+template<class Basis, class LocalInterpolationRule, class RT>
+class CosseratRodEnergy
+: public LocalEnergy<Basis, RigidBodyMotion<RT,3> >
+{
+ typedef RigidBodyMotion<RT,3> TargetSpace;
+
+ // grid types
+ using GridView = typename Basis::GridView;
+ typedef typename GridView::Grid::ctype DT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+
+ // some other sizes
+ enum {dim=GridView::dimension};
+
+ // Quadrature order used for the extension and shear energy
+ enum {shearQuadOrder = 2};
+
+ // Quadrature order used for the bending and torsion energy
+ enum {bendingQuadOrder = 2};
+
+public:
+
+ /** \brief The stress-free configuration
+
+ The number type cannot be RT, because RT become `adouble` when
+ using RodLocalStiffness together with an AD system.
+ The referenceConfiguration is not a variable, and we don't
+ want to use `adouble` for it.
+ */
+ std::vector<RigidBodyMotion<double,3> > referenceConfiguration_;
+
+public:
+
+ //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
+ static constexpr auto blocksize = TargetSpace::EmbeddedTangentVector::dimension;
+
+ // /////////////////////////////////
+ // The material parameters
+ // /////////////////////////////////
+
+ /** \brief Material constants */
+ std::array<double,3> K_;
+ std::array<double,3> A_;
+
+ GridView gridView_;
+
+ //! Constructor
+ CosseratRodEnergy(const GridView& gridView,
+ const std::array<double,3>& K, const std::array<double,3>& A)
+ : K_(K),
+ A_(A),
+ gridView_(gridView)
+ {}
+
+ /** \brief Constructor setting shape constants and material parameters
+ \param A The rod section area
+ \param J1, J2 The geometric moments (Flächenträgheitsmomente)
+ \param E Young's modulus
+ \param nu Poisson number
+ */
+ CosseratRodEnergy(const GridView& gridView,
+ double A, double J1, double J2, double E, double nu)
+ : gridView_(gridView)
+ {
+ // shear modulus
+ double G = E/(2+2*nu);
+
+ K_[0] = E * J1;
+ K_[1] = E * J2;
+ K_[2] = G * (J1 + J2);
+
+ A_[0] = G * A;
+ A_[1] = G * A;
+ A_[2] = E * A;
+ }
+
+ /** \brief Set the stress-free configuration
+ */
+ void setReferenceConfiguration(const std::vector<RigidBodyMotion<double,3> >& referenceConfiguration) {
+ referenceConfiguration_ = referenceConfiguration;
+ }
+
+ /** \brief Compute local element energy */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const std::vector<RigidBodyMotion<RT,3> >& localSolution) const override;
+
+ /** \brief Get the rod strain at one point in the rod
+ *
+ * \tparam Number This is a member template because the method has to work for double and adouble
+ */
+ template<class ReboundLocalInterpolationRule>
+ auto getStrain(const ReboundLocalInterpolationRule& localSolution,
+ const Entity& element,
+ const FieldVector<double,1>& pos) const;
+
+ /** \brief Get the rod stress at one point in the rod
+ *
+ * \tparam Number This is a member template because the method has to work for double and adouble
+ */
+ template<class Number>
+ auto getStress(const std::vector<RigidBodyMotion<Number,3> >& localSolution,
+ const Entity& element,
+ const FieldVector<double,1>& pos) const;
+
+ /** \brief Get average strain for each element */
+ void getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
+ BlockVector<FieldVector<double, blocksize> >& strain) const;
+
+ /** \brief Get average stress for each element */
+ void getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
+ BlockVector<FieldVector<double, blocksize> >& stress) const;
+
+ /** \brief Return resultant force across boundary in canonical coordinates
+
+ \note Linear run-time in the size of the grid */
+ template <class PatchGridView>
+ auto getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
+ const std::vector<RigidBodyMotion<double,3> >& sol) const;
+
+protected:
+
+ std::vector<RigidBodyMotion<double,3> > getLocalReferenceConfiguration(const typename Basis::LocalView& localView) const
+ {
+ unsigned int numOfBaseFct = localView.size();
+ std::vector<RigidBodyMotion<double,3> > localReferenceConfiguration(numOfBaseFct);
+
+ for (size_t i=0; i<numOfBaseFct; i++)
+ localReferenceConfiguration[i] = referenceConfiguration_[localView.index(i)];
+
+ return localReferenceConfiguration;
+ }
+
+ template <class T>
+ static FieldVector<T,3> darboux(const Rotation<T,3>& q, const FieldVector<T,4>& q_s)
+ {
+ FieldVector<T,3> u; // The Darboux vector
+
+ u[0] = 2 * (q.B(0) * q_s);
+ u[1] = 2 * (q.B(1) * q_s);
+ u[2] = 2 * (q.B(2) * q_s);
+
+ return u;
+ }
+
+};
+
+template<class Basis, class LocalInterpolationRule, class RT>
+RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+energy(const typename Basis::LocalView& localView,
+ const std::vector<RigidBodyMotion<RT,3> >& localCoefficients) const
+{
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
+
+ const auto& element = localView.element();
+
+ RT energy = 0;
+
+ std::vector<RigidBodyMotion<double,3> > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
+ using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<RigidBodyMotion<double,3> >::other;
+ InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
+
+ // ///////////////////////////////////////////////////////////////////////////////
+ // The following two loops are a reduced integration scheme. We integrate
+ // the transverse shear and extensional energy with a first-order quadrature
+ // formula, even though it should be second order. This prevents shear-locking.
+ // ///////////////////////////////////////////////////////////////////////////////
+
+ const auto& shearingQuad = QuadratureRules<double, 1>::rule(element.type(), shearQuadOrder);
+
+ for (size_t pt=0; pt<shearingQuad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const auto quadPos = shearingQuad[pt].position();
+
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+
+ double weight = shearingQuad[pt].weight() * integrationElement;
+
+ auto strain = getStrain(localConfiguration, element, quadPos);
+
+ // The reference strain
+ auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+
+ for (int i=0; i<3; i++)
+ energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
+
+ }
+
+ // Get quadrature rule
+ const auto& bendingQuad = QuadratureRules<double, 1>::rule(element.type(), bendingQuadOrder);
+
+ for (size_t pt=0; pt<bendingQuad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const FieldVector<double,1>& quadPos = bendingQuad[pt].position();
+
+ double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
+
+ auto strain = getStrain(localConfiguration, element, quadPos);
+
+ // The reference strain
+ auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+
+ // Part II: the bending and twisting energy
+ for (int i=0; i<3; i++)
+ energy += weight * 0.5 * K_[i] * (strain[i+3] - referenceStrain[i+3]) * (strain[i+3] - referenceStrain[i+3]);
+
+ }
+
+ return energy;
+}
+
+
+template<class Basis, class LocalInterpolationRule, class RT>
+template <class ReboundLocalInterpolationRule>
+auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+getStrain(const ReboundLocalInterpolationRule& localInterpolation,
+ const Entity& element,
+ const FieldVector<double,1>& pos) const
+{
+ const auto jit = element.geometry().jacobianInverseTransposed(pos);
+
+ auto value = localInterpolation.evaluate(pos);
+
+ auto referenceDerivative = localInterpolation.evaluateDerivative(pos);
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 8)
+ auto derivative = referenceDerivative * transpose(jit);
+#else
+ auto derivative = referenceDerivative;
+ derivative *= jit[0][0];
+#endif
+
+ using Number = std::decay_t<decltype(derivative[0][0])>;
+ FieldVector<Number,3> r_s = {derivative[0][0], derivative[1][0], derivative[2][0]};
+
+ // Transformation from the reference element
+ Quaternion<Number> q_s(derivative[3][0],
+ derivative[4][0],
+ derivative[5][0],
+ derivative[6][0]);
+
+ // /////////////////////////////////////////////
+ // Sum it all up
+ // /////////////////////////////////////////////
+
+ // Strain defined on each element
+ // Part I: the shearing and stretching strain
+ FieldVector<Number, 6> strain(0);
+ strain[0] = r_s * value.q.director(0); // shear strain
+ strain[1] = r_s * value.q.director(1); // shear strain
+ strain[2] = r_s * value.q.director(2); // stretching strain
+
+ // Part II: the Darboux vector
+
+ FieldVector<Number,3> u = darboux<Number>(value.q, q_s);
+ strain[3] = u[0];
+ strain[4] = u[1];
+ strain[5] = u[2];
+
+ return strain;
+}
+
+template<class Basis, class LocalInterpolationRule, class RT>
+template <class Number>
+auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+getStress(const std::vector<RigidBodyMotion<Number,3> >& localSolution,
+ const Entity& element,
+ const FieldVector<double, 1>& pos) const
+{
+ const auto& indexSet = gridView_.indexSet();
+ std::vector<TargetSpace> localRefConf = {referenceConfiguration_[indexSet.subIndex(element, 0, 1)],
+ referenceConfiguration_[indexSet.subIndex(element, 1, 1)]};
+
+ auto&& strain = getStrain(localSolution, element, pos);
+ auto&& referenceStrain = getStrain(localRefConf, element, pos);
+
+ FieldVector<RT, 6> stress;
+ for (int i=0; i < dim; i++)
+ stress[i] = (strain[i] - referenceStrain[i]) * A_[i];
+
+ for (int i=0; i < dim; i++)
+ stress[i+3] = (strain[i+3] - referenceStrain[i+3]) * K_[i];
+ return stress;
+}
+
+template<class Basis, class LocalInterpolationRule, class RT>
+void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
+ BlockVector<FieldVector<double, blocksize> >& strain) const
+{
+ const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
+
+ if (sol.size()!=this->basis_.size())
+ DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
+
+ // Strain defined on each element
+ strain.resize(indexSet.size(0));
+ strain = 0;
+
+ // Loop over all elements
+ for (const auto& element : elements(this->basis_.gridView()))
+ {
+ int elementIdx = indexSet.index(element);
+
+ // Extract local solution on this element
+ P1LocalFiniteElement<double,double,1> localFiniteElement;
+ int numOfBaseFct = localFiniteElement.localCoefficients().size();
+
+ std::vector<RigidBodyMotion<double,3> > localSolution(2);
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.subIndex(element,i,1)];
+
+ // Get quadrature rule
+ const int polOrd = 2;
+ const auto& quad = QuadratureRules<double, 1>::rule(element.type(), polOrd);
+
+ for (std::size_t pt=0; pt<quad.size(); pt++)
+ {
+ // Local position of the quadrature point
+ const auto quadPos = quad[pt].position();
+
+ double weight = quad[pt].weight() * element.geometry().integrationElement(quadPos);
+
+ auto localStrain = getStrain(localSolution, element, quad[pt].position());
+
+ // Sum it all up
+ strain[elementIdx].axpy(weight, localStrain);
+ }
+
+ // /////////////////////////////////////////////////////////////////////////
+ // We want the average strain per element. Therefore we have to divide
+ // the integral we just computed by the element volume.
+ // /////////////////////////////////////////////////////////////////////////
+ // we know the element is a line, therefore the integration element is the volume
+ FieldVector<double,1> dummyPos(0.5);
+ strain[elementIdx] /= element.geometry().integrationElement(dummyPos);
+ }
+}
+
+template<class Basis, class LocalInterpolationRule, class RT>
+void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
+ BlockVector<FieldVector<double, blocksize> >& stress) const
+{
+ // Get the strain
+ getStrain(sol,stress);
+
+ // Get reference strain
+ BlockVector<FieldVector<double, blocksize> > referenceStrain;
+ getStrain(referenceConfiguration_, referenceStrain);
+
+ // Linear diagonal constitutive law
+ for (size_t i=0; i<stress.size(); i++)
+ {
+ for (int j=0; j<3; j++)
+ {
+ stress[i][j] = (stress[i][j] - referenceStrain[i][j]) * A_[j];
+ stress[i][j+3] = (stress[i][j+3] - referenceStrain[i][j+3]) * K_[j];
+ }
+ }
+}
+
+template<class Basis, class LocalInterpolationRule, class RT>
+template <class PatchGridView>
+auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
+ const std::vector<RigidBodyMotion<double,3> >& sol) const
+{
+ const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
+
+ if (sol.size()!=indexSet.size(1))
+ DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
+
+ FieldVector<double,3> canonicalStress(0);
+ FieldVector<double,3> canonicalTorque(0);
+
+ // Loop over the given boundary
+ for (auto facet : boundary)
+ {
+ // //////////////////////////////////////////////
+ // Compute force across this boundary face
+ // //////////////////////////////////////////////
+
+ double pos = facet.geometryInInside().corner(0);
+
+ std::vector<RigidBodyMotion<double,3> > localSolution(2);
+ localSolution[0] = sol[indexSet.subIndex(*facet.inside(),0,1)];
+ localSolution[1] = sol[indexSet.subIndex(*facet.inside(),1,1)];
+
+ std::vector<RigidBodyMotion<double,3> > localRefConf(2);
+ localRefConf[0] = referenceConfiguration_[indexSet.subIndex(*facet.inside(),0,1)];
+ localRefConf[1] = referenceConfiguration_[indexSet.subIndex(*facet.inside(),1,1)];
+
+ auto strain = getStrain(localSolution, *facet.inside(), pos);
+ auto referenceStrain = getStrain(localRefConf, *facet.inside(), pos);
+
+ FieldVector<double,3> localStress;
+ for (int i=0; i<3; i++)
+ localStress[i] = (strain[i] - referenceStrain[i]) * A_[i];
+
+ FieldVector<double,3> localTorque;
+ for (int i=0; i<3; i++)
+ localTorque[i] = (strain[i+3] - referenceStrain[i+3]) * K_[i];
+
+ // Transform stress given with respect to the basis given by the three directors to
+ // the canonical basis of R^3
+
+ FieldMatrix<double,3,3> orientationMatrix;
+ sol[indexSet.subIndex(*facet.inside(),facet.indexInInside(),1)].q.matrix(orientationMatrix);
+
+ orientationMatrix.umv(localStress, canonicalStress);
+
+ orientationMatrix.umv(localTorque, canonicalTorque);
+
+ // Multiply force times boundary normal to get the transmitted force
+ canonicalStress *= facet.unitOuterNormal(FieldVector<double,0>(0))[0];
+ canonicalTorque *= facet.unitOuterNormal(FieldVector<double,0>(0))[0];
+ }
+
+ FieldVector<double,6> result;
+ for (int i=0; i<3; i++)
+ {
+ result[i] = canonicalStress[i];
+ result[i+3] = canonicalTorque[i];
+ }
+
+ return result;
+}
+
+} // namespace Dune::GFE
+
+#endif
+
diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index 7be9133aa8839008cef76454992ffd251fac41f4..cdd98f405361db054429f9288f69331cebbaa418 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -212,7 +212,7 @@ assembleGradient(const std::vector<TargetSpace>& sol,
{
localView.bind(element);
- // A 1d grid has two vertices
+ // The number of degrees of freedom of the current element
const auto nDofs = localView.tree().size();
// Extract local solution
diff --git a/dune/gfe/geodesicfeassemblerwrapper.hh b/dune/gfe/geodesicfeassemblerwrapper.hh
index 4570b4077798ce8fbe3e14a55444e91adf4ec1bc..c42a1e3ef63b001b01d231de9fc95a7947648096 100644
--- a/dune/gfe/geodesicfeassemblerwrapper.hh
+++ b/dune/gfe/geodesicfeassemblerwrapper.hh
@@ -84,7 +84,7 @@ splitVector(const std::vector<TargetSpace>& sol) const {
Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> solutionSplit;
solutionSplit[_0].resize(n);
solutionSplit[_1].resize(n);
- for (int i = 0; i < n; i++) {
+ for (std::size_t i = 0; i < n; i++) {
solutionSplit[_0][i] = sol[i].r; // Deformation part
solutionSplit[_1][i] = sol[i].q; // Rotational part
}
@@ -145,7 +145,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
hessian = 0;
gradient.resize(n);
gradient = 0;
- for (int i = 0; i < n; i++) {
+ for (std::size_t i = 0; i < n; i++) {
for (int j = 0; j < blocksize0 + blocksize1; j++)
gradient[i][j] = j < blocksize0 ? gradient0[i][j] : gradient1[i][j - blocksize0];
}
diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
index 9f5bc7ce47f47c43eddda15ff0fe52acd52fbc91..202fd21be8a8cba0865b3fa82bd1fc663e72babb 100644
--- a/dune/gfe/gramschmidtsolver.hh
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -27,7 +27,8 @@ class GramSchmidtSolver
static void normalize(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
Dune::FieldVector<field_type,embeddedDim>& v)
{
- v /= std::sqrt(matrix.energyScalarProduct(v,v));
+ using std::sqrt;
+ v /= sqrt(matrix.energyScalarProduct(v,v));
}
@@ -130,4 +131,4 @@ private:
};
-#endif
\ No newline at end of file
+#endif
diff --git a/dune/gfe/localfirstordermodel.hh b/dune/gfe/localfirstordermodel.hh
index a098c5bacc3c549d4caea909e6815fb3111c0fab..ad5c1903cf122a9caecc6a5a9179dc946eef13a0 100644
--- a/dune/gfe/localfirstordermodel.hh
+++ b/dune/gfe/localfirstordermodel.hh
@@ -13,9 +13,7 @@ class LocalFirstOrderModel
{
public:
- /** \brief Assemble the element gradient of the energy functional
-
- The default implementation in this class uses a finite difference approximation */
+ /** \brief Assemble the element gradient of the energy functional */
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<typename TargetSpace::TangentVector>& gradient) const = 0;
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 3502f371dc54ec4756c2a7701604a6eda35d4153..e13acbd8d4934b9b0fd6010175be22eb0a519e71 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -59,6 +59,13 @@ public:
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
+ /** \brief Rebind the FEFunction to another TargetSpace */
+ template<class U>
+ struct rebind
+ {
+ using other = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,U>;
+ };
+
/** \brief The number of Lagrange points */
unsigned int size() const
{
@@ -589,6 +596,13 @@ public:
}
+ /** \brief Rebind the FEFunction to another TargetSpace */
+ template<class U>
+ struct rebind
+ {
+ using other = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,U>;
+ };
+
/** \brief The number of Lagrange points */
unsigned int size() const
{
diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 259270ccb475e1b63b84227c543bba7391ff5761..ba538c55f28df4b9f0bbe37829e023d454d4e622 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -70,6 +70,13 @@ Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A,
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
+ /** \brief Rebind the FEFunction to another TargetSpace */
+ template<class U>
+ struct rebind
+ {
+ using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
+ };
+
/** \brief The number of Lagrange points */
unsigned int size() const
{
@@ -452,6 +459,13 @@ Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A,
orientationFunction_ = std::make_unique<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
}
+ /** \brief Rebind the FEFunction to another TargetSpace */
+ template<class U>
+ struct rebind
+ {
+ using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
+ };
+
/** \brief The number of Lagrange points */
unsigned int size() const
{
diff --git a/dune/gfe/parallel/globalp2mapper.hh b/dune/gfe/parallel/globalp2mapper.hh
index 437273b96e4b714dadfa1e9f3eb5258843735b74..82c0ac3401cfac2daeaf719f110af147137989e6 100644
--- a/dune/gfe/parallel/globalp2mapper.hh
+++ b/dune/gfe/parallel/globalp2mapper.hh
@@ -26,54 +26,96 @@ namespace Dune {
p2Mapper_(gridView)
{
#if !HAVE_DUNE_PARMG
- static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids");
+ static_assert(GridView::dimension<=2, "Only implemented for one- and two-dimensional grids");
if (gridView.size(GeometryTypes::triangle)>1)
DUNE_THROW(NotImplemented, "GlobalP2Mapper only works for quad grids!");
- GlobalIndexSet<GridView> globalVertexIndex(gridView,2);
- GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
+ GlobalIndexSet<GridView> globalVertexIndex(gridView,GridView::dimension);
GlobalIndexSet<GridView> globalElementIndex(gridView,0);
- // total number of degrees of freedom
- this->size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
-
auto localView = p2Mapper_.localView();
- // Determine
- for (const auto& element : elements(gridView))
+ if (GridView::dimension==1)
{
- localView.bind(element);
+ // total number of degrees of freedom
+ this->size_ = globalVertexIndex.size(GridView::dimension) + globalElementIndex.size(0);
- // Loop over all local degrees of freedom
- for (size_t i=0; i<localView.size(); i++)
+ // Determine
+ for (const auto& element : elements(gridView))
{
- int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
- int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
+ localView.bind(element);
- auto localIndex = localView.index(i);
- int globalIndex;
- switch (codim)
+ // Loop over all local degrees of freedom
+ for (size_t i=0; i<localView.size(); i++)
{
- case 2: // vertex dofs
- globalIndex = globalVertexIndex.index(element.template subEntity<2>(entity));
- break;
-
- case 1: // edge dofs
- globalIndex = globalEdgeIndex.index(element.template subEntity<1>(entity)) + globalVertexIndex.size(2);
- break;
-
- case 0: // element dofs
- globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
- + globalEdgeIndex.size(1)
- + globalVertexIndex.size(2);
- break;
-
- default:
- DUNE_THROW(Dune::Exception, "Impossible codimension!");
+ int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
+ int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
+
+ auto localIndex = localView.index(i);
+ int globalIndex;
+ switch (codim)
+ {
+ case 1: // vertex dofs
+ globalIndex = globalVertexIndex.index(element.template subEntity<1>(entity));
+ break;
+
+ case 0: // element dofs
+ globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
+ + globalVertexIndex.size(2);
+ break;
+
+ default:
+ DUNE_THROW(Dune::Exception, "Impossible codimension!");
+ }
+
+ localGlobalMap_[localIndex] = globalIndex;
}
+ }
+
+ }
+ else
+ {
+ GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
- localGlobalMap_[localIndex] = globalIndex;
+ // total number of degrees of freedom
+ this->size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
+
+ // Determine
+ for (const auto& element : elements(gridView))
+ {
+ localView.bind(element);
+
+ // Loop over all local degrees of freedom
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
+ int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
+
+ auto localIndex = localView.index(i);
+ int globalIndex;
+ switch (codim)
+ {
+ case 2: // vertex dofs
+ globalIndex = globalVertexIndex.index(element.template subEntity<GridView::dimension>(entity));
+ break;
+
+ case 1: // edge dofs
+ globalIndex = globalEdgeIndex.index(element.template subEntity<1>(entity)) + globalVertexIndex.size(2);
+ break;
+
+ case 0: // element dofs
+ globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
+ + globalEdgeIndex.size(1)
+ + globalVertexIndex.size(2);
+ break;
+
+ default:
+ DUNE_THROW(Dune::Exception, "Impossible codimension!");
+ }
+
+ localGlobalMap_[localIndex] = globalIndex;
+ }
}
}
#endif
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 4937720c4adabd0abb1ce46a66eed4ab1a2605be..1976ed3c5f9b630f27614fed39a24e2a3d92c332 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -243,7 +243,7 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
if (this->verbosity_ == Solver::FULL and rank==0)
- std::cout << "Oveall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+ std::cout << "Overall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
// Transfer matrix data
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 88d4e7713aa48b760478a93cbd3ddeb323979acc..89a0e3246dc7a94ccfe28096bd644b9ee30775ea 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -424,7 +424,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
std::cout << "Gradient norm: " << gradient.two_norm() << std::endl;
}
if (this->verbosity_ == Solver::FULL and rank==0)
- std::cout << "Oveall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+ std::cout << "Overall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
// Transfer matrix data
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
deleted file mode 100644
index 32300bd3fe334971b8bd850a69255938bb5afbfa..0000000000000000000000000000000000000000
--- a/dune/gfe/rodassembler.cc
+++ /dev/null
@@ -1,591 +0,0 @@
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/localfunctions/lagrange/p1.hh>
-
-#include <dune/gfe/rodlocalstiffness.hh>
-
-
-
-template <class Basis>
-void RodAssembler<Basis,3>::
-assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
-{
- using namespace Dune;
-
- if (sol.size()!=this->basis_.size())
- DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
-
- grad.resize(sol.size());
- grad = 0;
-
- // A view on the FE basis on a single element
- auto localView = this->basis_.localView();
-
- // Loop over all elements
- for (const auto& element : Dune::elements(this->basis_.gridView()))
- {
- localView.bind(element);
-
- // A 1d grid has two vertices
- static const int nDofs = 2;
-
- // Extract local solution
- std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
-
- for (int i=0; i<nDofs; i++)
- localSolution[i] = sol[localView.index(i)];
-
- // Assemble local gradient
- std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
-
- this->localStiffness_->assembleGradient(localView,
- localSolution,
- localGradient);
-
- // Add to global gradient
- for (int i=0; i<nDofs; i++)
- grad[localView.index(i)] += localGradient[i];
-
- }
-
-}
-
-
-template <class GridView>
-void RodAssembler<GridView,3>::
-getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const
-{
- using namespace Dune;
-
- const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
-
- if (sol.size()!=this->basis_.size())
- DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
-
- // Strain defined on each element
- strain.resize(indexSet.size(0));
- strain = 0;
-
- // Loop over all elements
- for (const auto& element : elements(this->basis_.gridView()))
- {
- int elementIdx = indexSet.index(element);
-
- // Extract local solution on this element
- Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
- int numOfBaseFct = localFiniteElement.localCoefficients().size();
-
- std::vector<RigidBodyMotion<double,3> > localSolution(2);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.subIndex(element,i,gridDim)];
-
- // Get quadrature rule
- const int polOrd = 2;
- const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), polOrd);
-
- for (std::size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- double weight = quad[pt].weight() * element.geometry().integrationElement(quadPos);
-
- FieldVector<double,blocksize> localStrain = std::dynamic_pointer_cast<RodLocalStiffness<GridView, double> >(this->localStiffness_)->getStrain(localSolution, element, quad[pt].position());
-
- // Sum it all up
- strain[elementIdx].axpy(weight, localStrain);
-
- }
-
- // /////////////////////////////////////////////////////////////////////////
- // We want the average strain per element. Therefore we have to divide
- // the integral we just computed by the element volume.
- // /////////////////////////////////////////////////////////////////////////
- // we know the element is a line, therefore the integration element is the volume
- FieldVector<double,1> dummyPos(0.5);
- strain[elementIdx] /= element.geometry().integrationElement(dummyPos);
-
- }
-
-}
-
-template <class GridView>
-void RodAssembler<GridView,3>::
-getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& stress) const
-{
- // Get the strain
- getStrain(sol,stress);
-
- // Get reference strain
- Dune::BlockVector<Dune::FieldVector<double, blocksize> > referenceStrain;
- getStrain(dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->referenceConfiguration_, referenceStrain);
-
- // Linear diagonal constitutive law
- for (size_t i=0; i<stress.size(); i++) {
- for (int j=0; j<3; j++) {
- stress[i][j] = (stress[i][j] - referenceStrain[i][j]) * dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->A_[j];
- stress[i][j+3] = (stress[i][j+3] - referenceStrain[i][j+3]) * dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->K_[j];
- }
- }
-}
-
-template <class GridView>
-template <class PatchGridView>
-Dune::FieldVector<double,6> RodAssembler<GridView,3>::
-getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
- const std::vector<RigidBodyMotion<double,3> >& sol) const
-{
- using namespace Dune;
-
- // if (gridView_ != &boundary.gridView())
- // DUNE_THROW(Dune::Exception, "The boundary patch has to match the grid view of the assembler!");
-
- const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
-
- if (sol.size()!=indexSet.size(gridDim))
- DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
-
- FieldVector<double,3> canonicalStress(0);
- FieldVector<double,3> canonicalTorque(0);
-
- // Loop over the given boundary
- typename BoundaryPatch<PatchGridView>::iterator it = boundary.begin();
- typename BoundaryPatch<PatchGridView>::iterator endIt = boundary.end();
-
- for (; it!=endIt; ++it) {
-
- // //////////////////////////////////////////////
- // Compute force across this boundary face
- // //////////////////////////////////////////////
-
- double pos = it->geometryInInside().corner(0);
-
- 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<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)];
-
- FieldVector<double, blocksize> strain = dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->getStrain(localSolution, *it->inside(), pos);
- FieldVector<double, blocksize> referenceStrain = dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->getStrain(localRefConf, *it->inside(), pos);
-
- FieldVector<double,3> localStress;
- for (int i=0; i<3; i++)
- localStress[i] = (strain[i] - referenceStrain[i]) * dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->A_[i];
-
- FieldVector<double,3> localTorque;
- for (int i=0; i<3; i++)
- localTorque[i] = (strain[i+3] - referenceStrain[i+3]) * dynamic_cast<RodLocalStiffness<GridView, double>* >(this->localStiffness_)->K_[i];
-
- // Transform stress given with respect to the basis given by the three directors to
- // the canonical basis of R^3
-
- FieldMatrix<double,3,3> orientationMatrix;
- sol[indexSet.subIndex(*it->inside(),it->indexInInside(),1)].q.matrix(orientationMatrix);
-
- orientationMatrix.umv(localStress, canonicalStress);
-
- orientationMatrix.umv(localTorque, canonicalTorque);
- // Reverse transformation to make sure we did the correct thing
-// assert( std::abs(localStress[0]-canonicalStress*sol[0].q.director(0)) < 1e-6 );
-// assert( std::abs(localStress[1]-canonicalStress*sol[0].q.director(1)) < 1e-6 );
-// assert( std::abs(localStress[2]-canonicalStress*sol[0].q.director(2)) < 1e-6 );
-
- // Multiply force times boundary normal to get the transmitted force
- canonicalStress *= it->unitOuterNormal(FieldVector<double,0>(0))[0];
- canonicalTorque *= it->unitOuterNormal(FieldVector<double,0>(0))[0];
-
- }
-
- Dune::FieldVector<double,6> result;
- for (int i=0; i<3; i++) {
- result[i] = canonicalStress[i];
- result[i+3] = canonicalTorque[i];
- }
-
- return result;
-}
-
-
-template <class GridView>
-void RodAssembler<GridView,2>::
-assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix)
-{
- Dune::MatrixIndexSet neighborsPerVertex;
- this->getNeighborsPerVertex(neighborsPerVertex);
-
- matrix = 0;
-
- Dune::Matrix<MatrixBlock> mat;
-
- for (const auto& element : Dune::elements(this->basis_.gridView()))
- {
- const int numOfBaseFct = 2;
-
- // Extract local solution
- std::vector<RigidBodyMotion<double,2> > localSolution(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[this->basis_.index(element,i)];
-
- // setup matrix
- getLocalMatrix(element, localSolution, mat);
-
- // Add element matrix to global stiffness matrix
- for(int i=0; i<numOfBaseFct; i++) {
-
- int row = this->basis_.index(element,i);
-
- for (int j=0; j<numOfBaseFct; j++ ) {
-
- int col = this->basis_.index(element,j);
- matrix[row][col] += mat[i][j];
-
- }
- }
-
- }
-
-}
-
-
-
-
-
-
-template <class GridView>
-void RodAssembler<GridView,2>::
-getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<double,2> >& localSolution,
- Dune::Matrix<MatrixBlock>& localMat) const
-{
- /* ndof is the number of vectors of the element */
- int ndof = localSolution.size();
-
- localMat.setSize(ndof,ndof);
- localMat = 0;
-
- Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
-
- // Get quadrature rule
- const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(entity.type(), 2);
-
- /* Loop over all integration points */
- for (int ip=0; ip<quad.size(); ip++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,gridDim>& quadPos = quad[ip].position();
-
- // calc Jacobian inverse before integration element is evaluated
- const Dune::FieldMatrix<double,gridDim,gridDim>& inv = entity.geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = entity.geometry().integrationElement(quadPos);
-
- /* Compute the weight of the current integration point */
- double weight = quad[ip].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(ndof);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(ndof);
-
- // multiply with jacobian inverse
- for (int dof=0; dof<ndof; dof++)
- inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
-
-
- std::vector<Dune::FieldVector<double,1> > shapeFunction;
- localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
- double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
-
- double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
-
- for (int i=0; i<localMat.N(); i++) {
-
- for (int j=0; j<localMat.M(); j++) {
-
- // \partial J^2 / \partial x_i \partial x_j
- localMat[i][j][0][0] += weight * shapeGrad[i][0] * shapeGrad[j][0]
- * (A1 * cos(theta) * cos(theta) + A3 * sin(theta) * sin(theta));
-
- // \partial J^2 / \partial x_i \partial y_j
- localMat[i][j][0][1] += weight * shapeGrad[i][0] * shapeGrad[j][0]
- * (-A1 + A3) * sin(theta)* cos(theta);
-
- // \partial J^2 / \partial x_i \partial theta_j
- localMat[i][j][0][2] += weight * shapeGrad[i][0] * shapeFunction[j]
- * (-A1 * (x_s*sin(theta) + y_s*cos(theta)) * cos(theta)
- - A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
- +A3 * (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
- +A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
-
- // \partial J^2 / \partial y_i \partial x_j
- localMat[i][j][1][0] += weight * shapeGrad[i][0] * shapeGrad[j][0]
- * (-A1 * sin(theta)* cos(theta) + A3 * cos(theta) * sin(theta));
-
- // \partial J^2 / \partial y_i \partial y_j
- localMat[i][j][1][1] += weight * shapeGrad[i][0] * shapeGrad[j][0]
- * (A1 * sin(theta)*sin(theta) + A3 * cos(theta)*cos(theta));
-
- // \partial J^2 / \partial y_i \partial theta_j
- localMat[i][j][1][2] += weight * shapeGrad[i][0] * shapeFunction[j]
- * (A1 * (x_s * sin(theta) + y_s * cos(theta)) * sin(theta)
- -A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
- +A3 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
- -A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
-
- // \partial J^2 / \partial theta_i \partial x_j
- localMat[i][j][2][0] += weight * shapeFunction[i] * shapeGrad[j][0]
- * (-A1 * (x_s*sin(theta) + y_s*cos(theta)) * cos(theta)
- - A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
- +A3 * (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
- +A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
-
- // \partial J^2 / \partial theta_i \partial y_j
- localMat[i][j][2][1] += weight * shapeFunction[i] * shapeGrad[j][0]
- * (A1 * (x_s * sin(theta) + y_s * cos(theta)) * sin(theta)
- -A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
- +A3 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
- -A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
-
- // \partial J^2 / \partial theta_i \partial theta_j
- localMat[i][j][2][2] += weight * B * shapeGrad[i][0] * shapeGrad[j][0];
- localMat[i][j][2][2] += weight * shapeFunction[i] * shapeFunction[j]
- * (+ A1 * (x_s*sin(theta) + y_s*cos(theta)) * (x_s*sin(theta) + y_s*cos(theta))
- + A1 * (x_s*cos(theta) - y_s*sin(theta)) * (-x_s*cos(theta)+ y_s*sin(theta))
- + A3 * (x_s*cos(theta) - y_s*sin(theta)) * (x_s*cos(theta) - y_s*sin(theta))
- - A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * (x_s*sin(theta) + y_s*cos(theta)));
-
-
-
- }
-
- }
-
-
- }
-
-#if 0
- static int eleme = 0;
- printf("********* Element %d **********\n", eleme++);
- for (int row=0; row<matSize; row++) {
-
- for (int rcomp=0; rcomp<dim; rcomp++) {
-
- for (int col=0; col<matSize; col++) {
-
- for (int ccomp=0; ccomp<dim; ccomp++)
- std::cout << mat[row][col][rcomp][ccomp] << " ";
-
- std::cout << " ";
-
- }
-
- std::cout << std::endl;
-
- }
-
- std::cout << std::endl;
-
- }
- exit(0);
-#endif
-
-}
-
-template <class GridView>
-void RodAssembler<GridView,2>::
-assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
-{
- if (sol.size()!=this->basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
-
- grad.resize(sol.size());
- grad = 0;
-
- // Loop over all elements
- for (const auto& element : Dune::elements(this->basis_gridView()))
- {
- // Extract local solution on this element
- Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
- const int numOfBaseFct = localFiniteElement.localBasis().size();
-
- RigidBodyMotion<double,2> localSolution[numOfBaseFct];
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[this->basis_.index(element,i)];
-
- // Get quadrature rule
- const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(element.type(), 2);
-
- for (int pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- const Dune::FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = element.geometry().integrationElement(quadPos);
-
- double weight = quad[pt].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
-
- // multiply with jacobian inverse
- for (int dof=0; dof<numOfBaseFct; dof++)
- inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
-
- // Get the values of the shape functions
- std::vector<Dune::FieldVector<double,1> > shapeFunction;
- localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
- double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
- double theta_s = localSolution[0].q.angle_*shapeGrad[0][0] + localSolution[1].q.angle_*shapeGrad[1][0];
-
- double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
-
- // /////////////////////////////////////////////
- // Sum it all up
- // /////////////////////////////////////////////
-
- double partA1 = A1 * (x_s * cos(theta) - y_s * sin(theta));
- double partA3 = A3 * (x_s * sin(theta) + y_s * cos(theta) - 1);
-
- for (int dof=0; dof<numOfBaseFct; dof++) {
-
- int globalDof = this->basis_.index(element,dof);
-
- //printf("globalDof: %d partA1: %g partA3: %g\n", globalDof, partA1, partA3);
-
- // \partial J / \partial x^i
- grad[globalDof][0] += weight * (partA1 * cos(theta) + partA3 * sin(theta)) * shapeGrad[dof][0];
-
- // \partial J / \partial y^i
- grad[globalDof][1] += weight * (-partA1 * sin(theta) + partA3 * cos(theta)) * shapeGrad[dof][0];
-
- // \partial J / \partial \theta^i
- grad[globalDof][2] += weight * (B * theta_s * shapeGrad[dof][0]
- + partA1 * (-x_s * sin(theta) - y_s * cos(theta)) * shapeFunction[dof]
- + partA3 * ( x_s * cos(theta) - y_s * sin(theta)) * shapeFunction[dof]);
-
- }
-
-
- }
-
- }
-
-}
-
-
-template <class GridView>
-double RodAssembler<GridView,2>::
-computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const
-{
- double energy = 0;
-
- if (sol.size()!=this->basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
-
- // Loop over all elements
- for (const auto& element : Dune::elements(this->basis_gridView()))
- {
- // Extract local solution on this element
- Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
-
- int numOfBaseFct = localFiniteElement.localBasis().size();
-
- std::vector<RigidBodyMotion<double,2> > localSolution(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[this->basis_.index(element,i)];
-
- // Get quadrature rule
- const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(element.type(), 2);
-
- for (int pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- const Dune::FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = element.geometry().integrationElement(quadPos);
-
- double weight = quad[pt].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
-
- // multiply with jacobian inverse
- for (int dof=0; dof<numOfBaseFct; dof++)
- inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
-
- // Get the value of the shape functions
- std::vector<Dune::FieldVector<double,1> > shapeFunction;
- localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
- double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
- double theta_s = localSolution[0].q.angle_*shapeGrad[0][0] + localSolution[1].q.angle_*shapeGrad[1][0];
-
- double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
-
- // /////////////////////////////////////////////
- // Sum it all up
- // /////////////////////////////////////////////
-
- double partA1 = x_s * cos(theta) - y_s * sin(theta);
- double partA3 = x_s * sin(theta) + y_s * cos(theta) - 1;
-
-
- energy += 0.5 * weight * (B * theta_s * theta_s
- + A1 * partA1 * partA1
- + A3 * partA3 * partA3);
-
- }
-
- }
-
- return energy;
-
-}
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
deleted file mode 100644
index 64c1327ecae83ff146d6675be04243003a8158ae..0000000000000000000000000000000000000000
--- a/dune/gfe/rodassembler.hh
+++ /dev/null
@@ -1,160 +0,0 @@
-#ifndef DUNE_EXTENSIBLE_ROD_ASSEMBLER_HH
-#define DUNE_EXTENSIBLE_ROD_ASSEMBLER_HH
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/gfe/localgeodesicfefdstiffness.hh>
-#include "rigidbodymotion.hh"
-#include "rodlocalstiffness.hh"
-#include "geodesicfeassembler.hh"
-
-/** \brief The FEM operator for an extensible, shearable rod in 3d
- */
-template <class Basis, int spaceDim>
-class RodAssembler
-{
- static_assert(spaceDim==2 || spaceDim==3,
- "You can only instantiate the class RodAssembler for 2d and 3d spaces");
-};
-
-/** \brief The FEM operator for an extensible, shearable rod in 3d
- */
-template <class Basis>
-class RodAssembler<Basis,3> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >
-{
- typedef typename Basis::GridView GridView;
-
- //! Dimension of the grid.
- enum { gridDim = GridView::dimension };
- static_assert(gridDim==1, "RodAssembler can only be used with one-dimensional grids!");
-
- enum { elementOrder = 1};
-
- //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
- enum { blocksize = 6 };
-
-public:
- //! ???
- RodAssembler(const Basis& basis,
- LocalGeodesicFEStiffness<Basis, RigidBodyMotion<double,3> >& localStiffness)
- : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
- {
- std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.size());
-
- for (const auto vertex : Dune::vertices(basis.gridView()))
- {
- auto idx = basis.gridView().indexSet().index(vertex);
-
- referenceConfiguration[idx].r[0] = 0;
- referenceConfiguration[idx].r[1] = 0;
- referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
- referenceConfiguration[idx].q = Rotation<double,3>::identity();
- }
-
- rodEnergy()->setReferenceConfiguration(referenceConfiguration);
- }
-
- auto rodEnergy()
- {
- // TODO: Does not work for other stiffness implementations
- auto localFDStiffness = std::dynamic_pointer_cast<LocalGeodesicFEFDStiffness<Basis, RigidBodyMotion<double,3> > >(this->localStiffness_);
- return const_cast<RodLocalStiffness<GridView,double>*>(dynamic_cast<const RodLocalStiffness<GridView,double>*>(localFDStiffness->localEnergy_));
- }
-
- std::vector<RigidBodyMotion<double,3> > getRefConfig()
- {
- return rodEnergy()->referenceConfiguration_;
- }
-
- virtual void assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const override;
-
- void getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const;
-
- 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
-
- \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<double,3> >& sol) const;
- }; // end class
-
-
-/** \brief The FEM operator for a 2D extensible, shearable rod
- */
-template <class Basis>
-class RodAssembler<Basis,2> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >
-{
-
- typedef typename Basis::GridView GridView;
- typedef typename GridView::template Codim<0>::Entity EntityType;
-
- //! Dimension of the grid.
- enum { gridDim = GridView::dimension };
- static_assert(gridDim==1, "RodAssembler can only be used with one-dimensional grids!");
-
- enum { elementOrder = 1};
-
- //! Each block is x, y, theta
- enum { blocksize = 3 };
-
- //!
- typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
- /** \brief Material constants */
- double B;
- double A1;
- double A3;
-
-public:
-
- //! ???
- RodAssembler(const GridView &gridView)
- : GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >(gridView,nullptr)
- {
- B = 1;
- A1 = 1;
- A3 = 1;
- }
-
- ~RodAssembler() {}
-
- void setParameters(double b, double a1, double a3) {
- B = b;
- A1 = a1;
- A3 = a3;
- }
-
- /** \brief Assemble the tangent stiffness matrix and the right hand side
- */
- void assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix);
-
- virtual void assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const override;
-
- /** \brief Compute the energy of a deformation state */
- virtual double computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const override;
-
-protected:
-
- /** \brief Compute the element tangent stiffness matrix */
- void getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<double,2> >& localSolution,
- Dune::Matrix<MatrixBlock>& mat) const;
-
-}; // end class
-
-#include "rodassembler.cc"
-
-#endif
-
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
deleted file mode 100644
index 5157a30e3f7a4e4d67f292c938a13ccc505b8f2d..0000000000000000000000000000000000000000
--- a/dune/gfe/rodlocalstiffness.hh
+++ /dev/null
@@ -1,720 +0,0 @@
-#ifndef ROD_LOCAL_STIFFNESS_HH
-#define ROD_LOCAL_STIFFNESS_HH
-
-#include <array>
-
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-
-#include <dune/gfe/localfirstordermodel.hh>
-#include "rigidbodymotion.hh"
-
-template<class GridView, class RT>
-class RodLocalStiffness
-: public Dune::GFE::LocalFirstOrderModel<Dune::Functions::LagrangeBasis<GridView,1>, RigidBodyMotion<RT,3> >
-{
- typedef RigidBodyMotion<RT,3> TargetSpace;
- typedef Dune::Functions::LagrangeBasis<GridView,1> Basis;
-
- // grid types
- typedef typename GridView::Grid::ctype DT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- enum {dim=GridView::dimension};
-
- // Quadrature order used for the extension and shear energy
- enum {shearQuadOrder = 2};
-
- // Quadrature order used for the bending and torsion energy
- enum {bendingQuadOrder = 2};
-
-public:
-
- /** \brief The stress-free configuration */
- std::vector<RigidBodyMotion<RT,3> > referenceConfiguration_;
-
-public:
-
- //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
- enum { blocksize = 6 };
-
- // define the number of components of your system, this is used outside
- // to allocate the correct size of (dense) blocks with a FieldMatrix
- enum {m=blocksize};
-
- // types for matrics, vectors and boundary conditions
- typedef Dune::FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
- typedef Dune::FieldVector<RT,m> VBlockType; // one entry in the global vectors
-
- // /////////////////////////////////
- // The material parameters
- // /////////////////////////////////
-
- /** \brief Material constants */
- std::array<double,3> K_;
- std::array<double,3> A_;
-
- GridView gridView_;
-
- //! Constructor
- RodLocalStiffness (const GridView& gridView,
- const std::array<double,3>& K, const std::array<double,3>& A)
- : K_(K),
- A_(A),
- gridView_(gridView)
- {}
-
- /** \brief Constructor setting shape constants and material parameters
- \param A The rod section area
- \param J1, J2 The geometric moments (Flächenträgheitsmomente)
- \param E Young's modulus
- \param nu Poisson number
- */
- RodLocalStiffness (const GridView& gridView,
- double A, double J1, double J2, double E, double nu)
- : gridView_(gridView)
- {
- // shear modulus
- double G = E/(2+2*nu);
-
- K_[0] = E * J1;
- K_[1] = E * J2;
- K_[2] = G * (J1 + J2);
-
- A_[0] = G * A;
- A_[1] = G * A;
- A_[2] = E * A;
- }
-
-
-
- void setReferenceConfiguration(const std::vector<RigidBodyMotion<RT,3> >& referenceConfiguration) {
- referenceConfiguration_ = referenceConfiguration;
- }
-
- /** \brief Compute local element energy */
- virtual RT energy (const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& localSolution) const override;
-
- /** \brief Assemble the element gradient of the energy functional */
- void assembleGradient(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& solution,
- std::vector<Dune::FieldVector<double,6> >& gradient) const override;
-
- Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
- const Entity& element,
- const Dune::FieldVector<double,1>& pos) const;
-
- Dune::FieldVector<RT, 6> getStress(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<RT,3> >& localReferenceConfiguration) const {
-
- unsigned int numOfBaseFct = element.subEntities(dim);
- localReferenceConfiguration.resize(numOfBaseFct);
-
- for (size_t i=0; i<numOfBaseFct; i++)
- localReferenceConfiguration[i] = referenceConfiguration_[gridView_.indexSet().subIndex(element,i,dim)];
- }
-
-public:
- static void interpolationDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
- std::array<Quaternion<double>,6>& grad);
-
- static void interpolationVelocityDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
- double intervalLength, std::array<Quaternion<double>,6>& grad);
-
-protected:
- template <class T>
- 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
-
- u[0] = 2 * (q.B(0) * q_s);
- u[1] = 2 * (q.B(1) * q_s);
- u[2] = 2 * (q.B(2) * q_s);
-
- return u;
- }
-
-};
-
-template <class GridType, class RT>
-RT RodLocalStiffness<GridType, RT>::
-energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& localSolution) const
-{
- assert(localSolution.size()==2);
- const auto& element = localView.element();
-
- RT energy = 0;
-
- std::vector<RigidBodyMotion<RT,3> > localReferenceConfiguration;
- getLocalReferenceConfiguration(element, localReferenceConfiguration);
-
- // ///////////////////////////////////////////////////////////////////////////////
- // The following two loops are a reduced integration scheme. We integrate
- // the transverse shear and extensional energy with a first-order quadrature
- // formula, even though it should be second order. This prevents shear-locking.
- // ///////////////////////////////////////////////////////////////////////////////
-
- const Dune::QuadratureRule<double, 1>& shearingQuad
- = Dune::QuadratureRules<double, 1>::rule(element.type(), shearQuadOrder);
-
- // hack: convert from std::array to std::vector
- std::vector<RigidBodyMotion<RT,3> > localSolutionVector(localSolution.begin(), localSolution.end());
-
- for (size_t pt=0; pt<shearingQuad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,1>& quadPos = shearingQuad[pt].position();
-
- const double integrationElement = element.geometry().integrationElement(quadPos);
-
- double weight = shearingQuad[pt].weight() * integrationElement;
-
- Dune::FieldVector<double,6> strain = getStrain(localSolutionVector, element, quadPos);
-
- // The reference strain
- Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
-
- for (int i=0; i<3; i++)
- energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
-
- }
-
- // Get quadrature rule
- const Dune::QuadratureRule<double, 1>& bendingQuad
- = Dune::QuadratureRules<double, 1>::rule(element.type(), bendingQuadOrder);
-
- for (size_t pt=0; pt<bendingQuad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,1>& quadPos = bendingQuad[pt].position();
-
- double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
-
- Dune::FieldVector<double,6> strain = getStrain(localSolutionVector, element, quadPos);
-
- // The reference strain
- Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
-
- // Part II: the bending and twisting energy
- for (int i=0; i<3; i++)
- energy += weight * 0.5 * K_[i] * (strain[i+3] - referenceStrain[i+3]) * (strain[i+3] - referenceStrain[i+3]);
-
- }
-
- return energy;
-}
-
-
-template <class GridType, class RT>
-void RodLocalStiffness<GridType, RT>::
-interpolationDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
- std::array<Quaternion<double>,6>& grad)
-{
- // Clear output array
- for (int i=0; i<6; i++)
- grad[i] = 0;
-
- // The derivatives with respect to w^0
-
- // Compute q_1^{-1}q_0
- 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<RT,3>::expInv(q1InvQ0);
- v *= (1-s);
-
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(SkewMatrix<RT,3>(v));
-
- 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++)
- 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<RT> dw;
- for (int j=0; j<4; j++)
- dw[j] = 0.5 * (i==j); // dExp[j][i] at v=0
-
- dw = q1InvQ0.Quaternion<double>::mult(dw);
-
- mat.umv(dw,grad[i]);
- grad[i] = q1.Quaternion<double>::mult(grad[i]);
-
- }
- }
-
-
- // The derivatives with respect to w^1
-
- // Compute q_0^{-1}
- 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<RT,3>::expInv(q0InvQ1);
- v *= s;
-
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(SkewMatrix<RT,3>(v));
-
- 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++)
- 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<RT> dw;
- for (int j=0; j<4; j++)
- dw[j] = 0.5 * ((i-3)==j); // dExp[j][i-3] at v=0
-
- dw = q0InvQ1.Quaternion<double>::mult(dw);
-
- mat.umv(dw,grad[i]);
- grad[i] = q0.Quaternion<double>::mult(grad[i]);
-
- }
- }
-}
-
-
-
-template <class GridType, class RT>
-void RodLocalStiffness<GridType, RT>::
-interpolationVelocityDerivative(const Rotation<RT,3>& q0, const Rotation<RT,3>& q1, double s,
- double intervalLength, std::array<Quaternion<double>,6>& grad)
-{
- // Clear output array
- for (int i=0; i<6; i++)
- grad[i] = 0;
-
- // Compute q_0^{-1}
- Rotation<RT,3> q0Inv = q0;
- q0Inv.invert();
-
-
- // Compute v = s \exp^{-1} ( q_0^{-1} q_1)
- Dune::FieldVector<RT,3> v = Rotation<RT,3>::expInv(q0Inv.mult(q1));
- v *= s/intervalLength;
-
- Dune::FieldMatrix<RT,4,3> dExp_v = Rotation<RT,3>::Dexp(SkewMatrix<RT,3>(v));
-
- std::array<Dune::FieldMatrix<RT,3,3>, 4> ddExp;
- Rotation<RT,3>::DDexp(v, ddExp);
-
- 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++)
- 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<RT> 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<RT,3> xi = Rotation<RT,3>::expInv(q0Inv.mult(q1));
-
- Quaternion<RT> addend0;
- addend0 = 0;
- dExp_v.umv(xi,addend0);
- addend0 = dw.mult(addend0);
- addend0 /= intervalLength;
-
- // \parder{\xi}{w^1_j} = ...
- Quaternion<RT> 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<RT,3> dxi(0);
- Rotation<RT,3>::DexpInv(q0Inv.mult(q1)).umv(dwConj, dxi);
-
- Quaternion<RT> 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];
- }
-
- vHv *= s/intervalLength/intervalLength;
-
- // Third addend
- mat.umv(dwConj,grad[i]);
-
- // add up
- grad[i] += addend0;
- grad[i] += vHv;
-
- grad[i] = q0.Quaternion<double>::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<RT> 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<RT,3> xi = Rotation<RT,3>::expInv(q0Inv.mult(q1));
-
- // \parder{\xi}{w^1_j} = ...
- Dune::FieldVector<RT,3> dxi(0);
- dExpInv.umv(q0Inv.Quaternion<double>::mult(q1.Quaternion<double>::mult(dw)), dxi);
-
- Quaternion<RT> 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
- // ///////////////////////////////////
-
-
- dw = q0Inv.Quaternion<double>::mult(q1.Quaternion<double>::mult(dw));
-
- mat.umv(dw,grad[i]);
- grad[i] += vHv;
-
- grad[i] = q0.Quaternion<double>::mult(grad[i]);
-
- }
-
-}
-
-template <class GridType, class RT>
-Dune::FieldVector<double, 6> RodLocalStiffness<GridType, RT>::
-getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
- const Entity& element,
- const Dune::FieldVector<double,1>& pos) const
-{
- if (!element.isLeaf())
- DUNE_THROW(Dune::NotImplemented, "Only for leaf elements");
-
- assert(localSolution.size() == 2);
-
- // Strain defined on each element
- Dune::FieldVector<double, 6> strain(0);
-
- // Extract local solution on this element
- Dune::P1LocalFiniteElement<double,double,1> localFiniteElement;
- int numOfBaseFct = localFiniteElement.localCoefficients().size();
-
- const Dune::FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(pos);
-
- // ///////////////////////////////////////
- // Compute deformation gradient
- // ///////////////////////////////////////
- std::vector<Dune::FieldMatrix<double,1,1> > shapeGrad;
-
- localFiniteElement.localBasis().evaluateJacobian(pos, shapeGrad);
-
- for (int dof=0; dof<numOfBaseFct; dof++) {
-
- // multiply with jacobian inverse
- Dune::FieldVector<double,1> tmp(0);
- inv.umv(shapeGrad[dof][0], tmp);
- shapeGrad[dof][0] = tmp;
-
- }
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- Dune::FieldVector<double,3> r_s;
- for (int i=0; i<3; i++)
- 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<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<RT,3>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
- pos);
- // Transformation from the reference element
- q_s *= inv[0][0];
-
- // /////////////////////////////////////////////
- // Sum it all up
- // /////////////////////////////////////////////
-
- // Part I: the shearing and stretching strain
- strain[0] = r_s * q.director(0); // shear strain
- strain[1] = r_s * q.director(1); // shear strain
- strain[2] = r_s * q.director(2); // stretching strain
-
- // Part II: the Darboux vector
-
- Dune::FieldVector<double,3> u = darboux(q, q_s);
-
- strain[3] = u[0];
- strain[4] = u[1];
- strain[5] = u[2];
-
- return strain;
-}
-template <class GridType, class RT>
-Dune::FieldVector<RT, 6> RodLocalStiffness<GridType, RT>::
-getStress(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
- const Entity& element,
- const Dune::FieldVector<double, 1>& pos) const
-{
- const auto& indexSet = gridView_.indexSet();
- std::vector<TargetSpace> localRefConf = {referenceConfiguration_[indexSet.subIndex(element, 0, 1)],
- referenceConfiguration_[indexSet.subIndex(element, 1, 1)]};
-
- auto&& strain = getStrain(localSolution, element, pos);
- auto&& referenceStrain = getStrain(localRefConf, element, pos);
-
- Dune::FieldVector<RT, 6> stress;
- for (int i=0; i < dim; i++)
- stress[i] = (strain[i] - referenceStrain[i]) * A_[i];
-
- for (int i=0; i < dim; i++)
- stress[i+3] = (strain[i+3] - referenceStrain[i+3]) * K_[i];
- return stress;
-}
-
-
-template <class GridType, class RT>
-void RodLocalStiffness<GridType, RT>::
-assembleGradient(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& solution,
- std::vector<Dune::FieldVector<double,6> >& gradient) const
-{
- using namespace Dune;
- const auto& element = localView.element();
-
- std::vector<RigidBodyMotion<RT,3> > localReferenceConfiguration;
- getLocalReferenceConfiguration(element, localReferenceConfiguration);
-
- // Extract local solution on this element
- Dune::P1LocalFiniteElement<double,double,1> localFiniteElement;
- int numOfBaseFct = localFiniteElement.localCoefficients().size();
-
- // init
- gradient.resize(numOfBaseFct);
- for (size_t i=0; i<gradient.size(); i++)
- gradient[i] = 0;
-
- double intervalLength = element.geometry().corner(1)[0] - element.geometry().corner(0)[0];
-
- // ///////////////////////////////////////////////////////////////////////////////////
- // Reduced integration to avoid locking: assembly is split into a shear part
- // and a bending part. Different quadrature rules are used for the two parts.
- // This avoids locking.
- // ///////////////////////////////////////////////////////////////////////////////////
-
- // Get quadrature rule
- const QuadratureRule<double, 1>& shearingQuad = QuadratureRules<double, 1>::rule(element.type(), shearQuadOrder);
-
- for (size_t pt=0; pt<shearingQuad.size(); pt++) {
-
- // Local position of the quadrature point
- const FieldVector<double,1>& quadPos = shearingQuad[pt].position();
-
- const FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = element.geometry().integrationElement(quadPos);
-
- double weight = shearingQuad[pt].weight() * integrationElement;
-
- // ///////////////////////////////////////
- // Compute deformation gradient
- // ///////////////////////////////////////
- std::vector<Dune::FieldMatrix<double,1,1> > shapeGrad;
- localFiniteElement.localBasis().evaluateJacobian(quadPos, shapeGrad);
-
- for (int dof=0; dof<numOfBaseFct; dof++) {
-
- // multiply with jacobian inverse
- FieldVector<double,1> tmp(0);
- inv.umv(shapeGrad[dof][0], tmp);
- shapeGrad[dof][0] = tmp;
-
- }
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- FieldVector<double,3> r_s;
- for (int i=0; i<3; i++)
- r_s[i] = solution[0].r[i]*shapeGrad[0] + solution[1].r[i]*shapeGrad[1];
-
- // Interpolate current rotation at this quadrature point
- 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);
-
- // The reference strain
- FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
-
-
- // dd_dvij[m][i][j] = \parder {(d_k)_i} {q}
- Tensor3<double,3 ,3, 4> dd_dq;
- q.getFirstDerivativesOfDirectors(dd_dq);
-
- // First derivatives of the position
- std::array<Quaternion<double>,6> dq_dwij;
- interpolationDerivative(solution[0].q, solution[1].q, quadPos, dq_dwij);
-
- // /////////////////////////////////////////////
- // Sum it all up
- // /////////////////////////////////////////////
-
- for (int i=0; i<numOfBaseFct; i++) {
-
- // /////////////////////////////////////////////
- // The translational part
- // /////////////////////////////////////////////
-
- // \partial \bar{W} / \partial r^i_j
- for (int j=0; j<3; j++) {
-
- for (int m=0; m<3; m++)
- gradient[i][j] += weight
- * ( A_[m] * (strain[m] - referenceStrain[m]) * shapeGrad[i] * q.director(m)[j]);
-
- }
-
- // \partial \bar{W}_v / \partial v^i_j
- for (int j=0; j<3; j++) {
-
- for (int m=0; m<3; m++) {
- FieldVector<double,3> tmp(0);
- dd_dq[m].umv(dq_dwij[3*i+j], tmp);
- gradient[i][3+j] += weight
- * A_[m] * (strain[m] - referenceStrain[m]) * (r_s*tmp);
- }
- }
-
- }
-
- }
-
- // /////////////////////////////////////////////////////
- // Now: the bending/torsion part
- // /////////////////////////////////////////////////////
-
- // Get quadrature rule
- const QuadratureRule<double, 1>& bendingQuad = QuadratureRules<double, 1>::rule(element.type(), bendingQuadOrder);
-
- for (std::size_t pt=0; pt<bendingQuad.size(); pt++) {
-
- // Local position of the quadrature point
- const FieldVector<double,1>& quadPos = bendingQuad[pt].position();
-
- const FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = element.geometry().integrationElement(quadPos);
-
- double weight = bendingQuad[pt].weight() * integrationElement;
-
- // Interpolate current rotation at this quadrature point
- 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<RT,3>::interpolateDerivative(solution[0].q, solution[1].q,
- quadPos);
- // Transformation from the reference element
- q_s *= inv[0][0];
-
-
- // The current strain
- FieldVector<double,blocksize> strain = getStrain(solution, element, quadPos);
-
- // The reference strain
- FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
-
- // First derivatives of the position
- std::array<Quaternion<double>,6> dq_dwij;
- interpolationDerivative(solution[0].q, solution[1].q, quadPos, dq_dwij);
-
- std::array<Quaternion<double>,6> dq_ds_dwij;
- interpolationVelocityDerivative(solution[0].q, solution[1].q, quadPos[0]*intervalLength, intervalLength,
- dq_ds_dwij);
-
- // /////////////////////////////////////////////
- // Sum it all up
- // /////////////////////////////////////////////
-
- for (int i=0; i<numOfBaseFct; i++) {
-
- // /////////////////////////////////////////////
- // The rotational part
- // /////////////////////////////////////////////
-
- // \partial \bar{W}_v / \partial v^i_j
- for (int j=0; j<3; j++) {
-
- for (int m=0; m<3; m++) {
-
- // Compute derivative of the strain
- /** \todo Is this formula correct? It seems strange to call
- B(m) for a _derivative_ of a rotation */
- double du_dvij_m = 2 * (dq_dwij[i*3+j].B(m) * q_s)
- + 2* ( q.B(m) * dq_ds_dwij[i*3+j]);
-
- // Sum it up
- gradient[i][3+j] += weight * K_[m]
- * (strain[m+3]-referenceStrain[m+3]) * du_dvij_m;
-
- }
-
- }
-
- }
-
- }
-}
-
-
-#endif
-
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 2ca609ed8c98312c084c5a0345824db4b9323134..96ea13db09cfbd34112410fdb530091c8ede82f0 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -145,7 +145,8 @@ class Rotation<T,3> : public Quaternion<T>
/** \brief Computes sin(x/2) / x without getting unstable for small x */
static T sincHalf(const T& x) {
- return (x < 1e-4) ? 0.5 - (x*x/48) + (x*x*x*x)/3840 : std::sin(x/2)/x;
+ using std::sin;
+ return (x < 1e-4) ? 0.5 - (x*x/48) + (x*x*x*x)/3840 : sin(x/2)/x;
}
/** \brief Computes sin(sqrt(x)/2) / sqrt(x) without getting unstable for small x
@@ -154,7 +155,9 @@ class Rotation<T,3> : public Quaternion<T>
* which ADOL-C doesn't like.
*/
static T sincHalfOfSquare(const T& x) {
- return (x < 1e-15) ? 0.5 - (x/48) + (x*x)/(120*32) : std::sin(std::sqrt(x)/2)/std::sqrt(x);
+ using std::sin;
+ using std::sqrt;
+ return (x < 1e-15) ? 0.5 - (x/48) + (x*x)/(120*32) : sin(sqrt(x)/2)/sqrt(x);
}
/** \brief Computes sin(sqrt(x)) / sqrt(x) without getting unstable for small x
@@ -163,7 +166,9 @@ class Rotation<T,3> : public Quaternion<T>
* which ADOL-C doesn't like.
*/
static T sincOfSquare(const T& x) {
- return (x < 1e-15) ? 1 - (x/6) + (x*x)/120 : std::sin(std::sqrt(x))/std::sqrt(x);
+ using std::sin;
+ using std::sqrt;
+ return (x < 1e-15) ? 1 - (x/6) + (x*x)/120 : sin(sqrt(x))/sqrt(x);
}
public:
@@ -262,6 +267,9 @@ public:
/** \brief The exponential map from \f$ \mathfrak{so}(3) \f$ to \f$ SO(3) \f$
*/
static Rotation<T,3> exp(const SkewMatrix<T,3>& v) {
+ using std::cos;
+ using std::sqrt;
+
Rotation<T,3> q;
Dune::FieldVector<T,3> vAxial = v.axial();
@@ -282,7 +290,7 @@ public:
// 1 - x*x/8 + x*x*x*x/384 - ...
q[3] = (normV2 < 1e-4)
? 1 - normV2/8 + normV2*normV2 / 384
- : std::cos(std::sqrt(normV2)/2);
+ : cos(sqrt(normV2)/2);
return q;
}
@@ -302,7 +310,10 @@ public:
*/
static Rotation<T,3> exp(const Rotation<T,3>& p, const Dune::FieldVector<T,4>& v) {
- assert( std::abs(p*v) < 1e-8 );
+ using std::abs;
+ using std::cos;
+ using std::sqrt;
+ assert( abs(p*v) < 1e-8 );
const T norm2 = v.two_norm2();
@@ -310,10 +321,10 @@ public:
// The series expansion of cos(x) at x=0 is
// 1 - x*x/2 + x*x*x*x/24 - ...
- T cos = (norm2 < 1e-4)
+ T cosValue = (norm2 < 1e-4)
? 1 - norm2/2 + norm2*norm2 / 24
- : std::cos(std::sqrt(norm2));
- result *= cos;
+ : cos(sqrt(norm2));
+ result *= cosValue;
result.axpy(sincOfSquare(norm2), v);
return result;
@@ -385,6 +396,8 @@ public:
static Dune::FieldMatrix<T,4,3> Dexp(const SkewMatrix<T,3>& v) {
+ using std::cos;
+
Dune::FieldMatrix<T,4,3> result(0);
Dune::FieldVector<T,3> vAxial = v.axial();
T norm = vAxial.two_norm();
@@ -395,8 +408,8 @@ public:
result[m][i] = (norm<1e-10)
/** \todo Isn't there a better way to implement this stably? */
- ? 0.5 * (i==m)
- : 0.5 * std::cos(norm/2) * vAxial[i] * vAxial[m] / (norm*norm) + sincHalf(norm) * ( (i==m) - vAxial[i]*vAxial[m]/(norm*norm));
+ ? T(0.5) * (i==m)
+ : 0.5 * cos(norm/2) * vAxial[i] * vAxial[m] / (norm*norm) + sincHalf(norm) * ( (i==m) - vAxial[i]*vAxial[m]/(norm*norm));
}
@@ -409,6 +422,9 @@ public:
static void DDexp(const Dune::FieldVector<T,3>& v,
std::array<Dune::FieldMatrix<T,3,3>, 4>& result) {
+ using std::cos;
+ using std::sin;
+
T norm = v.two_norm();
if (norm<=1e-10) {
@@ -427,15 +443,15 @@ public:
for (int m=0; m<3; m++) {
- result[m][i][j] = -0.25*std::sin(norm/2)*v[i]*v[j]*v[m]/(norm*norm*norm)
+ result[m][i][j] = -0.25*sin(norm/2)*v[i]*v[j]*v[m]/(norm*norm*norm)
+ ((i==j)*v[m] + (j==m)*v[i] + (i==m)*v[j] - 3*v[i]*v[j]*v[m]/(norm*norm))
- * (0.5*std::cos(norm/2) - sincHalf(norm)) / (norm*norm);
+ * (0.5*cos(norm/2) - sincHalf(norm)) / (norm*norm);
}
result[3][i][j] = -0.5/(norm*norm)
- * ( 0.5*std::cos(norm/2)*v[i]*v[j] + std::sin(norm/2) * ((i==j)*norm - v[i]*v[j]/norm));
+ * ( 0.5*cos(norm/2)*v[i]*v[j] + sin(norm/2) * ((i==j)*norm - v[i]*v[j]/norm));
}
@@ -457,7 +473,8 @@ public:
} else {
- T invSinc = 1/sincHalf(2*std::acos(q[3]));
+ using std::acos;
+ T invSinc = 1/sincHalf(2*acos(q[3]));
v[0] = q[0] * invSinc;
v[1] = q[1] * invSinc;
@@ -575,7 +592,9 @@ public:
T sp = a.globalCoordinates() * b.globalCoordinates();
// Scalar product may be larger than 1.0, due to numerical dirt
- T dist = 2*std::acos( std::min(sp,T(1.0)) );
+ using std::acos;
+ using std::min;
+ T dist = 2*acos( min(sp,T(1.0)) );
// Make sure we do the right thing if a and b are not in the same sheet
// of the double covering of the unit quaternions over SO(3)
@@ -601,7 +620,11 @@ public:
} else {
- T dist = 2*std::acos( diff[3] );
+ // TODO: ADOL-C does not like this part of the code,
+ // because arccos is not differentiable at -1 and 1.
+ // (Even though the overall 'difference' function is differentiable.)
+ using std::acos;
+ T dist = 2*acos( diff[3] );
// Make sure we do the right thing if a and b are not in the same sheet
// of the double covering of the unit quaternions over SO(3)
@@ -725,7 +748,8 @@ public:
result.axpy(-1*(q*p), q);
// The ternary operator comes from the derivative of the absolute value function
- result *= 4 * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp)) * ( (sp<0) ? -1 : 1 );
+ using std::abs;
+ result *= 4 * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp)) * ( (sp<0) ? -1 : 1 );
return result;
}
@@ -742,7 +766,8 @@ public:
for (int j=0; j<=i; j++)
A(i,j) = pProjected[i]*pProjected[j];
- A *= 4*UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp));
+ using std::abs;
+ A *= 4*UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp));
// Compute matrix B (see notes)
Dune::SymmetricMatrix<T,4> Pq;
@@ -751,7 +776,7 @@ public:
Pq(i,j) = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
// Bring it all together
- A.axpy(-4* ((sp<0) ? -1 : 1) * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp))*sp, Pq);
+ A.axpy(-4* ((sp<0) ? -1 : 1) * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp))*sp, Pq);
return A;
}
@@ -773,7 +798,8 @@ public:
Dune::FieldMatrix<T,4,4> A;
// A = row * column
Dune::FMatrixHelp::multMatrix(column,row,A);
- A *= 4 * UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp));
+ using std::abs;
+ A *= 4 * UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp));
// Compute matrix B (see notes)
Dune::FieldMatrix<T,4,4> Pp, Pq;
@@ -787,7 +813,7 @@ public:
Dune::FMatrixHelp::multMatrix(Pp,Pq,B);
// Bring it all together
- A.axpy(4 * ( (sp<0) ? -1 : 1) * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp)), B);
+ A.axpy(4 * ( (sp<0) ? -1 : 1) * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp)), B);
return A;
}
@@ -812,16 +838,17 @@ public:
double plusMinus = (sp < 0) ? -1 : 1;
+ using std::abs;
for (int i=0; i<4; i++)
for (int j=0; j<4; j++)
for (int k=0; k<4; k++) {
- result[i][j][k] = plusMinus * UnitVector<T,4>::thirdDerivativeOfArcCosSquared(std::abs(sp)) * pProjected[i] * pProjected[j] * pProjected[k]
- - UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp)) * ((i==j)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[j])*pProjected[k]
- - UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp)) * ((i==k)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[k])*pProjected[j]
- - UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp)) * pProjected[i] * Pq[j][k] * sp
- + plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp)) * ((i==j)*q.globalCoordinates()[k] + (i==k)*q.globalCoordinates()[j]) * sp
- - plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp)) * p.globalCoordinates()[i] * Pq[j][k];
+ result[i][j][k] = plusMinus * UnitVector<T,4>::thirdDerivativeOfArcCosSquared(abs(sp)) * pProjected[i] * pProjected[j] * pProjected[k]
+ - UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp)) * ((i==j)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[j])*pProjected[k]
+ - UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp)) * ((i==k)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[k])*pProjected[j]
+ - UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp)) * pProjected[i] * Pq[j][k] * sp
+ + plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp)) * ((i==j)*q.globalCoordinates()[k] + (i==k)*q.globalCoordinates()[j]) * sp
+ - plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp)) * p.globalCoordinates()[i] * Pq[j][k];
}
Pq *= 4;
@@ -862,10 +889,11 @@ public:
double plusMinus = (sp < 0) ? -1 : 1;
- result = plusMinus * UnitVector<T,4>::thirdDerivativeOfArcCosSquared(std::abs(sp)) * Tensor3<T,4,4,4>::product(qProjected,pProjected,pProjected)
- + UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp)) * derivativeOfPqOTimesPq
- - UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp)) * sp * Tensor3<T,4,4,4>::product(qProjected,Pq)
- - plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp)) * Tensor3<T,4,4,4>::product(qProjected,Pq);
+ using std::abs;
+ result = plusMinus * UnitVector<T,4>::thirdDerivativeOfArcCosSquared(abs(sp)) * Tensor3<T,4,4,4>::product(qProjected,pProjected,pProjected)
+ + UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp)) * derivativeOfPqOTimesPq
+ - UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp)) * sp * Tensor3<T,4,4,4>::product(qProjected,Pq)
+ - plusMinus * UnitVector<T,4>::derivativeOfArcCosSquared(abs(sp)) * Tensor3<T,4,4,4>::product(qProjected,Pq);
result *= 4;
@@ -955,6 +983,8 @@ public:
We tacitly assume that the matrix really is orthogonal */
void set(const Dune::FieldMatrix<T,3,3>& m) {
+ using std::sqrt;
+
// Easier writing
Dune::FieldVector<T,4>& p = (*this);
// The following equations for the derivation of a unit quaternion from a rotation
@@ -969,7 +999,7 @@ public:
// avoid rounding problems
if (p[0] >= p[1] && p[0] >= p[2] && p[0] >= p[3]) {
- p[0] = std::sqrt(p[0]);
+ p[0] = sqrt(p[0]);
// r_x r_y = (R_12 + R_21) / 4
p[1] = (m[0][1] + m[1][0]) / 4 / p[0];
@@ -982,7 +1012,7 @@ public:
} else if (p[1] >= p[0] && p[1] >= p[2] && p[1] >= p[3]) {
- p[1] = std::sqrt(p[1]);
+ p[1] = sqrt(p[1]);
// r_x r_y = (R_12 + R_21) / 4
p[0] = (m[0][1] + m[1][0]) / 4 / p[1];
@@ -995,7 +1025,7 @@ public:
} else if (p[2] >= p[0] && p[2] >= p[1] && p[2] >= p[3]) {
- p[2] = std::sqrt(p[2]);
+ p[2] = sqrt(p[2]);
// r_x r_z = (R_13 + R_31) / 4
p[0] = (m[0][2] + m[2][0]) / 4 / p[2];
@@ -1008,7 +1038,7 @@ public:
} else {
- p[3] = std::sqrt(p[3]);
+ p[3] = sqrt(p[3]);
// r_0 r_x = (R_32 - R_23) / 4
p[0] = (m[2][1] - m[1][2]) / 4 / p[3];
@@ -1028,6 +1058,9 @@ public:
We tacitly assume that the matrix really is orthogonal */
static Tensor3<T,4,3,3> derivativeOfMatrixToQuaternion(const Dune::FieldMatrix<T,3,3>& m)
{
+ using std::pow;
+ using std::sqrt;
+
Tensor3<T,4,3,3> result;
Dune::FieldVector<T,4> p;
@@ -1045,10 +1078,10 @@ public:
if (p[0] >= p[1] && p[0] >= p[2] && p[0] >= p[3])
{
result[0] = {{1,0,0},{0,-1,0},{0,0,-1}};
- result[0] *= 1.0/(8.0*std::sqrt(p[0]));
+ result[0] *= 1.0/(8.0*sqrt(p[0]));
- T denom = 32 * std::pow(p[0],1.5);
- T offDiag = 1.0/(4*std::sqrt(p[0]));
+ T denom = 32 * pow(p[0],1.5);
+ T offDiag = 1.0/(4*sqrt(p[0]));
result[1] = { {-(m[0][1]+m[1][0]) / denom, offDiag, 0},
{offDiag, (m[0][1]+m[1][0]) / denom, 0},
@@ -1065,10 +1098,10 @@ public:
else if (p[1] >= p[0] && p[1] >= p[2] && p[1] >= p[3])
{
result[1] = {{-1,0,0},{0,1,0},{0,0,-1}};
- result[1] *= 1.0/(8.0*std::sqrt(p[1]));
+ result[1] *= 1.0/(8.0*sqrt(p[1]));
- T denom = 32 * std::pow(p[1],1.5);
- T offDiag = 1.0/(4*std::sqrt(p[1]));
+ T denom = 32 * pow(p[1],1.5);
+ T offDiag = 1.0/(4*sqrt(p[1]));
result[0] = { {(m[0][1]+m[1][0]) / denom, offDiag, 0},
{offDiag, -(m[0][1]+m[1][0]) / denom, 0},
@@ -1085,10 +1118,10 @@ public:
else if (p[2] >= p[0] && p[2] >= p[1] && p[2] >= p[3])
{
result[2] = {{-1,0,0},{0,-1,0},{0,0,1}};
- result[2] *= 1.0/(8.0*std::sqrt(p[2]));
+ result[2] *= 1.0/(8.0*sqrt(p[2]));
- T denom = 32 * std::pow(p[2],1.5);
- T offDiag = 1.0/(4*std::sqrt(p[2]));
+ T denom = 32 * pow(p[2],1.5);
+ T offDiag = 1.0/(4*sqrt(p[2]));
result[0] = { {(m[0][2]+m[2][0]) / denom, 0, offDiag},
{0, (m[0][2]+m[2][0]) / denom, 0},
@@ -1105,10 +1138,10 @@ public:
else
{
result[3] = {{1,0,0},{0,1,0},{0,0,1}};
- result[3] *= 1.0/(8.0*std::sqrt(p[3]));
+ result[3] *= 1.0/(8.0*sqrt(p[3]));
- T denom = 32 * std::pow(p[3],1.5);
- T offDiag = 1.0/(4*std::sqrt(p[3]));
+ T denom = 32 * pow(p[3],1.5);
+ T offDiag = 1.0/(4*sqrt(p[3]));
result[0] = { {-(m[2][1]-m[1][2]) / denom, 0, 0},
{0, -(m[2][1]-m[1][2]) / denom, -offDiag},
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 734188488d5accf9b290a535ec518148507d1350..83f1813c53e1911611795c7358339102d175a167 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -3,7 +3,12 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
+#include <dune/common/math.hh>
+#else
#include <dune/common/power.hh>
+#endif
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/symmetricmatrix.hh>
@@ -24,42 +29,54 @@ class UnitVector
/** \brief Computes sin(x) / x without getting unstable for small x */
static T sinc(const T& x) {
- return (x < 1e-4) ? 1 - (x*x/6) : std::sin(x)/x;
+ using std::sin;
+ return (x < 1e-4) ? 1 - (x*x/6) : sin(x)/x;
}
/** \brief Compute arccos^2 without using the (at 1) nondifferentiable function acos x close to 1 */
static T arcCosSquared(const T& x) {
+ using std::acos;
const T eps = 1e-4;
if (x > 1-eps) { // acos is not differentiable, use the series expansion instead
return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
} else
- return Dune::Power<2>::eval(std::acos(x));
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
+ return Dune::power(acos(x),2);
+#else
+ return Dune::Power<2>::eval(acos(x));
+#endif
}
/** \brief Compute the derivative of arccos^2 without getting unstable for x close to 1 */
static T derivativeOfArcCosSquared(const T& x) {
+ using std::acos;
+ using std::sqrt;
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
- return -2*std::acos(x) / std::sqrt(1-x*x);
+ return -2*acos(x) / sqrt(1-x*x);
}
/** \brief Compute the second derivative of arccos^2 without getting unstable for x close to 1 */
static T secondDerivativeOfArcCosSquared(const T& x) {
+ using std::acos;
+ using std::pow;
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return 2.0/3 - 8*(x-1)/15;
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
- return 2/(1-x*x) - 2*x*std::acos(x) / std::pow(1-x*x,1.5);
+ return 2/(1-x*x) - 2*x*acos(x) / pow(1-x*x,1.5);
}
/** \brief Compute the third derivative of arccos^2 without getting unstable for x close to 1 */
static T thirdDerivativeOfArcCosSquared(const T& x) {
+ using std::acos;
+ using std::sqrt;
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return -8.0/15 + 24*(x-1)/35;
@@ -67,7 +84,7 @@ class UnitVector
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else {
T d = 1-x*x;
- return 6*x/(d*d) - 6*x*x*std::acos(x)/(d*d*std::sqrt(d)) - 2*std::acos(x)/(d*std::sqrt(d));
+ return 6*x/(d*d) - 6*x*x*acos(x)/(d*d*sqrt(d)) - 2*acos(x)/(d*sqrt(d));
}
}
@@ -156,11 +173,13 @@ public:
/** \brief The exponential map */
static UnitVector exp(const UnitVector& p, const EmbeddedTangentVector& v) {
- assert( std::abs(p.data_*v) < 1e-5 );
+ using std::abs;
+ using std::cos;
+ assert( abs(p.data_*v) < 1e-5 );
const T norm = v.two_norm();
UnitVector result = p;
- result.data_ *= std::cos(norm);
+ result.data_ *= cos(norm);
result.data_.axpy(sinc(norm), v);
return result;
}
@@ -176,15 +195,18 @@ public:
/** \brief Length of the great arc connecting the two points */
static T distance(const UnitVector& a, const UnitVector& b) {
+ using std::acos;
+ using std::min;
+
// Not nice: we are in a class for unit vectors, but the class is actually
// supposed to handle perturbations of unit vectors as well. Therefore
// we normalize here.
T x = a.data_ * b.data_/a.data_.two_norm()/b.data_.two_norm();
// paranoia: if the argument is just eps larger than 1 acos returns NaN
- x = std::min(x,1.0);
+ x = min(x,1.0);
- return std::acos(x);
+ return acos(x);
}
#if ADOLC_ADOUBLE_H
@@ -197,7 +219,8 @@ public:
adouble x = a.data_ * b.data_ / (a.data_.two_norm()*b.data_.two_norm());
// paranoia: if the argument is just eps larger than 1 acos returns NaN
- x = std::min(x,1.0);
+ using std::min;
+ x = min(x,1.0);
// Special implementation that remains AD-differentiable near x==1
return arcCosSquared(x);
@@ -219,7 +242,8 @@ public:
result = b.projectOntoTangentSpace(result);
// Gradient must be a tangent vector at b, in other words, orthogonal to it
- assert( std::abs(b.data_ * result) < 1e-5);
+ using std::abs;
+ assert(abs(b.data_ * result) < 1e-5);
return result;
}
diff --git a/problems/staticrod.parset b/problems/staticrod.parset
index 63d7044a8b016a38fa700b79777c062bc7c3e75f..2d001ecad0fbe15fff3188010c99d570ae8c0909 100644
--- a/problems/staticrod.parset
+++ b/problems/staticrod.parset
@@ -51,18 +51,17 @@ instrumented = no
# Problem specifications
############################
+# Interpolation method
+interpolationMethod = geodesic
+
A = 1
J1 = 1
J2 = 1
E = 2.5e5
nu = 0.3
-dirichletValueX = 1
-dirichletValueY = 0
-dirichletValueZ = 0
+dirichletValue = 1 0 0
-dirichletAxisX = 1
-dirichletAxisY = 0
-dirichletAxisZ = 0
+dirichletAxis = 1 0 0
dirichletAngle = 0
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index eb763ee627ae9684bef0a6436126649715778155..a34e177404fe0bb92d59a6ca1f35952bf9a08838 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -419,7 +419,7 @@ int main (int argc, char *argv[]) try
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary,
neumannFunction,
- nullptr);
+ volumeLoad);
MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
RealTuple<double,3>,
diff --git a/src/rod3d.cc b/src/rod3d.cc
index 8d4a9abcac5d91b5051fddb28b0debbe0e5a30ca..866addb252c4b8174963e3d3d6f52c1087b0ffa2 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -1,5 +1,12 @@
#include <config.h>
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/drivers/drivers.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <optional>
+
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
@@ -8,34 +15,53 @@
#include <dune/istl/io.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
+#if HAVE_DUNE_VTK
+#include <dune/vtk/vtkwriter.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
+#else
#include <dune/gfe/cosseratvtkwriter.hh>
+#endif
+
+#include <dune/gfe/cosseratrodenergy.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/geodesicdifference.hh>
#include <dune/gfe/rotation.hh>
-#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
typedef RigidBodyMotion<double,3> TargetSpace;
const int blocksize = TargetSpace::TangentVector::dimension;
+// Approximation order of the finite element space
+constexpr int order = 2;
+
using namespace Dune;
-using std::string;
int main (int argc, char *argv[]) try
{
+ MPIHelper::instance(argc, argv);
+
typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
// parse data file
ParameterTree parameterSet;
- if (argc==2)
- ParameterTreeParser::readINITree(argv[1], parameterSet);
- else
- ParameterTreeParser::readINITree("rod3d.parset", parameterSet);
+ if (argc < 2)
+ DUNE_THROW(Exception, "Usage: ./rod3d <parameter file>");
+
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
// read solver settings
const int numLevels = parameterSet.get<int>("numLevels");
@@ -71,33 +97,43 @@ int main (int argc, char *argv[]) try
using GridView = GridType::LeafGridView;
GridView gridView = grid.leafGridView();
- using FEBasis = Functions::LagrangeBasis<GridView,1>;
+ using FEBasis = Functions::LagrangeBasis<GridView,order>;
FEBasis feBasis(gridView);
SolutionType x(feBasis.size());
- // //////////////////////////
- // Initial solution
- // //////////////////////////
+ //////////////////////////////////////////////
+ // Create the stress-free configuration
+ //////////////////////////////////////////////
+
+ std::vector<double> referenceConfigurationX(feBasis.size());
+
+ auto identity = [](const FieldVector<double,1>& x) { return x; };
- for (size_t i=0; i<x.size(); i++) {
- x[i].r[0] = 0;
- x[i].r[1] = 0;
- x[i].r[2] = double(i)/(x.size()-1);
- x[i].q = Rotation<double,3>::identity();
+ Functions::interpolate(feBasis, referenceConfigurationX, identity);
+
+ std::vector<RigidBodyMotion<double,3> > referenceConfiguration(feBasis.size());
+
+ for (std::size_t i=0; i<referenceConfiguration.size(); i++)
+ {
+ referenceConfiguration[i].r[0] = 0;
+ referenceConfiguration[i].r[1] = 0;
+ referenceConfiguration[i].r[2] = referenceConfigurationX[i];
+ referenceConfiguration[i].q = Rotation<double,3>::identity();
}
+ /////////////////////////////////////////////////////////////////
+ // Select the reference configuration as initial iterate
+ /////////////////////////////////////////////////////////////////
+
+ x = referenceConfiguration;
+
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
- x.back().r[0] = parameterSet.get<double>("dirichletValueX");
- x.back().r[1] = parameterSet.get<double>("dirichletValueY");
- x.back().r[2] = parameterSet.get<double>("dirichletValueZ");
-
- FieldVector<double,3> axis;
- axis[0] = parameterSet.get<double>("dirichletAxisX");
- axis[1] = parameterSet.get<double>("dirichletAxisY");
- axis[2] = parameterSet.get<double>("dirichletAxisZ");
+ x.back().r = parameterSet.get<FieldVector<double,3> >("dirichletValue");
+
+ auto axis = parameterSet.get<FieldVector<double,3> >("dirichletAxis");
double angle = parameterSet.get<double>("dirichletAngle");
x.back().q = Rotation<double,3>(axis, M_PI*angle/180);
@@ -120,17 +156,43 @@ int main (int argc, char *argv[]) try
dirichletNodes[0] = true;
dirichletNodes.back() = true;
-
- // ///////////////////////////////////////////
- // Create a solver for the rod problem
- // ///////////////////////////////////////////
- RodLocalStiffness<GridView,double> localStiffness(gridView,
- A, J1, J2, E, nu);
+ //////////////////////////////////////////////
+ // Create the energy and assembler
+ //////////////////////////////////////////////
+
+ using ATargetSpace = TargetSpace::rebind<adouble>::other;
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<1, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<1, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+
+ // Assembler using ADOL-C
+ std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localRodEnergy;
- LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
+ if (parameterSet["interpolationMethod"] == "geodesic")
+ {
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<FEBasis, GeodesicInterpolationRule, adouble> >(gridView,
+ A, J1, J2, E, nu);
+ energy->setReferenceConfiguration(referenceConfiguration);
+ localRodEnergy = energy;
+ }
+ else if (parameterSet["interpolationMethod"] == "projected")
+ {
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<FEBasis, ProjectedInterpolationRule, adouble> >(gridView,
+ A, J1, J2, E, nu);
+ energy->setReferenceConfiguration(referenceConfiguration);
+ localRodEnergy = energy;
+ }
+ else
+ DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
- RodAssembler<FEBasis,3> rodAssembler(gridView, localFDStiffness);
+ LocalGeodesicFEADOLCStiffness<FEBasis,
+ TargetSpace> localStiffness(localRodEnergy.get());
+
+ GeodesicFEAssembler<FEBasis,TargetSpace> rodAssembler(gridView, localStiffness);
+
+ /////////////////////////////////////////////
+ // Create a solver for the rod problem
+ /////////////////////////////////////////////
RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
@@ -162,87 +224,56 @@ int main (int argc, char *argv[]) try
// //////////////////////////////
// Output result
// //////////////////////////////
- CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "rod3d-result");
-
- BlockVector<FieldVector<double, 6> > strain(x.size()-1);
- rodAssembler.getStrain(x,strain);
-
- // If convergence measurement is not desired stop here
- if (!instrumented)
- exit(0);
-
- // //////////////////////////////////////////////////////////
- // Recompute and compare against exact solution
- // //////////////////////////////////////////////////////////
-
- SolutionType exactSolution = x;
-
- // //////////////////////////////////////////////////////////
- // Compute hessian of the rod functional at the exact solution
- // for use of the energy norm it creates.
- // //////////////////////////////////////////////////////////
-
- BCRSMatrix<FieldMatrix<double, 6, 6> > hessian;
- MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
- rodAssembler.getNeighborsPerVertex(indices);
- indices.exportIdx(hessian);
- BlockVector<FieldVector<double,6> > dummyRhs(x.size());
- rodAssembler.assembleGradientAndHessian(exactSolution, dummyRhs, hessian);
-
-
- double error = std::numeric_limits<double>::max();
-
- SolutionType intermediateSolution(x.size());
-
- // Create statistics file
- std::ofstream statisticsFile((resultPath + "trStatistics").c_str());
-
- // Compute error of the initial iterate
- typedef BlockVector<FieldVector<double,6> > RodDifferenceType;
- RodDifferenceType rodDifference = computeGeodesicDifference(exactSolution, initialIterate);
- double oldError = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(rodDifference, hessian));
-
- int i;
- for (i=0; i<maxTrustRegionSteps; i++) {
-
- // /////////////////////////////////////////////////////
- // Read iteration from file
- // /////////////////////////////////////////////////////
- char iSolFilename[100];
- sprintf(iSolFilename, "tmp/intermediateSolution_%04d", i);
-
- FILE* fp = fopen(iSolFilename, "rb");
- if (!fp)
- DUNE_THROW(IOError, "Couldn't open intermediate solution '" << iSolFilename << "'");
- for (size_t j=0; j<intermediateSolution.size(); j++) {
- fread(&intermediateSolution[j].r, sizeof(double), 3, fp);
- fread(&intermediateSolution[j].q, sizeof(double), 4, fp);
- }
-
- fclose(fp);
-
- // /////////////////////////////////////////////////////
- // Compute error
- // /////////////////////////////////////////////////////
-
- rodDifference = computeGeodesicDifference(exactSolution, intermediateSolution);
-
- error = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(rodDifference, hessian));
-
-
- double convRate = error / oldError;
-
- // Output
- std::cout << "Trust-region iteration: " << i << " error : " << error << ", "
- << "convrate " << convRate << std::endl;
- statisticsFile << i << " " << error << " " << convRate << std::endl;
-
- if (error < 1e-12)
- break;
+#if HAVE_DUNE_VTK
+ using DataCollector = Vtk::LagrangeDataCollector<GridView,order>;
+ DataCollector dataCollector(gridView);
+ VtkUnstructuredGridWriter<GridView,DataCollector> vtkWriter(gridView, Vtk::ASCII);
+
+ // Make basis for R^3-valued data
+ using namespace Functions::BasisFactory;
+
+ auto worldBasis = makeBasis(
+ gridView,
+ power<3>(lagrange<order>())
+ );
+
+ // The rod displacement field
+ BlockVector<FieldVector<double,3> > displacement(worldBasis.size());
+ for (std::size_t i=0; i<x.size(); i++)
+ displacement[i] = x[i].r;
+
+ std::vector<double> xEmbedding;
+ Functions::interpolate(feBasis, xEmbedding, [](FieldVector<double,1> x){ return x; });
+
+ BlockVector<FieldVector<double,3> > gridEmbedding(xEmbedding.size());
+ for (std::size_t i=0; i<gridEmbedding.size(); i++)
+ gridEmbedding[i] = {xEmbedding[i], 0, 0};
+
+ displacement -= gridEmbedding;
+
+ auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(worldBasis, displacement);
+ vtkWriter.addPointData(displacementFunction, "displacement", 3);
+
+ // The three director fields
+ using FunctionType = decltype(displacementFunction);
+ std::array<std::optional<FunctionType>, 3> directorFunction;
+ std::array<BlockVector<FieldVector<double, 3> >, 3> director;
+ for (int i=0; i<3; i++)
+ {
+ director[i].resize(worldBasis.size());
+ for (std::size_t j=0; j<x.size(); j++)
+ director[i][j] = x[j].q.director(i);
+
+ directorFunction[i] = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(worldBasis, std::move(director[i]));
+ vtkWriter.addPointData(*directorFunction[i], "director " + std::to_string(i), 3);
+ }
- oldError = error;
-
- }
+ vtkWriter.write(resultPath + "rod3d-result");
+#else
+ std::cout << "Falling back to legacy file writing. Get dune-vtk for better results" << std::endl;
+ // Fall-back solution for users without dune-vtk
+ CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "rod3d-result");
+#endif
} catch (Exception& e)
{
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e1952d1af6c8993a4dcc161acec079f6a3436980..042be43e0174cf40caf74ceea13da7cbb8140169 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -2,14 +2,13 @@ set(TESTS
adolctest
averagedistanceassemblertest
cosseratenergytest
- frameinvariancetest
+ cosseratrodenergytest
harmonicenergytest
localgeodesicfefunctiontest
localgfetestfunctiontest
localprojectedfefunctiontest
orthogonalmatrixtest
rigidbodymotiontest
- rodassemblertest
rotationtest
svdtest
targetspacetest
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 256fba6fe2adfd69f9ff72019df1b60d9071c1f1..cdfc0f38c9d95afdf49e1aadd9cb380d69296214 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -15,8 +15,6 @@
#include "multiindex.hh"
#include "valuefactory.hh"
-const double eps = 1e-4;
-
typedef RigidBodyMotion<double,3> TargetSpace;
using namespace Dune;
diff --git a/test/frameinvariancetest.cc b/test/cosseratrodenergytest.cc
similarity index 60%
rename from test/frameinvariancetest.cc
rename to test/cosseratrodenergytest.cc
index dcf53d480975d78a4b2fce61765376dd61a5e121..ff65271b6c6ed816cded185956c0158e4b1a1fd9 100644
--- a/test/frameinvariancetest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -4,15 +4,10 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
-#include <dune/gfe/quaternion.hh>
-
-#include <dune/gfe/rodassembler.hh>
-
+#include <dune/gfe/cosseratrodenergy.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/rigidbodymotion.hh>
-// Number of degrees of freedom:
-// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
-const int blocksize = 6;
using namespace Dune;
@@ -74,14 +69,36 @@ int main (int argc, char *argv[]) try
rotatedX[i].q = rotation.mult(x[i].q);
}
- RodLocalStiffness<GridView,double> localRodFirstOrderModel(gridView,
- 1,1,1,1e6,0.3);
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<1, double,
+ FEBasis::LocalView::Tree::FiniteElement,
+ RigidBodyMotion<double,3> >;
+
+ GFE::CosseratRodEnergy<FEBasis,
+ GeodesicInterpolationRule,
+ double> localRodEnergy(gridView,
+ 1,1,1,1e6,0.3);
+
+ std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(1));
+
+ for (const auto vertex : vertices(gridView))
+ {
+ auto idx = gridView.indexSet().index(vertex);
+
+ referenceConfiguration[idx].r[0] = 0;
+ referenceConfiguration[idx].r[1] = 0;
+ referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
+ referenceConfiguration[idx].q = Rotation<double,3>::identity();
+ }
+
+ localRodEnergy.setReferenceConfiguration(referenceConfiguration);
- LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localRodFirstOrderModel);
+ auto localView = feBasis.localView();
+ localView.bind(*gridView.begin<0>());
- RodAssembler<FEBasis,3> assembler(feBasis, localFDStiffness);
+ SolutionType localX = {x[0], x[1]};
+ SolutionType localRotatedX = {rotatedX[0], rotatedX[1]};
- if (std::abs(assembler.computeEnergy(x) - assembler.computeEnergy(rotatedX)) > 1e-6)
+ if (std::abs(localRodEnergy.energy(localView, localX) - localRodEnergy.energy(localView, localRotatedX)) > 1e-6)
DUNE_THROW(Dune::Exception, "Rod energy not invariant under rigid body motions!");
} catch (Exception& e) {
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 3e78813537d2fc681af3723c304491920ccfc1f0..36250b3c9480635914866ebbd654feeb9e12ee25 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -114,7 +114,6 @@ int main (int argc, char *argv[])
));
using CompositeBasis = decltype(compositeBasis);
- using LocalView = typename CompositeBasis::LocalView;
/////////////////////////////////////////////////////////////////////////
// Create the energy functions with their parameters
@@ -170,7 +169,7 @@ int main (int argc, char *argv[])
x[_0].resize(compositeBasis.size({0}));
x[_1].resize(compositeBasis.size({1}));
std::vector<RBM> xRBM(compositeBasis.size({0}));
- for (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
for (int j = 0; j < gridDim; j++)
initialDeformation[i][j] = identity[i][j];
x[_0][i] = initialDeformation[i];
diff --git a/test/rodassemblertest.cc b/test/rodassemblertest.cc
deleted file mode 100644
index 554007834a7010f3f3a0eb1f35572d9ac4d88685..0000000000000000000000000000000000000000
--- a/test/rodassemblertest.cc
+++ /dev/null
@@ -1,587 +0,0 @@
-#include <config.h>
-
-#include <array>
-
-#include <dune/grid/onedgrid.hh>
-
-#include <dune/istl/io.hh>
-
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/rodassembler.hh>
-
-// Number of degrees of freedom:
-// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
-const int blocksize = 6;
-
-using namespace Dune;
-
-void infinitesimalVariation(RigidBodyMotion<double,3>& c, double eps, int i)
-{
- if (i<3)
- c.r[i] += eps;
- else {
- Dune::FieldVector<double,3> axial(0);
- axial[i-3] = eps;
- SkewMatrix<double,3> variation(axial);
- c.q = c.q.mult(Rotation<double,3>::exp(variation));
- }
-}
-
-template <class FEBasis>
-void strainFD(const std::vector<RigidBodyMotion<double,3> >& x,
- double pos,
- std::array<FieldMatrix<double,2,6>, 6>& fdStrainDerivatives,
- const RodAssembler<FEBasis,3>& assembler)
-{
- assert(x.size()==2);
- double eps = 1e-8;
-
- auto element = assembler.basis_.gridView_.template begin<0>();
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
-
- for (size_t i=0; i<x.size(); i++) {
-
- Dune::FieldVector<double,6> fdGradient;
-
- for (int j=0; j<6; j++) {
-
- infinitesimalVariation(forwardSolution[i], eps, j);
- infinitesimalVariation(backwardSolution[i], -eps, j);
-
-// fdGradient[j] = (assembler.computeEnergy(forwardSolution) - assembler.computeEnergy(backwardSolution))
-// / (2*eps);
- Dune::FieldVector<double,6> strain;
- strain = assembler.getStrain(forwardSolution, element, pos);
- strain -= assembler.getStrain(backwardSolution, element, pos);
- strain /= 2*eps;
-
- for (int m=0; m<6; m++)
- fdStrainDerivatives[m][i][j] = strain[m];
-
- forwardSolution[i] = x[i];
- backwardSolution[i] = x[i];
- }
-
- }
-
-}
-
-
-template <class FEBasis>
-void strain2ndOrderFD(const std::vector<RigidBodyMotion<double,3> >& x,
- double pos,
- std::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer,
- std::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer,
- const RodAssembler<FEBasis,3>& assembler)
-{
- assert(x.size()==2);
- double eps = 1e-3;
-
- auto element = assembler.basis_.gridView_.template begin<0>();
-
- for (int m=0; m<3; m++) {
-
- translationDer[m].setSize(2,2);
- translationDer[m] = 0;
-
- rotationDer[m].setSize(2,2);
- rotationDer[m] = 0;
-
- }
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
-
- std::vector<RigidBodyMotion<double,3> > forwardForwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > forwardBackwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardForwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardBackwardSolution = x;
-
- for (int i=0; i<2; i++) {
-
- for (int j=0; j<3; j++) {
-
- for (int k=0; k<2; k++) {
-
- for (int l=0; l<3; l++) {
-
- if (i==k && j==l) {
-
- infinitesimalVariation(forwardSolution[i], eps, j+3);
- infinitesimalVariation(backwardSolution[i], -eps, j+3);
-
- // Second derivative
-// fdHessian[j][k] = (assembler.computeEnergy(forwardSolution)
-// - 2*assembler.computeEnergy(x)
-// + assembler.computeEnergy(backwardSolution)) / (eps*eps);
-
- Dune::FieldVector<double,6> strain;
- strain = assembler.getStrain(forwardSolution, element, pos);
- strain += assembler.getStrain(backwardSolution, element, pos);
- strain.axpy(-2, assembler.getStrain(x, element, pos));
- strain /= eps*eps;
-
- for (int m=0; m<3; m++)
- rotationDer[m][i][k][j][l] = strain[3+m];
-
- forwardSolution = x;
- backwardSolution = x;
-
- } else {
-
- infinitesimalVariation(forwardForwardSolution[i], eps, j+3);
- infinitesimalVariation(forwardForwardSolution[k], eps, l+3);
- infinitesimalVariation(forwardBackwardSolution[i], eps, j+3);
- infinitesimalVariation(forwardBackwardSolution[k], -eps, l+3);
- infinitesimalVariation(backwardForwardSolution[i], -eps, j+3);
- infinitesimalVariation(backwardForwardSolution[k], eps, l+3);
- infinitesimalVariation(backwardBackwardSolution[i],-eps, j+3);
- infinitesimalVariation(backwardBackwardSolution[k],-eps, l+3);
-
-
-// fdHessian[j][k] = (assembler.computeEnergy(forwardForwardSolution)
-// + assembler.computeEnergy(backwardBackwardSolution)
-// - assembler.computeEnergy(forwardBackwardSolution)
-// - assembler.computeEnergy(backwardForwardSolution)) / (4*eps*eps);
-
- Dune::FieldVector<double,6> strain;
- strain = assembler.getStrain(forwardForwardSolution, element, pos);
- strain += assembler.getStrain(backwardBackwardSolution, element, pos);
- strain -= assembler.getStrain(forwardBackwardSolution, element, pos);
- strain -= assembler.getStrain(backwardForwardSolution, element, pos);
- strain /= 4*eps*eps;
-
-
- for (int m=0; m<3; m++)
- rotationDer[m][i][k][j][l] = strain[3+m];
-
- forwardForwardSolution = x;
- forwardBackwardSolution = x;
- backwardForwardSolution = x;
- backwardBackwardSolution = x;
-
-
- }
-
- }
-
- }
- }
-
- }
-
-}
-
-
-template <class GridType>
-void strain2ndOrderFD2(const std::vector<RigidBodyMotion<double,3> >& x,
- double pos,
- Dune::FieldVector<double,1> shapeGrad[2],
- Dune::FieldVector<double,1> shapeFunction[2],
- std::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer,
- std::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer,
- const RodAssembler<GridType,3>& assembler)
-{
- assert(x.size()==2);
- double eps = 1e-3;
-
- for (int m=0; m<3; m++) {
-
- translationDer[m].setSize(2,2);
- translationDer[m] = 0;
-
- rotationDer[m].setSize(2,2);
- rotationDer[m] = 0;
-
- }
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
-
- for (int k=0; k<2; k++) {
-
- for (int l=0; l<3; l++) {
-
- infinitesimalVariation(forwardSolution[k], eps, l+3);
- infinitesimalVariation(backwardSolution[k], -eps, l+3);
-
- std::array<FieldMatrix<double,2,6>, 6> forwardStrainDer;
- assembler.strainDerivative(forwardSolution, pos, shapeGrad, shapeFunction, forwardStrainDer);
-
- std::array<FieldMatrix<double,2,6>, 6> backwardStrainDer;
- assembler.strainDerivative(backwardSolution, pos, shapeGrad, shapeFunction, backwardStrainDer);
-
- for (int i=0; i<2; i++) {
- for (int j=0; j<3; j++) {
- for (int m=0; m<3; m++) {
- rotationDer[m][i][k][j][l] = (forwardStrainDer[m][i][j]-backwardStrainDer[m][i][j]) / (2*eps);
- }
- }
- }
-
- forwardSolution = x;
- backwardSolution = x;
-
- }
-
- }
-
-}
-
-
-
-
-
-template <class GridType>
-void expHessianFD()
-{
- using namespace Dune;
- double eps = 1e-3;
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- SkewMatrix<double,3> forward;
- SkewMatrix<double,3> backward;
-
- SkewMatrix<double,3> forwardForward;
- SkewMatrix<double,3> forwardBackward;
- SkewMatrix<double,3> backwardForward;
- SkewMatrix<double,3> backwardBackward;
-
- for (int i=0; i<3; i++) {
-
- for (int j=0; j<3; j++) {
-
- Quaternion<double> hessian;
-
- if (i==j) {
-
- forward = backward = 0;
- forward.axial()[i] += eps;
- backward.axial()[i] -= eps;
-
- // Second derivative
- // fdHessian[j][k] = (assembler.computeEnergy(forward)
- // - 2*assembler.computeEnergy(x)
- // + assembler.computeEnergy(backward)) / (eps*eps);
-
- hessian = Rotation<double,3>::exp(forward);
- hessian += Rotation<double,3>::exp(backward);
- SkewMatrix<double,3> zero(0);
- hessian.axpy(-2, Rotation<double,3>::exp(zero));
- hessian /= eps*eps;
-
- } else {
-
- forwardForward = forwardBackward = 0;
- backwardForward = backwardBackward = 0;
-
- forwardForward.axial()[i] += eps;
- forwardForward.axial()[j] += eps;
- forwardBackward.axial()[i] += eps;
- forwardBackward.axial()[j] -= eps;
- backwardForward.axial()[i] -= eps;
- backwardForward.axial()[j] += eps;
- backwardBackward.axial()[i] -= eps;
- backwardBackward.axial()[j] -= eps;
-
-
- hessian = Rotation<double,3>::exp(forwardForward);
- hessian += Rotation<double,3>::exp(backwardBackward);
- hessian -= Rotation<double,3>::exp(forwardBackward);
- hessian -= Rotation<double,3>::exp(backwardForward);
- hessian /= 4*eps*eps;
-
- }
-
- printf("i: %d, j: %d ", i, j);
- std::cout << hessian << std::endl;
-
- }
-
- }
-}
-
-
-template <class GridType>
-void gradientFDCheck(const std::vector<RigidBodyMotion<double,3> >& x,
- const Dune::BlockVector<Dune::FieldVector<double,6> >& gradient,
- const RodAssembler<GridType,3>& assembler)
-{
-#ifndef ABORT_ON_ERROR
- static int gradientError = 0;
- double maxError = -1;
-#endif
- double eps = 1e-6;
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
-
- std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
-
- for (size_t i=0; i<x.size(); i++) {
-
- Dune::FieldVector<double,6> fdGradient;
-
- for (int j=0; j<6; j++) {
-
- infinitesimalVariation(forwardSolution[i], eps, j);
- infinitesimalVariation(backwardSolution[i], -eps, j);
-
- fdGradient[j] = (assembler.computeEnergy(forwardSolution) - assembler.computeEnergy(backwardSolution))
- / (2*eps);
-
- forwardSolution[i] = x[i];
- backwardSolution[i] = x[i];
- }
-
- if ( (fdGradient-gradient[i]).two_norm() > 1e-6 ) {
-#ifdef ABORT_ON_ERROR
- std::cout << "Differing gradients at " << i
- << "! (" << (fdGradient-gradient[i]).two_norm() / gradient[i].two_norm()
- << ") Analytical: " << gradient[i] << ", fd: " << fdGradient << std::endl;
- //std::cout << "Current configuration " << std::endl;
- for (size_t j=0; j<x.size(); j++)
- std::cout << x[j] << std::endl;
- //abort();
-#else
- gradientError++;
- maxError = std::max(maxError, (fdGradient-gradient[i]).two_norm());
-#endif
- }
-
- }
-
-#ifndef ABORT_ON_ERROR
- std::cout << gradientError << " errors in the gradient corrected (max: " << maxError << ")!" << std::endl;
-#endif
-
-}
-
-
-template <class GridType>
-void hessianFDCheck(const std::vector<RigidBodyMotion<double,3> >& x,
- const Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >& hessian,
- const RodAssembler<GridType,3>& assembler)
-{
-#ifndef ABORT_ON_ERROR
- static int hessianError = 0;
-#endif
- double eps = 1e-2;
-
- typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >::row_type::const_iterator ColumnIterator;
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
-
- std::vector<RigidBodyMotion<double,3> > forwardForwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > forwardBackwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardForwardSolution = x;
- std::vector<RigidBodyMotion<double,3> > backwardBackwardSolution = x;
-
-
- // ///////////////////////////////////////////////////////////////
- // Loop over all blocks of the outer matrix
- // ///////////////////////////////////////////////////////////////
- for (size_t i=0; i<hessian.N(); i++) {
-
- ColumnIterator cIt = hessian[i].begin();
- ColumnIterator cEndIt = hessian[i].end();
-
- for (; cIt!=cEndIt; ++cIt) {
-
- // ////////////////////////////////////////////////////////////////////////////
- // Compute a finite-difference approximation of this hessian matrix block
- // ////////////////////////////////////////////////////////////////////////////
-
- Dune::FieldMatrix<double,6,6> fdHessian;
-
- for (int j=0; j<6; j++) {
-
- for (int k=0; k<6; k++) {
-
- if (i==cIt.index() && j==k) {
-
- infinitesimalVariation(forwardSolution[i], eps, j);
- infinitesimalVariation(backwardSolution[i], -eps, j);
-
- // Second derivative
- fdHessian[j][k] = (assembler.computeEnergy(forwardSolution)
- + assembler.computeEnergy(backwardSolution)
- - 2*assembler.computeEnergy(x)
- ) / (eps*eps);
-
- forwardSolution[i] = x[i];
- backwardSolution[i] = x[i];
-
- } else {
-
- infinitesimalVariation(forwardForwardSolution[i], eps, j);
- infinitesimalVariation(forwardForwardSolution[cIt.index()], eps, k);
- infinitesimalVariation(forwardBackwardSolution[i], eps, j);
- infinitesimalVariation(forwardBackwardSolution[cIt.index()], -eps, k);
- infinitesimalVariation(backwardForwardSolution[i], -eps, j);
- infinitesimalVariation(backwardForwardSolution[cIt.index()], eps, k);
- infinitesimalVariation(backwardBackwardSolution[i], -eps, j);
- infinitesimalVariation(backwardBackwardSolution[cIt.index()],-eps, k);
-
-
- fdHessian[j][k] = (assembler.computeEnergy(forwardForwardSolution)
- + assembler.computeEnergy(backwardBackwardSolution)
- - assembler.computeEnergy(forwardBackwardSolution)
- - assembler.computeEnergy(backwardForwardSolution)) / (4*eps*eps);
-
- forwardForwardSolution[i] = x[i];
- forwardForwardSolution[cIt.index()] = x[cIt.index()];
- forwardBackwardSolution[i] = x[i];
- forwardBackwardSolution[cIt.index()] = x[cIt.index()];
- backwardForwardSolution[i] = x[i];
- backwardForwardSolution[cIt.index()] = x[cIt.index()];
- backwardBackwardSolution[i] = x[i];
- backwardBackwardSolution[cIt.index()] = x[cIt.index()];
-
- }
-
- }
-
- }
- //exit(0);
- // /////////////////////////////////////////////////////////////////////////////
- // Compare analytical and fd Hessians
- // /////////////////////////////////////////////////////////////////////////////
-
- Dune::FieldMatrix<double,6,6> diff = fdHessian;
- diff -= *cIt;
-
- if ( diff.frobenius_norm() > 1e-5 ) {
-#ifdef ABORT_ON_ERROR
- std::cout << "Differing hessians at [(" << i << ", "<< cIt.index() << ")]!" << std::endl
- << "Analytical: " << std::endl << *cIt << std::endl
- << "fd: " << std::endl << fdHessian << std::endl;
- abort();
-#else
- hessianError++;
-#endif
- }
-
- }
-
- }
-
-#ifndef ABORT_ON_ERROR
- std::cout << hessianError << " errors in the Hessian corrected!" << std::endl;
-#endif
-
-}
-
-
-int main (int argc, char *argv[]) try
-{
- typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
-
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
- typedef OneDGrid GridType;
- GridType grid(1, 0, 1);
-
- using GridView = GridType::LeafGridView;
- GridView gridView = grid.leafGridView();
-
- SolutionType x(gridView.size(1));
-
- // //////////////////////////
- // Initial solution
- // //////////////////////////
- FieldVector<double,3> xAxis(0);
- xAxis[0] = 1;
- FieldVector<double,3> zAxis(0);
- zAxis[2] = 1;
-
- for (size_t i=0; i<x.size(); i++) {
- x[i].r[0] = 0; // x
- x[i].r[1] = 0; // y
- x[i].r[2] = double(i)/(x.size()-1); // z
- //x[i].r[2] = i+5;
- x[i].q = Quaternion<double>::identity();
- //x[i].q = Quaternion<double>(zAxis, M_PI/2 * double(i)/(x.size()-1));
- }
-
- //x.back().r[1] = 0.1;
- //x.back().r[2] = 2;
- //x.back().q = Quaternion<double>(zAxis, M_PI/4);
-
- std::cout << "Left boundary orientation:" << std::endl;
- std::cout << "director 0: " << x[0].q.director(0) << std::endl;
- std::cout << "director 1: " << x[0].q.director(1) << std::endl;
- std::cout << "director 2: " << x[0].q.director(2) << std::endl;
- std::cout << std::endl;
- std::cout << "Right boundary orientation:" << std::endl;
- std::cout << "director 0: " << x[x.size()-1].q.director(0) << std::endl;
- std::cout << "director 1: " << x[x.size()-1].q.director(1) << std::endl;
- std::cout << "director 2: " << x[x.size()-1].q.director(2) << std::endl;
-// exit(0);
-
- //x[0].r[2] = -1;
-
- // ///////////////////////////////////////////
- // Create a solver for the rod problem
- // ///////////////////////////////////////////
-
- using Basis = Functions::LagrangeBasis<GridView,1>;
- Basis basis(gridView);
-
- RodLocalStiffness<GridView,double> localStiffness(gridView,
- 0.01, 0.0001, 0.0001, 2.5e5, 0.3);
-
- LocalGeodesicFEFDStiffness<Basis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
-
- RodAssembler<Basis,3> rodAssembler(basis, localFDStiffness);
-
- std::cout << "Energy: " << rodAssembler.computeEnergy(x) << std::endl;
-
- double pos = (argc==2) ? atof(argv[1]) : 0.5;
-
- FieldVector<double,1> shapeGrad[2];
- shapeGrad[0] = -1;
- shapeGrad[1] = 1;
-
- FieldVector<double,1> shapeFunction[2];
- shapeFunction[0] = 1-pos;
- shapeFunction[1] = pos;
-
- exit(0);
- BlockVector<FieldVector<double,6> > rhs(x.size());
- BCRSMatrix<FieldMatrix<double,6,6> > hessianMatrix;
- MatrixIndexSet indices(grid.size(1), grid.size(1));
- rodAssembler.getNeighborsPerVertex(indices);
- indices.exportIdx(hessianMatrix);
-
- rodAssembler.assembleGradientAndHessian(x, rhs, hessianMatrix);
-
- gradientFDCheck(x, rhs, rodAssembler);
- hessianFDCheck(x, hessianMatrix, rodAssembler);
-
- // //////////////////////////////
- } catch (Exception& e) {
-
- std::cout << e.what() << std::endl;
-
- }
diff --git a/test/rotationtest.cc b/test/rotationtest.cc
index ac4f2f6a32aeceb42bd568c64157ad318e595458..7b51c0a049f6463c74aae36205c72e20989b574b 100644
--- a/test/rotationtest.cc
+++ b/test/rotationtest.cc
@@ -5,12 +5,7 @@
#include <dune/common/fmatrix.hh>
-#warning Do not include onedgrid.hh
-#include <dune/grid/onedgrid.hh>
-
#include <dune/gfe/rotation.hh>
-#warning Do not include rodlocalstiffness.hh
-#include <dune/gfe/rodlocalstiffness.hh>
#include "valuefactory.hh"
using namespace Dune;
@@ -94,123 +89,6 @@ void testDDExp()
}
}
-void testDerivativeOfInterpolatedPosition()
-{
- std::array<Rotation<double,3>, 6> q;
-
- FieldVector<double,3> xAxis = {1,0,0};
- FieldVector<double,3> yAxis = {0,1,0};
- FieldVector<double,3> zAxis = {0,0,1};
-
- q[0] = Rotation<double,3>(xAxis, 0);
- q[1] = Rotation<double,3>(xAxis, M_PI/2);
- q[2] = Rotation<double,3>(yAxis, 0);
- q[3] = Rotation<double,3>(yAxis, M_PI/2);
- q[4] = Rotation<double,3>(zAxis, 0);
- q[5] = Rotation<double,3>(zAxis, M_PI/2);
-
- double eps = 1e-7;
-
- for (int i=0; i<6; i++) {
-
- for (int j=0; j<6; j++) {
-
- for (int k=0; k<7; k++) {
-
- double s = k/6.0;
-
- std::array<Quaternion<double>,6> fdGrad;
-
- // ///////////////////////////////////////////////////////////
- // First: test the interpolated position
- // ///////////////////////////////////////////////////////////
- for (int l=0; l<3; l++) {
- SkewMatrix<double,3> forward(FieldVector<double,3>(0));
- SkewMatrix<double,3> backward(FieldVector<double,3>(0));
- forward.axial()[l] += eps;
- backward.axial()[l] -= eps;
- fdGrad[l] = Rotation<double,3>::interpolate(q[i].mult(Rotation<double,3>::exp(forward)), q[j], s);
- fdGrad[l] -= Rotation<double,3>::interpolate(q[i].mult(Rotation<double,3>::exp(backward)), q[j], s);
- fdGrad[l] /= 2*eps;
-
- fdGrad[3+l] = Rotation<double,3>::interpolate(q[i], q[j].mult(Rotation<double,3>::exp(forward)), s);
- fdGrad[3+l] -= Rotation<double,3>::interpolate(q[i], q[j].mult(Rotation<double,3>::exp(backward)), s);
- fdGrad[3+l] /= 2*eps;
- }
-
- // Compute analytical gradient
- std::array<Quaternion<double>,6> grad;
- RodLocalStiffness<OneDGrid,double>::interpolationDerivative(q[i], q[j], s, grad);
-
- for (int l=0; l<6; l++) {
- Quaternion<double> diff = fdGrad[l];
- diff -= grad[l];
- if (diff.two_norm() > 1e-6) {
- std::cout << "Error in position " << l << ": fd: " << fdGrad[l]
- << " analytical: " << grad[l] << std::endl;
- }
-
- }
-
- // ///////////////////////////////////////////////////////////
- // Second: test the interpolated velocity vector
- // ///////////////////////////////////////////////////////////
-
- for (const double intervalLength : {1.0/3, 2.0/3, 3.0/3, 4.0/3, 5.0/3, 6.0/3, 7.0/3})
- {
- Dune::FieldVector<double,3> variation;
-
- for (int m=0; m<3; m++) {
- variation = 0;
- variation[m] = eps;
- fdGrad[m] = Rotation<double,3>::interpolateDerivative(q[i].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))),
- q[j], s);
- variation = 0;
- variation[m] = -eps;
- fdGrad[m] -= Rotation<double,3>::interpolateDerivative(q[i].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))),
- q[j], s);
- fdGrad[m] /= 2*eps;
-
- variation = 0;
- variation[m] = eps;
- fdGrad[3+m] = Rotation<double,3>::interpolateDerivative(q[i], q[j].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))), s);
- variation = 0;
- variation[m] = -eps;
- fdGrad[3+m] -= Rotation<double,3>::interpolateDerivative(q[i], q[j].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))), s);
- fdGrad[3+m] /= 2*eps;
-
- }
-
- // Scale the finite difference gradient with the interval length
- for (auto& parDer : fdGrad)
- parDer /= intervalLength;
-
- // Compute analytical velocity vector gradient
- RodLocalStiffness<OneDGrid,double>::interpolationVelocityDerivative(q[i], q[j], s*intervalLength, intervalLength, grad);
-
- for (int m=0; m<6; m++) {
- Quaternion<double> diff = fdGrad[m];
- diff -= grad[m];
- if (diff.two_norm() > 1e-6) {
- std::cout << "Error in velocity " << m
- << ": s = " << s << " of (" << intervalLength << ")"
- << " fd: " << fdGrad[m] << " analytical: " << grad[m] << std::endl;
- }
-
- }
-
- }
-
- }
-
- }
-
- }
-
-}
-
-
-
void testRotation(Rotation<double,3> q)
{
// Make sure it really is a unit quaternion
@@ -388,11 +266,6 @@ int main (int argc, char *argv[]) try
// //////////////////////////////////////////////
testDDExp();
- // //////////////////////////////////////////////
- // Test derivative of interpolated position
- // //////////////////////////////////////////////
- testDerivativeOfInterpolatedPosition();
-
return not passed;
} catch (Exception& e) {