a package for elasticity problems of all sorts

[[Imported from SVN: r9230]]

Makefile.am

+# $Id$
+
+IPOPT_DIR = /home/haile/sander/COIN/Ipopt
+PARAM_DIR = /home/haile/sander/libpsurface
+
+# possible options
+LDADD = $(AMIRAMESH_LDFLAGS) $(AMIRAMESH_LIBS) -L$(IPOPT_DIR)/lib -lipopt -llapack -lblas -lg2c
+AM_CPPFLAGS += $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) -I$(IPOPT_DIR)/IPOPT/include
+
+# Include and linker paths to the contact handling library
+LDADD       += -L$(PARAM_DIR)/lib -lpsurface
+AM_CPPFLAGS += -I$(PARAM_DIR)/include
+
+LDADD       += $(MPI_LDFLAGS)
+AM_CPPFLAGS += $(MPI_CPPFLAGS)
+
+noinst_PROGRAMS = linelast nonlinelast threegrids
+
+linelast_SOURCES        = linelast.cc
+nonlinelast_SOURCES     = nonlinelast.cc
+
+threegrids_SOURCES      = threegrids.cc
+threegrids_CXXFLAGS     = $(ALBERTA_CPPFLAGS) $(ALUGRID_CPPFLAGS) $(AM_CPPFLAGS) -DNDEBUG
+threegrids_LDADD        = $(ALBERTA_LDFLAGS) $(ALBERTA_LIBS) $(ALUGRID_LDFLAGS) $(ALUGRID_LIBS)
+
+# don't follow the full GNU-standard
+# we need automake 1.5
+AUTOMAKE_OPTIONS = foreign 1.5

README

