localstiffness.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_LOCAL_STIFFNESS_HH
#define DUMUX_LOCAL_STIFFNESS_HH
 
#include <array>
#include<iostream>
#include<vector>
#include<set>
#include<map>
 
#include<dune/common/timer.hh>
#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/geometry/type.hh>
#include<dune/grid/common/grid.hh>
#include<dune/istl/operators.hh>
#include<dune/istl/bvector.hh>
#include<dune/istl/matrix.hh>
 
#include<dumux/common/boundaryconditions.hh>
#include<dumux/porousmediumflow/sequential/properties.hh>
 
namespace Dumux {
 
template<class TypeTag, int m>
class LocalStiffness
{
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    // grid types
    using Entity = typename GridView::template Codim<0>::Entity;
    enum {n=GridView::dimension};
 
public:
    // types for matrics, vectors and boundary conditions
    using MBlockType = Dune::FieldMatrix<Scalar, m, m>;                      // one entry in the stiffness matrix
    using VBlockType = Dune::FieldVector<Scalar, m>;                        // one entry in the global vectors
    using BCBlockType = std::array<BoundaryConditions::Flags, m>; // componentwise boundary conditions
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
 
    virtual ~LocalStiffness ()
    {
    }
 
    virtual void assemble (const Entity& e, int k=1) = 0;
 
    virtual void assemble (const Entity& e, const Dune::BlockVector<VBlockType>& localSolution, int k=1) = 0;
 
    virtual void assembleBoundaryCondition (const Entity& e, int k=1) = 0;
 
    void print (std::ostream& s, int width, int precision)
    {
      // set the output format
      s.setf(std::ios_base::scientific, std::ios_base::floatfield);
      int oldprec = s.precision();
      s.precision(precision);
 
      for (int i=0; i<currentsize(); i++)
        {
          s << "FEM";    // start a new row
          s << " ";      // space in front of each entry
          s.width(4);    // set width for counter
          s << i;        // number of first entry in a line
          for (int j=0; j<currentsize(); j++)
            {
              s << " ";         // space in front of each entry
              s.width(width);   // set width for each entry anew
              s << mat(i,j);     // yeah, the number !
            }
          s << " ";         // space in front of each entry
          s.width(width);   // set width for each entry anew
          s << rhs(i);
          s << " ";         // space in front of each entry
          s.width(width);   // set width for each entry anew
          s << bc(i)[0];
          s << std::endl;// start a new line
        }
 
 
      // reset the output format
      s.precision(oldprec);
      s.setf(std::ios_base::fixed, std::ios_base::floatfield);
    }
 
    const MBlockType& mat (int i, int j) const
    {
      return A[i][j];
    }
 
    const VBlockType& rhs (int i) const
    {
      return b[i];
    }
 
     const BoundaryTypes& bc (int i) const
     {
       return bctype[i];
     }
 
    void setcurrentsize (int s)
    {
          A.setSize(s,s);
          b.resize(s);
          bctype.resize(s);
    }
 
    int currentsize ()
    {
            return A.N();
    }
 
protected:
    // assembled data
    Dune::Matrix<MBlockType> A;
    std::vector<VBlockType> b;
    std::vector<BoundaryTypes> bctype;
};
 
template<class TypeTag, int m>
class LinearLocalStiffness : public LocalStiffness<TypeTag,m>
{
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    // grid types
    using Entity = typename GridView::template Codim<0>::Entity;
    enum {n=GridView::dimension};
 
public:
    // types for matrics, vectors and boundary conditions
    using MBlockType = Dune::FieldMatrix<Scalar, m, m>;                      // one entry in the stiffness matrix
    using VBlockType = Dune::FieldVector<Scalar, m>;                        // one entry in the global vectors
    using BCBlockType = std::array<BoundaryConditions::Flags, m>;    // componentwise boundary conditions
 
    LinearLocalStiffness ()
    {}
 
    virtual ~LinearLocalStiffness ()
    {
    }
 
    virtual void assemble (const Entity& e, int k=1) = 0;
 
    virtual void assemble (const Entity& e, const Dune::BlockVector<VBlockType>& localSolution, int k=1)
    {
        assemble(e,k);
    }
};
 
} // end namespace Dumux
#endif