asmhandler.hpp

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef ASMHANDLER_HPP_
#define ASMHANDLER_HPP_
 
#include <opm/common/utility/platform_dependent/disable_warnings.h>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/common/fvector.hh>
 
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 
#include <opm/elasticity/logutils.hpp>
#include <opm/elasticity/mpc.hh>
#include <opm/elasticity/matrixops.hpp>
 
namespace Opm {
namespace Elasticity {
 
 
  template<class GridType>
class ASMHandler {
  public:
    static const int dim = GridType::dimension;
 
    typedef typename GridType::LeafGridView::IndexSet LeafIndexSet;
 
    typedef typename GridType::LeafGridView::template Codim<0>::Iterator LeafIterator;
 
    explicit ASMHandler(const GridType& gv_) : gv(gv_), maxeqn(0)
    {
    }
 
   ~ASMHandler()
   {
     for (MPCMap::iterator it=mpcs.begin(); it != mpcs.end();++it)
       delete it->second;
   }
 
   ASMHandler(const ASMHandler&) = delete;
 
   ASMHandler& operator=(const ASMHandler&) = delete;
 
    typedef Dune::FieldVector<double,dim> NodeValue;
 
    size_t getEqns() const
    {
      return maxeqn;
    }
 
    int getEquationForDof(int node, int dof)
    {
      return meqn[node*dim+dof];
    }
 
    Matrix& getOperator()
    {
      return A;
    }
 
    Vector& getLoadVector()
    {
      return b;
    }
 
    void initForAssembly();
 
      template<int esize>
    void addElement(const Dune::FieldMatrix<double,esize,esize>* K,
                    const Dune::FieldVector<double,esize>* S,
                    const LeafIterator& cell,
                    Vector* b=NULL);
 
    void addMatElement(const int i, const int j, const double val)
    {
      A[i][j] += val;
    }
 
      template<int comp>
    void extractValues(Dune::FieldVector<double,comp>& v, 
                       const Vector& u, const LeafIterator& it);
 
    void expandSolution(Vector& result, const Vector& u);
 
    void addMPC(MPC* mpc);
 
    MPC* getMPC(int node, int dof);
 
    void updateFixedNode(int node, 
                         const std::pair<Direction,NodeValue>& entry);
 
    bool isFixed(int node)
    {
      return (fixedNodes.find(node) != fixedNodes.end());
    }
 
    void printOperator() const;
 
    void printLoadVector() const;
 
    AdjacencyPattern& getAdjacencyPattern()
    {
      return adjacencyPattern;
    }
  protected:
    void resolveMPCChains()
    {
      for (MPCMap::iterator it  = mpcs.begin();
                            it != mpcs.end();++it)
        resolveMPCChain(it->second);
    }
 
    void resolveMPCChain(MPC* mpc);
 
    void preprocess();
 
    void nodeAdjacency(const LeafIterator& it, int vertexsize, int row);
 
    void determineAdjacencyPattern();
 
      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);
 
    MPCMap mpcs;
 
    std::vector<int> meqn;
 
    typedef std::pair<Direction,NodeValue> fixEntry;
 
    typedef std::map<int,fixEntry> fixMap;
 
    typedef typename fixMap::iterator fixIt;
 
    fixMap fixedNodes;
 
    AdjacencyPattern adjacencyPattern;
 
    Matrix A;
 
    Vector b;
 
    const GridType& gv;
 
    size_t maxeqn;
};
 
}
}
 
#include "asmhandler_impl.hpp"
 
#endif