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