new project about nonlinear rods

[[Imported from SVN: r410]]

Makefile.am

+# $Id$
+
+# possible options
+#LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS)
+#AM_CPPFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS)
+
+noinst_PROGRAMS = staticrod
+
+staticrod_SOURCES = staticrod.cc
+
+# don't follow the full GNU-standard
+# we need automake 1.5
+AUTOMAKE_OPTIONS = foreign 1.5

README

+Preparing the Sources
+=========================
+
+Additional to the software mentioned in README you'll need the
+following programs installed on your system:
+
+  automake >= 1.5
+
+  autoconf >= 2.50
+
+  libtool
+
+Getting started
+---------------
+
+If these preliminaries are met, you should run the script
+
+  ./autogen.sh
+
+which calls the GNU autoconf/automake to create a ./configure-script
+and the Makefiles. Most probably you'll have to provide where to find
+the DUNE-files by
+
+   ./autogen.sh --with-dune=PATH
+
+where PATH is a directory with a dune/-subdirectory inside (this
+convention is needed to keep the #include-syntax consistent even when
+the headers are installed into /usr/include/dune later).
+
+
+Passing options to ./configure
+------------------------------
+
+autogen.sh also calls the newly created configure-script to
+conveniently pass on options about the used compiler. Thus you'll have
+to provide autogen.sh any options you want configure to get, e.g.
+
+  ./autogen.sh --with-dune=... --with-albert=... --without-x
+
+
+Choosing the compiler and the options
+-------------------------------------
+
+The selection of the compiler works as follows: if --gnu or --intel is
+passed to autogen it reads the content of gcc.opts or icc.opts to get
+the default compiler flags. With the option --optim you can switch the
+compiler-specific optimization parameters on.
+
+If you want to change the compiler options to your favourites you can
+either
+
+ - adapt the appropriate .opts-file and rerun autogen.sh. Please don't
+   commit this changed file to CVS if you're not sure if the options
+   work for everybody.
+
+ - copy an existing .opts-file to a new name, change the options and
+   use
+
+      ./autogen.sh --opts=my.opts
+
+
+More info
+---------
+
+See
+
+     ./autogen.sh --help
+   
+and (if it exists)
+ 
+     ./configure --help
+
+for further options.
+
+
+The full build-system is described in the dune/doc/Buildsystem (not in
+duneapps/doc!)
+
+$Id$
+

autogen.sh

+#!/bin/sh
+# $Id$
+
+#### barf on errors
+set -e
+
+# may be used to force a certain automake-version e.g. 1.7
+AMVERS=1.8
+
+# everybody who checks out the CVS wants the maintainer-mode to be enabled
+# (should be off for source distributions, this should happen automatically)
+#
+DEFAULTCONFOPT="--enable-maintainer-mode"
+
+# default values
+DEBUG=1
+OPTIM=0
+
+usage () {
+    echo "Usage: ./autogen.sh [options]"
+    echo "  -i, --intel        use intel compiler"
+    echo "  -g, --gnu          use gnu compiler (default)"
+    echo "  --opts=FILE        use compiler-options from FILE"
+    echo "  -d, --debug        switch debug-opts on"
+    echo "  -n, --nodebug      switch debug-opts off"
+    echo "  -o, --optim        switch optimization on"
+    echo "  --with-dune=PATH   directory with dune/ inside"
+    echo "  -h, --help         you already found this :)"
+    echo
+    echo "Parameters not in the list above are directly passed to configure. See"
+    echo
+    echo "    ./configure --help"
+    echo
+    echo "for a list of additional options"
+}
+
+# no compiler set yet
+COMPSET=0
+for OPT in $* ; do
+
+    set +e
+    # stolen from configure...
+    # when no option is set, this returns an error code
+    arg=`expr "x$OPT" : 'x[^=]*=\(.*\)'`
+    set -e
+
+    case "$OPT" in
+	-i|--intel)   . ./icc.opts ; COMPSET=1 ;;
+	-g|--gnu)     . ./gcc.opts ; COMPSET=1 ;;
+	--opts=*)
+	    if [ -r $arg ] ; then
+	      echo "reading options from $arg..."
+	      . ./$arg ;
+	      COMPSET=1;
+	    else
+	      echo "Cannot open compiler options file $arg!" ;
+	      exit 1;
+	    fi ;;
+	-d|--debug)   DEBUG=1 ;;
+	-n|--nodebug) DEBUG=0 ;;
+	-o|--optim)   OPTIM=1 ;;
+	-h|--help) usage ; exit 0 ;;
+	# special hack: use the with-dune-dir for aclocal-includes
+	--with-dune=*)
+	    eval DUNEDIR=$arg
+	    # add the option anyway
+	    CONFOPT="$CONFOPT $OPT" ;;
+	# pass unknown opts to ./configure
+	*) CONFOPT="$CONFOPT $OPT" ;;
+    esac
+done
+
+# set special m4-path if --with-dune is set
+if [ x$DUNEDIR != x ] ; then
+    # aclocal from automake 1.8 seems to need an absolute path for inclusion
+    FULLDIR=`cd $DUNEDIR && pwd`
+
+    # automagically use directory above if complete Dune-dir was supplied
+    if test `basename $FULLDIR` = "dune" ; then
+      FULLDIR=`cd $FULLDIR/.. && pwd`
+    fi
+
+    ACLOCALOPT="-I $FULLDIR/dune/m4/"
+fi
+
+# use the free compiler as default :-)
+if [ "$COMPSET" != "1" ] ; then
+    echo "No compiler set, using GNU compiler as default"
+    . ./gcc.opts
+fi
+
+# create flags
+COMPFLAGS="$FLAGS"
+
+# maybe add debug flag
+if [ "$DEBUG" = "1" ] ; then	
+    COMPFLAGS="$COMPFLAGS $DEBUGFLAGS"
+fi
+
+# maybe add optimization flag
+if [ "$OPTIM" = "1" ] ; then	
+    COMPFLAGS="$COMPFLAGS $OPTIMFLAGS"
+fi
+
+# check if automake-version was set
+if test "x$AMVERS" != x ; then
+  echo Warning: explicitly using automake version $AMVERS
+  # binaries are called automake-$AMVERS
+  AMVERS="-$AMVERS"
+fi
+
+#### create all autotools-files
+
+echo "--> libtoolize..."
+# force to write new versions of files, otherwise upgrading libtools
+# doesn't do anything...
+libtoolize --force
+
+echo "--> aclocal..."
+aclocal$AMVERS $ACLOCALOPT
+
+# sanity check to catch missing --with-dune
+if ! grep DUNE aclocal.m4 > /dev/null ; then
+    echo "aclocal.m4 doesn't contain any DUNE-macros, this would crash autoconf"
+    echo "or automake later. Maybe you should provide a --with-dune=PATH parameter"
+    exit 1
+fi
+
+echo "--> autoheader..."
+autoheader
+
+echo "--> automake..."
+automake$AMVERS --add-missing
+
+echo "--> autoconf..."
+autoconf
+
+#### start configure with special environment
+
+export CC="$COMP"
+export CXX="$CXXCOMP"
+export CPP="$COMP -E"
+
+export CFLAGS="$COMPFLAGS"
+export CXXFLAGS="$COMPFLAGS"
+
+./configure $DEFAULTCONFOPT $CONFOPT
+

configure.ac

+#                                               -*- Autoconf -*-
+# Process this file with autoconf to produce a configure script.
+
+AC_PREREQ(2.50)
+AC_INIT(rods, 1.0, sander@mi.fu-berlin.de)
+AM_INIT_AUTOMAKE(rods, 1.0, sander@mi.fu-berlin.de)
+AC_CONFIG_SRCDIR([staticrod.cc])
+AM_CONFIG_HEADER([config.h])
+
+# we need no more than the standard DUNE-stuff
+DUNE_CHECK_ALL
+
+# implicitly set the Dune-flags everywhere
+AC_SUBST(AM_CPPFLAGS, $DUNE_CPPFLAGS)
+AC_SUBST(AM_LDFLAGS, $DUNE_LDFLAGS)
+LIBS="$DUNE_LIBS"
+
+AC_CONFIG_FILES([Makefile])
+AC_OUTPUT
+
+# finally print the summary information
+DUNE_SUMMARY_ALL

