Merge pull request #395 from bska/geomech

Add Logic for Geomechanical Use Cases

CMakeLists_files.cmake

@@ -212,6 +212,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/elasticity/material.hh
   opm/elasticity/materials.hh
   opm/elasticity/matrixops.hpp
+  opm/elasticity/meshcolorizer.hpp
   opm/elasticity/mortar_evaluator.hpp
   opm/elasticity/mortar_schur.hpp
   opm/elasticity/mortar_schur_precond.hpp

opm/elasticity/asmhandler.hpp

@@ -111,6 +111,11 @@ class ASMHandler {
                     const LeafIterator& cell,
                     Vector* b=NULL);
 
+    void addMatElement(const int i, const int j, const double val)
+    {
+      A[i][j] += val;
+    }
+
     //! \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

opm/elasticity/elasticity.hpp

@@ -35,6 +35,10 @@ class Elasticity {
     //! \param[in] point (Reference) coordinates of quadrature point 
     //! \param[in] Jinv Jacobian matrix in quadrature point
     //! \param[out] B The B matrix
+      template<int funcdim>
+    void getBVector(Dune::FieldVector<ctype,funcdim>& BVector,
+                    const Dune::FieldVector<ctype,dim>& point);
+  
       template<int components, int funcdim>
     void getBmatrix(Dune::FieldMatrix<ctype,components,funcdim>& B,
                     const Dune::FieldVector<ctype,dim>& point,

opm/elasticity/elasticity_impl.hpp

@@ -71,6 +71,21 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
       B[i][j] = temp[k%rows][k/rows];
 }
 
+
+    template<class GridType>
+    template<int funcdims>
+void Elasticity<GridType>::getBVector(Dune::FieldVector<ctype,funcdims>& Bvector,
+                                      const Dune::FieldVector<ctype,dim>& point)
+{
+    P1ShapeFunctionSet<ctype,ctype,dim> basis 
+        = P1ShapeFunctionSet<ctype,ctype,dim>::instance();
+
+    Dune::FieldMatrix<ctype,funcdims,dim> N;
+    for (int i = 0; i < funcdims; i++) {
+        Bvector[i] = basis[i].evaluateFunction(point);
+    }
+}
+    
   template <class GridType>
     template<int comp, int funcdim>
 void Elasticity<GridType>::getStiffnessMatrix(

opm/elasticity/meshcolorizer.hpp

@@ -44,7 +44,16 @@ class MeshColorizer {
     }
 
     //! \brief Calculate the coloring
-    void calcGroups();
+    void calcGroups()
+    {
+      tg[0].resize(1);
+      tg[0].reserve(grid.size(0));
+      for (int i = 0; i < grid.size(0); ++i) {
+        tg[0][0].push_back(i);
+      }
+      tg[1].clear();
+    }
+
   private:
     IntMat tg[2]; //!< The color groups
     const GridType& grid; //!< Reference to grid being colored