added: integrate dune_elasticity_upscale

this merges the elasticity upscale module with the
opm upscale module. summary of changes:
 - add a new subfolder dune/elasticity with the helper classes
 - add the main program as examples/upscale_elasticity

configure.ac

@@ -37,6 +37,7 @@ AC_CONFIG_FILES([
   doc/doxygen/Doxyfile
   dune/Makefile
   dune/upscaling/Makefile
+  dune/elasticity/Makefile
   dune/upscaling/test/Makefile
   opm-upscaling.pc
   examples/Makefile

dune/Makefile.am

@@ -1,3 +1,3 @@
-SUBDIRS = upscaling
+SUBDIRS = upscaling elasticity
 
 include $(top_srcdir)/am/global-rules

dune/elasticity/Makefile.am

@@ -0,0 +1,9 @@
+noinst_LTLIBRARIES = libelasticityupscale_noinst.la
+
+libelasticityupscale_noinst_la_SOURCES = \
+        boundarygrid.cpp \
+        material.cpp \
+        materials.cpp \
+        mpc.cpp
+
+include $(top_srcdir)/am/global-rules

dune/elasticity/asmhandler.hpp

@@ -0,0 +1,223 @@
+//==============================================================================
+//!
+//! \file asmhandler.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Class handling finite element assembly
+//!
+//==============================================================================
+#pragma once
+
+//! \brief A vector holding our RHS
+typedef Dune::BlockVector<Dune::FieldVector<double,1> > Vector;
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/common/fvector.hh>
+#include <dune/elasticity/mpc.hh>
+#include <dune/elasticity/matrixops.hpp>
+
+#include <iostream>
+
+//!\brief Class handling finite element assembly
+  template<class GridType>
+class ASMHandler {
+  public:
+    //! \brief The dimension of the grid
+    static const int dim = GridType::dimension;
+
+    //! \brief A set of indices
+    typedef typename GridType::LeafGridView::IndexSet LeafIndexSet;
+
+    //! \brief An iterator over grid cells
+    typedef typename GridType::LeafGridView::template Codim<0>::Iterator LeafIterator;
+
+    //! \brief The default constructor
+    //! \param[in] gv_ The grid the operator is assembled over
+    ASMHandler(const GridType& gv_) : gv(gv_), maxeqn(0)
+    {
+    }
+
+    //! \brief Destructor
+   ~ASMHandler()
+   {
+     for (MPCMap::iterator it=mpcs.begin(); it != mpcs.end();++it)
+       delete it->second;
+   }
+
+    //! \brief A vectorial node value
+    typedef Dune::FieldVector<double,dim> NodeValue;
+
+    //! \brief Get the number of equations in the system
+    //! \returns The number of equations in the system
+    int getEqns() const
+    {
+      return maxeqn;
+    }
+
+    //! \brief Get the equation number for a given dof on a given node
+    //! \param[in] node The node number
+    //! \param[in] dof The DOF
+    //! \returns The equation number if active, -1 otherwise
+    int getEquationForDof(int node, int dof)
+    {
+      return meqn[node*dim+dof];
+    }
+
+    //! \brief Obtain a reference to the linear operator
+    //! \returns Reference to linear operator
+    Matrix& getOperator()
+    {
+      return A;
+    }
+
+    //! \brief Obtain a reference to the load vector
+    //! \returns Reference to load vector
+    Vector& getLoadVector()
+    {
+      return b;
+    }
+
+    //! \brief This function needs to be called before starting
+    //!        the element assembly. 
+    void initForAssembly();
+
+    //! \brief Add an element matrix/vector to the system
+    //! \param[in] K Pointer to the element matrix. Can be NULL
+    //! \param[in] S Pointer to the element load vector. Can be NULL
+    //! \param[in] cell An iterator pointing to the cell we're assembling for
+    //! \param[in] b Vector to add contributions to. If not given, 
+    //!              we use the internal vector
+      template<int esize>
+    void addElement(const Dune::FieldMatrix<double,esize,esize>* K,
+                    const Dune::FieldVector<double,esize>* S,
+                    const LeafIterator& cell,
+                    Vector* b=NULL);
+
+    //! \brief Extract values corresponding to cell
+    //! \param[in] u The global load vector
+    //! \param[in] it An iterator to the cell we want to extract values for
+    //! \param[out] v Vector holding the values requested
+      template<int comp>
+    void extractValues(Dune::FieldVector<double,comp>& v, 
+                       const Vector& u, const LeafIterator& it);
+
+    //! \brief Expand a system vector to a solution vector
+    void expandSolution(Vector& result, const Vector& u);
+
+    //! \brief Add a MPC
+    //! \param[in] mpc Pointer to the MPC to add.
+    //! \note This class handles destruction
+    void addMPC(MPC* mpc);
+
+    //! \brief Look for and return a MPC for a specified node+dof
+    //! \param[in] node The requested node
+    //! \param[in] dof The requested DOF at given node
+    //! \returns The MPC for the node/dof if found, else NULL
+    MPC* getMPC(int node, int dof);
+
+    //! \brief Update/add a fixed node
+    //! \param[in] node The node number
+    //! \param[in] entry The fixed values
+    void updateFixedNode(int node, 
+                         const std::pair<Direction,NodeValue>& entry);
+
+    //! \brief Check if a node is marked as fixed (in any direction)
+    //! \param[in] node The node to query for
+    bool isFixed(int node)
+    {
+      return (fixedNodes.find(node) != fixedNodes.end());
+    }
+
+    //! \brief Print the current operator
+    void printOperator() const;
+
+    //! \brief Print the current load vector
+    void printLoadVector() const;
+  protected:
+    //! \brief Resolve chained MPCs
+    void resolveMPCChains()
+    {
+      for (MPCMap::iterator it  = mpcs.begin();
+                            it != mpcs.end();++it)
+        resolveMPCChain(it->second);
+    }
+
+    //! \brief Internal function. Handles a single MPC
+    //! \param[in] mpc Pointer to the MPC to resolve
+    void resolveMPCChain(MPC* mpc);
+
+    //! \brief Internal function. Generate meqn for registered MPC/fixed nodes
+    void preprocess();
+
+    //! \brief Internal function. Generate adjacency pattern for a given node
+    //! \param[in] it Iterator pointing to the cell in of the node
+    //! \param[in] vertexsize Number of vertices in the cell
+    //! \param[in] row The equation number/row in matrix
+    void nodeAdjacency(const LeafIterator& it, int vertexsize, int row);
+
+    //! \brief Internal function. Calculate adjacency pattern 
+    void determineAdjacencyPattern();
+
+    //! \brief Internal function. Assemble entries for a single DOF
+    //! \param[in] row The row in the global matrix
+    //! \param[in] erow The row in the element matrix
+    //! \param[in] K Pointer to the element matrix. Can be NULL
+    //! \param[in] S Pointer to the element load vector. Can be NULL
+    //! \param[in] set The index set
+    //! \param[in] cell An iterator pointing to the cell we're assembling for
+    //! \param[in] b Vector to add contributions to
+    //! \param[in] scale Scale for elements. Used with MPC couplings
+      template<int esize>
+    void addDOF(int row, int erow,
+                const Dune::FieldMatrix<double,esize,esize>* K,
+                const Dune::FieldVector<double,esize>* S,
+                const LeafIndexSet& set,
+                const LeafIterator& cell,
+                Vector* b,
+                double scale=1.f);
+
+    //! \brief The set of MPC
+    MPCMap mpcs;
+
+    //! \brief Vector of (interleaved) dof<->eqn mapping
+    std::vector<int> meqn;
+
+    //! \brief Fixed nodes
+    typedef std::pair<Direction,NodeValue> fixEntry;
+
+    //! \brief A mapping from dof to fix value info
+    typedef std::map<int,fixEntry> fixMap;
+
+    //! \brief Iterator over a fixmap
+    typedef typename fixMap::iterator fixIt;
+
+    //! \brief The map holding information about our fixed nodes
+    fixMap fixedNodes;
+
+    //! \brief Holds the adjacency pattern of the sparse matrix
+    AdjacencyPattern adjacencyPattern;
+
+    //! \brief The linear operator
+    Matrix A;
+
+    //! \brief The load vector
+    Vector b;
+
+    //! \brief A reference to the grid in use
+    const GridType& gv;
+
+    //! \brief The number of equations in the system
+    size_t maxeqn;
+  private:
+    //! \brief No copying of this class
+    ASMHandler(const ASMHandler&) {}
+    //! \brief No copying of this class
+    ASMHandler& operator=(const ASMHandler&) {}
+};
+
+#include "asmhandler_impl.hpp"

dune/elasticity/asmhandler_impl.hpp

@@ -0,0 +1,399 @@
+//==============================================================================
+//!
+//! \file asmhandler_impl.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Class handling finite element assembly - template implementations
+//!
+//==============================================================================
+
+  template<class GridType>
+void ASMHandler<GridType>::initForAssembly()
+{
+  resolveMPCChains();
+  preprocess();
+  determineAdjacencyPattern();
+
+  // workaround what looks to be a bug in bcrs matrix
+  A.setBuildMode(Matrix::random);
+  A.endrowsizes();
+
+  MatrixOps::fromAdjacency(A,adjacencyPattern,maxeqn,maxeqn);
+  b.resize(maxeqn);
+  b = 0;
+
+  // print some information
+  std::cout << "Number of nodes: " << gv.size(dim) << std::endl;
+  std::cout << "Number of elements: " << gv.size(0) << std::endl;
+  std::cout << "Number of constraints: " << mpcs.size() << std::endl;
+  int fixedDofs=0;
+  for (fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) {
+    if (it->second.first & X)
+      fixedDofs++;
+    if (it->second.first & Y)
+      fixedDofs++;
+    if (it->second.first & Z)
+      fixedDofs++;
+  }
+  std::cout << "Number of fixed dofs: " << fixedDofs << std::endl;
+}
+
+  template<class GridType>
+    template<int esize>
+void ASMHandler<GridType>::addDOF(int row, int erow,
+                              const Dune::FieldMatrix<double,esize,esize>* K,
+                              const Dune::FieldVector<double,esize>* S,
+                              const LeafIndexSet& set,
+                              const LeafIterator& cell, 
+                              Vector* b,
+                              double scale)
+{
+  if (row == -1)
+    return;
+  if (K) {
+    for (int j=0;j<esize/dim;++j) {
+      int index2 = set.subIndex(*cell,j,dim);
+      for (int l=0;l<dim;++l) {
+        MPC* mpc = getMPC(index2,l);
+        if (mpc) {
+          for (int n=0;n<mpc->getNoMaster();++n) {
+            int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
+            if (idx != -1)
+              A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l];
+          }
+        } else if (meqn[index2*dim+l] != -1) {
+          A[row][meqn[index2*dim+l]] += scale*(*K)[erow][j*dim+l];
+        }
+      }
+    }
+  }
+  if (S)
+    (*b)[row] += scale*(*S)[erow];
+}
+
+  template<class GridType>
+    template<int esize>
+void ASMHandler<GridType>::addElement(
+                    const Dune::FieldMatrix<double,esize,esize>* K,
+                    const Dune::FieldVector<double,esize>* S,
+                    const LeafIterator& cell,
+                    Vector* b2)
+{
+  if (!b2)
+    b2 = &b;
+  const LeafIndexSet& set = gv.leafView().indexSet();
+  for (int i=0;i<esize/dim;++i) {
+    int index1 = set.subIndex(*cell,i,dim);
+    fixIt it = fixedNodes.find(index1);
+    if (it != fixedNodes.end() && it->second.first == XYZ)
+      continue;
+    for (int k=0;k<dim;++k) {
+      MPC* mpc = getMPC(index1,k);
+      if (mpc) {
+        for (int n=0;n<mpc->getNoMaster();++n) {
+          int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
+          addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff);
+        }
+      } else
+        addDOF(meqn[index1*dim+k],i*dim+k,K,S,set,cell,b2);
+    }
+  }
+}
+
+  template<class GridType>
+    template<int comp>
+void ASMHandler<GridType>::extractValues(Dune::FieldVector<double,comp>& v,
+                                         const Vector& u,
+                                         const LeafIterator& it)
+{
+  v = 0;
+  const LeafIndexSet& set = gv.leafView().indexSet();
+  Dune::GeometryType gt = it->type();
+  const Dune::template GenericReferenceElement<double,dim> &ref =
+                      Dune::GenericReferenceElements<double,dim>::general(gt);
+  int vertexsize = ref.size(dim);
+  int l=0;
+  for (int i=0;i<vertexsize;++i) {
+    int indexi = set.subIndex(*it,i,dim);
+    fixIt it2 = fixedNodes.find(indexi);
+    for (int n=0;n<dim;++n) {
+      MPC* mpc = getMPC(indexi,n);
+      if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
+        v[l++] = it2->second.second[n];
+      else if (mpc) {
+        for (int m=0;m<mpc->getNoMaster();++m) {
+          int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1];
+          if (idx != -1)
+            v[l] += u[idx]*mpc->getMaster(m).coeff;
+        }
+        l++;
+      } else
+        v[l++] = u[meqn[indexi*dim+n]];
+    }
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
+{
+  int nodes = gv.size(dim);
+  result.resize(nodes*dim);
+  result = 0;
+  int l=0;
+  for (size_t i=0;i<nodes;++i) {
+    fixIt it = fixedNodes.find(i);
+    Direction dir;
+    if (it == fixedNodes.end())
+      dir = NONE;
+    else
+      dir = it->second.first;
+
+    int flag=1;
+    for (int j=0;j<dim;++j) {
+      if (dir & flag)
+        result[l] = it->second.second[j];
+      else if (meqn[l] != -1)
+        result[l] = u[meqn[l]];
+      l++;
+      flag *= 2;
+    }
+  }
+  // second loop - handle MPC couplings
+  l = 0;
+  for (size_t i=0;i<nodes;++i) {
+    for (int j=0;j<dim;++j) {
+      MPC* mpc = getMPC(i,j);
+      if (mpc) {
+        for (int n=0;n<mpc->getNoMaster();++n) {
+          int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1;
+          if (meqn[idx] != -1)
+            result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff;
+        }
+      }
+      ++l;
+    }
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::addMPC(MPC* mpc)
+{
+  if (!mpc)
+    return;
+
+  int slaveNode = mpc->getSlave().node*dim+mpc->getSlave().dof-1;
+  fixIt it = fixedNodes.find(mpc->getSlave().node);
+  int flag = 1 << (mpc->getSlave().dof-1);
+  if (it == fixedNodes.end() || 
+      !(it->second.first & (1 << (mpc->getSlave().dof-1)))) {
+    mpcs.insert(std::make_pair(slaveNode,mpc));
+    return;
+  }
+
+  delete mpc;
+}
+
+  template<class GridType>
+MPC* ASMHandler<GridType>::getMPC(int node, int dof)
+{
+  if (mpcs.find(node*dim+dof) != mpcs.end())
+    return mpcs[node*dim+dof];
+
+  return NULL;
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::updateFixedNode(int node,
+                                  const std::pair<Direction,NodeValue>& entry)
+{
+  fixIt it = fixedNodes.find(node);
+  // same type or new - update values/add and return
+  if (it == fixedNodes.end() || it->second.first == entry.first) {
+    fixedNodes[node] = entry;
+    return;
+  }
+  int temp = it->second.first;
+  temp |= entry.first;
+  it->second.first = (Direction)temp;
+  int flag = 1;
+  for (int i=0;i<dim;++i) {
+    if (entry.first & flag)
+      it->second.second[i] = entry.second[i];
+    flag *= 2;
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::printOperator() const
+{
+  MatrixOps::print(A);
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::printLoadVector() const
+{
+  for (int i=0;i<b.size();++i) {
+    double val = b[i];
+    std::cout << (fabs(val)>1.e-12?val:0.f) << std::endl;
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::resolveMPCChain(MPC* mpc)
+{
+  size_t nMaster = mpc->getNoMaster();
+  if (nMaster == 0) return; // no masters, prescribed displacement only
+
+  for (size_t i = 0; i < nMaster; i++)
+  {
+    // Check if the master node has a local coordinate system attached. If yes,
+    // the slave DOF might be coupled to all (local) DOFs of the master node.
+    const MPC::DOF& master = mpc->getMaster(i);
+    Dune::FieldVector<double,dim> coeff;
+    coeff = 0;
+    coeff[master.dof-1] = master.coeff;
+
+    int removeOld = 0;
+    for (size_t d = 1; d <= dim; d++)
+      if (fabs(coeff[d-1]) > 1.0e-8)
+      {
+        MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1);
+        if (mpc2)
+        {
+          // We have a master DOF which is a slave in another MPC.
+          // Invoke resolveMPCchain recursively to ensure that all master DOFs
+          // of that equation are not slaves themselves.
+          resolveMPCChain(mpc2);
+
+          // Remove current master specification, unless it has been updated
+          if (!removeOld) removeOld = 1;
+
+          // Add constant offset from the other equation
+          mpc->addOffset(mpc2->getSlave().coeff*coeff[d-1]);
+
+          // Add masters from the other equations
+          for (size_t j = 0; j < mpc2->getNoMaster(); j++)
+            mpc->addMaster(mpc2->getMaster(j).node,
+                           mpc2->getMaster(j).dof,
+                           mpc2->getMaster(j).coeff*coeff[d-1]);
+        }
+        else
+          // The master node is free, but has a local coordinate system
+          if (d != mpc->getMaster(i).dof)
+            // Add coupling to the other local DOFs of this master node.
+            mpc->addMaster(mpc->getMaster(i).node,d,coeff[d-1]);
+          else if (coeff[d-1] != mpc->getMaster(i).coeff)
+          {
+            // Update the coupling coefficient of this master DOF
+            // due to the local-to-global transformation
+            mpc->updateMaster(i,coeff[d-1]);
+            removeOld = -1;
+          }
+      }
+      else if (d == mpc->getMaster(i).dof && !removeOld)
+        // The coefficient of the current master DOF is zero,
+        removeOld = 1; // so remove it from the contraint equation
+
+    if (removeOld == 1)
+    {
+      // Remove the old master DOF specification when it has been replaced
+      mpc->removeMaster(i--);
+      nMaster--; // we don't need to check the added masters
+    }
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::preprocess()
+{
+  int nodes = gv.size(dim);
+  meqn.resize(nodes*dim);
+
+  // iterate over nodes
+  for (int indexi=0;indexi<nodes;++indexi) {
+    fixIt it2 = fixedNodes.find(indexi);
+    if (it2 == fixedNodes.end()) {
+      for (int i=0;i<dim;++i) {
+        MPC* mpc = getMPC(indexi,i);
+        if (!mpc)
+          meqn[indexi*dim+i] = maxeqn++;
+        else
+          meqn[indexi*dim+i] = -1;
+      }
+    } else {
+      int flag=1;
+      for (int i=0;i<dim;++i) {
+        if (it2->second.first & flag)
+          meqn[indexi*dim+i] = -1;
+        else {
+          MPC* mpc = getMPC(indexi,i);
+          if (!mpc)
+            meqn[indexi*dim+i] = maxeqn++;
+          else
+            meqn[indexi*dim+i] = -1;
+        }
+        flag *= 2;
+      }
+    }
+  }
+  std::cout << "number of equations: " << maxeqn << std::endl;
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::nodeAdjacency(const LeafIterator& it,
+                                         int vertexsize, int row)
+{
+  if (row == -1)
+    return;
+  const LeafIndexSet& set = gv.leafView().indexSet();
+  for (int j=0;j<vertexsize;++j) {
+    int indexj = set.subIndex(*it,j,dim);
+    for (int l=0;l<dim;++l) {
+      MPC* mpc = getMPC(indexj,l);
+      if (mpc) {
+        for (int i=0;i<mpc->getNoMaster();++i) {
+          int idx = meqn[mpc->getMaster(i).node*dim+
+            mpc->getMaster(i).dof-1];
+          if (idx != -1)
+            adjacencyPattern[row].insert(idx);
+        }
+      } else if (meqn[indexj*dim+l] != -1)
+        adjacencyPattern[row].insert(meqn[indexj*dim+l]);
+    }
+  }
+}
+
+  template<class GridType>
+void ASMHandler<GridType>::determineAdjacencyPattern()
+{
+  adjacencyPattern.resize(maxeqn);
+
+  const LeafIndexSet& set = gv.leafView().indexSet();
+  LeafIterator itend = gv.leafView().template end<0>();
+
+  // iterate over cells
+  for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it) {
+    Dune::GeometryType gt = it->type();
+    const Dune::template GenericReferenceElement<double,dim>& ref =
+      Dune::GenericReferenceElements<double,dim>::general(gt);
+
+    int vertexsize = ref.size(dim);
+    for (int i=0; i < vertexsize; i++) {
+      int indexi = set.subIndex(*it,i,dim);
+      for (int k=0;k<dim;++k) {
+        MPC* mpc = getMPC(indexi,k);
+        if (mpc) {
+          for (int l=0;l<mpc->getNoMaster();++l) {
+            nodeAdjacency(it,vertexsize,
+                          meqn[mpc->getMaster(l).node*dim+
+                               mpc->getMaster(l).dof-1]);
+          }
+        } else
+          nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
+      }
+    }
+  }
+}

dune/elasticity/boundarygrid.cpp

@@ -0,0 +1,379 @@
+//==============================================================================
+//!
+//! \file boundarygrid.cpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Class describing 2D quadrilateral grids
+//!
+//==============================================================================
+#include "boundarygrid.hh"
+
+BoundaryGrid BoundaryGrid::uniform(const FaceCoord& min, const FaceCoord& max,
+                                   int k1, int k2, bool dc)
+{
+  double dx = (max[0]-min[0])/k1;
+  double dy = (max[1]-min[1])/k2;
+  BoundaryGrid result;
+  result.fixNodes.resize((k1+1)*(k2+1));
+  for (int i=0;i<k1;++i) {
+    for (int j=0;j<k2;++j) {
+      Quad q;
+      int k = j*(k1+1)+i;
+      result.fixNodes[k] = false;
+      q.v[0].i = k;
+      q.v[0].c[0] = min[0]+i*dx;
+      q.v[0].c[1] = min[1]+j*dy;
+      q.v[1].i = k+1;
+      result.fixNodes[k+1] = false;
+      q.v[1].c[0] = q.v[0].c[0]+dx;
+      int y=dc?3:2;
+      q.v[y].i = k+k1+2;
+      result.fixNodes[k+k1+2] = false;
+      q.v[y].c[0] = q.v[0].c[0]+dx;
+      q.v[y].c[1] = q.v[0].c[1]+dy;
+      y=dc?2:3;
+      q.v[y].i = k+k1+1;
+      result.fixNodes[k+k1+1] = false;
+      q.v[y].c[0] = q.v[0].c[0];
+      q.v[y].c[1] = q.v[0].c[1]+dy;
+      if (!dc) {
+        if (i == 0 && j == 0)
+          result.fixNodes[q.v[0].i] = true;
+        if (i == k1-1 && j == 0)
+          result.fixNodes[q.v[1].i] = true;
+        if (i == 0 && j == k2-1)
+          result.fixNodes[q.v[3].i] = true;
+        if (i == k1-1 && j == k2-1)
+          result.fixNodes[q.v[2].i] = true;
+      }
+      result.add(q);
+    }
+  }
+
+  result.nodes = (k1+1)*(k2+1);
+  return result;
+}
+
+void BoundaryGrid::extract(FaceCoord& res,
+                           const GlobalCoordinate& coord, int dir)
+{
+  int j=0;
+  for (int i=0;i<3;++i) {
+    if (i != dir)
+      res[j++] = coord[i];
+  }
+}
+
+void BoundaryGrid::extract(Vertex& res, const Dune::FieldVector<double,3>& coord, int dir)
+{
+  extract(res.c,coord,dir);
+}
+
+void BoundaryGrid::add(const Quad& quad)
+{
+  grid.push_back(quad);
+  Quad& q = grid.back();
+
+  // establish bounding box
+  q.bb[0] = q.v[0].c;
+  q.bb[1] = q.v[0].c;
+  for (int k=1;k<4;++k) {
+    if (q.v[k].c[0] < q.bb[0][0])
+      q.bb[0][0] = q.v[k].c[0];
+    if (q.v[k].c[0] > q.bb[1][0])
+      q.bb[1][0] = q.v[k].c[0];
+    if (q.v[k].c[1] < q.bb[0][1])
+      q.bb[0][1] = q.v[k].c[1];
+    if (q.v[k].c[1] > q.bb[1][1])
+      q.bb[1][1] = q.v[k].c[1];
+  }
+  if ((q.bb[1][0]-q.bb[0][0])*(q.bb[1][1]-q.bb[0][1]) < 1.e-7)
+    grid.pop_back();
+}
+
+//! \brief Check that two points are sufficiently equal
+//! \param[in] x First point
+//! \param[in] y Second point
+//! \param[in] tol Tolerance of comparison
+inline bool EQUAL2(const BoundaryGrid::FaceCoord& x,
+                   const BoundaryGrid::FaceCoord& y, double tol)
+{
+  return hypot(x[0]-y[0],x[1]-y[1]) < tol;
+}
+
+bool BoundaryGrid::find(Vertex& res, const Vertex& coord) const
+{
+  // find first quad with coord within bounding box 
+  std::vector<Quad>::const_iterator it = std::find_if(grid.begin(),grid.end(),
+                                                      BoundedPredicate(coord.c));
+
+  res.i = -1;
+  while (it != grid.end()) {
+    res.q = const_cast<Quad*>(&(*it));
+    // check if we have exactly a node
+    bool ok=false;
+    for (int i=0;i<4;++i) {
+      if (EQUAL2(coord.c,it->v[i].c,1.e-8)) {
+        res.i = it->v[i].i;
+        ok = true;
+        break;
+      }
+    } 
+    if (!ok && Q4inv(res.c,*it,coord.c,1.e-8,1.e-8) > 0) {
+      ok = true;
+    }
+
+    if (ok)
+      break;
+    it = std::find_if(it+1,grid.end(),BoundedPredicate(coord.c));
+  }
+  if (it == grid.end()) {
+    std::cout << " failed to locate " << coord.c << std::endl;
+    assert(0);
+  }
+  return it != grid.end();
+}
+
+int BoundaryGrid::Q4inv(FaceCoord& res, const Quad& q,
+                        const FaceCoord& coord,
+                        double epsZero, double epsOut) const
+{
+  double A[4];
+  double B[4];
+
+  /* Find coefficients of the bi-linear equation */
+  A[0] =        ( q.v[0].c[0] - q.v[1].c[0] +
+                  q.v[2].c[0] - q.v[3].c[0]);
+  A[1] =        ( -q.v[0].c[0]+q.v[1].c[0]);
+  A[2] =        ( -q.v[0].c[0]+q.v[3].c[0]);
+  A[3] = coord[0]-q.v[0].c[0];
+
+  B[0] =        ( q.v[0].c[1] - q.v[1].c[1] +
+                  q.v[2].c[1] - q.v[3].c[1]);
+  B[1] =        ( -q.v[0].c[1]+q.v[1].c[1]);
+  B[2] =        ( -q.v[0].c[1]+q.v[3].c[1]);
+  B[3] = coord[1]-q.v[0].c[1];
+
+  // We have to solve the following set of equations:
+
+  //  A1*XI*ETA + A2*XI + A3*ETA = A4
+  //  B1*XI*ETA + B2*XI + B3*ETA = B4
+
+  //  The way that we may solve this nonlinear (in XI and ETA)
+  //  set of equations depends on the coefficients Ai,Bi.
+  //  The solution is unique for proper input.
+  std::vector<double> xi;
+  std::vector<double> eta;
+  bilinearSolve(epsZero,1,A,B,xi,eta);
+
+  // check that obtained solutions are inside element
+  double tol = 1+epsOut;
+  int nInside=0;
+  for (int i=0;i<xi.size();++i) {
+    if (xi[i] < tol && eta[i] < tol) {
+      if (++nInside > 1) {
+        std::cout << "multiple solutions" << std::endl;
+        FaceCoord old = q.pos(res[0],res[1]);
+        FaceCoord ny  = q.pos(xi[i],eta[i]);
+        double d1 = hypot(coord[0]-ny[0],coord[1]-ny[1]);
+        double d2 = hypot(coord[0]-old[0],coord[1]-old[1]);
+        if (d2 < d1)
+          continue;
+      }
+    } else if (nInside == 0) {
+      if (i > 0) {
+        FaceCoord old = q.pos(res[0],res[1]);
+        FaceCoord ny  = q.pos(xi[i],eta[i]);
+        double d1 = hypot(coord[0]-ny[0],coord[1]-ny[1]);
+        double d2 = hypot(coord[0]-old[0],coord[1]-old[1]);
+        if (d2 < d1)
+          continue;
+      }
+    }
+    res[0] = xi[i];
+    res[1] = eta[i];
+  }
+
+  return nInside;
+}
+
+bool BoundaryGrid::bilinearSolve(double epsilon, double order,
+                                 const double* A, const double* B,
+                                 std::vector<double>& X,
+                                 std::vector<double>& Y) const
+{
+  double tol = 0;
+  // geometric tolerance ?
+  for (int i=0;i<4;++i) {
+    double det = fabs(A[i]);
+    if (det > tol) tol = det;
+    det = fabs(B[i]);
+    if (det > tol) tol = det;
+  }
+  tol *= epsilon;
+
+  double det = A[1]*B[2]-B[1]*A[2];
+  if (fabs(A[0]) < tol && fabs(B[0]) < tol) {
+    // linear eqs
+    if (fabs(det) < tol*tol) return false;
+    X.push_back((B[2]*A[3]-A[2]*B[3])/det);
+    Y.push_back((-B[1]*A[3]+A[1]*B[3])/det);
+    return true;
+  }
+  // second order
+  double Q0 = B[2]*A[3]-A[2]*B[3];
+  double Q1 = B[0]*A[3]-A[0]*B[3]-det;
+  double Q2 = A[0]*B[1]-B[0]*A[1];
+  std::vector<double> Z;
+  cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
+  for (int i=0;i<Z.size();++i) {
+    Q0 = A[0]*Z[i]+A[2];
+    if (fabs(Q0) > tol) {
+      X.push_back(Z[i]);
+      Y.push_back((A[3]-A[1]*Z[i])/Q0);
+    }
+  }
+  Q0 = B[1]*A[3]-A[1]*B[3];
+  Q1 = B[0]*A[3]-A[0]*B[3]+det;
+  Q2 = A[0]*B[2]-B[0]*A[2];
+  Z.clear();
+  cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
+  for (int i=0;i<Z.size();++i) {
+    Q0 = A[0]*Z[i]+A[1];
+    if (fabs(Q0) > tol) {
+      int j=0;
+      for (j=0;j<Y.size();++j)
+        if (fabs(Y[j]-Z[i]) <= epsilon*order) break;
+      if (j == Y.size()) {
+        X.push_back((A[3]-A[2]*Z[i])/Q0);
+        Y.push_back(Z[i]);
+      }
+    }
+  }
+
+  return X.size() > 0;
+}
+
+bool BoundaryGrid::cubicSolve(double eps, double A, double B, double C,
+                              double D, std::vector<double>& X) const
+{
+  if (fabs(A) > eps) { // cubic
+    double epsmall = pow(eps,6.f);
+    double P = (C-B*B/(3*A))/(3*A);
+    double Q = ((2*B*B/(27*A)-C/3)*B/A+D)/(2*A);
+    double W = Q*Q+P*P*P;
+    if (W <= -epsmall && P < 0) {
+      double FI = acos(-Q/sqrt(-P*P*P));
+      X.push_back( 2*sqrt(-P)*cos(FI/3));
+      X.push_back(-2*sqrt(-P)*cos((FI+M_PI)/3));
+      X.push_back(-2*sqrt(-P)*cos((FI-M_PI)/3));
+    } else if (fabs(W) < epsmall && Q < 0) {
+      X.push_back(2*pow(-Q,1.f/3));
+      X.push_back(-.5f*X[0]);
+      X.push_back(X[1]);
+    } else if (W > -epsmall && Q+sqrt(W) < 0 && Q-sqrt(W) < 0) {
+      X.push_back(pow(-Q+sqrt(W),1.f/3)+pow(-Q-sqrt(W),1.f/3));
+      X.push_back(-.5f*X[0]);
+      X.push_back(X[1]);
+    } else if (W >= epsmall && fabs(Q) > epsmall && P > 0) {
+      double FI = atan(sqrt(P*P*P)/Q);
+      double KI = atan(pow(tan(.5f*FI),1.f/3));
+      X.push_back(-2*sqrt(P)/tan(KI+KI));
+      X.push_back(-.5f*X[0]);
+      X.push_back(X[1]);
+    } else if (W > -epsmall && fabs(Q) > epsmall && P < 0)
+      return false;
+    else
+      return false;
+
+    W = B/(3*A);
+    X[0] -= W;
+    X[1] -= W;
+    X[2] -= W;
+
+    return true;
+  } else if (fabs(B) > eps) {
+    double epsmall = pow(eps,4.f);
+    double P = C*C-4*B*D;
+    if (P > 0) {
+      double Q = sqrt(P);
+      X.push_back((-C+Q)/(B+B));
+      X.push_back((-C-Q)/(B+B));
+    } else if (P > -epsmall) {
+      X.push_back(-C/(B+B));
+      X.push_back(X[0]);
+    }
+    else
+      return false;
+  } else if (fabs(C) > eps) {
+    X.push_back(-D/C);
+  } else
+    return false;
+
+  return true;
+}
+
+BoundaryGrid::FaceCoord BoundaryGrid::Quad::pos(double xi, double eta) const
+{
+  BoundaryGrid::FaceCoord res;
+  std::vector<double> N = evalBasis(xi,eta);
+
+  res[0] = N[0]*v[0].c[0]+N[1]*v[1].c[0]+N[2]*v[2].c[0]+N[3]*v[3].c[0];
+  res[1] = N[0]*v[0].c[1]+N[1]*v[1].c[1]+N[2]*v[2].c[1]+N[3]*v[3].c[1];
+
+  return res;
+}
+
+std::vector<double> BoundaryGrid::Quad::evalBasis(double xi, double eta) const
+{
+  std::vector<double> res;
+  res.push_back((1-xi)*(1-eta));
+  res.push_back(xi*(1-eta));
+  res.push_back(xi*eta);
+  res.push_back((1-xi)*eta);
+
+  return res;
+}
+
+BoundaryGrid::Vertex minXminY(std::vector<BoundaryGrid::Vertex>& in)
+{
+  // find the nodes with minimal X
+  // then find the minimum Y among these 
+  std::vector<BoundaryGrid::Vertex> s(in);
+  std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0));
+  std::sort(s.begin(),s.begin()+2,BoundaryGrid::VertexLess(1));
+  return *s.begin();
+}
+
+BoundaryGrid::Vertex maxXminY(std::vector<BoundaryGrid::Vertex>& in)
+{
+  // find the nodes with maximum X
+  // then find the minimum Y among these 
+  std::vector<BoundaryGrid::Vertex> s(in);
+  std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0));
+  std::sort(s.begin()+2,s.end(),BoundaryGrid::VertexLess(1));
+  return *(s.end()-2);
+}
+
+BoundaryGrid::Vertex maxXmaxY(std::vector<BoundaryGrid::Vertex>& in)
+{
+  // find the nodes with maximum X
+  // then find the maximum Y among these 
+  std::vector<BoundaryGrid::Vertex> s(in);
+  std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0));
+  std::sort(s.begin()+2,s.end(),BoundaryGrid::VertexLess(1));
+  return *(s.end()-1);
+}
+
+BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in)
+{
+  // find the nodes with minimum X
+  // then find the maximum Y among these 
+  std::vector<BoundaryGrid::Vertex> s(in);
+  std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0));
+  std::sort(s.begin(),s.begin()+2,BoundaryGrid::VertexLess(1));
+  return *(s.begin()+1);
+}

dune/elasticity/boundarygrid.hh

@@ -0,0 +1,479 @@
+//==============================================================================
+//!
+//! \file boundarygrid.hh
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Class describing 2D quadrilateral grids
+//!
+//==============================================================================
+#pragma once
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/genericgeometry/matrixhelper.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <vector>
+
+//! \brief A class describing a quad grid
+class BoundaryGrid {
+  public:
+    //! \brief A coordinate on the underlying 3D grid
+    typedef Dune::FieldVector<double,3> GlobalCoordinate;
+    //! \brief A coordinate on the quad grid
+    typedef Dune::FieldVector<double,2> FaceCoord;
+
+    //! \brief Establish an uniform quad grid
+    //! \param[in] min Lower left corner
+    //! \param[in] max Upper right corner
+    //! \param[in] k1 Number of elements in the first direction
+    //! \param[in] k2 Number of elements in the second direction
+    //! \param[in] dc If true, order quads according to dune conventions
+    //! \returns A quad grid spanning the given area. Nodes are numbered
+    //              in natural order along the first direction
+    static BoundaryGrid uniform(const FaceCoord& min, const FaceCoord& max,
+                                int k1, int k2, bool dc=false);
+
+    //! \brief Default (empty) constructor
+    BoundaryGrid() {}
+
+    //! \brief Default (empty) destructor
+    virtual ~BoundaryGrid() {}
+
+    //! \brief Holds the indices and relevant coordinates of the vertices 
+    //         on a boundary
+    class Quad;
+
+    //! \brief A class describing a 2D vertex
+    class Vertex {
+      public:
+        //! \brief Default constructor
+        Vertex() : i(-1), c(0), q(0), fixed(false) {}
+
+        //! \brief Index of the vertex
+        int i;
+
+        //! \brief Coordinates of the vertex (2D)
+        FaceCoord c;
+
+        //! \brief The quad containing the vertex (if search has been done)
+        Quad* q;
+
+        //! \brief Whether or not this node is fixed
+        bool fixed;
+
+        //! \brief Comparison operator
+        bool operator==(const Vertex& v2)
+        {
+          return hypot(v2.c[0]-c[0],v2.c[1]-c[1]) < 1.e-8;
+        }
+    };
+    //! \brief A class describing a linear, quadrilateral element
+    class Quad {
+      public:
+        //! \brief Default constructor
+        Quad()
+        {
+          v[0].i = v[1].i = v[2].i = v[3].i = 0;
+          v[0].c = v[1].c = v[2].c = v[3].c = 0.f;
+        }
+        //! \brief Return the physical coordinates corresponding to the
+        //!        given local coordinates
+        //! \param[in] xi The local coordinate in the first direction
+        //! \param[in] eta The local coordinate in the second direction
+        FaceCoord pos(double xi, double eta) const;
+
+        //! \brief Evaluate the basis functions
+        //! \param[in] xi The local coordinate in the first direction
+        //! \param[in] eta The local coordinate in the second direction
+        std::vector<double> evalBasis(double xi, double eta) const;
+
+        //! \brief Vertices
+        Vertex v[4];
+        //! \brief Bounding box
+        FaceCoord bb[2];
+      protected:
+        //! \brief Print to a stream
+        friend std::ostream& operator <<(std::ostream& os, const Quad& q)
+        {
+          os << "Nodes: " << q.v[0].i << "," << q.v[1].i << "," << q.v[2].i << "," << q.v[3].i << std::endl;
+          os << "(" << q.v[0].c << ")(" << q.v[1].c << ")(" << q.v[2].c << ")(" << q.v[3].c << ")";
+          return os;
+        }
+    };
+
+    //! \brief Add a quad to the grid
+    //! \param[in] quad The quad to add
+    void add(const Quad& quad);
+
+    //! \brief Obtain a reference to a quad
+    //! \param[in] index The index of the requested quad
+    //! \returns A reference to the requested quad
+    Quad& operator[](int index)
+    {
+      return grid[index];
+    }
+
+    //! \brief Obtain a const reference to a quad
+    //! \param[in] index The index of the requested quad
+    //! \returns A const reference to the requested quad
+    const Quad& operator[](int index) const
+    {
+      return grid[index];
+    }
+
+    //! \brief Obtain the number of quads in the grid
+    size_t size() const
+    {
+      return grid.size();
+    }
+
+    //! \brief Return the total number of nodes on the grid when known
+    //! \sa uniform
+    size_t totalNodes() const
+    {
+      return nodes;
+    }
+
+    //! \brief Check if a given node is marked as fixed
+    //! \param[in] node The requested node
+    //! \returns Whether or not the node is marked as fixed
+    bool isFixed(int node) const
+    {
+      return fixNodes[node];
+    }
+
+    //! \brief Locate the position of a vertex on the grid
+    //! \param[in] coord The coordinate of the vertex
+    //! \param[out] res The resulting coordinates
+    bool find(Vertex& res, const Vertex& coord) const;
+
+    //! \brief Helper function for extracting given 2D coordinates from a 3D coordinate
+    //! \param[in] coord The 3D coordinates of the vertex
+    //! \param[in] dir The direction to ignore
+    //! \param[out] res The resulting coordinates
+    static void extract(FaceCoord& res,
+                        const GlobalCoordinate& coord, int dir);
+
+    //! \brief Helper function for extracting given 2D coordinates from a 3D vector
+    //! \param[in] coord The 3D coordinates of the vertex
+    //! \param[in] dir The direction to ignore
+    //! \param[out] res The resulting coordinates
+    static void extract(Vertex& res,
+                        const Dune::FieldVector<double,3>& coord, int dir);
+
+    //! \brief Predicate for sorting vertices
+    struct VertexLess {
+        //! \brief Default constructor.
+        //! \param[in] comp Direction to use for comparison. -1 to use index
+        VertexLess(int comp) : dir(comp) {}
+
+        //! \brief The comparison operator
+        bool operator()(const Vertex& q1, const Vertex& q2)
+        {
+          if (dir >= 0)
+            return q1.c[dir] < q2.c[dir];
+          return q1.i < q2.i;
+        }
+
+        //! \brief Compare using this direction, if -1 compare using index values
+        int dir;
+    };
+
+    //! \brief Predicate for checking if a vertex falls within a quads bounding box
+    struct BoundedPredicate {
+      //! \brief Default constructor
+      //! \param[in] coord_ The coordinates to check
+      BoundedPredicate(const FaceCoord& coord_) : coord(coord_) {}
+
+      //! \brief The comparison operator
+      bool operator()(const Quad& q)
+      {
+        double eps = 1.e-8;
+        return (coord[0] >= q.bb[0][0]-eps && 
+                coord[0] <= q.bb[1][0]+eps &&
+                coord[1] >= q.bb[0][1]-eps &&
+                coord[1] <= q.bb[1][1]+eps);
+      }
+
+      //! \brief The coordinates to check
+      FaceCoord coord;
+    };
+  protected:
+    //! \brief Our quadrilateral elements
+    std::vector<Quad> grid;
+
+    //! \brief Whether or not a given node is marked as fixed
+    std::vector<bool> fixNodes;
+
+    //! \brief Total number of nodes on grid
+    size_t nodes;
+
+    //! \brief Print to a stream
+    friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
+    {
+      for (int i=0;i<g.size();++i)
+        os << g[i] << std::endl;
+      return os;
+    }
+
+    //! \brief Solves a bi-linear set of equations in x and y.
+    //!        A1 * x*y  +  A2 * x  +  A3 * y  =  A4
+    //!        B1 * x*y  +  B2 * x  +  B3 * y  =  B4
+    //! \param[in] epsilon The tolerance for equality checks with zero
+    //! \param[in] order The expected order of the solution (used for unique checks)
+    //! \param[in] A The coefficients of the first equation
+    //! \param[in] B The coefficients of the second equation
+    //! \param[out] X The first component of the solutions
+    //! \param[out] Y The second component of the solutions
+    bool bilinearSolve(double epsilon, double order,
+                       const double* A, const double* B,
+                       std::vector<double>& X,
+                       std::vector<double>& Y) const;
+
+    //! \brief Solves the cubic equation A*x^3+B*x^2+C*x+D=0
+    //! \param[in] eps The tolerance for equality checks with zero
+    //! \param[in] A Equation coefficient
+    //! \param[in] B Equation coefficient
+    //! \param[in] C Equation coefficient
+    //! \param[in] D Equation coefficient
+    //! \param[out] X The obtained solutions
+    bool cubicSolve(double eps, double A, double B, double C,
+                    double D, std::vector<double>& X) const;
+
+    //! \brief Find the local coordinates of a given point within a given quad
+    //! \param[in] q The quad to search within
+    //! \param[in] coord The coordinates to search for
+    //! \param[in] epsZero The tolerance for equality checks with zero
+    //! \param[in] epsOut The tolerance check for outside checks
+    //! \param[out] res The obtained result
+    int Q4inv(FaceCoord& res, const Quad& q, const FaceCoord& coord,
+              double epsZero, double epsOut) const;
+};
+
+//! \brief Geometry class for general hexagons
+  template<int mydim, int cdim, class GridImp>
+class HexGeometry
+{
+};
+
+//! \brief Specialization for 2D quadrilaterals
+  template<int cdim, class GridImp>
+class HexGeometry<2, cdim, GridImp>
+{
+  public:
+    //! \brief The dimension of the grid.
+    enum {   dimension = 2};
+
+    //! \brief Dimension of the domain space
+    enum { mydimension = 2};
+
+    //! \brief Dimension of the range space
+    enum { coorddimension = cdim };
+
+    //! \brief World dimension of underlying grid
+    enum {dimensionworld = 2 };
+
+    //! \brief Coordinate element type
+    typedef double ctype;
+
+    //! \brief Domain type
+    typedef Dune::FieldVector<ctype,mydimension> LocalCoordinate;
+
+    //! \brief Range type
+    typedef Dune::FieldVector<ctype,coorddimension> GlobalCoordinate;
+
+    //! \brief Type of Jacobian matrix
+    typedef Dune::FieldMatrix<ctype,coorddimension,mydimension> Jacobian;
+
+    //! \brief Type of transposed Jacobian matrix
+    typedef Dune::FieldMatrix<ctype,mydimension,coorddimension> JacobianTransposed;
+
+    //! \brief Construct integration element extracted from a 3D grid
+    //! \param[in] q Quad describing element
+    //! \param[in] gv Underlying 3D grid quads are extracted from
+    //! \param[in] dir The direction of the normal vector on the face
+    HexGeometry(const BoundaryGrid::Quad& q, const GridImp& gv, int dir)
+    {
+      Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridImp,
+                                           Dune::MCMGVertexLayout> mapper(gv);
+      typename GridImp::LeafGridView::template Codim<3>::Iterator start=gv.leafView().template begin<3>();
+      const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafView().template end<3>();
+      for (int i=0;i<4;++i) {
+        typename GridImp::LeafGridView::template Codim<3>::Iterator it=start;
+        for (it ; it != itend; ++it) {
+          if (mapper.map(*it) == q.v[i].i)
+            break;
+        }
+        BoundaryGrid::extract(c[i],it->geometry().corner(0),dir);
+      }
+
+      m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
+    }
+
+    //! \brief Construct integration element
+    //! \param[in] q Quad describing element
+    HexGeometry(const BoundaryGrid::Quad& q)
+    {
+      for (int i=0;i<4;++i)
+        c[i] = q.v[i].c;
+      m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
+    }
+
+    //! \brief Returns entity type (a 2D cube)
+    Dune::GeometryType type() const
+    {
+      Dune::GeometryType t;
+      t.makeCube(mydimension);
+      return t;
+    }
+
+    //! \brief Returns number of corners
+    int corners() const
+    {
+      return 4;
+    }
+
+    //! \brief Returns volume (area) of quadrilateral
+    ctype volume() const
+    {
+      return m_volume;
+    }
+
+    //! \brief Returns center of quadrilateral
+    GlobalCoordinate center() const
+    {
+      LocalCoordinate local;
+      local = .5f;
+      return Global(local);
+    }
+
+    //! \brief Returns coordinates to requested corner
+    //! \param[in] cor The requested corner (0..3)
+    GlobalCoordinate corner(int cor) const
+    {
+      return c[cor];
+    }
+
+    //! \brief Map from local coordinates to global coordinates
+    //! \param[in] local The local coordinates
+    GlobalCoordinate global(const LocalCoordinate& local) const
+    {
+      // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
+      LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
+      uvw[0] -= local;
+      // Access pattern for uvw matching ordering of corners.
+      const int pat[4][2] = {{ 0, 0 },
+                             { 1, 0 },
+                             { 0, 1 },
+                             { 1, 1 }};
+      GlobalCoordinate xyz(0.0);
+      for (int i = 0; i < 4; ++i) {
+        GlobalCoordinate corner_contrib = corner(i);
+        double factor = 1.0;
+        for (int j = 0; j < 2; ++j) {
+          factor *= uvw[pat[i][j]][j];
+        }
+        corner_contrib *= factor;
+        xyz += corner_contrib;
+      }
+      return xyz;
+    }
+
+    //! \brief Map from global coordinates to local coordinates
+    //! \param[in] y The global coordinates
+    LocalCoordinate local(const GlobalCoordinate& y) const
+    {
+      const ctype epsilon = 1e-10;
+      const Dune::GenericReferenceElement< ctype , 2 > & refElement =
+        Dune::GenericReferenceElements< ctype, 2 >::general(type());
+      LocalCoordinate x = refElement.position(0,0);
+      LocalCoordinate dx;
+      int i=0;
+      do {
+        using namespace Dune::GenericGeometry;
+        // DF^n dx^n = F^n, x^{n+1} -= dx^n
+        JacobianTransposed JT = jacobianTransposed(x);
+        GlobalCoordinate z = global(x);
+        z -= y;
+        MatrixHelper<DuneCoordTraits<double> >::template xTRightInvA<2, 2>(JT, z, dx );
+        x -= dx;
+      } while (dx.two_norm2() > epsilon*epsilon);
+      return x;
+    }
+
+    //! \brief Return the transposed jacobian
+    //! \param[in] local The local coordinates
+    const Dune::FieldMatrix<ctype, mydimension, coorddimension>
+      jacobianTransposed(const LocalCoordinate& local) const
+    {
+      // uvw = { (1-u, 1-v, (u, v) }
+      LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
+      uvw[0] -= local;
+      // Access pattern for uvw matching ordering of corners.
+      const int pat[4][3] = {{ 0, 0 },
+                             { 1, 0 },
+                             { 0, 1 },
+                             { 1, 1 }};
+      Dune::FieldMatrix<ctype, mydimension, coorddimension> Jt(0.0);
+      for (int i = 0; i < 4; ++i) {
+        for (int deriv = 0; deriv < 2; ++deriv) {
+          // This part contributing to dg/du_{deriv}
+          double factor = 1.0;
+          for (int j = 0; j < 2; ++j) {
+            factor *= (j != deriv) ? uvw[pat[i][j]][j]
+              : (pat[i][j] == 0 ? -1.0 : 1.0);
+          }
+          GlobalCoordinate corner_contrib = corner(i);
+          corner_contrib *= factor;
+          Jt[deriv] += corner_contrib; // using FieldMatrix row access.
+        }
+      }
+      return Jt;
+    }
+
+    //! \brief Returns the inverse, transposed Jacobian
+    //! \param[in] local The local coordinates
+    const Dune::FieldMatrix<ctype, coorddimension, mydimension> 
+      jacobianInverseTransposed(const LocalCoordinate& local) const
+    {
+      Dune::FieldMatrix<ctype, coorddimension, mydimension> Jti = jacobianTransposed(local);
+      Jti.invert();
+      return Jti;
+    }
+
+    //! \brief Returns the integration element (|J'*J|)^(1/2)
+    //! \param[in] local The local coordinates
+    ctype integrationElement(const LocalCoordinate& local) const
+    {
+      Dune::FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local);
+      using namespace Dune::GenericGeometry;
+      return MatrixHelper<DuneCoordTraits<double> >::template sqrtDetAAT<2, 2>(Jt);
+    }
+  private:
+    //! \brief The coordinates of the corners
+    GlobalCoordinate c[4];
+
+    //! \brief The volume (area) of the quadrilateral
+    ctype m_volume;
+};
+
+//! \brief Find the vertex in the vector with minimum X and minimum Y
+//! \returns The requested vertex
+BoundaryGrid::Vertex minXminY(std::vector<BoundaryGrid::Vertex>& in);
+
+//! \brief Find the vertex in the vector with maximum X and minimum Y
+//! \returns The requested vertex
+BoundaryGrid::Vertex maxXminY(std::vector<BoundaryGrid::Vertex>& in);
+
+//! \brief Find the vertex in the vector with maximum X and maximum Y
+//! \returns The requested vertex
+BoundaryGrid::Vertex maxXmaxY(std::vector<BoundaryGrid::Vertex>& in);
+
+//! \brief Find the vertex in the vector with minimum X and maximum Y
+//! \returns The requested vertex
+BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in);

