diff --git a/src/Makefile.am b/src/Makefile.am
index da2fc0566462e5eb42f5633b38ca3adf10c98175..38eb13a84681803ea6f508bcb048b66888d68c9e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/..
srcincludedir = $(includedir)/dune/common
srcinclude_HEADERS = averagedistanceassembler.hh localgeodesicfefunction.hh \
- quaternion.hh rodrefine.hh averageinterface.hh \
+ quaternion.hh rodrefine.hh averageinterface.hh localstiffness.hh \
localgeodesicfestiffness.hh riemanniantrsolver.hh rodwriter.hh \
geodesicdifference.hh makestraightrod.hh rigidbodymotion.hh \
rotation.hh geodesicfeassembler.hh maxnormtrustregion.hh \
diff --git a/src/localstiffness.hh b/src/localstiffness.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3a049a91d47744cc5145fc460e685d93007c2d11
--- /dev/null
+++ b/src/localstiffness.hh
@@ -0,0 +1,184 @@
+// $Id: localstiffness.hh 560 2009-06-10 09:32:22Z sander $
+
+#ifndef DUNE_LOCALSTIFFNESS_HH
+#define DUNE_LOCALSTIFFNESS_HH
+
+#include<iostream>
+#include<vector>
+#include<set>
+#include<map>
+#include<stdio.h>
+#include<stdlib.h>
+
+#include<dune/common/timer.hh>
+#include<dune/common/fvector.hh>
+#include<dune/common/fmatrix.hh>
+#include<dune/common/array.hh>
+#include<dune/common/exceptions.hh>
+#include<dune/common/geometrytype.hh>
+#include<dune/grid/common/grid.hh>
+#include<dune/istl/operators.hh>
+#include<dune/istl/bvector.hh>
+#include<dune/istl/matrix.hh>
+
+/**
+ * @file
+ * @brief defines a class for piecewise linear finite element functions
+ * @author Peter Bastian
+ */
+
+/*! @defgroup DISC_Operators Operators
+ @ingroup DISC
+ @brief
+
+ @section D1 Introduction
+ <!--=================-->
+
+ To be written
+*/
+
+namespace Dune
+{
+ /** @addtogroup DISC_Operators
+ *
+ * @{
+ */
+ /**
+ * @brief base class for assembling local stiffness matrices
+ *
+ */
+
+
+ /*! @brief Base class for local assemblers
+
+ This class serves as a base class for local assemblers. It provides
+ space and access to the local stiffness matrix. The actual assembling is done
+ in a derived class via the virtual assemble method.
+
+ \tparam GV A grid view type
+ \tparam RT The field type used in the elements of the stiffness matrix
+ \tparam m number of degrees of freedom per node (system size)
+ */
+ template<class GV, class RT, int m>
+ class LocalStiffness
+ {
+ // grid types
+ typedef typename GV::Grid::ctype DT;
+ typedef typename GV::template Codim<0>::Entity Entity;
+ enum {n=GV::dimension};
+
+ public:
+ // types for matrics, vectors and boundary conditions
+ typedef FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
+ typedef FieldVector<RT,m> VBlockType; // one entry in the global vectors
+
+ virtual ~LocalStiffness ()
+ {
+ }
+
+ //! assemble local stiffness matrix including boundary conditions for given element and order
+ /*! On exit the following things have been done:
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method.
+ The boundary conditions are either neumann, process or dirichlet. Neumann indicates
+ that the corresponding node (assuming a nodal basis) is at the Neumann boundary, process
+ indicates that the node is at a process boundary (arising from the parallel decomposition of the mesh).
+ Process boundaries are treated as homogeneous Dirichlet conditions, i.e. the corresponding value
+ in the right hand side is set to 0. Finally, Dirichlet indicates that the node is at the Dirichlet
+ boundary.
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition.
+ It is accessible via the rhs() method.
+
+ @param[in] e a codim 0 entity reference
+ @param[in] k order of Lagrange basis (default is 1)
+ */
+ virtual void assemble (const Entity& e, int k=1) = 0;
+
+ /** \brief assemble local stiffness matrix including boundary conditions for given element and order
+
+ Unlike the method with only two arguments, this one additionally takes the local solution in order
+ to allow assembly of nonlinear operators.
+
+ On exit the following things have been done:
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method.
+ The boundary conditions are either neumann, process or dirichlet. Neumann indicates
+ that the corresponding node (assuming a nodal basis) is at the Neumann boundary, process
+ indicates that the node is at a process boundary (arising from the parallel decomposition of the mesh).
+ Process boundaries are treated as homogeneous Dirichlet conditions, i.e. the corresponding value
+ in the right hand side is set to 0. Finally, Dirichlet indicates that the node is at the Dirichlet
+ boundary.
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition.
+ It is accessible via the rhs() method.
+
+ @param[in] e a codim 0 entity reference
+ @param[in] localSolution The current solution on the entity, which is needed by nonlinear assemblers
+ @param[in] k order of Lagrange basis (default is 1)
+ */
+ virtual void assemble (const Entity& e, const BlockVector<VBlockType>& localSolution, int k=1) = 0;
+
+
+ //! print contents of local stiffness matrix
+ void print (std::ostream& s, int width, int precision)
+ {
+ // set the output format
+ s.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ int oldprec = s.precision();
+ s.precision(precision);
+
+ for (int i=0; i<currentsize(); i++)
+ {
+ s << "FEM"; // start a new row
+ s << " "; // space in front of each entry
+ s.width(4); // set width for counter
+ s << i; // number of first entry in a line
+ for (int j=0; j<currentsize(); j++)
+ {
+ s << " "; // space in front of each entry
+ s.width(width); // set width for each entry anew
+ s << mat(i,j); // yeah, the number !
+ }
+ s << std::endl;// start a new line
+ }
+
+
+ // reset the output format
+ s.precision(oldprec);
+ s.setf(std::ios_base::fixed, std::ios_base::floatfield);
+ }
+
+ //! access local stiffness matrix
+ /*! Access elements of the local stiffness matrix. Elements are
+ undefined without prior call to the assemble method.
+ */
+ const MBlockType& mat (int i, int j) const
+ {
+ return A[i][j];
+ }
+
+ //! set the current size of the local stiffness matrix
+ void setcurrentsize (int s)
+ {
+ A.setSize(s,s);
+ }
+
+ //! get the current size of the local stiffness matrix
+ int currentsize ()
+ {
+ return A.N();
+ }
+
+ protected:
+ // assembled data
+ Matrix<MBlockType> A;
+
+ };
+
+ /** @} */
+
+}
+#endif