boundarygrid.hh

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef BOUNDARYGRID_HH_
#define BOUNDARYGRID_HH_
 
#include <opm/common/utility/platform_dependent/disable_warnings.h>
 
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/defaultgridview.hh>
 
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 
#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)
    {
      using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridImp>>;
      Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
      typename LeafGridView::template Codim<3>::Iterator start=gv.leafGridView().template begin<3>();
      const typename LeafGridView::template Codim<3>::Iterator itend = gv.leafGridView().template end<3>();
      for (int i=0;i<4;++i) {
        typename LeafGridView::template Codim<3>::Iterator it=start;
        for (; it != itend; ++it) {
          if (mapper.index(*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 = Dune::GeometryTypes::cube(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;
      auto refElement = Dune::ReferenceElements<ctype,2>::cube();
      LocalCoordinate x = refElement.position(0,0);
      LocalCoordinate dx;
      do {
        // DF^n dx^n = F^n, x^{n+1} -= dx^n
        JacobianTransposed JT = jacobianTransposed(x);
        GlobalCoordinate z = global(x);
        z -= y;
        Dune::Impl::FieldMatrixHelper<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);
      return Dune::Impl::FieldMatrixHelper<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