boundarygrid.hh

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef BOUNDARYGRID_HH_
#define BOUNDARYGRID_HH_

#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/genericgeometry/matrixhelper.hh>
#include <dune/grid/common/mcmgmapper.hh>

#include <vector>

namespace Opm {
namespace Elasticity {

class BoundaryGrid {
  public:
    typedef Dune::FieldVector<double,3> GlobalCoordinate;
    typedef Dune::FieldVector<double,2> FaceCoord;

    //              in natural order along the first direction
    static BoundaryGrid uniform(const FaceCoord& min, const FaceCoord& max,
                                int k1, int k2, bool dc=false);

    BoundaryGrid() {}

    virtual ~BoundaryGrid() {}

    //         on a boundary
    class Quad;

    class Vertex {
      public:
        Vertex() : i(-1), c(0), q(0), fixed(false) {}

        int i;

        FaceCoord c;

        Quad* q;

        bool fixed;

        bool operator==(const Vertex& v2)
        {
          return hypot(v2.c[0]-c[0],v2.c[1]-c[1]) < 1.e-8;
        }
    };
    class Quad {
      public:
        Quad()
        {
          v[0].i = v[1].i = v[2].i = v[3].i = 0;
          v[0].c = v[1].c = v[2].c = v[3].c = 0.f;
        }
        FaceCoord pos(double xi, double eta) const;

        std::vector<double> evalBasis(double xi, double eta) const;

        Vertex v[4];
        FaceCoord bb[2];
      protected:
        friend std::ostream& operator <<(std::ostream& os, const Quad& q)
        {
          os << "Nodes: " << q.v[0].i << "," << q.v[1].i << "," << q.v[2].i << "," << q.v[3].i << std::endl;
          os << "(" << q.v[0].c << ")(" << q.v[1].c << ")(" << q.v[2].c << ")(" << q.v[3].c << ")";
          return os;
        }
    };

    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);
    }

    Quad& operator[](int index)
    {
      return grid[index];
    }

    const Quad& operator[](int index) const
    {
      return grid[index];
    }

    const Quad& getQuad(int col, int index) const
    {
      return colGrids[col][index];
    }

    size_t size() const
    {
      return grid.size();
    }

    size_t colSize(int i) const
    {
      return colGrids[i].size();
    }

    size_t totalNodes() const
    {
      return nodes;
    }

    bool isFixed(int node) const
    {
      return fixNodes[node];
    }

    bool find(Vertex& res, const Vertex& coord) const;

    static void extract(FaceCoord& res,
                        const GlobalCoordinate& coord, int dir);

    static void extract(Vertex& res,
                        const Dune::FieldVector<double,3>& coord, int dir);

    struct VertexLess {
        VertexLess(int comp) : dir(comp) {}

        bool operator()(const Vertex& q1, const Vertex& q2)
        {
          if (dir >= 0)
            return q1.c[dir] < q2.c[dir];
          return q1.i < q2.i;
        }

        int dir;
    };

    struct BoundedPredicate {
      BoundedPredicate(const FaceCoord& coord_) : coord(coord_) {}

      bool operator()(const Quad& q)
      {
        double eps = 1.e-8;
        return (coord[0] >= q.bb[0][0]-eps && 
                coord[0] <= q.bb[1][0]+eps &&
                coord[1] >= q.bb[0][1]-eps &&
                coord[1] <= q.bb[1][1]+eps);
      }

      FaceCoord coord;
    };
  protected:
    std::vector<Quad> grid;
    std::vector<std::vector<Quad> > colGrids;

    std::vector<bool> fixNodes;

    size_t nodes;

    friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
    {
      for (size_t i=0;i<g.size();++i)
        os << g[i] << std::endl;
      return os;
    }

    bool bilinearSolve(double epsilon, double order,
                       const double* A, const double* B,
                       std::vector<double>& X,
                       std::vector<double>& Y) const;

    bool cubicSolve(double eps, double A, double B, double C,
                    double D, std::vector<double>& X) const;

    int Q4inv(FaceCoord& res, const Quad& q, const FaceCoord& coord,
              double epsZero, double epsOut) const;
};

  template<int mydim, int cdim, class GridImp>
class HexGeometry
{
};

  template<int cdim, class GridImp>
class HexGeometry<2, cdim, GridImp>
{
  public:
    enum {   dimension = 2};

    enum { mydimension = 2};

    enum { coorddimension = cdim };

    enum {dimensionworld = 2 };

    typedef double ctype;

    typedef Dune::FieldVector<ctype,mydimension> LocalCoordinate;

    typedef Dune::FieldVector<ctype,coorddimension> GlobalCoordinate;

    typedef Dune::FieldMatrix<ctype,coorddimension,mydimension> Jacobian;

    typedef Dune::FieldMatrix<ctype,mydimension,coorddimension> JacobianTransposed;