gcc.opts

+# $Id$
+
+# options for gcc/g++
+# remember to run ./autogen.sh after changing these values!
+
+# name of compiler binaries
+COMP="gcc-4.0"
+CXXCOMP="g++-4.0"
+
+# flags set in any case
+FLAGS=""
+
+# additional flags for debugging
+DEBUGFLAGS="-g -DDUNE_ISTL_WITH_CHECKING"
+
+# additional flags for optimization
+OPTIMFLAGS="-O3"
+

icc.opts

+# $Id$
+
+# options for icc
+# remember to run ./autogen.sh after changing these values!
+
+# name of compiler binaries
+COMP="icc"
+CXXCOMP="icc"
+
+# flags set in any case
+FLAGS="-Wall"
+
+# additional flags for debugging
+DEBUGFLAGS="-O0 -g"
+
+# additional flags for optimization
+OPTIMFLAGS="-O3 -Ob2 -unroll"
+

src/rodassembler.cc

+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/common/matrix.hh>
+
+#include <dune/quadrature/quadraturerules.hh>
+
+
+
+template <class FunctionSpaceType, int polOrd>
+void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+getNeighborsPerVertex(MatrixIndexSet& nb) const
+{
+    const int gridDim = GridType::dimension;
+    
+    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 = it->template subIndex<gridDim>(i);
+                int jIdx = it->template subIndex<gridDim>(j);
+                
+                nb.add(iIdx, jIdx);
+                
+            }
+            
+        }
+        
+    }
+    
+}
+
+template <class FunctionSpaceType, int polOrd>
+void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
+               BCRSMatrix<MatrixBlock>& matrix)
+{
+
+    //int n = grid_->size(grid_->maxlevel(), dim);
+    
+    MatrixIndexSet neighborsPerVertex;
+    getNeighborsPerVertex(neighborsPerVertex);
+    
+    //neighborsPerVertex.exportIdx(*matrix_);
+    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 BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+        const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();  
+        
+        mat.resize(numOfBaseFct, numOfBaseFct);
+
+        // Extract local solution
+        BlockVector<FieldVector<double, blocksize> > localSolution(numOfBaseFct);
+        //BlockVector<FieldVector<double, dim> > localRhs(numOfBaseFct);
+        
+        for (int i=0; i<numOfBaseFct; i++)
+            localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+
+        // setup matrix 
+        getLocalMatrix( *it, localSolution, numOfBaseFct, mat);
+        
+        // Add element matrix to global stiffness matrix
+        for(int i=0; i<numOfBaseFct; i++) { 
+            
+            int row = functionSpace_.mapToGlobal( *it , i );
+
+            for (int j=0; j<numOfBaseFct; j++ ) {
+                
+                int col = functionSpace_.mapToGlobal( *it , j );    
+                matrix[row][col] += mat[i][j];
+                
+            }
+        }
+
+    }
+    
+}
+
+
+
+
+
+
+template <class FunctionSpaceType, int polOrd>
+template <class MatrixType>
+void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+getLocalMatrix( EntityType &entity, 
+                const BlockVector<FieldVector<double, blocksize> >& localSolution,
+                const int matSize, MatrixType& localMat) const
+{
+    
+    /* 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;
+    
+    typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+    const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
+    
+    // Get quadrature rule
+    const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity.geometry().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().jacobianInverse(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   */
+        /**********************************************/
+        Array<JacobianRange> shapeGrad(ndof);
+        
+        for (int dof=0; dof<ndof; dof++) {
+            
+            baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+            
+            // multiply with jacobian inverse 
+            JacobianRange tmp(0);
+            inv.umv(shapeGrad[dof][0], tmp[0]);
+            shapeGrad[dof][0] = tmp[0];
+
+        }
+        
+        
+        RangeType shapeFunction[matSize];
+        for(int i=0; i<matSize; i++) 
+            baseSet.eval(i,quadPos,shapeFunction[i]);
+        
+        // //////////////////////////////////
+        //   Interpolate
+        // //////////////////////////////////
+        
+        double x_s     = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
+        double y_s     = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
+        double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+
+        double theta   = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*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][0] * shapeGrad[j][0][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][0] * shapeGrad[j][0][0]
+                    * (-A1 * + A3) * sin(theta)* cos(theta);
+
+                // \partial J^2 / \partial x_i \partial theta_j
+                localMat[i][j][0][2] += weight * shapeGrad[i] * 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][0] * shapeGrad[j][0][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] * shapeGrad[j]
+                    * (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] * 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]
+                    * (-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]
+                    * (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] * shapeGrad[j];
+                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 FunctionSpaceType, int polOrd>
+void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
+                 BlockVector<FieldVector<double, blocksize> >& grad) const
+{
+    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
+        const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+        const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();  
+        
+        FieldVector<double, blocksize> localSolution[numOfBaseFct];
+        
+        for (int i=0; i<numOfBaseFct; i++)
+            localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+
+        // Get quadrature rule
+        const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().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().jacobianInverse(quadPos);
+            const double integrationElement = it->geometry().integrationElement(quadPos);
+        
+            double weight = quad[pt].weight() * integrationElement;
+            
+            // ///////////////////////////////////////
+            //   Compute deformation gradient
+            // ///////////////////////////////////////
+            Array<JacobianRange> shapeGrad(numOfBaseFct);
+            
+            for (int dof=0; dof<numOfBaseFct; dof++) {
+                
+                baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+                
+                // multiply with jacobian inverse 
+                JacobianRange tmp(0);
+                inv.umv(shapeGrad[dof][0], tmp[0]);
+                shapeGrad[dof][0] = tmp[0];
+                
+            }
+
+            // Get the value of the shape functions
+            RangeType shapeFunction[2];
+            for(int i=0; i<2; i++) 
+                baseSet.eval(i,quadPos,shapeFunction[i]);
+
+            // //////////////////////////////////
+            //   Interpolate
+            // //////////////////////////////////
+
+            double x_s     = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
+            double y_s     = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
+            double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+
+            double theta   = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*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 = functionSpace_.mapToGlobal(*it,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][0];
+
+                // \partial J / \partial y^i
+                grad[globalDof][1] += weight * (-partA1 * sin(theta) + partA3 * cos(theta)) * shapeGrad[dof][0][0];
+
+                // \partial J / \partial \theta^i
+                grad[globalDof][2] += weight * (B * theta_s * shapeGrad[dof][0][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 FunctionSpaceType, int polOrd>
+double Dune::RodAssembler<FunctionSpaceType, polOrd>::
+computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
+{
+    const int maxlevel = grid_->maxlevel();
+    double energy = 0;
+
+    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
+        const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+        const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();  
+        
+        FieldVector<double, blocksize> localSolution[numOfBaseFct];
+        
+        for (int i=0; i<numOfBaseFct; i++)
+            localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+
+        // Get quadrature rule
+        const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().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().jacobianInverse(quadPos);
+            const double integrationElement = it->geometry().integrationElement(quadPos);
+        
+            double weight = quad[pt].weight() * integrationElement;
+            
+            // ///////////////////////////////////////
+            //   Compute deformation gradient
+            // ///////////////////////////////////////
+            Array<JacobianRange> shapeGrad(numOfBaseFct);
+            
+            for (int dof=0; dof<numOfBaseFct; dof++) {
+                
+                baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+                //std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
+                
+                // multiply with jacobian inverse 
+                JacobianRange tmp(0);
+                inv.umv(shapeGrad[dof][0], tmp[0]);
+                shapeGrad[dof][0] = tmp[0];
+                //std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
+
+                
+                
+            }
+
+            // Get the value of the shape functions
+            RangeType shapeFunction[2];
+            for(int i=0; i<2; i++) 
+                baseSet.eval(i,quadPos,shapeFunction[i]);
+
+            // //////////////////////////////////
+            //   Interpolate
+            // //////////////////////////////////
+
+            double x_s     = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
+            double y_s     = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
+            double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+
+            double theta   = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*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);
+
+//             printf("energy %g\n", energy);
+        }
+
+    }
+
+    return energy;
+
+}

