added: support for higher order lagrangian multipliers

- add PN shape function support
- make the way we group boundary elements into pillars more robust
- use the new shape functions in mortar matrix assembly
- add options to specify multiplier order (defaults to linears for now)

dune/elasticity/boundarygrid.hh

@@ -113,6 +113,13 @@ class BoundaryGrid {
     //! \param[in] quad The quad to add
     void add(const Quad& quad);
 
+    void addToColumn(size_t col, const Quad& quad)
+    {
+      if (col >= colGrids.size())
+        colGrids.resize(col+1);
+      colGrids[col].push_back(quad);
+    }
+
     //! \brief Obtain a reference to a quad
     //! \param[in] index The index of the requested quad
     //! \returns A reference to the requested quad
@@ -129,12 +136,22 @@ class BoundaryGrid {
       return grid[index];
     }
 
+    const Quad& getQuad(int col, int index) const
+    {
+      return colGrids[col][index];
+    }
+
     //! \brief Obtain the number of quads in the grid
     size_t size() const
     {
       return grid.size();
     }
 
+    size_t colSize(int i) const
+    {
+      return colGrids[i].size();
+    }
+
     //! \brief Return the total number of nodes on the grid when known
     //! \sa uniform
     size_t totalNodes() const
@@ -209,6 +226,7 @@ class BoundaryGrid {
   protected:
     //! \brief Our quadrilateral elements
     std::vector<Quad> grid;
+    std::vector<std::vector<Quad> > colGrids;
 
     //! \brief Whether or not a given node is marked as fixed
     std::vector<bool> fixNodes;

dune/elasticity/elasticity_upscale.hpp

@@ -151,8 +151,11 @@ class ElasticityUpscale
     //! \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
+    //! \param[in] p1 The order of multipliers in the X direction
+    //! \param[in] p2 The order of multipliers in the Y direction
     void periodicBCsMortar(const double* min,
-                           const double* max, int n1, int n2);
+                           const double* max, int n1, int n2,
+                           int p1, int p2);
 
     //! \brief Assemble (optionally) stiffness matrix A and load vector
     //! \param[in] loadcase The strain load case. Set to -1 to skip
@@ -283,8 +286,11 @@ class ElasticityUpscale
     //! \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)
+    //! \param[in] n1 The multipler order in the first direction
+    //! \param[in] n2 The multipler order in the second direction
     Matrix findLMatrixMortar(const BoundaryGrid& b1,
-                             const BoundaryGrid& interface, int dir);
+                             const BoundaryGrid& interface, int dir,
+                             int n1, int n2);
 
     //! \brief This function loads and maps materials to active grid cells
     //! \param[in] file The eclipse grid to read materials from

dune/elasticity/elasticity_upscale_impl.hpp

@@ -56,6 +56,9 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
   int c = 0;
   int i = log2(dir);
   BoundaryGrid result;
+  // we first group nodes into this map through the coordinate of lower left 
+  // vertex. we then split this up into pillars for easy processing later
+  std::map<double, std::vector<BoundaryGrid::Quad> > nodeMap;
   for (LeafIterator cell  = gv.leafView().template begin<0>(); 
                     cell != cellend; ++cell, ++c) {
     std::vector<BoundaryGrid::Vertex> verts;
@@ -99,10 +102,27 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
         q.v[2] = maxXmaxY(verts);
         q.v[3] = minXmaxY(verts);
       }
+      std::map<double, std::vector<BoundaryGrid::Quad> >::iterator it;
+      for (it  = nodeMap.begin(); it != nodeMap.end(); ++it) {
+        if (fabs(it->first-q.v[0].c[0]) < 1.e-7) {
+          it->second.push_back(q);
+          break;
+        }
+      }
+      if (it == nodeMap.end())
+        nodeMap[q.v[0].c[0]].push_back(q);
+
       result.add(q);
     }
   }
 