    HexGeometry(const BoundaryGrid::Quad& q, const GridImp& gv, int dir)
    {
      Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridImp,
                                           Dune::MCMGVertexLayout> mapper(gv);
      typename GridImp::LeafGridView::template Codim<3>::Iterator start=gv.leafGridView().template begin<3>();
      const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafGridView().template end<3>();
      for (int i=0;i<4;++i) {
        typename GridImp::LeafGridView::template Codim<3>::Iterator it=start;
        for (; it != itend; ++it) {
          if (mapper.map(*it) == q.v[i].i)
            break;
        }
        BoundaryGrid::extract(c[i],it->geometry().corner(0),dir);
      }

      m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
    }

    HexGeometry(const BoundaryGrid::Quad& q)
    {
      for (int i=0;i<4;++i)
        c[i] = q.v[i].c;
      m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
    }

    Dune::GeometryType type() const
    {
      Dune::GeometryType t;
      t.makeCube(mydimension);
      return t;
    }

    int corners() const
    {
      return 4;
    }

    ctype volume() const
    {
      return m_volume;
    }

    GlobalCoordinate center() const
    {
      LocalCoordinate local;
      local = .5f;
      return Global(local);
    }

    GlobalCoordinate corner(int cor) const
    {
      return c[cor];
    }

    GlobalCoordinate global(const LocalCoordinate& local) const
    {
      // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
      LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
      uvw[0] -= local;
      // Access pattern for uvw matching ordering of corners.
      const int pat[4][2] = {{ 0, 0 },
                             { 1, 0 },
                             { 0, 1 },
                             { 1, 1 }};
      GlobalCoordinate xyz(0.0);
      for (int i = 0; i < 4; ++i) {
        GlobalCoordinate corner_contrib = corner(i);
        double factor = 1.0;
        for (int j = 0; j < 2; ++j) {
          factor *= uvw[pat[i][j]][j];
        }
        corner_contrib *= factor;
        xyz += corner_contrib;
      }
      return xyz;
    }

    LocalCoordinate local(const GlobalCoordinate& y) const
    {
      const ctype epsilon = 1e-10;
      const Dune::ReferenceElement< ctype , 2 > & refElement =
        Dune::ReferenceElements< ctype, 2 >::general(type());
      LocalCoordinate x = refElement.position(0,0);
      LocalCoordinate dx;
      do {
        using namespace Dune::GenericGeometry;
        // DF^n dx^n = F^n, x^{n+1} -= dx^n
        JacobianTransposed JT = jacobianTransposed(x);
        GlobalCoordinate z = global(x);
        z -= y;
        MatrixHelper<DuneCoordTraits<double> >::template xTRightInvA<2, 2>(JT, z, dx );
        x -= dx;
      } while (dx.two_norm2() > epsilon*epsilon);
      return x;
    }

    const Dune::FieldMatrix<ctype, mydimension, coorddimension>
      jacobianTransposed(const LocalCoordinate& local) const
    {
      // uvw = { (1-u, 1-v, (u, v) }
      LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
      uvw[0] -= local;
      // Access pattern for uvw matching ordering of corners.
      const int pat[4][3] = {{ 0, 0 },
                             { 1, 0 },
                             { 0, 1 },
                             { 1, 1 }};
      Dune::FieldMatrix<ctype, mydimension, coorddimension> Jt(0.0);
      for (int i = 0; i < 4; ++i) {
        for (int deriv = 0; deriv < 2; ++deriv) {
          // This part contributing to dg/du_{deriv}
          double factor = 1.0;
          for (int j = 0; j < 2; ++j) {
            factor *= (j != deriv) ? uvw[pat[i][j]][j]
              : (pat[i][j] == 0 ? -1.0 : 1.0);
          }
          GlobalCoordinate corner_contrib = corner(i);
          corner_contrib *= factor;
          Jt[deriv] += corner_contrib; // using FieldMatrix row access.
        }
      }
      return Jt;
    }

    const Dune::FieldMatrix<ctype, coorddimension, mydimension> 
      jacobianInverseTransposed(const LocalCoordinate& local) const
    {
      Dune::FieldMatrix<ctype, coorddimension, mydimension> Jti = jacobianTransposed(local);
      Jti.invert();
      return Jti;
    }

    ctype integrationElement(const LocalCoordinate& local) const
    {
      Dune::FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local);
      using namespace Dune::GenericGeometry;
      return MatrixHelper<DuneCoordTraits<double> >::template sqrtDetAAT<2, 2>(Jt);
    }
  private:
    GlobalCoordinate c[4];

    ctype m_volume;
};

BoundaryGrid::Vertex minXminY(std::vector<BoundaryGrid::Vertex>& in);

BoundaryGrid::Vertex maxXminY(std::vector<BoundaryGrid::Vertex>& in);

BoundaryGrid::Vertex maxXmaxY(std::vector<BoundaryGrid::Vertex>& in);

BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in);

}
}

#endif