dune/elasticity/elasticity.hpp

@@ -0,0 +1,65 @@
+//==============================================================================
+//!
+//! \file elasticity.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Elasticity helper class
+//!
+//==============================================================================
+#pragma once
+
+//! \brief Elasticity helper class
+  template<class GridType>
+class Elasticity {
+  public:
+    //! \brief The dimension of our grid
+    static const int dim = GridType::dimension;
+
+    //! \brief A basic number
+    typedef typename GridType::LeafGridView::ctype ctype;
+
+    //! \brief Default constructor
+    //! \param[in] gv_ The grid we are doing the calculations on
+    Elasticity(const GridType& gv_) : gv(gv_) {}
+
+    //! \brief Returns the B matrix in a quadrature point
+    //! \param[in] point (Reference) coordinates of quadrature point 
+    //! \param[in] Jinv Jacobian matrix in quadrature point
+    //! \param[out] B The B matrix
+      template<int components, int funcdim>
+    void getBmatrix(Dune::FieldMatrix<ctype,components,funcdim>& B,
+                    const Dune::FieldVector<ctype,dim>& point,
+                    const Dune::FieldMatrix<ctype,dim,dim>& Jinv);
+
+    //! \brief Return the stiffness matrix contributions in a quadrature point
+    //! \param[in] B The B matrix in the quadrature point
+    //! \param[in] C The constitutive matrix for the cell material
+    //! \param[in] detJW Det J times quadrature weight
+    //! \param[out] A The stiffness matrix contributions in the quadrature point
+      template<int comp, int funcdim>
+    void getStiffnessMatrix(Dune::FieldMatrix<ctype,funcdim,funcdim>& A,
+                            const Dune::FieldMatrix<ctype,comp,funcdim>& B,
+                            const Dune::FieldMatrix<ctype,comp,comp>& C,
+                            ctype detJW);
+
+    //! \brief Return the stress vector in a quadrature point
+    //! \param[in] v The displacements in the quadrature point
+    //! \param[in] eps0 The load case vector
+    //! \param[in] B The B matrix in the quadrature point
+    //! \param[in] C The constitutive matrix for the cell material
+    //! \param[out] sigma The stress vector in the given quadrature point
+      template<int comp, int funcdim>
+    void getStressVector(Dune::FieldVector<ctype,comp>& sigma,
+                         const Dune::FieldVector<ctype,funcdim>& v, 
+                         const Dune::FieldVector<ctype,comp>& eps0,
+                         const Dune::FieldMatrix<ctype,comp,funcdim>& B,
+                         const Dune::FieldMatrix<ctype,comp,comp>& C);
+  protected:
+     //! \brief Const reference to our grid
+    const GridType& gv;
+};
+
+#include "elasticity_impl.hpp"

dune/elasticity/elasticity_impl.hpp

@@ -0,0 +1,92 @@
+//==============================================================================
+//!
+//! \file elasticity_impl.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Elasticity helper class - template implementations
+//!
+//==============================================================================
+
+  template<class GridType>
+    template<int components, int funcdim>
+void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim>& B,
+                               const Dune::FieldVector<ctype,dim>& point,
+                               const Dune::FieldMatrix<ctype,dim,dim>& Jinv)
+{
+  P1ShapeFunctionSet<ctype,ctype,dim> basis 
+                             = P1ShapeFunctionSet<ctype,ctype,dim>::instance();
+  int funcs = funcdim/dim;
+
+  Dune::FieldMatrix<ctype,funcdim/dim,dim> dNdX;
+  for (int i=0;i < funcs; i++)
+    Jinv.mv(basis[i].evaluateGradient(point),dNdX[i]);
+  static const int rows = 6+(dim-2)*12;
+
+  Dune::FieldMatrix<ctype,rows,funcdim/dim> temp;
+  temp = 0;
+
+  if (dim == 3) {
+#define INDEX(i,j) i+6*j
+    for (size_t i=0; i < funcs; ++i) {
+      // normal strain part
+      temp[INDEX(0,0)][i] = dNdX[i][0];
+      temp[INDEX(1,1)][i] = dNdX[i][1];
+      temp[INDEX(2,2)][i] = dNdX[i][2];
+
+      // shear strain part
+      temp[INDEX(3,0)][i] = dNdX[i][1];
+      temp[INDEX(3,1)][i] = dNdX[i][0];
+
+      temp[INDEX(4,0)][i] = dNdX[i][2];
+      temp[INDEX(4,2)][i] = dNdX[i][0];
+      temp[INDEX(5,1)][i] = dNdX[i][2];
+      temp[INDEX(5,2)][i] = dNdX[i][1];
+    }
+  } else if (dim == 2) {
+#undef INDEX
+#define INDEX(i,j) i+3*j
+    for (size_t i=0; i < funcs; ++i) {
+      // normal strain part
+      temp[INDEX(0,0)][i] = dNdX[i][0];
+      temp[INDEX(1,1)][i] = dNdX[i][1];
+
+      // shear strain part
+      temp[INDEX(2,0)][i] = dNdX[i][1];
+      temp[INDEX(2,1)][i] = dNdX[i][0];
+    }
+  }
+  size_t k=0;
+  for (size_t j=0;j<funcs*dim;++j)
+    for (size_t i=0;i<components;++i,++k)
+      B[i][j] = temp[k%rows][k/rows];
+}
+
+  template <class GridType>
+    template<int comp, int funcdim>
+void Elasticity<GridType>::getStiffnessMatrix(
+                                Dune::FieldMatrix<ctype,funcdim,funcdim>& A,
+                                const Dune::FieldMatrix<ctype,comp,funcdim>& B,
+                                const Dune::FieldMatrix<ctype,comp,comp>& C,
+                                ctype detJW)
+{
+  Dune::FieldMatrix<ctype,funcdim,comp> B2;
+  for (int i=0;i<comp;++i)
+    for (int j=0;j<funcdim;++j)
+      B2[j][i] = B[i][j];
+  A = B.leftmultiplyany(C).leftmultiplyany(B2);
+  A *= detJW;
+}
+
+  template<class GridType>
+    template<int comp, int funcdim>
+void Elasticity<GridType>::getStressVector(Dune::FieldVector<ctype,comp>& sigma,
+                             const Dune::FieldVector<ctype,funcdim>& v,
+                             const Dune::FieldVector<ctype,comp>& eps0,
+                             const Dune::FieldMatrix<ctype,comp,funcdim>& B,
+                             const Dune::FieldMatrix<ctype,comp,comp>& C)
+{
+  sigma = Dune::FMatrixHelp::mult(C,Dune::FMatrixHelp::mult(B,v)+eps0);
+}

dune/elasticity/elasticity_upscale.hpp

@@ -0,0 +1,334 @@
+//==============================================================================
+//!
+//! \file elasticity_upscale.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Elasticity upscale class
+//!
+//==============================================================================
+#pragma once
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <dune/common/fmatrix.hh>
+#include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/istl/ilu.hh>
+#include <dune/istl/solvers.hh>
+#if HAVE_SUPERLU
+#include <dune/istl/superlu.hh>
+#endif
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/overlappingschwarz.hh>
+#include <dune/grid/CpGrid.hpp>
+#include <dune/elasticity/shapefunctions.hpp>
+
+#include <dune/elasticity/asmhandler.hpp>
+#include <dune/elasticity/boundarygrid.hh>
+#include <dune/elasticity/elasticity.hpp>
+#include <dune/elasticity/materials.hh>
+#include <dune/elasticity/mpc.hh>
+
+//! \brief An enumeration of available linear solvers
+enum Solver {
+  SLU,
+   CG 
+};
+
+typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::SymmetricMultiplicativeSchwarzMode> OverlappingSchwarz;
+
+//! \brief The main driver class
+  template<class GridType>
+class ElasticityUpscale
+{
+  public:
+    //! \brief Dimension of our grid
+    static const int dim = GridType::dimension;
+
+    //! \brief A basic number
+    typedef typename GridType::LeafGridView::ctype ctype;
+
+    //! \brief A vectorial node value
+    typedef Dune::FieldVector<double,dim> NodeValue;
+
+    //! \brief A global coordinate
+    typedef typename GridType::LeafGridView::template Codim<1>::Geometry::GlobalCoordinate GlobalCoordinate;
+
+    //! \brief A set of indices
+    typedef typename GridType::LeafGridView::IndexSet LeafIndexSet;
+
+    //! \brief An iterator over grid cells
+    typedef typename GridType::LeafGridView::template Codim<0>::Iterator LeafIterator;
+
+    //! \brief The linear operator
+    ASMHandler<GridType> A;
+
+    //! \brief The solution vectors
+    Vector u[6];
+    //! \brief The load vectors
+    Vector b[6];
+
+    //! \brief Vector holding the volume fractions for materials
+    std::vector<double> volumeFractions;
+
+    //! \brief Main constructor
+    //! \param[in] gv_ The grid to operate on
+    //! \param[in] tol_ The tolerance to use when deciding whether or not a coordinate falls on a plane/line/point. \sa tol
+    //! \param[in] Escale_ A scale value for E-moduluses to avoid numerical issues
+    //! \param[in] file The eclipse grid file
+    //! \param[in] rocklist If true, file is a rocklist
+    ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_, 
+                      const std::string& file, const std::string& rocklist)
+      : gv(gv_), tol(tol_), A(gv_), E(gv_), Escale(Escale_)
+    {
+      if (rocklist.empty())
+        loadMaterialsFromGrid(file);
+      else
+        loadMaterialsFromRocklist(file,rocklist);
+#if HAVE_SUPERLU
+      slu = 0;
+#endif
+      cgsolver = 0;
+      op = 0;
+      ilu = 0;
+      ovl = 0;
+    }
+
+    //! \brief The destructor
+    ~ElasticityUpscale()
+    {
+      // sort the pointers so unique can do its job
+      std::sort(materials.begin(),materials.end());
+      // this reorders the vector so we only get one entry per pointer
+      std::vector<Material*>::iterator itend = std::unique(materials.begin(),materials.end());
+      // now delete the pointers
+      for (std::vector<Material*>::iterator it  = materials.begin();
+                                            it != itend; ++it)
+        delete *it;
+
+#if HAVE_SUPERLU
+      delete slu;
+#endif
+      delete cgsolver;
+      delete ilu;
+      delete op;
+      delete ovl;
+    }
+
+    //! \brief Find boundary coordinates
+    //! \param[out] min The miminum coordinates of the grid
+    //! \param[out] max The maximum coordinates of the grid
+    void findBoundaries(double* min, double* max);
+
+    //! \brief Add a MPC equation
+    //! \param[in] dir The direction of the MPC
+    //! \param[in] slave The slave node index
+    //! \param[in] m The vertices on the master grid
+    void addMPC(Direction dir, int slave, 
+                const BoundaryGrid::Vertex& m);
+
+    //! \brief Establish periodic boundaries using the MPC approach
+    //! \param[in] min The minimum coordinates of the grid
+    //! \param[in] max The maximum coordinates of the grid
+    void periodicBCs(const double* min, const double* max);
+
+    //! \brief Establish periodic boundaries using the LLM approach
+    //! \param[in] min The minimum coordinates of the grid
+    //! \param[in] max The maximum coordinates of the grid
+    //! \param[in] n1 The number of elements on the interface grid in the X direction
+    //! \param[in] n2 The number of elements on the interface grid in the Y direction
+    void periodicBCsLLM(const double* min,
+                        const double* max,
+                        int n1, int n2);
+
+    //! \brief Establish periodic boundaries using the mortar approach
+    //! \param[in] min The minimum coordinates of the grid
+    //! \param[in] max The maximum coordinates of the grid
+    //! \param[in] n1 The number of elements on the lambda grid in the X direction
+    //! \param[in] n2 The number of elements on the lambda grid in the Y direction
+    void periodicBCsMortar(const double* min,
+                           const double* max, int n1, int n2);
+
+    //! \brief Assemble (optionally) stiffness matrix A and load vector
+    //! \param[in] loadcase The strain load case. Set to -1 to skip
+    //! \param[in] matrix Whether or not to assemble the matrix
+    void assemble(int loadcase, bool matrix);
+
+    //! \brief Calculate the average stress vector for the given loadcase
+    //! \param[out] sigma The stress vector
+    //! \param[in] u The displacement vector
+    //! \param[in] loadcase The strain load case considered
+      template <int comp>
+    void averageStress(Dune::FieldVector<ctype,comp>& sigma,
+                       const Vector& u, int loadcase);
+
+    //! \brief Solve Au = b for u
+    //! \param[in] solver The linear equation solver to employ
+    //! \param[in] tol The tolerance for iterative solvers
+    void solve(Solver solver, double tol, int loadcase);
+
+    void setupSolvers(Solver solver);
+  private:
+    //! \brief An iterator over grid vertices
+    typedef typename GridType::LeafGridView::template Codim<dim>::Iterator LeafVertexIterator;
+    
+    //! \brief A reference to our grid
+    const GridType& gv;
+
+    //! \brief Tolerance used to decide whether or not a coordinate falls on a plane/line/point.
+    ctype tol;
+
+    //! \brief Minimum E-modulus (scaling factor)
+    ctype Escale;
+
+    //! \brief Minimum real E for materials
+    ctype Emin;
+
+    //! \brief Check if the given coordinate falls on a given plane
+    //! \param[in] plane The plane of interest
+    //! \param[in] coord The coordinates to check
+    //! \param[in] value The constant coordinate describing the plane
+    bool isOnPlane(Direction plane, GlobalCoordinate coord, ctype value);
+
+    //! \brief Check if the given coordinate falls on a given line 
+    //! \param[in] dir The line direction 
+    //! \param[in] coord The coordinates to check
+    //! \param[in] x The first coordinate of the line
+    //! \param[in] y The second coordinate of the line
+    bool isOnLine(Direction dir, GlobalCoordinate coord, ctype x, ctype y);
+
+    //! \brief Check if the given coordinate is the given point
+    //! \param[in] coord The coordinates to check
+    //! \param[in] point The point
+    bool isOnPoint(GlobalCoordinate coord, GlobalCoordinate point);
+
+    //! \brief Vector holding material parameters for each active grid cell
+    std::vector<Material*> materials;
+
+    //! \brief Fix corner nodes
+    //! \param[in] min The minimum coordinates on the grid
+    //! \param[in] max The maximum coordinates on the grid
+    void fixCorners(const double* min, const double* max);
+
+    //! \brief Extract the vertices on a given face
+    //! \param[in] dir The direction of the face normal
+    //! \param[in] coord The coordinate of the face plane
+    //! \returns A vector holding the matching vertices
+    std::vector<BoundaryGrid::Vertex> extractFace(Direction dir, ctype coord);
+
+    //! \brief An enumeration used to indicate which side to extract from a cube
+    enum SIDE {
+      LEFT,
+      RIGHT
+    };
+
+    //! \brief Extract a quad grid over a given face
+    //! \param[in] dir The direction of the face normal
+    //! \param[in] coord the coordinate of the face plance
+    //! \param[in] side Extract left or right side
+    //! \param[in] dc If true, order vertices in dune convention
+    //! \returns A quad grid spanning the face
+    BoundaryGrid extractMasterFace(Direction dir, ctype coord,
+                                   SIDE side=LEFT, bool dc=false);
+
+    //! \brief Find and establish master/slave grids (MPC)
+    //! \param[in] min The minimum coordinates of the grid
+    //! \param[in] max The maximum coordinates of the grid
+    void determineSideFaces(const double* min, const double* max);
+
+    //! \brief Establish periodic boundary conditions on a given plane using MPC couplings
+    //! \param[in] plane The direction of the plane normal
+    //! \param[in] dir The coordinate direction to enforce periodicity in
+    //! \param[in] slave The slave point grid
+    //! \param[in] master The master quad grid
+    void periodicPlane(Direction plane, Direction dir,
+                       const std::vector<BoundaryGrid::Vertex>& slave,
+                       const BoundaryGrid& master);
+
+    //! \brief Fix the DOFs in a given point on the grid
+    //! \param[in] dir The coordinate direction to fix in
+    //! \param[in] coord The coordinates of the node to fix
+    //! \param[in] value The values to fix the given DOFs to
+    void fixPoint(Direction dir, GlobalCoordinate coord,
+                  const NodeValue& value = NodeValue(0));
+
+    //! \brief Fix the DOFs in a given line on the grid
+    //! \param[in] dir The coordinate direction to fix in
+    //! \param[in] x The first coordinate of the line
+    //! \param[in] y The second coordinate of the line
+    //! \param[in] value The values to fix the given DOFs to
+    void fixLine(Direction dir, ctype x, ctype y,
+                 const NodeValue& value = NodeValue(0));
+
+    //! \brief A point grid
+    typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
+
+    //! \brief Find the B matrix associated with a particular slave grid (LLM)
+    //! \param[in] slave The slave grid to extract
+    //! \returns The B matrix associated with the given slave grid
+    Matrix findBMatrixLLM(const SlaveGrid& slave);
+
+    //! \brief Find the L matrix associated with a particular slave grid (LLM)
+    //! \param[in] slave The slave grid we want to couple
+    //! \param[in] master The interface grid we want to couple to
+    //! \returns The L matrix describing the coupling between slave and master 
+    Matrix findLMatrixLLM(const SlaveGrid& slave, const BoundaryGrid& master);
+
+    //! \brief Find the L matrix associated with mortar couplings
+    //! \param[in] b1 The primal boundary to match
+    //! \param[in] interface The interface/multiplier grid
+    //! \param[in] dir The normal direction on the boundary (0/1)
+    Matrix findLMatrixMortar(const BoundaryGrid& b1,
+                             const BoundaryGrid& interface, int dir);
+
+    //! \brief This function loads and maps materials to active grid cells
+    //! \param[in] file The eclipse grid to read materials from
+    void loadMaterialsFromGrid(const std::string& file);
+
+    //! \brief This function loads and maps materials to active grid cells
+    //! \param[in] file The grid file to read SATNUM from
+    //! \param[in] rocklist The rock list to read materials from
+    void loadMaterialsFromRocklist(const std::string& file,
+                                   const std::string& rocklist);
+
+    //! \brief Setup an overlapping schwarz preconditioner with one 
+    //!        subdomain per pillar.
+    void setupPreconditioner();
+
+    //! \brief Master grids
+    std::vector<BoundaryGrid> master;
+
+    //! \brief Slave point grids
+    std::vector< std::vector<BoundaryGrid::Vertex> > slave;
+
+    //! \brief Vector of matrices holding boolean coupling matrices (LLM)
+    std::vector<Matrix> B;
+    
+    //! \brief Vector of matrices holding Lagrangian multipliers
+    std::vector<Matrix> L;
+
+#if HAVE_SUPERLU
+    //! \brief SuperLU solver
+    Dune::SuperLU<Matrix>* slu;
+#endif
+    //! \brief CG solver
+    Dune::CGSolver<Vector>* cgsolver;
+    //! \brief The linear operator
+    Dune::MatrixAdapter<Matrix,Vector,Vector>* op;
+
+    //! \brief ILU preconditioner
+    Dune::SeqILU0<Matrix,Vector,Vector>* ilu;
+    //! \brief Overlapping Schwarz preconditioner
+    OverlappingSchwarz* ovl;
+
+    //! \brief Elasticity helper class
+    Elasticity<GridType> E;
+};
+
+#include "elasticity_upscale_impl.hpp"

dune/elasticity/material.cpp

@@ -0,0 +1,137 @@
+//==============================================================================
+//!
+//! \file material.cpp
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Material interface.
+//!
+//==============================================================================
+
+#include "material.hh"
+#include "materials.hh"
+#include <fstream>
+
+
+Material* Material::create (int ID, const Dune::DynamicVector<double>& params)
+{
+  switch (params.size())
+    {
+      // Isotropic material, E and nu given
+    case 1: return new Isotropic(ID,params[0],double(0));
+    case 2: return new Isotropic(ID,params[0],params[1]);
+    case 3: return new Isotropic(ID,params[0],params[1],params[2]);
+
+      // Diagonal-orthotropic material, Ex, Ey, Ez, Gxy, Gxz, Gyz given
+    case 4: return new OrthotropicD(ID,
+                                    params[0],params[1],params[2],
+                                    params[3]);
+    case 5: return new OrthotropicD(ID,
+                                    params[0],params[1],params[2],
+                                    params[3],params[4]);
+    case 6: return new OrthotropicD(ID,
+                                    params[0],params[1],params[2],
+                                    params[3],params[4],params[5]);
+
+      // General symmetric orthotropic material, upper triangle of C given
+    case 9:
+      {
+        Dune::DynamicVector<double>mat(21,double(0));
+	mat[0]  = params[0]; mat[1] = params[1]; mat[2] = params[2];
+	mat[6]  = params[3]; mat[7] = params[4];
+	mat[11] = params[5];
+	mat[15] = params[6];
+	mat[18] = params[7];
+	mat[20] = params[8];
+	return new OrthotropicSym(ID,mat);
+      }
+      break;
+
+    case 21: return new OrthotropicSym(ID,params);
+    }
+
+  std::cerr <<"Material::create: Invalid number of material parameters, "
+            << params.size() << std::endl;
+  return 0;
+}
+
+Material* Material::create(int ID, const std::string& file)
+{
+  std::ifstream f;
+  f.open(file.c_str());
+  if (!f.good()) {
+    std::cerr << "Rock file " << file << " not found, using default isotropic material" << std::endl;
+    return new Isotropic(ID,100,0.38,3);
+  }
+  std::string str,str2;
+  f >> str;
+  if (str == "ti") { // transverse isotropic material with unit axis in z-direction. 
+      Dune::DynamicVector<double>mat(21,double(0));
+      double a,b,c,d,e,rho;
+      f >> a >> b >> c >> d >> e;
+      f >> str2;
+      if (str2 != "density") {
+          std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"ti\", " << std::endl 
+                    << "the next line need 5 a,b,c,d and e for C matrix with C11=C22=a, C33=c, " << std::endl
+                    << "C44=C55=d, C66=e, C12=a-2e and C13=b. Then keyword density on next " << std::endl
+                    << "line followed by material density (double) on last line." << std::endl;
+          exit(1);
+          assert(0);
+      }
+      f >> rho;
+      mat[0] = a; mat[1] = a-2.f*e; mat[2] = b; 
+      mat[6] = a; mat[7] = b; mat[11] = c;
+      mat[15] = d; mat[18] = d; mat[20] = e;
+      return new  OrthotropicSym(ID,mat);
+  }
+  else if (str == "anisotropic") { // full symmetric matrix with 21 elements (upper triangle, indexed c11 c12 ... c16 c22 c23 ... c66)
+      Dune::DynamicVector<double>mat(21,double(0));
+      double rho;
+      for (int cuidx=0; cuidx<21; cuidx++) {
+          f >> mat[cuidx];
+      } 
+      f >> str2;
+      if (str2 != "density") {
+          std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"anisotropic\", " << std::endl 
+                    << "the next line need 21 entries for upper C matrix. Then keyword density on next " << std::endl
+                    << "line followed by material density (double) on last line." << std::endl;
+          exit(1);
+          assert(0);
+      }
+      
+      f >> rho;
+      return new  OrthotropicSym(ID,mat);
+  }
+  else {
+      double p1, p2, rho, E, nu;
+      f >> p1 >> p2;
+      f >> str2;
+      if (str2 != "density") {
+          std::cerr << "Rock file " << file << " has wrong format. when isotropycase is either \"km\", " << std::endl
+                    << "\"lm\", \"en\" or \"vpvs\", the next line need 21 entries for upper C matrix. " << std::endl
+                    << "Then keyword density on next line followed by material density (double) on last line." << std::endl;
+          exit(1);
+          assert(0);
+      }
+      f >> rho;
+      if (str == "vpvs") {
+          p1 = rho*p1*p1-4.f/3*rho*p2*p2;
+          p2 = rho*p2*p2;
+          str = "km";
+      }
+      if (str == "km") { // bulk modulus and shear modulus
+          nu = (3*p1-2*p2)/(6*p1+2*p2);
+          E = 2*p2*(1+nu);
+      } else if (str == "lm") { // lame's lambda and shear modulus
+          nu = p1/(2*(p2+p1));
+          E = (1+nu)*(1-2*nu)*p1/nu;
+      } else if (str == "en") { // young's modulus and poisson's ratio
+          E = p1;
+          nu = p2;
+      }
+      
+      return new Isotropic(ID,E,nu,rho);
+  }
+}

dune/elasticity/material.hh

@@ -0,0 +1,93 @@
+//==============================================================================
+//!
+//! \file material.hh
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Material interface.
+//!
+//==============================================================================
+#pragma once
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+
+/*!
+  \brief This is a base class for linear elastic materials.
+  \details It is an abstract class since some of the member functions are
+  purely virtual.
+*/
+
+class Material
+{
+protected:
+  //! \brief Default constructor creating an empty material.
+  Material(int ID = 0, double density = 0.0)
+    : id(ID), rho(density)
+  {
+  }
+
+  //! \brief Prints the material properties to a stream.
+  virtual std::ostream& write(std::ostream& os) const
+  { 
+    return os;
+  }
+public:
+  //! \brief Empty virtual destructor.
+  virtual ~Material() {}
+
+  //! \brief Returns the external material id.
+  int num() const
+  {
+    return id;
+  }
+
+  //! \brief Returns the number of parameters describing this material.
+  virtual int numPar() const = 0;
+  //! \brief Returns the \a ipar'th parameter describing this material.
+  virtual double getPar(int ipar = 1) const
+  {
+    return double(0);
+  }
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,6,6>& C,
+				     bool invers = false) const = 0;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,3,3>& C,
+				     bool invers = false) const = 0;
+
+  //! \brief Returns the mass density of this material.
+  double getMassDensity() const
+  {
+    return rho;
+  }
+
+  //! \brief Global stream operator printing a material properties object.
+  friend std::ostream& operator<<(std::ostream& os, const Material& m)
+  {
+    return m.write(os);
+  }
+
+  //! \brief Creates a material object of a given type.
+  //! \details The material type depends on the number of parameters provided.
+  //! \param[in] ID External number for this material
+  //! \param[in] params Array of material parameters
+  static Material* create(int ID, const Dune::DynamicVector<double>& params);
+
+  //! \brief Creates a material object from a rocklist
+  //! \param[in] ID ID of the material
+  //! \param[in] file The URL to the rocklist
+  static Material* create(int ID, const std::string& file);
+private:
+  int  id;  //!< External material number
+  double rho; //!< Mass density
+};

dune/elasticity/materials.cpp

@@ -0,0 +1,216 @@
+//==============================================================================
+//!
+//! \file materials.cpp
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Material properties.
+//!
+//==============================================================================
+
+#include "materials.hh"
+#include <string.h>
+
+
+std::ostream& Isotropic::write (std::ostream& os) const
+{
+  os <<"Isotropic material "<< this->num()
+     <<": E = "<< E <<" nu = "<< nu;
+  return os << std::endl;
+}
+
+
+/*!
+  For edim = 3: \f[
+  [C] = \frac{E}{(1+\nu)(1-2\nu)} \left[\begin{array}{cccccc}
+  1-\nu & \nu & \nu \\
+  \nu & 1-\nu & \nu \\
+  \nu & \nu & 1-\nu \\
+  & & & \frac{1}{2}-\nu \\
+  & & & & \frac{1}{2}-\nu \\
+  & & & & & \frac{1}{2}-\nu
+  \end{array}\right] \f]
+
+  For edim = 2 (plain strain): \f[
+  [C] = \frac{E}{(1+\nu)(1-2\nu)} \left[\begin{array}{ccc}
+  1-\nu & \nu \\
+  \nu & 1-\nu \\
+  & & \frac{1}{2}-\nu
+  \end{array}\right] \f]
+*/
+
+bool Isotropic::getConstitutiveMatrix (Dune::FieldMatrix<double,6,6>& C,
+                                       bool invers) const
+{
+  C = 0;
+  const double one  = 1.f;
+  const double two  = 2.f;
+  const double half = 0.5f;
+  if (nu < 0 || nu >= half) return false;
+
+  const double fact = E / ((one + nu) * (one - nu - nu));
+
+  C[0][0] = invers ? one / E : (one - nu) * fact;
+  C[1][0] = invers ? -nu / E : nu * fact;
+  C[2][0] = C[1][0];
+
+  C[0][1] = C[1][0];
+  C[1][1] = C[0][0];
+  C[2][1] = C[1][0];
+
+  C[0][2] = C[1][0];
+  C[1][2] = C[1][0];
+  C[2][2] = C[0][0];
+
+  C[3][3] = invers ? (two + nu + nu) / E : (half - nu) * fact;
+  C[4][4] = C[3][3];
+  C[5][5] = C[3][3];
+
+  return true;
+}
+
+bool Isotropic::getConstitutiveMatrix (Dune::FieldMatrix<double,3,3>& C,
+                                       bool invers) const
+{
+  const double one  = 1.f;
+  const double two  = 2.f;
+  const double half = 0.5f;
+  if (nu < 0 || nu >= half) return false;
+
+  const double fact = E / ((one + nu) * (one - nu - nu));
+
+  C[0][0] = invers ? (one - nu*nu) / E : (one - nu) * fact;
+  C[1][0] = invers ? (-nu - nu*nu) / E : nu * fact;
+
+  C[0][1] = C[1][0];
+  C[1][1] = C[0][0];
+
+  C[2][2] = invers ? (two + nu + nu) / E : (half - nu) * fact;
+
+  return true;
+}
+
+
+OrthotropicD::OrthotropicD (int ID, double Ex, double Ey, double Ez,
+                            double Gxy, double Gxz, double Gyz) : Material(ID)
+{
+  E[0] = Ex;
+  E[1] = Ey;
+  E[2] = Ez;
+  E[3] = Gxy;
+  E[4] = Gxz > double(0) ? Gxz : E[3];
+  E[5] = Gyz > double(0) ? Gyz : E[4];
+}
+
+
+std::ostream& OrthotropicD::write (std::ostream& os) const
+{
+  os <<"Orthotropic material "<< this->num()
+     <<": Ex = "<< E[0] <<" Ey = "<< E[1] <<" Ez = "<< E[2];
+  if (E[4] == E[3] && E[5] == E[3])
+    os <<" G = "<< E[3];
+  else
+    os <<" Gxy = "<< E[3] <<" Gxz = "<< E[4] <<" Gyz = "<< E[5];
+  return os << std::endl;
+}
+
+
+double OrthotropicD::getPar (int ipar) const
+{
+  if (ipar > 0 && ipar <= 6)
+    return E[ipar-1];
+  else
+    return double(0);
+}
+
+
+/*!
+  \f[ [C] = \left[\begin{array}{cccccc}
+  E_x \\ & E_y \\ & & E_z \\
+  & & & G_{xy} \\ & & & & G_{xz} \\ & & & & & G_{yz}
+  \end{array}\right] \f]
+*/
+
+bool OrthotropicD::getConstitutiveMatrix (Dune::FieldMatrix<double,6,6>& C,
+                                          bool invers) const
+{
+  for (size_t i = 0; i < 6; i++)
+    C[i][i] = invers ? double(1)/E[i] : E[i];
+
+  return true;
+}
+
+bool OrthotropicD::getConstitutiveMatrix (Dune::FieldMatrix<double,3,3>& C,
+                                          bool invers) const
+{
+  return false;
+}
+
+
+OrthotropicSym::OrthotropicSym (int ID,
+                        const Dune::DynamicVector<double>& Cu) : Material(ID)
+{
+  // The input matrix is assumed with Voigt ordering (xx,yy,zz,yz,zx,xy)
+  // Have to reorder to the internal order required  (xx,yy,zz,xy,xz,yz)
+  memcpy(Cupper,&Cu[0],21*sizeof(double));
+  std::swap(Cupper[ 3],Cupper[ 5]);
+  std::swap(Cupper[ 8],Cupper[10]);
+  std::swap(Cupper[12],Cupper[14]);
+  std::swap(Cupper[15],Cupper[20]);
+  std::swap(Cupper[16],Cupper[19]);
+}
+
+
+std::ostream& OrthotropicSym::write (std::ostream& os) const
+{
+  //TODO: Print matrix
+  return os;
+}
+
+
+double OrthotropicSym::getPar (int ipar) const
+{
+  if (ipar > 0 && ipar <= 21)
+    return Cupper[ipar-1];
+  else
+    return double(0);
+}
+
+
+/*!
+  \f[ [C] = \left[\begin{array}{cccccc}
+    C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16} \\
+           & C_{22} & C_{23} & C_{24} & C_{25} & C_{26} \\
+           &        & C_{33} & C_{34} & C_{35} & C_{36} \\
+           &        &        & C_{44} & C_{45} & C_{46} \\
+    \multicolumn{3}{c}{\rm symm.} &   & C_{55} & C_{56} \\
+           &        &        &        &        & C_{66}
+  \end{array}\right] \f]
+*/
+
+bool OrthotropicSym::getConstitutiveMatrix (Dune::FieldMatrix<double,6,6>& C,
+                                            bool invers) const
+{
+  size_t i, j, k = 0;
+  for (i = 0; i < 6; i++)
+  {
+    C[i][i] = Cupper[k++];
+    for (j = i+1; j < 6; j++)
+      C[i][j] = C[j][i] = Cupper[k++];
+  }
+
+  if (invers) {
+    //TODO!
+    assert(0);
+  }
+
+  return true;
+}
+
+bool OrthotropicSym::getConstitutiveMatrix (Dune::FieldMatrix<double,3,3>& C,
+                                            bool invers) const
+{
+  return false;
+}

dune/elasticity/materials.hh

@@ -0,0 +1,169 @@
+//==============================================================================
+//!
+//! \file materials.hh
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Material properties.
+//!
+//==============================================================================
+#pragma once
+
+#include "material.hh"
+
+
+/*!
+  \brief Isotropic linear elastic material.
+*/
+
+class Isotropic : public Material
+{
+public:
+  //! \brief Constructor creating a new isotropic material.
+  //! \param[in] ID      External material number
+  //! \param[in] Emod    Young's modulus
+  //! \param[in] Poisson Poisson's ratio
+  //! \param[in] rho     Mass density
+  Isotropic(int ID, double Emod, double Poisson, double rho = 0.0)
+    : Material(ID,rho)
+  {
+    E = Emod;
+    nu = Poisson;
+  }
+
+  //! \brief Empty virtual destructor.
+  virtual ~Isotropic() {}
+
+  //! \brief Returns the number of parameters describing this material.
+  virtual int numPar() const
+  {
+    return 2;
+  }
+
+  //! \brief Returns the \a ipar'th parameter describing this material.
+  virtual double getPar(int ipar = 1) const
+  {
+    return ipar == 1 ? E : nu;
+  }
+
+  //! \brief Set the E modulus of the material
+  //! \param[in] E_ The E modulus to be set
+  void setE(double E_)
+  {
+    E = E_;
+  }
+
+  //! \brief Returns the E modulus of the material
+  double getE() const
+  {
+    return E;
+  }
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,6,6>& C,
+                                     bool invers = false) const;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,3,3>& C,
+                                     bool invers = false) const;
+protected:
+  //! \brief Prints the material properties to a stream.
+  virtual std::ostream& write(std::ostream& os) const;
+
+private:
+  double E;  //!< Young's modulus
+  double nu; //!< Poisson's ratio
+};
+
+//! \brief Orthotropic linear elastic material with diagonal constitutive matrix.
+class OrthotropicD : public Material
+{
+public:
+  //! \brief Constructor creating a new material.
+  //! \param[in] ID External material number
+  //! \param[in] Ex Elasticity modulus in local x-direction
+  //! \param[in] Ey Elasticity modulus in local y-direction
+  //! \param[in] Ez Elasticity modulus in local z-direction
+  //! \param[in] Gxy Shear modulus in the local xy-plane
+  //! \param[in] Gxz Shear modulus in the local xz-plane, default = Gxy
+  //! \param[in] Gyz Shear modulus in the local yz-plane, default = Gxz
+  OrthotropicD(int ID, double Ex, double Ey, double Ez,
+               double Gxy, double Gxz = double(-1), double Gyz = double(-1));
+
+  //! \brief Empty virtual destructor.
+  virtual ~OrthotropicD() {}
+
+  //! \brief Returns the number of parameters describing this material.
+  virtual int numPar() const
+  {
+    return 6;
+  }
+
+  //! \brief Returns the \a ipar'th parameter describing this material.
+  virtual double getPar(int ipar = 1) const;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,6,6>& C,
+				     bool invers = false) const;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,3,3>& C,
+				     bool invers = false) const;
+
+protected:
+  //! \brief Prints the material properties to a stream.
+  virtual std::ostream& write(std::ostream& os) const;
+
+private:
+  double E[6]; //!< The diagonal of the constitutive matrix
+};
+
+//! \brief Orthotropic linear elastic material with symmetric constitutive matrix.
+class OrthotropicSym : public Material
+{
+public:
+  //! \brief Constructor creating a new material.
+  //! \param[in] ID External material number
+  //! \param[in] Cu Upper triangle of the symmetric constitutive matrix
+  OrthotropicSym(int ID, const Dune::DynamicVector<double>& Cu);
+
+  //! \brief Empty virtual destructor.
+  virtual ~OrthotropicSym() {}
+
+  //! \brief Returns the number of parameters describing this material.
+  virtual int numPar() const
+  {
+    return 21;
+  }
+
+  //! \brief Returns the \a ipar'th parameter describing this material.
+  virtual double getPar(int ipar = 1) const;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,6,6>& C,
+				     bool invers = false) const;
+
+  //! \brief Establishes the full constitutive matrix for this material.
+  //! \param[out] C The constitutive matrix
+  //! \param[in] invers If \e true, set up the inverse matrix instead
+  virtual bool getConstitutiveMatrix(Dune::FieldMatrix<double,3,3>& C,
+				     bool invers = false) const;
+protected:
+  //! \brief Prints the material properties to a stream.
+  virtual std::ostream& write(std::ostream& os) const;
+
+private:
+  double Cupper[21]; //!< Upper triangle of the symmetric constitutive matrix
+};

dune/elasticity/matrixops.hpp

@@ -0,0 +1,54 @@
+//==============================================================================
+//!
+//! \file matrixops.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Helper class with some matrix operations
+//!
+//==============================================================================
+#pragma once
+
+#include <dune/istl/bcrsmatrix.hh>
+
+//! \brief A sparse matrix holding our operator
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > Matrix;
+
+//! \brief For storing matrix adjacency/sparsity patterns
+typedef std::vector< std::set<int> > AdjacencyPattern;
+
+//! \brief Helper class with some matrix operations
+class MatrixOps {
+  public:
+    //! \brief Create a sparse matrix from a given adjacency pattern
+    //! \param[in] adj The adjacency pattern
+    //! \param[in] rows The number of rows in the matrix
+    //! \param[in] cols The number of columns in the matrix
+    //! \param[out] A The created matrix
+    static void fromAdjacency(Matrix& A, const AdjacencyPattern& adj,
+                              int rows, int cols);
+
+    //! \brief Print a matrix to stdout
+    //! \param[in] A The matrix to print
+    static void print(const Matrix& A);
+
+    //! \brief axpy like operation - returns A+alpha*B
+    //! \param[in] A The matrix to subtract from
+    //! \param[in] B The matrix to subtract
+    //! \param[in] alpha The constant in front of B
+    //! \returns A+alpha*B
+    static Matrix Axpy(const Matrix& A, const Matrix& B, double alpha);
+
+    //! \brief Augment a matrix with another
+    //! \param[in] A The matrix to be augmented
+    //! \param[in] B The matrix to augment with
+    //! \param[in] r0 The starting row of the augment matrix
+    //! \param[in] c0 The starting column of the augment matrix
+    //! \param[in] symmetric If true, augment symmetrically
+    static Matrix augment(const Matrix& A, const Matrix& B,
+                          size_t r0, size_t c0, bool symmetric);
+};
+
+#include "matrixops_impl.hpp"

dune/elasticity/matrixops_impl.hpp

@@ -0,0 +1,146 @@
+//==============================================================================
+//!
+//! \file matrixops_impl.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Helper class with some matrix operations - template implementations
+//!
+//==============================================================================
+
+void MatrixOps::fromAdjacency(Matrix& A, const std::vector< std::set<int> >& adj,
+                              int rows, int cols)
+{  
+  size_t sum=0;
+  for (size_t i=0;i<adj.size();++i)
+    sum += adj[i].size();
+  A.setSize(rows, cols, sum);
+  A.setBuildMode(Matrix::random);
+
+  for (int i = 0; i < adj.size(); i++)
+    A.setrowsize(i,adj[i].size());
+  A.endrowsizes();
+
+  for (size_t i = 0; i < adj.size(); i++) {
+    std::set<int>::iterator setend = adj[i].end();
+    for (std::set<int>::iterator setit = adj[i].begin();
+        setit != setend; ++setit) {
+      A.addindex(i,*setit);
+    }
+  }
+  A.endindices();
+  A = 0;
+}
+
+void MatrixOps::print(const Matrix& A)
+{
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end(); ++it) {
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2) {
+      double val = *it2; 
+      if (fabs(val) < 1.e-14)
+        continue;
+      std::cout << it.index() << " " << it2.index() << " : " << val << std::endl;
+    }
+  }
+}
+
+Matrix MatrixOps::Axpy(const Matrix& A, const Matrix& B, double alpha)
+{
+  assert(A.M() == B.M() && A.N() == B.N());
+
+  // establish union adjacency pattern
+  std::vector<std::set<int> > adj;
+  adj.resize(A.N());
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end(); ++it) {
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2)
+      adj[it.index()].insert(it2.index());
+  }
+  for (Matrix::ConstIterator it  = B.begin();
+                             it != B.end(); ++it) {
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2)
+      adj[it.index()].insert(it2.index());
+  }
+  Matrix result;
+  fromAdjacency(result,adj,A.N(),A.M());
+  // now insert elements from A
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end(); ++it) {
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2)
+      result[it.index()][it2.index()] = *it2;
+  }
+  // and subtract elements from B
+  for (Matrix::ConstIterator it  = B.begin();
+                             it != B.end(); ++it) {
+    for (Matrix::ConstColIterator it2  = it->begin();
+                                  it2 != it->end();++it2)
+      result[it.index()][it2.index()] += alpha*(*it2);
+  }
+
+  return result;
+}
+
+Matrix MatrixOps::augment(const Matrix& A, const Matrix& B,
+                     size_t r0, size_t c0, bool symmetric)
+{
+  std::cout << "Augmenting matrix of dimension " << A.N() << "x" << A.M()
+            << " with matrix of dimension " << B.N() << "x" << B.M() << std::endl;
+  size_t nrow = A.N();
+  size_t ncol = A.M();
+  if (r0+B.N() > nrow) nrow = r0+B.N();
+  if (symmetric && r0+B.N() > ncol) ncol = r0+B.N();
+  if (c0+B.M() > ncol) ncol = c0+B.M();
+  if (symmetric && c0+B.M() > nrow) nrow = c0+B.M();
+  std::cout << "Resulting size: " << nrow << "x" << ncol << std::endl;
+
+  AdjacencyPattern adj;
+  adj.resize(nrow);
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end();++it) {
+    for (Matrix::ConstColIterator it2  = it->begin(); 
+                                  it2 != it->end();++it2) {
+      adj[it.index()].insert(it2.index());
+    }
+  }
+  for (Matrix::ConstIterator it  = B.begin();
+                             it != B.end();++it) {
+    for (Matrix::ConstColIterator it2  = it->begin(); 
+                                  it2 != it->end();++it2) {
+      adj[it.index()+r0].insert(it2.index()+c0);
+      if (symmetric)
+        adj[it2.index()+c0].insert(it.index()+r0);
+    }
+  }
+  if (symmetric) {
+    // always establish diagonal elements or superLU crashes
+    for (int i=0;i<nrow;++i)
+      adj[i].insert(i);
+  }
+  Matrix result;
+  fromAdjacency(result,adj,nrow,ncol);
+  for (Matrix::ConstIterator it  = A.begin();
+                             it != A.end();++it) {
+    for (Matrix::ConstColIterator it2  = it->begin(); 
+                                  it2 != it->end();++it2) {
+      result[it.index()][it2.index()] = *it2;
+    }
+  }
+  for (Matrix::ConstIterator it  = B.begin();
+                             it != B.end();++it) {
+    for (Matrix::ConstColIterator it2  = it->begin(); 
+                                  it2 != it->end();++it2) {
+      result[it.index()+r0][it2.index()+c0] = *it2;
+      if (symmetric)
+        result[it2.index()+c0][it.index()+r0] = *it2;
+    }
+  }
+  
+  return result;
+}

dune/elasticity/mpc.cpp

@@ -0,0 +1,79 @@
+//==============================================================================
+//!
+//! \file mpc.cpp
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Representation of multi-point constraint (MPC) equations.
+//!
+//==============================================================================
+
+#include "mpc.hh"
+
+bool MPC::Less::compareSlaveDofOnly = false;
+
+
+//! \brief Global operator for comparing two MPC::DOF objects.
+
+bool operator< (const MPC::DOF& lhs, const MPC::DOF& rhs)
+{
+  if (lhs.node < rhs.node) return true;
+  if (lhs.node > rhs.node) return false;
+  if (lhs.dof  < rhs.dof)  return true;
+  if (lhs.dof  > rhs.dof)  return false;
+
+  if (MPC::Less::compareSlaveDofOnly)
+    return false; // ignore coefficient differences, if any
+  else
+    return lhs.coeff < rhs.coeff ? true : false;
+}
+
+
+bool MPC::Less::operator() (const MPC* lhs, const MPC* rhs) const
+{
+  if (!rhs) return false;
+  if (!lhs) return true;
+
+  if (lhs->getSlave() < rhs->getSlave()) return true;
+  if (rhs->getSlave() < lhs->getSlave()) return false;
+
+  if (compareSlaveDofOnly) return false; // ignore master differences, if any
+
+  size_t lMaster = lhs->getNoMaster();
+  size_t rMaster = rhs->getNoMaster();
+  for (size_t i = 0; i < lMaster && i < rMaster; i++)
+  {
+    if (lhs->getMaster(i) < rhs->getMaster(i)) return true;
+    if (rhs->getMaster(i) < lhs->getMaster(i)) return false;
+  }
+
+  return lMaster < rMaster ? true : false;
+}
+
+
+//! \brief Global stream operator printing a constraint equation.
+
+std::ostream& operator<< (std::ostream& s, const MPC& mpc)
+{
+  s <<"Slave "<< mpc.slave.node <<","<< mpc.slave.dof;
+  if (mpc.slave.coeff != double(0))
+    s <<" = "<< mpc.slave.coeff;
+  else if (mpc.master.empty())
+    return s <<" = 0"<< std::endl;
+  else
+    s <<" = ";
+
+  for (size_t i = 0; i < mpc.master.size(); i++)
+  {
+    if (i == 0 && mpc.slave.coeff == double(0))
+      s << mpc.master[i].coeff;
+    else if (mpc.master[i].coeff >= double(0))
+      s <<" + "<< mpc.master[i].coeff;
+    else
+      s <<" - "<< -mpc.master[i].coeff;
+    s <<"*("<< mpc.master[i].node <<","<< mpc.master[i].dof <<")";
+  }
+  return s << std::endl;
+}

dune/elasticity/mpc.hh

@@ -0,0 +1,155 @@
+//==============================================================================
+//!
+//! \file mpc.hh
+//!
+//! \date Oct 1 2007
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Representation of multi-point constraint (MPC) equations.
+//!
+//==============================================================================
+#pragma once
+
+#include <map>
+#include <vector>
+#include <iostream>
+#include <set>
+
+//! \brief An enum for specification of global coordinate directions
+enum Direction { NONE = 0, X = 1, Y = 2, Z = 4,
+                 XY = 1+2, XZ = 1+4, YZ = 2+4, 
+                 XYZ = 1+2+4 };
+
+/*!
+  \brief A class for representing a general multi-point constraint equation.
+
+  \details A multi-point constraint (MPC) equation is used to introduce
+  additional coupling between the degrees of freedom (DOF) in a FE grid.
+
+  An MPC equation is in general a linear coupling of one (slave) dof to a set
+  of master dofs, and may be expressed as
+
+  \f[ s = c_0 + c_1 m_1 + c_2 m_2 + \ldots + c_n m_n \f]
+
+  where
+  <ul>
+  <li> \f$s\f$ denotes the slave dof
+  <li> \f$c_0\f$ is the value of the slave dof when all master dofs are zero,
+       or when there are no masters in the equation (prescribed dof).
+  <li> \f$m_i\f$ is a master dof
+  <li> \f$c_i, i > 0\f$ is a constant scaling coefficient
+       associated with the \a i'th master
+  <li> \a n denotes the total number of master dofs in the equation
+  </ul>
+  When \a n = 0, the above equation represents a fixed (\f$c_0=0\f$) or
+  prescribed (\f$c_0\ne0\f$) dof.
+
+  One or more of the master dofs may also be a slave in another constraint
+  equation (chained constrains). In this case the two equations are combined
+  to eliminate the master dof that is constrained. This is done while
+  preprocessing the model by the resolveMPCchains function.
+*/
+
+class MPC
+{
+public:
+
+  /*!
+    \brief A struct for representing one term (DOF number and associated
+    coefficient) in a MPC equation.
+  */
+  struct DOF
+  {
+    //! \brief Default constructor.
+    DOF() : node(0), dof(0), coeff(double(0)) {}
+
+    //! \brief Convenience constructor creating a valid DOF object.
+    //! \param[in] n Node number (1...NNOD)
+    //! \param[in] d The local DOF number (1...3)
+    //! \param[in] c Associated coefficient or constrained value
+    DOF(int n, int d, double c = double(0)) : node(n), dof(d), coeff(c) {}
+
+    //! \brief Global stream operator printing a DOF instance.
+    friend std::ostream& operator<<(std::ostream& s, const DOF& dof)
+    {
+      return s <<"u_"<< char('w'+dof.dof) <<" in node "<< dof.node;
+    }
+
+    int  node;  //!< Node number identifying this DOF
+    int  dof;   //!< Local DOF number within \a node
+    double coeff; //!< The constrained value, or master DOF scaling coefficient
+  };
+
+  //! \brief Comparison predicate for MPCs
+  class Less {
+  public:
+    //! \brief Comparison operator used when inserting an MPC-pointer into a
+    //! \a set<MPC*,MPC::Less> object.
+    bool operator()(const MPC* lhs, const MPC* rhs) const;
+    //! \brief Indicates whether only the slave dof number should affect sorting.
+    //! \details The default is to also compare the associated coefficients.
+    static bool compareSlaveDofOnly;
+  };
+
+  //! \brief Constructor creating a constraint for a specified slave DOF
+  //! with no master DOFs.
+  //! \param[in] n The node number of the slave DOF (1...NNOD)
+  //! \param[in] d The local DOF number of the slave DOF (1...3)
+  //! \param[in] c The actual value that this slave DOF is constrained to
+  //! when there are no master DOFs, or all master DOFs are zero
+  MPC(int n, int d, double c = double(0)) : slave(n,d,c) { iceq = -1; }
+
+  //! \brief Adds a master DOF to the constraint equation.
+  //! \param[in] n The node number of the master DOF (1...NNOD)
+  //! \param[in] d The local DOF number of the master DOF (1...3)
+  //! \param[in] c The coefficient that this master should be scaled with
+  //! \param[in] tol Tolerance for comparison with zero,
+  //! if the coefficient \a c is zero, the master DOF is not added
+  void addMaster(int n, int d, double c = double(1), double tol = double(1.0e-8))
+  {
+    if (c < -tol || c > tol) master.push_back(DOF(n,d,c));
+  }
+
+  //! \brief Updates the coefficient of the \a pos'th master DOF.
+  void updateMaster(size_t pos, double c)
+  {
+    if (pos < master.size())
+      master[pos].coeff = c;
+  }
+
+  //! \brief Removes the \a pos'th master DOF from the constraint equation.
+  void removeMaster(size_t pos)
+  {
+    if (pos < master.size())
+      master.erase(master.begin()+pos);
+  }
+
+  //! \brief Increments the \a c0 coefficient by a given \a offset.
+  void addOffset(double offset) { slave.coeff += offset; }
+
+  //! \brief Assigns a new \a c0 coefficient to the constraint equation.
+  void setSlaveCoeff(double c0) { slave.coeff = c0; }
+
+  //! \brief Returns a reference to the slave DOF.
+  const DOF& getSlave() const          { return slave; }
+  //! \brief Returns a reference to the \a i'th master DOF.
+  const DOF& getMaster(size_t i) const { return master[i]; }
+
+  //! \brief Returns the number of master DOFs.
+  size_t getNoMaster() const { return master.size(); }
+
+  //! \brief Global stream operator printing a constraint equation.
+  friend std::ostream& operator<<(std::ostream& s, const MPC& mpc);
+
+  int              iceq;   //!< Global constraint equation identifier
+private:
+  DOF              slave;  //!< The slave DOF of this constraint equation
+  std::vector<DOF> master; //!< The master DOFs of this constraint equation
+};
+
+//! \brief A set of MPCs
+typedef std::set<MPC*,MPC::Less> MPCSet;
+
+//! \brief A mapping from dof to MPCs
+typedef std::map<int,MPC*> MPCMap;

dune/elasticity/shapefunctions.hpp

@@ -0,0 +1,171 @@
+//==============================================================================
+//!
+//! \file shapefunctions.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Class for linear shape functions. Loosely based on code in dune-grid-howto 
+//!
+//==============================================================================
+#pragma once
+
+#include <dune/common/fvector.hh>
+
+//! \brief Represents a linear shape function on a Q4/Q8 element
+  template<class ctype, class rtype, int dim>
+class LinearShapeFunction
+{
+  public:
+    //! \brief The dimension of the shape function
+    enum { dimension = dim };
+
+    //! \brief Default constructor
+    LinearShapeFunction() : coeff0(0.0), coeff1(0.0) {}
+
+    //! \brief Construct a shape function with the given coefficients
+    //! \param[in] coeff0_ The constant coefficients
+    //! \param[in] coeff1_ The linear coefficients
+    LinearShapeFunction(const Dune::FieldVector<rtype,dim>& coeff0_,
+                        const Dune::FieldVector<rtype,dim>& coeff1_)
+      : coeff0(coeff0_), coeff1(coeff1_) {}
+
+    //! \brief Set the given conefficients
+    //! \param[in] coeff0_ The constant coefficients
+    //! \param[in] coeff1_ The linear coefficients
+    void setCoeff(const Dune::FieldVector<rtype,dim>& coeff0_, const Dune::FieldVector<rtype,dim>& coeff1_)
+    {
+      coeff0 = coeff0_;
+      coeff1 = coeff1_;
+    }
+
+    //! \brief Evaluate the shape function
+    //! \param[in] local The local coordinates
+    rtype evaluateFunction(const Dune::FieldVector<ctype,dim>& local) const
+    {
+      rtype result = 1;
+      for (int i = 0; i < dim; ++i)
+        result *= (coeff0[i]+coeff1[i] * local[i]);
+      return result;
+    }
+
+    //! \brief Evaluate the gradient of the shape function
+    //! \param[in] local The local coordinates
+    Dune::FieldVector<rtype,dim>
+    evaluateGradient(const Dune::FieldVector<ctype,dim>& local) const
+    {
+      Dune::FieldVector<rtype, dim> result;
+      for (int i=0;i<dim;++i) {
+        result[i] = 1;
+        for (int j=0;j<dim;++j) {
+          if (i == j)
+            result[i] *= coeff1[j];
+          else
+            result[i] *= (coeff0[j]+coeff1[j]*local[j]);
+        }
+      }
+
+      return result;
+    }
+
+  private:
+    //! \brief The constant coefficients
+    Dune::FieldVector<rtype,dim> coeff0;
+
+    //! \brief The linear coefficients
+    Dune::FieldVector<rtype,dim> coeff1;
+};
+
+//! \brief Singleton handler for the set of LinearShapeFunction
+  template<class ctype, class rtype, int dim>
+class P1ShapeFunctionSet
+{
+public:
+    //! \brief The number of shape functions in the set
+    enum { n = dim==2?4:8 };
+
+    //! \brief A single shape function
+    typedef LinearShapeFunction<ctype,rtype,dim> ShapeFunction;
+
+    //! \brief The type of the return value from a shape function
+    typedef rtype resulttype;
+
+    //! \brief Get the only instance of this class
+    static const P1ShapeFunctionSet& instance()
+    {
+      static const P1ShapeFunctionSet sfs;
+      return sfs;
+    }
+
+    //! \brief Obtain a given shape function
+    //! \param[in] i The requested shape function
+    const ShapeFunction& operator[](int i) const
+    {
+      return f[i];
+    }
+
+private:
+    //! \brief Private constructor prevents additional instances
+    P1ShapeFunctionSet()
+    {
+        static rtype coeffs11[] = {0,
+                                   1};
+        static rtype coeffs12[] = {1,
+                                   -1};
+
+        static rtype coeffs21[] = { 1, 1,
+                                    0, 1,
+                                    1, 0,
+                                    0, 0};
+        static rtype coeffs22[] = {-1,-1, 
+                                    1,-1, 
+                                   -1, 1,
+                                    1, 1};
+
+        static rtype coeffs31[] = { 1, 1, 1,
+                                    0, 1, 1,
+                                    1, 0, 1,
+                                    0, 0, 1,
+                                    1, 1, 0,
+                                    0, 1, 0,
+                                    1, 0, 0,
+                                    0, 0, 0};
+        static rtype coeffs32[] = {-1,-1,-1, 
+                                    1,-1,-1, 
+                                   -1, 1,-1,
+                                    1, 1,-1,
+                                   -1,-1, 1,
+                                    1,-1, 1,
+                                   -1, 1, 1,
+                                    1, 1, 1};
+
+        rtype* coeffs1;
+        rtype* coeffs2;
+        if (dim == 1) {
+          coeffs1 = coeffs11;
+          coeffs2 = coeffs12;
+        }
+        if (dim == 2) {
+          coeffs1 = coeffs21;
+          coeffs2 = coeffs22;
+        }
+        if (dim == 3) {
+          coeffs1 = coeffs31;
+          coeffs2 = coeffs32;
+        }
+
+        Dune::FieldVector<rtype,dim> c1;
+        Dune::FieldVector<rtype,dim> c2;
+        for (int i=0;i<n;++i) {
+          for (int j=0;j<dim;++j) {
+            c1[j] = coeffs1[i*dim+j];
+            c2[j] = coeffs2[i*dim+j];
+          }
+          f[i].setCoeff(c1,c2);
+        }
+    }
+
+    //! \brief Our shape functions
+    ShapeFunction f[n];
+};