+  int p=0;
+  std::map<double, std::vector<BoundaryGrid::Quad> >::const_iterator it;
+  for (it = nodeMap.begin(); it != nodeMap.end(); ++it, ++p) {
+    for (size_t i=0;i<it->second.size();++i)
+      result.addToColumn(p,it->second[i]);
+  }
+
   return result;
 }
 
@@ -301,34 +321,65 @@ Matrix ElasticityUpscale<GridType>::findLMatrixLLM(const SlaveGrid& slave,
   return L;
 }
 
+static std::vector< std::vector<int> > renumber(const BoundaryGrid& b,
+                                                int n1, int n2)
+{
+  std::vector<std::vector<int> > nodes;
+  nodes.resize(b.size());
+  // loop over elements
+  int ofs = 0;
+  for (size_t e=0; e < b.size(); ++e) {
+    // first direction major ordered nodes within each element
+    for (int i2=0; i2 < n2; ++i2) {
+      if (e != 0)
+        nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
+      for (int i1=(e==0?0:1); i1 < n1; ++i1)
+        nodes[e].push_back(ofs++);
+    }
+  }
+
+  return nodes;
+}
+
   template<class GridType>
 Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
-                                                  const BoundaryGrid& interface,
-                                                  int dir)
+                                                      const BoundaryGrid& interface,
+                                                      int dir, int n1, int n2)
 {
   std::vector< std::set<int> > adj;
   adj.resize(A.getEqns());
 
+#define QINDEX(p,q,dir) (q)
+
+  // get a set of P1 shape functions for the displacements
+  P1ShapeFunctionSet<ctype,ctype,2> ubasis = 
+                P1ShapeFunctionSet<ctype,ctype,2>::instance();
+
+  // get a set of PN shape functions for the multipliers
+  PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
+
+  std::vector<std::vector<int> > lnodes = renumber(interface,n1,n2);
+  int totalNodes = lnodes.back().back()+1;
+
   // process pillar by pillar
-  size_t per_pillar = b1.size()/interface.size();
   for (size_t p=0;p<interface.size();++p) {
-    for (size_t q=0;q<per_pillar;++q) {
-      for (size_t i=0;i<4;++i) {
+    for (size_t q=0;q<b1.colSize(p);++q) {
+      for (size_t i=0;i<ubasis.n;++i) {
         for (size_t d=0;d<3;++d) {
-          MPC* mpc = A.getMPC(b1[p*per_pillar+q].v[i].i,d);
+          MPC* mpc = A.getMPC(b1.getQuad(p,QINDEX(p,q,dir)).v[i].i,d);
           if (mpc) {
-            for (int n=0;n<mpc->getNoMaster();++n) {
+            for (size_t n=0;n<mpc->getNoMaster();++n) {
               int dof = A.getEquationForDof(mpc->getMaster(n).node,d);
               if (dof > -1) {
-                for (int j=0;j<4;++j)
-                  adj[dof].insert(3*interface[p].v[j].i+d);
+                for (size_t j=0;j<lnodes[p].size();++j)
+                  adj[dof].insert(3*lnodes[p][j]+d);
               }
             }
           } else {
-            int dof = A.getEquationForDof(b1[p*per_pillar+q].v[i].i,d);
+            int dof = A.getEquationForDof(b1.getQuad(p,QINDEX(p,q,dir)).v[i].i,d);
             if (dof > -1) {
-              for (int j=0;j<4;++j)
-                adj[dof].insert(3*interface[p].v[j].i+d);
+              for (size_t j=0;j<lnodes[p].size();++j)
+                adj[dof].insert(3*lnodes[p][j]+d);
             }
           }
         }
@@ -337,41 +388,40 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
   }
 
   Matrix B;
-  MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*interface.totalNodes());
+  MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*totalNodes);
 
-  // get a set of P1 shape functions for the velocity
-  P1ShapeFunctionSet<ctype,ctype,2> ubasis = P1ShapeFunctionSet<ctype,ctype,2>::instance();
-  // get a set of P1 shape functions for multipliers
-  P1ShapeFunctionSet<ctype,ctype,2> lbasis = P1ShapeFunctionSet<ctype,ctype,2>::instance();
   // get a reference element
   Dune::GeometryType gt;
   gt.makeCube(2);
-  const Dune::template GenericReferenceElement<ctype,2> &ref =
-    Dune::GenericReferenceElements<ctype,2>::general(gt);
   // get a quadrature rule
+  int quadorder = std::max((1.0+n1+0.5)/2.0,(1.0+n2+0.5)/2.0);
+  quadorder = std::max(quadorder, 2);
   const Dune::QuadratureRule<ctype,2>& rule = 
-                        Dune::QuadratureRules<ctype,2>::rule(gt,2);
+                  Dune::QuadratureRules<ctype,2>::rule(gt,quadorder);
 
   // do the assembly loop
   typename Dune::QuadratureRule<ctype,2>::const_iterator r;
+  Dune::DynamicMatrix<ctype> E(ubasis.n,(n1+1)*(n2+1),0.0);
   for (size_t p=0;p<interface.size();++p) {
     const BoundaryGrid::Quad& qi(interface[p]);
     HexGeometry<2,2,GridType> lg(qi);
-    for (size_t q=0;q<per_pillar;++q) {
-      const BoundaryGrid::Quad& qu(b1[p*per_pillar+q]);
+    for (size_t q=0;q<b1.colSize(p);++q) {
+      const BoundaryGrid::Quad& qu = b1.getQuad(p,q);
       HexGeometry<2,2,GridType> hex(qu,gv,dir);
-      Dune::FieldMatrix<ctype,4,4> E;
       E = 0;
       for (r = rule.begin(); r != rule.end();++r) {
         ctype detJ = hex.integrationElement(r->position());
-        if (detJ < 0)
-          continue;
+        if (detJ < 1.e-4 && r == rule.begin())
+          std::cerr << "warning: interface cell (close to) degenerated, |J|=" << detJ << std::endl;
+
         typename HexGeometry<2,2,GridType>::LocalCoordinate loc = 
                                         lg.local(hex.global(r->position()));
-        for (int i=0;i<4;++i) {
-          for (int j=0;j<4;++j)
+        assert(loc[0] <= 1.0+1.e-4 && loc[0] >= 0.0 && loc[1] <= 1.0+1.e-4 && loc[1] >= 0.0);
+        for (int i=0;i<ubasis.n;++i) {
+          for (int j=0;j<lbasis.size();++j) {
             E[i][j] += ubasis[i].evaluateFunction(r->position())*
                        lbasis[j].evaluateFunction(loc)*detJ*r->weight();
+          }
         }
       }
       // and assemble element contributions
@@ -379,11 +429,11 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
         for (int i=0;i<4;++i) {
           MPC* mpc = A.getMPC(qu.v[i].i,d);
           if (mpc) {
-            for (int n=0;n<mpc->getNoMaster();++n) {
+            for (size_t n=0;n<mpc->getNoMaster();++n) {
               int indexi = A.getEquationForDof(mpc->getMaster(n).node,d);
               if (indexi > -1) {
-                for (int j=0;j<4;++j) {
-                  int indexj = qi.v[j].i*3+d;
+                for (size_t j=0;j<lnodes[p].size();++j) {
+                  int indexj = lnodes[p][j]*3+d;
                   B[indexi][indexj] += E[i][j];
                 }
               }
@@ -391,8 +441,8 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
           } else {
             int indexi = A.getEquationForDof(qu.v[i].i,d);
             if (indexi > -1) {
-              for (int j=0;j<4;++j) {
-                int indexj = qi.v[j].i*3+d;
+              for (size_t j=0;j<lnodes[p].size();++j) {
+                int indexj = lnodes[p][j]*3+d;
                 B[indexi][indexj] += E[i][j];
               }
             }
@@ -802,7 +852,8 @@ void ElasticityUpscale<GridType>::periodicBCs(const double* min,
   template<class GridType>
 void ElasticityUpscale<GridType>::periodicBCsMortar(const double* min, 
                                                     const double* max,
-                                                    int n1, int n2)
+                                                    int n1, int n2,
+                                                    int p1, int p2)
 {
   // this method
   // 1. fixes the primal corner dofs
@@ -844,13 +895,13 @@ void ElasticityUpscale<GridType>::periodicBCsMortar(const double* min,
   BoundaryGrid lambday = BoundaryGrid::uniform(fmin,fmax,n1,1,true);
 
   // step 5
-  Matrix L1 = findLMatrixMortar(master[0],lambdax,0);
-  Matrix L2 = findLMatrixMortar(master[1],lambdax,0);
+  Matrix L1 = findLMatrixMortar(master[0],lambdax,0, p2, 1);
+  Matrix L2 = findLMatrixMortar(master[1],lambdax,0, p2, 1);
   L.push_back(MatrixOps::Axpy(L1,L2,-1));
 
   // step 6
-  Matrix L3 = findLMatrixMortar(master[2],lambday,1);
-  Matrix L4 = findLMatrixMortar(master[3],lambday,1);
+  Matrix L3 = findLMatrixMortar(master[2],lambday,1, p1, 1);
+  Matrix L4 = findLMatrixMortar(master[3],lambday,1, p1, 1);
   L.push_back(MatrixOps::Axpy(L3,L4,-1));
 }
 

dune/elasticity/shapefunctions.hpp

@@ -6,12 +6,15 @@
 //!
 //! \author Arne Morten Kvarving / SINTEF
 //!
-//! \brief Class for linear shape functions. Loosely based on code in dune-grid-howto 
+//! \brief Classes for shape functions. Loosely based on code in dune-grid-howto
 //!
 //==============================================================================
 #pragma once
 
 #include <dune/common/fvector.hh>
+#include "dynmatrixev.hh"
+
+#include <complex>
 
 namespace Opm {
 namespace Elasticity {
@@ -80,6 +83,96 @@ class LinearShapeFunction
     Dune::FieldVector<rtype,dim> coeff1;
 };
 
+//! \brief Represents a cardinal function on a line
+  template<class ctype, class rtype>
+class LagrangeCardinalFunction
+{
+  public:
+    //! \brief Empty default constructor
+    LagrangeCardinalFunction() {}
+
+    //! \brief Construct a cardinal function with the given nodes
+    //! \param[in] nodes_ The nodes
+    //! \param[in] i The node this function is associated with
+    LagrangeCardinalFunction(const std::vector<rtype>& nodes_,
+                             size_t i)
+      : nodes(nodes_), node(i) {}
+
+    //! \brief Evaluate the shape function
+    //! \param[in] local The local coordinates
+    rtype evaluateFunction(const ctype& local) const
+    {
+      rtype result = 1;
+      for (size_t i=0; i < nodes.size(); ++i) {
+        if (i != node)
+          result *= (local-nodes[i])/(nodes[node]-nodes[i]);
+      }
+
+      return result;
+    }
+
+    //! \brief Evaluate the derivative of the cardinal function
+    //! \param[in] local The local coordinates
+    rtype evaluateGradient(const ctype& local) const
+    {
+      rtype result = 0;
+      for (size_t i=0; i < nodes.size(); ++i) {
+        rtype f = 1;
+        for (int j=0; j < nodes.size(); ++j) {
+          if (i != j && j != node)
+            f *= (local-nodes[j])/(nodes[node]-nodes[j]);
+        }
+        result += f/(nodes[node]-nodes[i]);
+      }
+
+      return result;
+    }
+
+  private:
+    //! \brief The nodes
+    std::vector<rtype> nodes;
+
+    size_t node;
+};
+
+//! \brief Represents a tensor-product of 1D functions
+  template<class rtype, class ctype, class ftype, int dim>
+class TensorProductFunction
+{
+  public:
+    //! \brief The dimension of the function
+    enum { dimension = dim };
+
+    //! \brief Empty default constructor
+    TensorProductFunction() {}
+
+    //! \brief Construct a tensor-product function
+    //! \param[in] funcs_ The functions
+    TensorProductFunction(const Dune::FieldVector<ftype, dim>& funcs_)
+      : funcs(funcs_) {}
+
+    //! \brief Evaluate the 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 *= funcs[i].evaluateFunction(local[i]);
+
+      return result;
+    }
+
+    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] = funcs[i].evaluateGradient(local[i]);
+    }
+  private:
+    Dune::FieldVector<ftype, dim> funcs;
+};
+
 //! \brief Singleton handler for the set of LinearShapeFunction
   template<class ctype, class rtype, int dim>
 class P1ShapeFunctionSet
@@ -121,8 +214,8 @@ private:
                                     0, 1,
                                     1, 0,
                                     0, 0};
-        static rtype coeffs22[] = {-1,-1, 
-                                    1,-1, 
+        static rtype coeffs22[] = {-1,-1,
+                                    1,-1,
                                    -1, 1,
                                     1, 1};
 
@@ -134,8 +227,8 @@ private:
                                     0, 1, 0,
                                     1, 0, 0,
                                     0, 0, 0};
-        static rtype coeffs32[] = {-1,-1,-1, 
-                                    1,-1,-1, 
+        static rtype coeffs32[] = {-1,-1,-1,
+                                    1,-1,-1,
                                    -1, 1,-1,
                                     1, 1,-1,
                                    -1,-1, 1,
@@ -173,5 +266,147 @@ private:
     ShapeFunction f[n];
 };
 
+  template<int dim>
+class PNShapeFunctionSet
+{
+public:
+    typedef LagrangeCardinalFunction<double, double> CardinalFunction;
+
+    typedef TensorProductFunction<double, double, CardinalFunction, dim>
+            ShapeFunction;
+
+    PNShapeFunctionSet(int n1, int n2, int n3=0)
+    {
+      int dims[3] = {n1, n2, n3};
+      cfuncs.resize(dim);
+      for (int i=0; i < dim; ++i) {
+        std::vector<double> grid;
+        grid = gaussLobattoLegendreGrid(dims[i]);
+        for (int j=0;j<dims[i];++j)
+          cfuncs[i].push_back(CardinalFunction(grid,j));
+      }
+      int l=0;
+      Dune::FieldVector<CardinalFunction,dim> fs;
+      if (dim == 3) {
+        f.resize(n1*n2*n3);
+        for (int k=0; k < n3; ++k) {
+          for (int j=0; j < n2; ++j)
+            for (int i=0; i< n1; ++i) {
+              fs[0] = cfuncs[0][i];
+              fs[1] = cfuncs[1][j];
+              fs[2] = cfuncs[2][k];
+              f[l++] = ShapeFunction(fs);
+            }
+        }
+      } else {
+        f.resize(n1*n2);
+        for (int j=0; j < n2; ++j) {
+          for (int i=0; i< n1; ++i) {
+            fs[0] = cfuncs[0][i];
+            fs[1] = cfuncs[1][j];
+            f[l++] = ShapeFunction(fs);
+          }
+        }
+      }
+    }
+
+    //! \brief Obtain a given shape function
+    //! \param[in] i The requested shape function
+    const ShapeFunction& operator[](int i) const
+    {
+      return f[i];
+    }
+
+    int size()
+    {
+      return f.size();
+    }
+protected:
+    std::vector< std::vector<CardinalFunction> > cfuncs;
+    std::vector<ShapeFunction> f;
+
+    double legendre(double x, int n)
+    {
+      std::vector<double> Ln;
+      Ln.resize(n+1);
+      Ln[0] = 1.f;
+      Ln[1] = x;
+      if( n > 1 ) {
+        for( int i=1;i<n;i++ )
+          Ln[i+1] = (2*i+1.0)/(i+1.0)*x*Ln[i]-i/(i+1.0)*Ln[i-1];
+      }
+
+      return Ln[n];
+    }
+
+    double legendreDerivative(double x, int n)
+    {
+      std::vector<double> Ln;
+      Ln.resize(n+1);
+
+      Ln[0] = 1.0; Ln[1] = x;
+
+      if( (x == 1.0) || (x == -1.0) )
+        return( pow(x,n-1)*n*(n+1.0)/2.0 );
+      else {
+        for( int i=1;i<n;i++ )
+          Ln[i+1] = (2.0*i+1.0)/(i+1.0)*x*Ln[i]-(double)i/(i+1.0)*Ln[i-1];
+        return( (double)n/(1.0-x*x)*Ln[n-1]-n*x/(1-x*x)*Ln[n] );
+      }
+    }
+
+    std::vector<double> gaussLegendreGrid(int n)
+    {
+      Dune::DynamicMatrix<double> A(n,n,0.0);
+
+      A[0][1] = 1.f;
+      for (int i=1;i<n-1;++i) {
+        A[i][i-1] = (double)i/(2.0*(i+1.0)-1.0);
+        A[i][i+1] = (double)(i+1.0)/(2*(i+1.0)-1.0);
+      }
+      A[n-1][n-2] = (n-1.0)/(2.0*n-1.0);
+
+      Dune::DynamicVector<std::complex<double> > eigenValues(n);
+      Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues);
+
+      std::vector<double> result(n);
+      for (int i=0; i < n; ++i)
+        result[i] = std::real(eigenValues[i]);
+      std::sort(result.begin(),result.begin()+n);
+
+      return result;
+    }
+
+    std::vector<double> gaussLobattoLegendreGrid(int n)
+    {
+      assert(n > 1);
+      const double tolerance = 1.e-15;
+
+      std::vector<double> result(n);
+      result[0] = 0.0;
+      result[n-1] = 1.0;
+      if (n == 3)
+        result[1] = 0.5;
+
+      if (n < 4)
+        return result;
+
+      std::vector<double> glgrid = gaussLegendreGrid(n-1);
+      for (int i=1;i<n-1;++i) {
+        result[i] = (glgrid[i-1]+glgrid[i])/2.0;
+        double old = 0.0;
+        while (std::abs(old-result[i]) > tolerance) {
+          old = result[i];
+          double L = legendre(old, n-1);
+          double Ld = legendreDerivative(old, n-1);
+          result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L);
+        }
+        result[i] = (result[i]+1.0)/2.0;
+      }
+
+      return result;
+    }
+};
+
 }
 }

examples/upscale_elasticity.cpp

@@ -91,6 +91,8 @@ struct Params {
   int cellsy;
   //! \brief cell in z
   int cellsz;
+  //! \brief Polynomial order of multipliers in first / second direction
+  int lambda[2];
 };
 
 //! \brief Parse the command line arguments
@@ -104,6 +106,8 @@ void parseCommandLine(int argc, char** argv, Params& p)
   p.min[1]   = param.getDefault("ymin",-1);
   p.min[2]   = param.getDefault("zmin",-1);
   //std::string method = param.getDefault<std::string>("method","mpc");
+  p.lambda[0] = param.getDefault("lambdax", 1);
+  p.lambda[1] = param.getDefault("lambday", 1);
   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);
@@ -254,8 +258,8 @@ int main(int argc, char** argv)
       upscale.A.initForAssembly();
     } else {
       std::cout << "using Mortar couplings.." << std::endl;
-      upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2);
-    } 
+      upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]);
+    }
     Dune::FieldMatrix<double,6,6> C;
     Dune::VTKWriter<GridType::LeafGridView> vtkwriter(grid.leafView());
     Opm::Elasticity::Vector field[6];