+Getting started
+===============
+
+You have to create the configure-script on your own!
+
+First, you'll need the followings programs installed:
+
+  automake >= 1.5
+
+  autoconf >= 2.50
+
+  libtool
+
+Then run
+
+  ./autogen.sh DUNEDIR
+
+where DUNEDIR is the (absolute or relative) path of the dune/-directory.
+
+The directory is needed because autogen.sh needs access to the tests
+stored in DUNEDIR/m4
+
+Now call
+
+  ./configure --with-dune=DIR
+
+where DIR is the directory _above_ dune/. If you want to include third
+party packages check
+
+  ./configure --help
+
+for options on how to include Albert, Grape, UG, ...
+
+If configure checked your system without problems you can use
+
+  make
+
+to build all examples.
+

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(elasticity, 1.0, sander)
+AM_INIT_AUTOMAKE(elasticity, 1.0, sander)
+AC_CONFIG_SRCDIR([linelast.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

linelast.cc

+#include <config.h>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/io/file/amirameshwriter.hh>
+#include <dune/io/file/amirameshreader.hh>
+
+#include "src/linelast.hh"
+
+#include <dune/istl/io.hh>
+
+#include "../common/boundarytreatment.hh"
+#include <dune/common/bitfield.hh>
+#include "../common/readbitfield.hh"
+
+#include "../common/linearipopt.hh"
+
+#include "../common/blockgsstep.hh"
+#include "../common/multigridstep.hh"
+#include <dune/solver/iterativesolver.hh>
+
+#include "../common/energynorm.hh"
+
+// Choose a solver
+//#define IPOPT
+//#define GAUSS_SEIDEL
+#define MULTIGRID
+
+#define IPOPT_BASE
+
+// The grid dimension
+const int dim = 3;
+
+// Number of grid levels
+const int numLevels = 1;
+
+// Number of multigrid iterations per time step
+const int numIt = 100;
+
+// Number of presmoothing steps
+const int nu1 = 3;
+
+// Number of postsmoothing steps
+const int nu2 = 3;
+
+// Number of coarse grid corrections
+const int mu = 1;
+
+// Number of base solver iterations
+const int baseIt = 100;
+
+// Tolerance of the solver
+const double tolerance = -1;//1e-5;
+
+// Tolerance of the base grid solver
+const double baseTolerance = -1;
+
+using namespace Dune;
+
+// ////////////////////////////////////
+//   Problem specifications
+// ////////////////////////////////////
+
+// std::string gridFile = "/home/haile/sander/data/elasticity/doublehexa_z/doublehexa.grid";
+// std::string dnFile   = "/home/haile/sander/data/elasticity/doublehexa_z/doublehexa.dn";
+// std::string dvFile   = "/home/haile/sander/data/elasticity/doublehexa_z/doublehexa.dv";
+
+// std::string gridFile = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.grid";
+// std::string dnFile   = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.dn";
+// std::string dvFile   = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.dv";
+
+// std::string gridFile = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.grid";
+// std::string dnFile   = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.dn";
+// std::string dvFile   = "/home/haile/sander/data/elasticity/doubletetra/doubletetra.dv";
+
+std::string gridFile = "/home/haile/sander/data/contact/vishum_stump/femur.grid";
+std::string dnFile   = "/home/haile/sander/data/contact/vishum_stump/femur.few.dn";
+std::string dvFile   = "/home/haile/sander/data/contact/vishum_stump/femur.few.dv";
+
+
+int main (int argc, char *argv[]) try
+{
+    // Some types that I need
+    typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
+    typedef BlockVector<FieldVector<double, dim> >     VectorType;
+
+    // /////////////////////////////
+    //   Generate the grid
+    // /////////////////////////////
+    typedef UGGrid<dim,dim> GridType;
+    GridType grid;
+
+    Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
+    
+    // Read Dirichlet boundary values
+    Array<BitField>  dirichletNodes(numLevels);
+    readBitField(dirichletNodes[0], grid.size(0, dim), dim, dnFile);
+
+    Array<VectorType> dirichletValues(numLevels);
+    dirichletValues[0].resize(grid.size(0, dim));
+    AmiraMeshReader<int>::readFunction(dirichletValues[0], dvFile);
+
+    for (int i=0; i<numLevels-1; i++)
+        grid.globalRefine(1);
+
+    int maxlevel = grid.maxLevel();
+
+  // Determine Dirichlet dofs
+    //prolongDirichletInformation(grid, dirichletNodes, dim);
+#warning Dirichlet handling not done yet!
+    sampleOnBitField(grid, dirichletValues, dirichletNodes);
+
+    VectorType rhs(grid.size(grid.maxLevel(), dim));
+    VectorType x(grid.size(grid.maxLevel(), dim));
+    
+    // Initial solution
+    x = 0;
+
+    rhs = 0;
+#if 1
+    for (int i=0; i<rhs.size(); i++)
+        for (int j=0; j<dim; j++) {
+            if (dirichletNodes[maxlevel][i*dim+j])
+                x[i][j] = dirichletValues[maxlevel][i][j];
+
+        }
+    //std::cout << "rhs\n " << rhs << std::endl;
+#endif
+
+    // Assemble elasticity problem
+    typedef Dune::LinElasticityOp<GridType, 3 > LinearElasticityType;
+    LinearElasticityType linElastOp0(grid);
+
+    linElastOp0.assembleMatrix();
+    const OperatorType* bilinearForm = linElastOp0.getMatrix();
+
+    // Create a solver
+#if defined IPOPT
+
+    typedef LinearIPOptSolver SolverType;
+    
+    SolverType solver(*bilinearForm, x, rhs);
+    solver.dirichletNodes_ = &dirichletNodes[maxlevel];
+
+#elif defined GAUSS_SEIDEL
+
+    L2Norm<VectorType> l2Norm;
+
+    typedef BlockGSStep<OperatorType, VectorType> SmootherType;
+    SmootherType blockGSStep(*bilinearForm, x, rhs);
+    blockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
+
+    IterativeSolver<OperatorType, VectorType> solver;
+    solver.iterationStep = &blockGSStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::FULL;
+    solver.errorNorm_ = &l2Norm;
+    solver.tolerance_ = tolerance;
+    
+#elif defined MULTIGRID
+
+    // First create a base solver
+#ifdef IPOPT_BASE
+
+    LinearIPOptSolver<VectorType> baseSolver;
+
+#else // Gauss-Seidel is the base solver
+
+    BlockGSStep<OperatorType, VectorType> baseSolverStep;
+
+    IterativeSolver<OperatorType, VectorType> baseSolver;
+    baseSolver.iterationStep = &baseSolverStep;
+    baseSolver.numIt = baseIt;
+    baseSolver.verbosity_ = Solver::QUIET;
+    baseSolver.errorNorm_ = &l2Norm;
+    baseSolver.tolerance_ = baseTolerance;
+#endif
+
+    // Make pre and postsmoothers
+    BlockGSStep<OperatorType, VectorType> presmoother;
+    BlockGSStep<OperatorType, VectorType> postsmoother;
+    
+
+    MultigridStep<OperatorType, VectorType, GridType> multigridStep(*bilinearForm, x, rhs, numLevels);
+
+    multigridStep.setMGType(1, nu1, nu2);
+    multigridStep.dirichletNodes_ = &dirichletNodes;
+    multigridStep.basesolver_     = &baseSolver;
+    multigridStep.presmoother_    = &presmoother;
+    multigridStep.postsmoother_   = &postsmoother;    
+    multigridStep.grid_           = &grid;
+
+    EnergyNorm<OperatorType, VectorType> energyNorm(multigridStep);
+
+    IterativeSolver<OperatorType, VectorType> solver;
+    solver.iterationStep = &multigridStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::FULL;
+    solver.errorNorm_ = &energyNorm;
+    solver.tolerance_ = tolerance;
+
+#else
+    #warning You have to specify a solver!
+#endif
+
+
+
+
+
+    
+  solver.preprocess();
+#ifdef MULTIGRID
+  multigridStep.preprocess();
+#endif
+
+  // Compute solution
+  solver.solve();
+
+#ifdef MULTIGRID
+  x = multigridStep.getSol();
+#endif
+
+  //std::cout << x << std::endl;
+
+  // Output result
+  AmiraMeshWriter<GridType>::writeGrid(grid, "resultGrid");
+  AmiraMeshWriter<GridType>::writeBlockVector(grid, x, "resultGridFunc");
+
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }
+
+
+

nonlinelast.cc

+#include <config.h>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/io/file/amirameshwriter.hh>
+#include <dune/io/file/amirameshreader.hh>
+
+#include "src/ogdenassembler.hh"
+
+#include <dune/istl/io.hh>
+
+#include "../common/prolongboundarypatch.hh"
+#include <dune/common/bitfield.hh>
+#include "../common/readbitfield.hh"
+
+#include "../common/linearipopt.hh"
+
+#include "../contact/src/contactmmgstep.hh"
+#include <dune/solver/iterativesolver.hh>
+
+#include "../common/energynorm.hh"
+
+// Choose a solver
+#define IPOPT
+//#define GAUSS_SEIDEL
+//#define MULTIGRID
+
+//#define IPOPT_BASE
+
+// The grid dimension
+const int dim = 3;
+
+// Number of grid levels
+const int numLevels = 1;
+
+// Number of steps for the boundary data homotopy
+const int numHomotopySteps = 1;
+
+// Number of Newton steps
+const int numNewtonSteps = 1;
+
+// Number of multigrid iterations per time step
+const int numIt = 1000;
+
+// Number of presmoothing steps
+const int nu1 = 3;
+
+// Number of postsmoothing steps
+const int nu2 = 3;
+
+// Number of coarse grid corrections
+const int mu = 1;
+
+// Number of base solver iterations
+const int baseIt = 100;
+
+// Tolerance of the solver
+const double tolerance = 1e-5;
+
+// Tolerance of the base grid solver
+const double baseTolerance = -1;
+
+using namespace Dune;
+
+// ////////////////////////////////////
+//   Problem specifications
+// ////////////////////////////////////
+
+// std::string gridFile = "data_1b/unithexacube.am";
+// std::string dnFile   = "data_1b/unithexacube.x.dn";
+// std::string dvFile   = "data_1b/unithexacube.x.vectors";
+// double obs(const FieldVector<double, dim>& x) {return -(x[0]-1.005);}
+// double obsDistance = 0.01;
+
+std::string gridFile = "/home/haile/sander/data/elasticity/cubus27/Cubus27.grid";
+std::string dnFile   = "/home/haile/sander/data/elasticity/cubus27/cubus27.dn";
+std::string dvFile   = "/home/haile/sander/data/elasticity/cubus27/cubus27.dv";
+
+// std::string gridFile = "data_1b/sphere.coarse.grid";
+// std::string dnFile   = "data_1b/sphere.coarse.dn";
+// std::string dvFile   = "data_1b/sphere.coarse.vectors";
+// double obs(const FieldVector<double, dim>& x) {return -(x[0]-1);}
+// double obsDistance = 0.2;
+
+int main (int argc, char *argv[]) try
+{
+    // Some types that I need
+    typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
+    typedef BlockVector<FieldVector<double, dim> >     VectorType;
+
+    // /////////////////////////////
+    //   Generate the grid
+    // /////////////////////////////
+    typedef UGGrid<dim,dim> GridType;
+    GridType grid;
+
+    Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
+    
+    // Read Dirichlet boundary values
+    Array<BoundaryPatch<GridType> >  dirichletBoundary(numLevels);
+    readBoundaryPatch(dirichletBoundary[0], dnFile);
+
+    Array<VectorType> dirichletValues(numLevels);
+    dirichletValues[0].resize(grid.size(0, dim));
+    AmiraMeshReader<int>::readFunction(dirichletValues[0], dvFile);
+
+    for (int i=0; i<numLevels-1; i++)
+        grid.globalRefine(1);
+
+    prolongBoundaryPatch(dirichletBoundary);
+
+    Array<BitField>  dirichletNodes(numLevels);
+    for (int i=0; i<numLevels; i++) {
+        dirichletNodes[i].resize(grid.size(i, dim)*dim);
+        for (int j=0; j<grid.size(i, dim); j++) {
+
+            for (int k=0; k<dim; k++)
+                dirichletNodes[i][j*dim+k] = dirichletBoundary[i].containsVertex(j);
+        }
+
+    }
+
+    int maxlevel = grid.maxLevel();
+
+  // Determine Dirichlet dofs
+    VectorType rhs(grid.size(grid.maxLevel(), dim));
+    VectorType x(grid.size(grid.maxLevel(), dim));
+    VectorType corr(grid.size(grid.maxLevel(), dim));
+    
+    // Initial solution
+    x = 0;
+    corr = 0;
+    rhs = 0;
+
+    // Create a solver
+#if defined IPOPT
+
+    LinearIPOptSolver<VectorType> solver;
+
+    solver.dirichletNodes_ = &dirichletNodes[maxlevel];
+
+#elif defined GAUSS_SEIDEL
+
+    L2Norm<ContactVector<dim> > l2Norm;
+
+    typedef ProjectedBlockGSStep<OperatorType, VectorType> SmootherType;
+    SmootherType projectedBlockGSStep(*bilinearForm, x, rhs);
+    projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
+
+    IterativeSolver<OperatorType, ContactVector<dim> > solver;
+    solver.iterationStep = &projectedBlockGSStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::FULL;
+    solver.errorNorm_ = &l2Norm;
+    solver.tolerance_ = tolerance;
+    
+#elif defined MULTIGRID
+
+    L2Norm<ContactVector<dim> > l2Norm;
+
+    // First create a base solver
+#ifdef IPOPT_BASE
+
+    LinearIPOptSolver baseSolver;
+
+#else // Gauss-Seidel is the base solver
+
+    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > baseSolverStep;
+
+    IterativeSolver<OperatorType, ContactVector<dim> > baseSolver;
+    baseSolver.iterationStep = &baseSolverStep;
+    baseSolver.numIt = baseIt;
+    baseSolver.verbosity_ = Solver::QUIET;
+    baseSolver.errorNorm_ = &l2Norm;
+    baseSolver.tolerance_ = baseTolerance;
+#endif
+
+    // Make pre and postsmoothers
+    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > presmoother;
+    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > postsmoother;
+    
+
+    ContactMMGStep<OperatorType, ContactVector<dim>, FuncSpaceType > contactMMGStep(*bilinearForm, x, rhs, numLevels);
+
+    contactMMGStep.setMGType(1, nu1, nu2);
+    contactMMGStep.dirichletNodes_ = &dirichletNodes;
+    contactMMGStep.basesolver_     = &baseSolver;
+    contactMMGStep.presmoother_    = &presmoother;
+    contactMMGStep.postsmoother_   = &postsmoother;    
+    contactMMGStep.funcSpace_      = funcSpace;
+
+    IterativeSolver<OperatorType, ContactVector<dim> > solver;
+    solver.iterationStep = &contactMMGStep;
+    solver.numIt = numIt;
+    solver.verbosity_ = Solver::FULL;
+    solver.errorNorm_ = &l2Norm;
+    solver.tolerance_ = tolerance;
+
+#else
+    #warning You have to specify a solver!
+#endif
+
+    typedef OgdenAssembler<GridType, 1> Assembler;
+    Assembler ogdenAssembler(grid);
+    
+    for (int i=0; i<numHomotopySteps; i++) {
+
+        // scale Dirichlet values
+        VectorType scaledDirichlet = dirichletValues[maxlevel];
+        scaledDirichlet *= (double(i)+1)/numHomotopySteps;
+
+        // Set new Dirichlet values in solution
+        for (int j=0; j<x.size(); j++)
+            for (int k=0; k<dim; k++)
+                if (dirichletNodes[maxlevel][j*dim+k])
+                    x[j][k] = scaledDirichlet[j][k];
+
+
+
+        for (int j=0; j<numNewtonSteps; j++) {
+
+          // Assemble the Jacobi matrix at the current solution
+            rhs = 0;
+          ogdenAssembler.assembleMatrix(x, rhs);
+
+          corr = 0;
+          
+          // Set new Dirichlet values in the right hand side
+          for (int k=0; k<rhs.size(); k++)
+              for (int l=0; l<dim; l++)
+                  if (dirichletNodes[maxlevel][k*dim+l])
+                      rhs[k][l] = 0;
+
+          solver.setProblem(*ogdenAssembler.matrix_, corr, rhs);
+
+          solver.preprocess();
+#ifdef MULTIGRID
+          contactMMGStep.preprocess();
+#endif
+
+          // Solve correction problem
+          solver.solve();
+
+          // Add damped correction to the current solution
+          x += corr;
+
+      }
+
+  }
+
+#ifdef MULTIGRID
+  x = contactMMGStep.getSol();
+#endif
+
+  //contactAssembler.postprocess(x);
+
+  //std::cout << x << std::endl;
+
+  // Output result
+  AmiraMeshWriter<GridType>::writeGrid(grid, "resultGrid");
+  AmiraMeshWriter<GridType>::writeBlockVector(grid, x, "resultGridFunc");
+
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }
+
+
+

threegrids.cc

+#include "config.h"
+
+#include <dune/grid/albertagrid.hh>
+#include <dune/grid/alu3dgrid.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/solvers.hh>
+
+#include <dune/common/bitfield.hh>
+
+#include <dune/io/file/amirameshreader.hh>
+#include <dune/io/file/amirameshwriter.hh>
+
+#include "../common/readbitfield.hh"
+#include "src/linelast.hh"
+
+using namespace Dune;
+
+const int dim = 3;
+
+typedef BlockVector<FieldVector<double, dim> > VectorType;
+typedef BCRSMatrix<FieldMatrix<double,dim,dim> > MatrixType;
+
+std::string path = "/home/haile/sander/data/kurbelwelle/";
+std::string dnFile = "kurbelwelle.dn";
+std::string dvFile = "kurbelwelle.dv";
+
+template <class GridType>
+void compute(GridType& grid)
+{
+
+    std::cout << "Computing on a " << grid.type() << std::endl;
+
+    // temporary
+    int numLevels = 1;
+
+    // Read Dirichlet boundary values
+    Array<BitField>  dirichletNodes(numLevels);
+    readBitField(dirichletNodes[0], grid.size(0, dim), dim, path + dnFile);
+
+    Array<VectorType> dirichletValues(numLevels);
+    dirichletValues[0].resize(grid.size(0, dim));
+    AmiraMeshReader<int>::readFunction(dirichletValues[0], path + dvFile);
+
+    // ///////////////////////////////////////
+    //   Assemble elasticity problem
+    // ///////////////////////////////////////
+
+    typedef Dune::LinElasticityOp<GridType, 3 > LinearElasticityType;
+    LinearElasticityType linElastOp0(grid);
+
+    linElastOp0.assembleMatrix();
+    MatrixType& stiffnessMatrix = *const_cast<MatrixType*>(linElastOp0.getMatrix());
+
+    // /////////////////////////////////////////
+    //   Modify Dirichlet rows
+    // /////////////////////////////////////////
+    typedef MatrixType::row_type::Iterator ColumnIterator;
+    
+    for (int i=0; i<stiffnessMatrix.N(); i++) {
+
+      if (dirichletNodes[0][i]) {
+
+          ColumnIterator cIt    = stiffnessMatrix[i].begin();
+          ColumnIterator cEndIt = stiffnessMatrix[i].end();
+
+          for (; cIt!=cEndIt; ++cIt) {
+              *cIt = 0;
+              if (i==cIt.index()) {
+                  for (int j=0; j<dim; j++)
+                      (*cIt)[j][j] = 1;
+              }
+          }
+
+      }
+
+    }
+
+    VectorType x(grid.size(0,dim));
+    VectorType rhs = dirichletValues[0];
+
+    //std::cout << rhs << std::endl;
+    // /////////////////////////
+    //   Compute solution
+    // /////////////////////////
+    
+    // Technicality:  turn the matrix into a linear operator
+    MatrixAdapter<MatrixType,VectorType,VectorType> op(stiffnessMatrix);
+    
+    // A preconditioner
+    SeqILU0<MatrixType,VectorType,VectorType> ilu0(stiffnessMatrix,1.0);
+    
+    // A preconditioned conjugate-gradient solver
+    int numIt = 1000;
+    CGSolver<VectorType> cg(op,ilu0,1E-4,numIt,2); 
+    
+    // Object storing some statistics about the solving process
+    InverseOperatorResult statistics;
+    
+    // Solve!
+    cg.apply(x, rhs, statistics);
+    
+
+    // ////////////////////////
+    //   Output result
+    // ////////////////////////
+    std::string gridName = transformToGridName(grid.type());
+    AmiraMeshWriter<GridType>::writeGrid(grid, "grid_" + gridName + ".am");
+    AmiraMeshWriter<GridType>::writeBlockVector(grid, x, "grid_" + gridName + ".sol");
+}
+
+
+int main() try
+{
+#if 0
+    // ///////////////////////////////////////////////////
+    //   Do the computation using an AlbertaGrid
+    // ///////////////////////////////////////////////////
+    AlbertaGrid<3,3> albertaGrid(path + "kurbelwelle.al");
+
+    compute(albertaGrid);
+
+    // ///////////////////////////////////////////////////
+    //   Do the computation using an Alu3dGrid
+    // ///////////////////////////////////////////////////
+    ALU3dGrid<3,3,tetra> alu3dGrid(path + "kurbelwelle.alu");
+
+    compute(alu3dGrid);
+#endif    
+
+    // ///////////////////////////////////////////////////
+    //   Do the computation using a UGGrid
+    // ///////////////////////////////////////////////////
+    UGGrid<3,3> uggrid;
+    AmiraMeshReader<UGGrid<3,3> >::read(uggrid, path + "kurbelwelle.am");
+
+    compute(uggrid);
+
+} catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }