From 663c54cdcf25d98d2501ea3d9dac071229328e1a Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 8 Nov 2020 21:19:36 +0100
Subject: [PATCH] Remove the RodAssembler class
GeodesicFEAssembler can replace it without problems.
---
dune/gfe/rodassembler.cc | 430 --------------------------
dune/gfe/rodassembler.hh | 148 ---------
src/rod3d.cc | 36 ++-
test/CMakeLists.txt | 1 -
test/frameinvariancetest.cc | 33 +-
test/rodassemblertest.cc | 587 ------------------------------------
test/rotationtest.cc | 127 --------
7 files changed, 49 insertions(+), 1313 deletions(-)
delete mode 100644 dune/gfe/rodassembler.cc
delete mode 100644 dune/gfe/rodassembler.hh
delete mode 100644 test/rodassemblertest.cc
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
deleted file mode 100644
index 01b5e2d7..00000000
--- a/dune/gfe/rodassembler.cc
+++ /dev/null
@@ -1,430 +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,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 70ab003e..00000000
--- a/dune/gfe/rodassembler.hh
+++ /dev/null
@@ -1,148 +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;
-
- }; // 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/src/rod3d.cc b/src/rod3d.cc
index e803181a..26054037 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -13,9 +13,11 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/localgeodesicfefdstiffness.hh>
#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/rodlocalstiffness.hh>
#include <dune/gfe/rotation.hh>
-#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
typedef RigidBodyMotion<double,3> TargetSpace;
@@ -118,16 +120,35 @@ int main (int argc, char *argv[]) try
dirichletNodes[0] = true;
dirichletNodes.back() = true;
+ //////////////////////////////////////////////
+ // Create the stress-free configuration
+ //////////////////////////////////////////////
+
+ auto localRodEnergy = std::make_shared<RodLocalStiffness<GridView, double> >(gridView,
+ A, J1, J2, E, nu);
+
+ 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);
+
// ///////////////////////////////////////////
// Create a solver for the rod problem
// ///////////////////////////////////////////
- RodLocalStiffness<GridView,double> localStiffness(gridView,
- A, J1, J2, E, nu);
+ LocalGeodesicFEFDStiffness<FEBasis,
+ TargetSpace> localStiffness(localRodEnergy.get());
- LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
-
- RodAssembler<FEBasis,3> rodAssembler(gridView, localFDStiffness);
+ GeodesicFEAssembler<FEBasis,TargetSpace> rodAssembler(gridView, localStiffness);
RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
@@ -161,9 +182,6 @@ int main (int argc, char *argv[]) try
// //////////////////////////////
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "rod3d-result");
- BlockVector<FieldVector<double, 6> > strain(x.size()-1);
- rodAssembler.getStrain(x,strain);
-
} catch (Exception& e)
{
std::cout << e.what() << std::endl;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e1952d1a..6b82573a 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,7 +9,6 @@ set(TESTS
localprojectedfefunctiontest
orthogonalmatrixtest
rigidbodymotiontest
- rodassemblertest
rotationtest
svdtest
targetspacetest
diff --git a/test/frameinvariancetest.cc b/test/frameinvariancetest.cc
index dcf53d48..fe1c1ccc 100644
--- a/test/frameinvariancetest.cc
+++ b/test/frameinvariancetest.cc
@@ -5,14 +5,9 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/gfe/quaternion.hh>
-
-#include <dune/gfe/rodassembler.hh>
-
#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/rodlocalstiffness.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,30 @@ 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);
+ RodLocalStiffness<GridView,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/rodassemblertest.cc b/test/rodassemblertest.cc
deleted file mode 100644
index 55400783..00000000
--- 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 ac4f2f6a..7b51c0a0 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) {
--
GitLab