src/rodassembler.hh

+#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/common/matrix.hh>
+
+namespace Dune 
+{
+
+    /** \brief The FEM operator for an extensible, shearable rod
+     */
+    template <class FunctionSpaceType, int polOrd>
+    class RodAssembler {
+        
+        //! The grid
+        typedef typename FunctionSpaceType::GridType GridType;
+        
+        typedef typename GridType::template Codim<0>::Entity EntityType;
+        typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
+        typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
+    
+
+        //! Dimension of the grid.  This needs to be one!
+        enum { gridDim = GridType::dimension };
+
+        //! Each block is x, y, theta
+        enum { blocksize = 3 };
+        
+        //!
+        typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+        
+        //! ???
+        typedef typename FunctionSpaceType::JacobianRange JacobianRange;
+        
+        //! ???
+        typedef typename FunctionSpaceType::RangeField RangeFieldType;
+        typedef typename FunctionSpaceType::Range       RangeType;
+        
+    public:
+        
+        /** \todo Does actually belong into the base class */
+        const GridType* grid_; 
+        
+        /** \todo Does actually belong into the base class */
+        const FunctionSpaceType& functionSpace_;
+
+        /** \brief Material constants */
+        double B;
+        double A1;
+        double A3;
+
+        //! ???
+        RodAssembler(const FunctionSpaceType &f) : 
+            functionSpace_(f)
+        { 
+            grid_ = &f.getGrid();
+            B = 1;
+            A1 = 1;
+            A3 = 1;
+        }
+
+        ~RodAssembler() {}
+
+        void getNeighborsPerVertex(MatrixIndexSet& nb) const;
+        
+        /** \brief Assemble the tangent stiffness matrix and the right hand side
+         */
+        void assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
+                            BCRSMatrix<MatrixBlock>& matrix);
+        
+        /** \brief Compute the element tangent stiffness matrix  */
+        template <class MatrixType>
+        void getLocalMatrix( EntityType &entity, 
+                             const BlockVector<FieldVector<double, blocksize> >& localSolution, 
+                             const int matSize, MatrixType& mat) const;
+
+        void assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
+                              BlockVector<FieldVector<double, blocksize> >& grad) const;
+
+        /** \brief Compute the energy of a deformation state */
+        double computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const;
+        
+        
+        
+        
+        
+        
+    }; // end class
+    
+} // end namespace 
+
+#include "rodassembler.cc"
+
+#endif
+