examples/Makefile.am

@@ -11,6 +11,7 @@ grdecldips                   			\
 upscale_avg					\
 upscale_cap					\
 upscale_cond					\
+upscale_elasticity				\
 upscale_perm					\
 upscale_relperm					\
 upscale_relpermvisc                             \
@@ -33,6 +34,10 @@ upscale_cap_SOURCES = upscale_cap.cpp
 
 upscale_cond_SOURCES = upscale_cond.cpp
 
+upscale_elasticity_SOURCES = upscale_elasticity.cpp
+upscale_elasticity_CXXFLAGS = ${OPENMP_CXXFLAGS}
+upscale_elasticity_LDADD = $(top_srcdir)/dune/elasticity/libelasticityupscale_noinst.la
+
 upscale_perm_SOURCES = upscale_perm.cpp
 
 upscale_relperm_SOURCES = upscale_relperm.cpp

examples/upscale_elasticity.cpp

@@ -0,0 +1,305 @@
+//==============================================================================
+//!
+//! \file upscale_elasticity.cpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Elasticity upscaling on cornerpoint grids
+//!
+//==============================================================================
+#ifdef HAVE_CONFIG_H
+# include "config.h"     
+#endif
+
+#include <iostream>
+#include <unistd.h>
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#if HAVE_OPENMP
+#include <omp.h>
+#endif
+
+#include <dune/elasticity/elasticity_upscale.hpp>
+
+
+//! \brief Display the available command line parameters
+void syntax(char** argv)
+{
+  std::cerr << "Usage: " << argv[0] << " gridfilename=filename.grdecl [method=]" << std::endl
+            << "\t[xmax=] [ymax=] [zmax=] [xmin=] [ymin=] [zmin=] [linsolver_type=]" << std::endl
+            <<" \t[Emin=] [ctol=] [ltol=] [rock_list=] [vtufilename=] [output=]" << std::endl << std::endl
+            << "\t gridfilename          - the grid file. can be 'uniform'" << std::endl
+            << "\t vtufilename           - save results to vtu file" << std::endl
+            << "\t rock_list             - use material information from given rocklist" << std::endl
+            << "\t output                - output results in text report format" << std::endl
+            << "\t method                - can be MPC, LLM or Mortar" << std::endl
+            << "\t\t if not specified, Mortar couplings are used" << std::endl
+            << "\t (x|y|z)max/(x|y|z)min - maximum/minimum grid coordinates." << std::endl
+            << "\t\t if not specified, coordinates are found by grid traversal" << std::endl
+            << "\t cells(x|y|z)          - number of cells if gridfilename=uniform." << std::endl
+            << "\t Emin                  - minimum E modulus (to avoid numerical issues)" << std::endl
+            << "\t ctol                  - collapse tolerance in grid parsing" << std::endl
+            << "\t ltol                  - tolerance in iterative linear solvers" << std::endl
+            << "\t linsolver_type=cg     - use the conjugate gradient method" << std::endl
+            << "\t linsolver_type=slu    - use the SuperLU sparse direct solver" << std::endl;
+}
+
+//! \brief Structure holding parameters configurable from command line
+struct Params {
+  //! \brief The eclipse grid file
+  std::string file;
+  //! \brief Rocklist overriding material params in .grdecl
+  std::string rocklist;
+  //! \brief VTU output file
+  std::string vtufile;
+  //! \brief Text output file
+  std::string output;
+  //! \brief Whether or not to use LLM couplings
+  bool LLM; 
+  //! \brief Whether or not to use Mortar couplings
+  bool mortar;
+  //! \brief Whether or not to use MPC couplings
+  bool mpc;
+  //! \brief A scaling parameter for the E-modulus to avoid numerical issues
+  double Emin;
+  //! \brief The tolerance for collapsing nodes in the Z direction
+  double ctol;
+  //! \brief The tolerance for the iterative linear solver
+  double ltol;
+  //! \brief Minimum grid coordinates
+  double min[3];
+  //! \brief Maximum grid coordinates
+  double max[3];
+  //! \brief The linear solver to employ
+  Solver solver;
+
+  // Debugging options
+  //! \brief Number of elements on interface grid in the x direction (LLM)
+  int n1;
+  //! \brief Number of elements on interface grid in the y direction (LLM)
+  int n2;
+  //! \brief cell in x
+  int cellsx;
+  //! \brief cell in y
+  int cellsy;
+  //! \brief cell in z
+  int cellsz;
+};
+
+//! \brief Parse the command line arguments
+void parseCommandLine(int argc, char** argv, Params& p)
+{
+  Opm::parameter::ParameterGroup param(argc, argv);
+  p.max[0]   = param.getDefault("xmax",-1);
+  p.max[1]   = param.getDefault("ymax",-1);
+  p.max[2]   = param.getDefault("zmax",-1);
+  p.min[0]   = param.getDefault("xmin",-1);
+  p.min[1]   = param.getDefault("ymin",-1);
+  p.min[2]   = param.getDefault("zmin",-1);
+  //std::string method = param.getDefault<std::string>("method","mpc");
+  std::string method = param.getDefault<std::string>("method","mortar");
+  p.LLM      = (strcasecmp(method.c_str(),"llm") == 0);
+  //p.mortar   = (strcasecmp(method.c_str(),"mortar") == 0);
+  p.mpc   = (strcasecmp(method.c_str(),"mpc") == 0);
+  p.Emin     = param.getDefault<double>("Emin",0.f);
+  p.ctol     = param.getDefault<double>("ctol",1.e-8);
+  p.ltol     = param.getDefault<double>("ltol",1.e-10);
+  p.file     = param.get<std::string>("gridfilename");
+  p.rocklist = param.getDefault<std::string>("rock_list","");
+  p.vtufile  = param.getDefault<std::string>("vtufilename","");
+  p.output = param.getDefault<std::string>("output","");
+  size_t i;
+  if ((i=p.vtufile.find(".vtu")) != std::string::npos)
+    p.vtufile = p.vtufile.substr(0,i);
+  //std::string solver = param.getDefault<std::string>("linsolver_type","slu");
+  std::string solver = param.getDefault<std::string>("linsolver_type","cg");
+  if (solver == "cg")
+    p.solver = CG;
+  else
+    p.solver = SLU;
+  if (p.file == "uniform") {
+    p.cellsx   = param.getDefault("cellsx",3);
+    p.cellsy   = param.getDefault("cellsy",3);
+    p.cellsz   = param.getDefault("cellsz",3);
+  }
+  p.n1       = -1;
+  p.n2       = -1;
+}
+
+//! \brief Write a log of the simulation to a text file
+void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
+                 const std::vector<double>& volume, 
+                 const Dune::FieldMatrix<double,6,6>& C)
+{
+  // get current time
+  time_t rawtime;
+  struct tm* timeinfo;
+  time(&rawtime);
+  timeinfo = localtime(&rawtime);
+
+  // get hostname
+  char hostname[1024];
+  gethostname(hostname,1024);
+
+  std::string method = "mortar";
+  if (p.LLM)
+    method = "llm";
+  if (p.mpc)
+    method = "mpc";
+
+  // write log
+  std::ofstream f;
+  f.open(p.output.c_str());
+  f << "######################################################################" << std::endl
+    << "# Results from upscaling elastic moduli." << std::endl
+    << "#" << std::endl
+    << "# Finished: " << asctime(timeinfo)
+    << "# Hostname: " << hostname << std::endl
+    << "#" << std::endl
+    << "# Upscaling time: " << watch.secsSinceStart() << " secs" << std::endl
+    << "#" << std::endl;
+  if (p.file == "uniform") {
+    f << "# Uniform grid used" << std::endl
+      << "#\t cells: " << p.cellsx*p.cellsy*p.cellsz << std::endl;
+  }
+  else {
+    f  << "# Eclipse file: " << p.file << std::endl
+       << "#\t cells: " << cells << std::endl;
+  }
+  f << "#" << std::endl;
+  if (!p.rocklist.empty()) {
+    f << "# Rock list: " << p.rocklist << std::endl
+      << "#" << std::endl;
+  }
+  f << "# Options used:" << std::endl
+    << "#\t         method: " << method << std::endl
+    << "#\t linsolver_type: " << (p.solver==SLU?"slu":"cg") << std::endl;
+  if (p.solver == CG)
+    f << "#\t           ltol: " << p.ltol << std::endl;
+  if (p.file == "uniform") {
+    f << "#\t          cellsx: " << p.cellsx << std::endl
+      << "#\t          cellsy: " << p.cellsy << std::endl
+      << "#\t          cellsz: " << p.cellsz << std::endl;
+  }
+  f << "#" << std::endl
+    <<"# Materials: " << volume.size() << std::endl;
+  for (int i=0;i<volume.size();++i)
+    f << "#\t Material" << i+1 << ": " << volume[i]*100 << "%" << std::endl;
+  f << "#" << std::endl
+    << "######################################################################" << std::endl
+    << C << std::endl;
+}
+
+//! \brief Main driver
+int main(int argc, char** argv)
+{
+  try {
+    static const int dim = 3;
+    static const int bfunc = 4+(dim-2)*4;
+
+    typedef Dune::CpGrid GridType;
+
+    if (argc < 2 || strcasecmp(argv[1],"-h") == 0 
+                 || strcasecmp(argv[1],"--help") == 0
+                 || strcasecmp(argv[1],"-?") == 0) {
+      syntax(argv);
+      exit(1);
+    }
+    Params p;
+    parseCommandLine(argc,argv,p);
+
+    Opm::time::StopWatch watch;
+    watch.start();
+
+    GridType grid;
+    if (p.file == "uniform") {
+      Dune::array<int,3> cells;
+      cells[0] = p.cellsx;
+      cells[1] = p.cellsy;
+      cells[2] = p.cellsz;
+      Dune::array<double,3> cellsize;
+      cellsize[0] = cellsize[1] = cellsize[2] = 1.f; 
+      grid.createCartesian(cells,cellsize);
+    } else
+      grid.readEclipseFormat(p.file,p.ctol,false);
+
+    typedef GridType::ctype ctype;
+    ElasticityUpscale<GridType> upscale(grid, p.ctol, p.Emin, p.file, p.rocklist);
+    if (p.max[0] < 0 || p.min[0] < 0) {
+      std::cout << "determine side coordinates..." << std::endl;
+      upscale.findBoundaries(p.min,p.max);
+      std::cout << "  min " << p.min[0] << " " << p.min[1] << " " << p.min[2] << std::endl;
+      std::cout << "  max " << p.max[0] << " " << p.max[1] << " " << p.max[2] << std::endl;
+    }
+    if (p.n1 == -1 || p.n2 == -1) {
+      p.n1 = grid.logicalCartesianSize()[0];
+      p.n2 = grid.logicalCartesianSize()[1];
+    }
+
+    if (p.LLM) {
+      std::cout << "using LLM couplings..." << std::endl;
+      upscale.periodicBCsLLM(p.min,p.max,p.n1,p.n2);
+    } else if (p.mpc) {
+      std::cout << "using MPC couplings in all directions..." << std::endl;
+      upscale.periodicBCs(p.min,p.max);
+      std::cout << "preprocessing grid..." << std::endl;
+      upscale.A.initForAssembly();
+    } else {
+      std::cout << "using Mortar couplings.." << std::endl;
+      upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2);
+    } 
+    Dune::FieldMatrix<double,6,6> C;
+    Dune::VTKWriter<GridType::LeafGridView> vtkwriter(grid.leafView());
+    Vector field[6];
+    std::cout << "assembling..." << "\n";
+    upscale.assemble(-1,true);
+    upscale.setupSolvers(p.solver);
+#pragma omp parallel for schedule(static)
+    for (int i=0;i<6;++i) {
+      upscale.assemble(i,false);
+      std::cout << "solving case " << i << "..." << "\n";
+      upscale.solve(p.solver,p.ltol,i);
+      upscale.A.expandSolution(field[i],upscale.u[i]);
+#define CLAMP(x) (fabs(x)<1.e-5?0.f:x)
+      for (int j=0;j<field[i].size();++j) {
+        double val = field[i][j];
+        field[i][j] = CLAMP(val);
+      }
+      Dune::FieldVector<double,6> v;
+      upscale.averageStress(v,upscale.u[i],i);
+      for (int j=0;j<6;++j)
+        C[i][j] = CLAMP(v[j]);
+    }
+    for (int i=0;i<6;++i) {
+      std::stringstream str;
+      str << "sol " << i+1;
+      vtkwriter.addVertexData(field[i], str.str().c_str(), dim);
+    }
+    if (!p.vtufile.empty())
+      vtkwriter.write(p.vtufile);
+    // voigt notation
+    for (int j=0;j<6;++j)
+      std::swap(C[3][j],C[5][j]);
+    for (int j=0;j<6;++j)
+      std::swap(C[j][3],C[j][5]);
+    std::cout << "---------" << std::endl;
+    std::cout << C << std::endl;
+    if (!p.output.empty()) {
+      writeOutput(p,watch,grid.size(0),upscale.volumeFractions,C);
+    }
+
+    return 0;
+  }
+  catch (Dune::Exception &e) {
+    std::cerr << "Dune reported error: " << e << std::endl;
+  }
+  catch (...) {
+    std::cerr << "Unknown exception thrown!" << std::endl;
+  }
+  return 1;
+}