diff --git a/src/Makefile.am b/src/Makefile.am
index 38eb13a84681803ea6f508bcb048b66888d68c9e..3a4555856fae0c9d5d9c2f1808d446f0d5d37f46 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -11,6 +11,6 @@ srcinclude_HEADERS = averagedistanceassembler.hh localgeodesicfefunction.hh \
geodesicdifference.hh makestraightrod.hh rigidbodymotion.hh \
rotation.hh geodesicfeassembler.hh maxnormtrustregion.hh \
rodassembler.hh svd.hh harmonicenergystiffness.hh \
- planarrodassembler.hh rodlocalstiffness.hh targetspacertrsolver.hh
+ rodlocalstiffness.hh targetspacertrsolver.hh
include $(top_srcdir)/am/global-rules
diff --git a/src/planarrodassembler.cc b/src/planarrodassembler.cc
deleted file mode 100644
index a3ee6a92377eff6e1af14e5d18e38fd199ff2483..0000000000000000000000000000000000000000
--- a/src/planarrodassembler.cc
+++ /dev/null
@@ -1,437 +0,0 @@
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-
-#include <dune/grid/common/quadraturerules.hh>
-
-#include <dune/localfunctions/lagrange/p1.hh>
-
-template <class GridType, int polOrd>
-void Dune::PlanarRodAssembler<GridType, polOrd>::
-getNeighborsPerVertex(MatrixIndexSet& nb) const
-{
- const int gridDim = GridType::dimension;
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
-
- int i, j;
- int n = grid_->size(grid_->maxLevel(), gridDim);
-
- nb.resize(n, n);
-
- ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
- ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
-
- for (; it!=endit; ++it) {
-
- for (i=0; i<it->template count<gridDim>(); i++) {
-
- for (j=0; j<it->template count<gridDim>(); j++) {
-
- int iIdx = indexSet.subIndex(*it,i,gridDim);
- int jIdx = indexSet.subIndex(*it,j,gridDim);
-
- nb.add(iIdx, jIdx);
-
- }
-
- }
-
- }
-
-}
-
-template <class GridType, int polOrd>
-void Dune::PlanarRodAssembler<GridType, polOrd>::
-assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
- BCRSMatrix<MatrixBlock>& matrix)
-{
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
-
- MatrixIndexSet neighborsPerVertex;
- getNeighborsPerVertex(neighborsPerVertex);
-
- matrix = 0;
-
- ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
- ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
-
- Matrix<MatrixBlock> mat;
-
- for( ; it != endit; ++it ) {
-
- const int numOfBaseFct = 2;
-
- mat.setSize(numOfBaseFct, numOfBaseFct);
-
- // Extract local solution
- std::vector<RigidBodyMotion<2> > localSolution(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
-
- // setup matrix
- getLocalMatrix( *it, localSolution, numOfBaseFct, mat);
-
- // Add element matrix to global stiffness matrix
- for(int i=0; i<numOfBaseFct; i++) {
-
- int row = indexSet.subIndex(*it,i,gridDim);
-
- for (int j=0; j<numOfBaseFct; j++ ) {
-
- int col = indexSet.subIndex(*it,j,gridDim);
- matrix[row][col] += mat[i][j];
-
- }
- }
-
- }
-
-}
-
-
-
-
-
-
-template <class GridType, int polOrd>
-template <class MatrixType>
-void Dune::PlanarRodAssembler<GridType, polOrd>::
-getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<2> >& localSolution,
- const int matSize, MatrixType& localMat) const
-{
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
-
- /* ndof is the number of vectors of the element */
- int ndof = matSize;
-
- for (int i=0; i<matSize; i++)
- for (int j=0; j<matSize; j++)
- localMat[i][j] = 0;
-
- P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
-
- // Get quadrature rule
- const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity.type(), polOrd);
-
- /* Loop over all integration points */
- for (int ip=0; ip<quad.size(); ip++) {
-
- // Local position of the quadrature point
- const FieldVector<double,gridDim>& quadPos = quad[ip].position();
-
- // calc Jacobian inverse before integration element is evaluated
- const 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<FieldMatrix<double,1,gridDim> > referenceElementGradients(ndof);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<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<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<matSize; i++) {
-
- for (int j=0; j<matSize; 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 GridType, int polOrd>
-void Dune::PlanarRodAssembler<GridType, polOrd>::
-assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
- BlockVector<FieldVector<double, blocksize> >& grad) const
-{
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
- const int maxlevel = grid_->maxLevel();
-
- if (sol.size()!=grid_->size(maxlevel, gridDim))
- DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
-
- grad.resize(sol.size());
- grad = 0;
-
- ElementIterator it = grid_->template lbegin<0>(maxlevel);
- ElementIterator endIt = grid_->template lend<0>(maxlevel);
-
- // Loop over all elements
- for (; it!=endIt; ++it) {
-
- // Extract local solution on this element
- P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
- const int numOfBaseFct = localFiniteElement.localBasis().size();
-
- RigidBodyMotion<2> localSolution[numOfBaseFct];
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
-
- // Get quadrature rule
- const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->type(), polOrd);
-
- for (int pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = it->geometry().integrationElement(quadPos);
-
- double weight = quad[pt].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- std::vector<FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<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<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 = indexSet.subIndex(*it,dof,gridDim);
-
- //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 GridType, int polOrd>
-double Dune::PlanarRodAssembler<GridType, polOrd>::
-computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
-{
- const int maxlevel = grid_->maxLevel();
- double energy = 0;
-
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(maxlevel);
-
- if (sol.size()!=grid_->size(maxlevel, gridDim))
- DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
-
- ElementIterator it = grid_->template lbegin<0>(maxlevel);
- ElementIterator endIt = grid_->template lend<0>(maxlevel);
-
- // Loop over all elements
- for (; it!=endIt; ++it) {
-
- // Extract local solution on this element
- P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
-
- int numOfBaseFct = localFiniteElement.localBasis().size();
-
- RigidBodyMotion<2> localSolution[numOfBaseFct];
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
-
- // Get quadrature rule
- const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->type(), polOrd);
-
- for (int pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = it->geometry().integrationElement(quadPos);
-
- double weight = quad[pt].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- std::vector<FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
- localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
-
- std::vector<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<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/src/planarrodassembler.hh b/src/planarrodassembler.hh
deleted file mode 100644
index 20bf1a199aec3b04878a1936a41d013666717cf0..0000000000000000000000000000000000000000
--- a/src/planarrodassembler.hh
+++ /dev/null
@@ -1,94 +0,0 @@
-#ifndef DUNE_PLANAR_ROD_ASSEMBLER_HH
-#define DUNE_PLANAR_ROD_ASSEMBLER_HH
-
-#include <dune/common/fmatrix.hh>
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-
-#include "rigidbodymotion.hh"
-
-namespace Dune
-{
-
- /** \brief The FEM operator for an extensible, shearable rod
- */
- template <class GridType, int polOrd>
- class PlanarRodAssembler {
-
- typedef typename GridType::template Codim<0>::Entity EntityType;
- typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
-
- //! Dimension of the grid. This needs to be one!
- enum { gridDim = GridType::dimension };
-
- enum { elementOrder = 1};
-
- //! Each block is x, y, theta
- enum { blocksize = 3 };
-
- //!
- typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
- const GridType* grid_;
-
- /** \brief Material constants */
- double B;
- double A1;
- double A3;
-
- public:
-
- //! ???
- PlanarRodAssembler(const GridType &grid) :
- grid_(&grid)
- {
- B = 1;
- A1 = 1;
- A3 = 1;
- }
-
- ~PlanarRodAssembler() {}
-
- 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<2> >& sol,
- BCRSMatrix<MatrixBlock>& matrix);
-
- void assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
- BlockVector<FieldVector<double, blocksize> >& grad) const;
-
- /** \brief Compute the energy of a deformation state */
- double computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const;
-
- void getNeighborsPerVertex(MatrixIndexSet& nb) const;
-
- protected:
-
- /** \brief Compute the element tangent stiffness matrix */
- template <class MatrixType>
- void getLocalMatrix( EntityType &entity,
- const std::vector<RigidBodyMotion<2> >& localSolution,
- const int matSize, MatrixType& mat) const;
-
-
-
-
-
-
-
- }; // end class
-
-} // end namespace
-
-#include "planarrodassembler.cc"
-
-#endif
-
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 4687812c31b94f189b6eea94e9068b770a18f035..47c1d20b563efc5587d9dee07707c95b57c0e622 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -6,6 +6,7 @@
#include <dune/grid/common/quadraturerules.hh>
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
+#include <dune/localfunctions/lagrange/p1.hh>
#include "src/rodlocalstiffness.hh"
@@ -251,3 +252,431 @@ getResultantForce(const BoundaryPatchBase<PatchGridView>& boundary,
return canonicalStress;
}
+
+template <class GridType, int polOrd>
+void PlanarRodAssembler<GridType, polOrd>::getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
+{
+ const int gridDim = GridType::dimension;
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+
+ int i, j;
+ int n = grid_->size(grid_->maxLevel(), gridDim);
+
+ nb.resize(n, n);
+
+ ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
+ ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
+
+ for (; it!=endit; ++it) {
+
+ for (i=0; i<it->template count<gridDim>(); i++) {
+
+ for (j=0; j<it->template count<gridDim>(); j++) {
+
+ int iIdx = indexSet.subIndex(*it,i,gridDim);
+ int jIdx = indexSet.subIndex(*it,j,gridDim);
+
+ nb.add(iIdx, jIdx);
+
+ }
+
+ }
+
+ }
+
+}
+
+template <class GridType, int polOrd>
+void PlanarRodAssembler<GridType, polOrd>::
+assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
+ Dune::BCRSMatrix<MatrixBlock>& matrix)
+{
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+
+ Dune::MatrixIndexSet neighborsPerVertex;
+ getNeighborsPerVertex(neighborsPerVertex);
+
+ matrix = 0;
+
+ ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
+ ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
+
+ Dune::Matrix<MatrixBlock> mat;
+
+ for( ; it != endit; ++it ) {
+
+ const int numOfBaseFct = 2;
+
+ mat.setSize(numOfBaseFct, numOfBaseFct);
+
+ // Extract local solution
+ std::vector<RigidBodyMotion<2> > localSolution(numOfBaseFct);
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ // setup matrix
+ getLocalMatrix( *it, localSolution, numOfBaseFct, mat);
+
+ // Add element matrix to global stiffness matrix
+ for(int i=0; i<numOfBaseFct; i++) {
+
+ int row = indexSet.subIndex(*it,i,gridDim);
+
+ for (int j=0; j<numOfBaseFct; j++ ) {
+
+ int col = indexSet.subIndex(*it,j,gridDim);
+ matrix[row][col] += mat[i][j];
+
+ }
+ }
+
+ }
+
+}
+
+
+
+
+
+
+template <class GridType, int polOrd>
+template <class MatrixType>
+void PlanarRodAssembler<GridType, polOrd>::
+getLocalMatrix( EntityType &entity,
+ const std::vector<RigidBodyMotion<2> >& localSolution,
+ const int matSize, MatrixType& localMat) const
+{
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+
+ /* ndof is the number of vectors of the element */
+ int ndof = matSize;
+
+ for (int i=0; i<matSize; i++)
+ for (int j=0; j<matSize; j++)
+ localMat[i][j] = 0;
+
+ Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
+
+ // Get quadrature rule
+ const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(entity.type(), polOrd);
+
+ /* 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<matSize; i++) {
+
+ for (int j=0; j<matSize; 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 GridType, int polOrd>
+void PlanarRodAssembler<GridType, polOrd>::
+assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
+{
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+ const int maxlevel = grid_->maxLevel();
+
+ if (sol.size()!=grid_->size(maxlevel, gridDim))
+ DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
+
+ grad.resize(sol.size());
+ grad = 0;
+
+ ElementIterator it = grid_->template lbegin<0>(maxlevel);
+ ElementIterator endIt = grid_->template lend<0>(maxlevel);
+
+ // Loop over all elements
+ for (; it!=endIt; ++it) {
+
+ // Extract local solution on this element
+ Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
+ const int numOfBaseFct = localFiniteElement.localBasis().size();
+
+ RigidBodyMotion<2> localSolution[numOfBaseFct];
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ // Get quadrature rule
+ const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), polOrd);
+
+ 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 = it->geometry().jacobianInverseTransposed(quadPos);
+ const double integrationElement = it->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 = indexSet.subIndex(*it,dof,gridDim);
+
+ //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 GridType, int polOrd>
+double PlanarRodAssembler<GridType, polOrd>::
+computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
+{
+ const int maxlevel = grid_->maxLevel();
+ double energy = 0;
+
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(maxlevel);
+
+ if (sol.size()!=grid_->size(maxlevel, gridDim))
+ DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
+
+ ElementIterator it = grid_->template lbegin<0>(maxlevel);
+ ElementIterator endIt = grid_->template lend<0>(maxlevel);
+
+ // Loop over all elements
+ for (; it!=endIt; ++it) {
+
+ // Extract local solution on this element
+ Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
+
+ int numOfBaseFct = localFiniteElement.localBasis().size();
+
+ RigidBodyMotion<2> localSolution[numOfBaseFct];
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ // Get quadrature rule
+ const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), polOrd);
+
+ 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 = it->geometry().jacobianInverseTransposed(quadPos);
+ const double integrationElement = it->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/src/rodassembler.hh b/src/rodassembler.hh
index 8814c87c0af30d830732030db99b8f3be8fd208e..6f275383558f7c845f9172c9ed9c36961457cc67 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -103,6 +103,75 @@ class RodAssembler : public GeodesicFEAssembler<GridView, RigidBodyMotion<3> >
}; // end class
+
+/** \brief The FEM operator for a 2D extensible, shearable rod
+ */
+template <class GridType, int polOrd>
+class PlanarRodAssembler {
+
+ typedef typename GridType::template Codim<0>::Entity EntityType;
+ typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
+
+ //! Dimension of the grid. This needs to be one!
+ enum { gridDim = GridType::dimension };
+
+ enum { elementOrder = 1};
+
+ //! Each block is x, y, theta
+ enum { blocksize = 3 };
+
+ //!
+ typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+
+ const GridType* grid_;
+
+ /** \brief Material constants */
+ double B;
+ double A1;
+ double A3;
+
+public:
+
+ //! ???
+ PlanarRodAssembler(const GridType &grid) :
+ grid_(&grid)
+ {
+ B = 1;
+ A1 = 1;
+ A3 = 1;
+ }
+
+ ~PlanarRodAssembler() {}
+
+ 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<2> >& sol,
+ Dune::BCRSMatrix<MatrixBlock>& matrix);
+
+ void assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+
+ /** \brief Compute the energy of a deformation state */
+ double computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const;
+
+ void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
+
+protected:
+
+ /** \brief Compute the element tangent stiffness matrix */
+ template <class MatrixType>
+ void getLocalMatrix( EntityType &entity,
+ const std::vector<RigidBodyMotion<2> >& localSolution,
+ const int matSize, MatrixType& mat) const;
+
+}; // end class
+
#include "rodassembler.cc"
#endif