staticrod.cc

+#include <config.h>
+
+#include <dune/grid/onedgrid.hh>
+
+#include <dune/fem/lagrangebase.hh>
+
+#include <dune/istl/io.hh>
+
+//#include "../common/boundarytreatment.hh"
+#include "../common/boundarypatch.hh"
+#include <dune/common/bitfield.hh>
+//#include "../common/readbitfield.hh"
+
+#include "src/rodassembler.hh"
+//#include "../common/linearipopt.hh"
+
+#include "../common/projectedblockgsstep.hh"
+#include <dune/solver/iterativesolver.hh>
+
+#include "../common/geomestimator.hh"
+#include "../common/energynorm.hh"
+#include <dune/common/configparser.hh>
+
+// Choose a solver
+//#define IPOPT
+#define GAUSS_SEIDEL
+//#define MULTIGRID
+
+//#define IPOPT_BASE
+
+// Number of degrees of freedom: 
+// 3 (x, y, theta) for a planar rod
+const int blocksize = 3;
+
+using namespace Dune;
+using std::string;
+
+int main (int argc, char *argv[]) try
+{
+    // Some types that I need
+    typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > MatrixType;
+    typedef BlockVector<FieldVector<double, blocksize> >     VectorType;
+
+    // parse data file
+    ConfigParser parameterSet;
+    parameterSet.parseFile("staticrod.parset");
+
+    // read solver settings
+    const int minLevel         = parameterSet.get("minLevel", int(0));
+    const int maxLevel         = parameterSet.get("maxLevel", int(0));
+    double loadIncrement       = parameterSet.get("loadIncrement", double(0));
+    const int maxNewtonSteps   = parameterSet.get("maxNewtonSteps", int(0));
+    const int numIt            = parameterSet.get("numIt", int(0));
+    const int nu1              = parameterSet.get("nu1", int(0));
+    const int nu2              = parameterSet.get("nu2", int(0));
+    const int mu               = parameterSet.get("mu", int(0));
+    const int baseIt           = parameterSet.get("baseIt", int(0));
+    const double tolerance     = parameterSet.get("tolerance", double(0));
+    const double baseTolerance = parameterSet.get("baseTolerance", double(0));
+    
+    // Problem settings
+    const int numRodElements = parameterSet.get("numRodElements", int(0));
+
+    // ///////////////////////////////////////
+    //    Create the two grids
+    // ///////////////////////////////////////
+    typedef OneDGrid<1,1> RodGridType;
+    RodGridType rod(numRodElements, 0, 1);
+
+    
+
+    Array<BitField> dirichletNodes;
+    dirichletNodes.resize(maxLevel+1);
+    dirichletNodes[0].resize( blocksize * (numRodElements+1) );
+
+    dirichletNodes[0].unsetAll();
+    dirichletNodes[0][0] = dirichletNodes[0][1] = dirichletNodes[0][2] = true;
+    dirichletNodes[0][blocksize*numRodElements+0] = true;
+    dirichletNodes[0][blocksize*numRodElements+1] = true;
+    dirichletNodes[0][blocksize*numRodElements+2] = true;
+
+    // refine uniformly until minlevel
+    for (int i=0; i<minLevel; i++)
+        rod.globalRefine(1);
+
+    int maxlevel = rod.maxlevel();
+
+    // //////////////////////////////////////////////////////////
+    //   Create obstacles
+    // //////////////////////////////////////////////////////////
+
+    Array<BitField> hasObstacle;
+    hasObstacle.resize(maxLevel+1);
+    hasObstacle[0].resize(numRodElements+1);
+    hasObstacle[0].unsetAll();
+
+
+
+    // //////////////////////////////////////////////////////////
+    //    Create discrete function spaces
+    // //////////////////////////////////////////////////////////
+
+    typedef FunctionSpace < double , double, 1, 1 > RodFuncSpace;
+    typedef DefaultGridIndexSet<RodGridType,LevelIndex> RodIndexSet;
+    typedef LagrangeDiscreteFunctionSpace < RodFuncSpace, RodGridType,RodIndexSet,  1> RodFuncSpaceType;
+
+    Array<RodIndexSet*> rodIndexSet(maxlevel+1);
+    Array<const RodFuncSpaceType*> rodFuncSpace(maxlevel+1);
+
+    for (int i=0; i<maxlevel+1; i++) {
+        rodIndexSet[i]  = new RodIndexSet(rod, i);
+        rodFuncSpace[i] = new RodFuncSpaceType(rod, *rodIndexSet[i], i);
+    }
+
+
+    // ////////////////////////////////////////////////////////////
+    //    Create solution and rhs vectors
+    // ////////////////////////////////////////////////////////////
+
+    VectorType rhs;
+    VectorType x;
+    VectorType corr;
+
+    MatrixType hessianMatrix;
+
+    rhs.resize(rodFuncSpace[maxlevel]->size());
+    x.resize(rodFuncSpace[maxlevel]->size());
+    corr.resize(rodFuncSpace[maxlevel]->size());
+    
+    // Initial solution
+    x = 0;
+
+    for (int i=0; i<numRodElements+1; i++) {
+        x[i][0] = i/((double)numRodElements);
+        x[i][1] = 0;
+        x[i][2] = M_PI/2;
+    }
+
+    x[0][1] = x[numRodElements][1] = 1;
+
+    RodAssembler<RodFuncSpaceType,2> test(*rodFuncSpace[0]);
+    test.assembleGradient(x, rhs);
+    //std::cout << "Solution: " << std::endl << x << std::endl;
+    //std::cout << "Gradient: " << std::endl << rhs << std::endl;
+    std::cout << "Energy: " << test.computeEnergy(x) << std::endl;
+
+    MatrixIndexSet indices(numRodElements+1, numRodElements+1);
+    test.getNeighborsPerVertex(indices);
+    indices.exportIdx(hessianMatrix);
+    test.assembleMatrix(x,hessianMatrix);
+
+    //printmatrix(std::cout, hessianMatrix, "hessianMatrix", "--");
+    //exit(0);
+
+    // Create a solver
+#if defined IPOPT
+
+    typedef LinearIPOptSolver<VectorType> SolverType;
+    
+    SolverType solver;
+    solver.dirichletNodes_ = &totalDirichletNodes[maxlevel];
+    solver.hasObstacle_    = &contactAssembler.hasObstacle_[maxlevel];
+    solver.obstacles_      = &contactAssembler.obstacles_[maxlevel];
+    solver.verbosity_      = Solver::FULL;
+
+#elif defined GAUSS_SEIDEL
+
+    typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType;
+    SmootherType projectedBlockGSStep(hessianMatrix, corr, rhs);
+    projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
+    projectedBlockGSStep.hasObstacle_    = &hasObstacle[maxlevel];
+    projectedBlockGSStep.obstacles_      = NULL;//&contactAssembler.obstacles_[maxlevel];
+
+    EnergyNorm<MatrixType, VectorType> energyNorm(projectedBlockGSStep);
+
+    IterativeSolver<MatrixType, VectorType> solver;
+    solver.iterationStep = &projectedBlockGSStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::QUIET;
+    solver.errorNorm_ = &energyNorm;
+    solver.tolerance_ = tolerance;
+    
+#elif defined MULTIGRID
+
+    // First create a base solver
+#ifdef IPOPT_BASE
+
+    LinearIPOptSolver<BlockVector<FieldVector<double,dim> >  > baseSolver;
+    baseSolver.verbosity_ = Solver::FULL;
+
+#else // Gauss-Seidel is the base solver
+
+    ProjectedBlockGSStep<MatrixType, BlockVector<FieldVector<double,dim> >  > baseSolverStep;
+
+    EnergyNorm<MatrixType, BlockVector<FieldVector<double,dim> >  > baseEnergyNorm(baseSolverStep);
+
+    IterativeSolver<MatrixType, BlockVector<FieldVector<double,dim> >  > baseSolver;
+    baseSolver.iterationStep = &baseSolverStep;
+    baseSolver.numIt = baseIt;
+    baseSolver.verbosity_ = Solver::QUIET;
+    baseSolver.errorNorm_ = &baseEnergyNorm;
+    baseSolver.tolerance_ = baseTolerance;
+#endif
+
+    // Make pre and postsmoothers
+    ProjectedBlockGSStep<MatrixType, BlockVector<FieldVector<double,dim> >  > presmoother;
+    ProjectedBlockGSStep<MatrixType, BlockVector<FieldVector<double,dim> >  > postsmoother;
+    
+
+    ContactMMGStep<MatrixType, BlockVector<FieldVector<double,dim> > , FuncSpaceType > contactMMGStep(maxlevel+1);
+
+    contactMMGStep.setMGType(1, nu1, nu2);
+    contactMMGStep.dirichletNodes_    = &totalDirichletNodes;
+    contactMMGStep.basesolver_        = &baseSolver;
+    contactMMGStep.presmoother_       = &presmoother;
+    contactMMGStep.postsmoother_      = &postsmoother;    
+    contactMMGStep.hasObstacle_       = &hasObstacle;
+    contactMMGStep.obstacles_         = &contactAssembler.obstacles_;
+
+    // Create the transfer operators
+    contactMMGStep.mgTransfer_.resize(maxlevel);
+    for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++)
+        contactMMGStep.mgTransfer_[i] = NULL;
+
+    EnergyNorm<MatrixType, VectorType> energyNorm(contactMMGStep);
+
+    IterativeSolver<MatrixType, BlockVector<FieldVector<double,dim> > > solver;
+    solver.iterationStep = &contactMMGStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::FULL;
+    solver.errorNorm_ = &energyNorm;
+    solver.tolerance_ = tolerance;
+
+#else
+    #warning You have to specify a solver!
+#endif
+
+    // ///////////////////////////////////////////////////
+    //   Do a homotopy of the Dirichlet boundary data
+    // ///////////////////////////////////////////////////
+    double loadFactor = 0;
+
+    do {
+
+        RodAssembler<RodFuncSpaceType, 1> rodAssembler(*rodFuncSpace[maxlevel]);
+
+        loadFactor += loadIncrement;
+
+        std::cout << "####################################################" << std::endl;
+        std::cout << "New load factor: " << loadFactor 
+                  << "    new load increment: " << loadIncrement << std::endl;
+        std::cout << "####################################################" << std::endl;
+
+        // /////////////////////////////////////////////////////
+        //   Newton Solver
+        // /////////////////////////////////////////////////////
+
+        for (int j=0; j<maxNewtonSteps; j++) {
+
+            rhs = 0;
+            corr = 0;
+
+            rodAssembler.assembleGradient(x, rhs);
+            rodAssembler.assembleMatrix(x, hessianMatrix);
+
+            rhs *= -1;
+
+            std::cout << "rhs: " << std::endl << rhs << std::endl;
+
+#ifndef IPOPT
+            solver.iterationStep->setProblem(hessianMatrix, corr, rhs);
+#else
+            solver.setProblem(hessianMatrix, corr, rhs);
+#endif
+
+            solver.preprocess();
+#ifdef MULTIGRID
+
+            contactMMGStep.preprocess();
+#endif
+
+            // /////////////////////////////
+            //    Solve !
+            // /////////////////////////////
+             solver.solve();
+
+#ifdef MULTIGRID
+             totalCorr = contactMMGStep.getSol();
+#endif
+
+             std::cout << "Correction: \n" << corr << std::endl;
+
+             // line search
+             printf("------  Line Search ---------\n");
+             int lSSteps = 10;
+             double smallestEnergy = std::numeric_limits<double>::max();
+             double smallestFactor = 1;
+             for (int k=0; k<lSSteps; k++) {
+
+                 double factor = double(k)/(lSSteps-1);
+                 VectorType sCorr = corr;
+                 sCorr *= factor;
+                 sCorr += x;
+
+                 double energy = rodAssembler.computeEnergy(sCorr);
+
+                 if (energy < smallestEnergy) {
+                     smallestEnergy = energy;
+                     smallestFactor = factor;
+                 }
+                 //printf("factor: %g,  energy: %g\n", factor, energy);
+             }
+
+             std::cout << "Damping factor: " << smallestFactor << std::endl;
+
+             //  Add correction to the current solution
+             x.axpy(smallestFactor, corr);
+
+             // Output result
+             std::cout << "Solution:" << std::endl << x << std::endl;
+
+             printf("infinity norm of the correction: %g\n", corr[0].infinity_norm());
+             if (corr.infinity_norm() < 1e-8)
+                 break;
+
+        }
+        
+        
+    } while (loadFactor < 1);
+
+
+
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }

staticrod.parset

+# Number of grid levels
+minLevel = 0
+
+maxLevel = 0
+
+# Initial load increment
+loadIncrement = 1
+
+# Max number of steps of the Newton solver
+maxNewtonSteps = 10
+
+# Number of multigrid iterations per time step
+numIt = 100
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the solver
+tolerance = -1
+
+# Tolerance of the base grid solver
+baseTolerance = -1
+
+############################
+#   Problem specifications
+############################
+
+# Number of elements of the rod grid
+numRodElements = 40