boundarygrid.hh
Go to the documentation of this file.
1 //==============================================================================
11 //==============================================================================
12 #ifndef BOUNDARYGRID_HH_
13 #define BOUNDARYGRID_HH_
14 
15 #include <opm/common/utility/platform_dependent/disable_warnings.h>
16 
17 #include <dune/common/version.hh>
18 #include <dune/common/fvector.hh>
19 #include <dune/common/fmatrix.hh>
20 #include <dune/geometry/referenceelements.hh>
21 #if ! DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5)
22 #include <dune/geometry/genericgeometry/matrixhelper.hh>
23 #endif
24 #include <dune/grid/common/mcmgmapper.hh>
25 
26 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
27 
28 
29 #include <vector>
30 
31 namespace Opm {
32 namespace Elasticity {
33 
35 class BoundaryGrid {
36  public:
38  typedef Dune::FieldVector<double,3> GlobalCoordinate;
40  typedef Dune::FieldVector<double,2> FaceCoord;
41 
49  // in natural order along the first direction
50  static BoundaryGrid uniform(const FaceCoord& min, const FaceCoord& max,
51  int k1, int k2, bool dc=false);
52 
55 
57  virtual ~BoundaryGrid() {}
58 
60  // on a boundary
61  class Quad;
62 
64  class Vertex {
65  public:
67  Vertex() : i(-1), c(0), q(0), fixed(false) {}
68 
70  int i;
71 
74 
76  Quad* q;
77 
79  bool fixed;
80 
82  bool operator==(const Vertex& v2)
83  {
84  return hypot(v2.c[0]-c[0],v2.c[1]-c[1]) < 1.e-8;
85  }
86  };
88  class Quad {
89  public:
91  Quad()
92  {
93  v[0].i = v[1].i = v[2].i = v[3].i = 0;
94  v[0].c = v[1].c = v[2].c = v[3].c = 0.f;
95  }
100  FaceCoord pos(double xi, double eta) const;
101 
105  std::vector<double> evalBasis(double xi, double eta) const;
106 
108  Vertex v[4];
111  protected:
113  friend std::ostream& operator <<(std::ostream& os, const Quad& q)
114  {
115  os << "Nodes: " << q.v[0].i << "," << q.v[1].i << "," << q.v[2].i << "," << q.v[3].i << std::endl;
116  os << "(" << q.v[0].c << ")(" << q.v[1].c << ")(" << q.v[2].c << ")(" << q.v[3].c << ")";
117  return os;
118  }
119  };
120 
123  void add(const Quad& quad);
124 
125  void addToColumn(size_t col, const Quad& quad)
126  {
127  if (col >= colGrids.size())
128  colGrids.resize(col+1);
129  colGrids[col].push_back(quad);
130  }
131 
135  Quad& operator[](int index)
136  {
137  return grid[index];
138  }
139 
143  const Quad& operator[](int index) const
144  {
145  return grid[index];
146  }
147 
148  const Quad& getQuad(int col, int index) const
149  {
150  return colGrids[col][index];
151  }
152 
154  size_t size() const
155  {
156  return grid.size();
157  }
158 
159  size_t colSize(int i) const
160  {
161  return colGrids[i].size();
162  }
163 
166  size_t totalNodes() const
167  {
168  return nodes;
169  }
170 
174  bool isFixed(int node) const
175  {
176  return fixNodes[node];
177  }
178 
182  bool find(Vertex& res, const Vertex& coord) const;
183 
188  static void extract(FaceCoord& res,
189  const GlobalCoordinate& coord, int dir);
190 
195  static void extract(Vertex& res,
196  const Dune::FieldVector<double,3>& coord, int dir);
197 
199  struct VertexLess {
202  VertexLess(int comp) : dir(comp) {}
203 
205  bool operator()(const Vertex& q1, const Vertex& q2)
206  {
207  if (dir >= 0)
208  return q1.c[dir] < q2.c[dir];
209  return q1.i < q2.i;
210  }
211 
213  int dir;
214  };
215 
220  BoundedPredicate(const FaceCoord& coord_) : coord(coord_) {}
221 
223  bool operator()(const Quad& q)
224  {
225  double eps = 1.e-8;
226  return (coord[0] >= q.bb[0][0]-eps &&
227  coord[0] <= q.bb[1][0]+eps &&
228  coord[1] >= q.bb[0][1]-eps &&
229  coord[1] <= q.bb[1][1]+eps);
230  }
231 
234  };
235  protected:
237  std::vector<Quad> grid;
238  std::vector<std::vector<Quad> > colGrids;
239 
241  std::vector<bool> fixNodes;
242 
244  size_t nodes;
245 
247  friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
248  {
249  for (size_t i=0;i<g.size();++i)
250  os << g[i] << std::endl;
251  return os;
252  }
253 
263  bool bilinearSolve(double epsilon, double order,
264  const double* A, const double* B,
265  std::vector<double>& X,
266  std::vector<double>& Y) const;
267 
275  bool cubicSolve(double eps, double A, double B, double C,
276  double D, std::vector<double>& X) const;
277 
284  int Q4inv(FaceCoord& res, const Quad& q, const FaceCoord& coord,
285  double epsZero, double epsOut) const;
286 };
287 
289  template<int mydim, int cdim, class GridImp>
291 {
292 };
293 
295  template<int cdim, class GridImp>
296 class HexGeometry<2, cdim, GridImp>
297 {
298  public:
300  enum { dimension = 2};
301 
303  enum { mydimension = 2};
304 
306  enum { coorddimension = cdim };
307 
309  enum {dimensionworld = 2 };
310 
312  typedef double ctype;
313 
315  typedef Dune::FieldVector<ctype,mydimension> LocalCoordinate;
316 
318  typedef Dune::FieldVector<ctype,coorddimension> GlobalCoordinate;
319 
321  typedef Dune::FieldMatrix<ctype,coorddimension,mydimension> Jacobian;
322 
324  typedef Dune::FieldMatrix<ctype,mydimension,coorddimension> JacobianTransposed;
325 
330  HexGeometry(const BoundaryGrid::Quad& q, const GridImp& gv, int dir)
331  {
332  Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridImp,
333  Dune::MCMGVertexLayout> mapper(gv);
334  typename GridImp::LeafGridView::template Codim<3>::Iterator start=gv.leafGridView().template begin<3>();
335  const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafGridView().template end<3>();
336  for (int i=0;i<4;++i) {
337  typename GridImp::LeafGridView::template Codim<3>::Iterator it=start;
338  for (; it != itend; ++it) {
339  if (mapper.index(*it) == q.v[i].i)
340  break;
341  }
342  BoundaryGrid::extract(c[i],it->geometry().corner(0),dir);
343  }
344 
345  m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
346  }
347 
351  {
352  for (int i=0;i<4;++i)
353  c[i] = q.v[i].c;
354  m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
355  }
356 
358  Dune::GeometryType type() const
359  {
360  Dune::GeometryType t;
361  t.makeCube(mydimension);
362  return t;
363  }
364 
366  int corners() const
367  {
368  return 4;
369  }
370 
372  ctype volume() const
373  {
374  return m_volume;
375  }
376 
379  {
380  LocalCoordinate local;
381  local = .5f;
382  return Global(local);
383  }
384 
387  GlobalCoordinate corner(int cor) const
388  {
389  return c[cor];
390  }
391 
395  {
396  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
397  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
398  uvw[0] -= local;
399  // Access pattern for uvw matching ordering of corners.
400  const int pat[4][2] = {{ 0, 0 },
401  { 1, 0 },
402  { 0, 1 },
403  { 1, 1 }};
404  GlobalCoordinate xyz(0.0);
405  for (int i = 0; i < 4; ++i) {
406  GlobalCoordinate corner_contrib = corner(i);
407  double factor = 1.0;
408  for (int j = 0; j < 2; ++j) {
409  factor *= uvw[pat[i][j]][j];
410  }
411  corner_contrib *= factor;
412  xyz += corner_contrib;
413  }
414  return xyz;
415  }
416 
420  {
421  const ctype epsilon = 1e-10;
422  const Dune::ReferenceElement< ctype , 2 > & refElement =
423  Dune::ReferenceElements< ctype, 2 >::general(type());
424  LocalCoordinate x = refElement.position(0,0);
425  LocalCoordinate dx;
426  do {
427  // DF^n dx^n = F^n, x^{n+1} -= dx^n
428  JacobianTransposed JT = jacobianTransposed(x);
429  GlobalCoordinate z = global(x);
430  z -= y;
431 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5)
432  Dune::Impl::FieldMatrixHelper<double>::template xTRightInvA<2, 2>(JT, z, dx );
433 #else
434  using namespace Dune::GenericGeometry;
435  MatrixHelper<DuneCoordTraits<double> >::template xTRightInvA<2, 2>(JT, z, dx );
436 #endif
437  x -= dx;
438  } while (dx.two_norm2() > epsilon*epsilon);
439  return x;
440  }
441 
444  const Dune::FieldMatrix<ctype, mydimension, coorddimension>
446  {
447  // uvw = { (1-u, 1-v, (u, v) }
448  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
449  uvw[0] -= local;
450  // Access pattern for uvw matching ordering of corners.
451  const int pat[4][3] = {{ 0, 0 },
452  { 1, 0 },
453  { 0, 1 },
454  { 1, 1 }};
455  Dune::FieldMatrix<ctype, mydimension, coorddimension> Jt(0.0);
456  for (int i = 0; i < 4; ++i) {
457  for (int deriv = 0; deriv < 2; ++deriv) {
458  // This part contributing to dg/du_{deriv}
459  double factor = 1.0;
460  for (int j = 0; j < 2; ++j) {
461  factor *= (j != deriv) ? uvw[pat[i][j]][j]
462  : (pat[i][j] == 0 ? -1.0 : 1.0);
463  }
464  GlobalCoordinate corner_contrib = corner(i);
465  corner_contrib *= factor;
466  Jt[deriv] += corner_contrib; // using FieldMatrix row access.
467  }
468  }
469  return Jt;
470  }
471 
474  const Dune::FieldMatrix<ctype, coorddimension, mydimension>
476  {
477  Dune::FieldMatrix<ctype, coorddimension, mydimension> Jti = jacobianTransposed(local);
478  Jti.invert();
479  return Jti;
480  }
481 
485  {
486  Dune::FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local);
487 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5)
488  return Dune::Impl::FieldMatrixHelper<double>::template sqrtDetAAT<2, 2>(Jt);
489 #else
490  using namespace Dune::GenericGeometry;
491  return MatrixHelper<DuneCoordTraits<double> >::template sqrtDetAAT<2, 2>(Jt);
492 #endif
493  }
494  private:
496  GlobalCoordinate c[4];
497 
499  ctype m_volume;
500 };
501 
504 BoundaryGrid::Vertex minXminY(std::vector<BoundaryGrid::Vertex>& in);
505 
508 BoundaryGrid::Vertex maxXminY(std::vector<BoundaryGrid::Vertex>& in);
509 
512 BoundaryGrid::Vertex maxXmaxY(std::vector<BoundaryGrid::Vertex>& in);
513 
516 BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in);
517 
518 }
519 }
520 
521 #endif
int dir
Compare using this direction, if -1 compare using index values.
Definition: boundarygrid.hh:213
bool isFixed(int node) const
Check if a given node is marked as fixed.
Definition: boundarygrid.hh:174
GlobalCoordinate corner(int cor) const
Returns coordinates to requested corner.
Definition: boundarygrid.hh:387
ctype integrationElement(const LocalCoordinate &local) const
Returns the integration element (|J&#39;*J|)^(1/2)
Definition: boundarygrid.hh:484
Dune::FieldVector< double, 3 > GlobalCoordinate
A coordinate on the underlying 3D grid.
Definition: boundarygrid.hh:38
LocalCoordinate local(const GlobalCoordinate &y) const
Map from global coordinates to local coordinates.
Definition: boundarygrid.hh:419
FaceCoord c
Coordinates of the vertex (2D)
Definition: boundarygrid.hh:73
std::vector< double > evalBasis(double xi, double eta) const
Evaluate the basis functions.
Definition: boundarygrid.cpp:340
FaceCoord bb[2]
Bounding box.
Definition: boundarygrid.hh:110
Dune::FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: boundarygrid.hh:324
size_t totalNodes() const
Return the total number of nodes on the grid when known.
Definition: boundarygrid.hh:166
HexGeometry(const BoundaryGrid::Quad &q, const GridImp &gv, int dir)
Construct integration element extracted from a 3D grid.
Definition: boundarygrid.hh:330
Dune::GeometryType type() const
Returns entity type (a 2D cube)
Definition: boundarygrid.hh:358
bool find(Vertex &res, const Vertex &coord) const
Locate the position of a vertex on the grid.
Definition: boundarygrid.cpp:117
Dune::FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: boundarygrid.hh:321
bool fixed
Whether or not this node is fixed.
Definition: boundarygrid.hh:79
std::vector< Quad > grid
Our quadrilateral elements.
Definition: boundarygrid.hh:237
virtual ~BoundaryGrid()
Default (empty) destructor.
Definition: boundarygrid.hh:57
bool operator==(const Vertex &v2)
Comparison operator.
Definition: boundarygrid.hh:82
Dune::FieldVector< ctype, coorddimension > GlobalCoordinate
Range type.
Definition: boundarygrid.hh:318
Quad * q
The quad containing the vertex (if search has been done)
Definition: boundarygrid.hh:76
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
const Dune::FieldMatrix< ctype, mydimension, coorddimension > jacobianTransposed(const LocalCoordinate &local) const
Return the transposed jacobian.
Definition: boundarygrid.hh:445
const Quad & operator[](int index) const
Obtain a const reference to a quad.
Definition: boundarygrid.hh:143
FaceCoord pos(double xi, double eta) const
Return the physical coordinates corresponding to the given local coordinates.
Definition: boundarygrid.cpp:329
VertexLess(int comp)
Default constructor.
Definition: boundarygrid.hh:202
bool cubicSolve(double eps, double A, double B, double C, double D, std::vector< double > &X) const
Solves the cubic equation A*x^3+B*x^2+C*x+D=0.
Definition: boundarygrid.cpp:270
A class describing a 2D vertex.
Definition: boundarygrid.hh:64
A class describing a linear, quadrilateral element.
Definition: boundarygrid.hh:88
BoundedPredicate(const FaceCoord &coord_)
Default constructor.
Definition: boundarygrid.hh:220
A class describing a quad grid.
Definition: boundarygrid.hh:35
FaceCoord coord
The coordinates to check.
Definition: boundarygrid.hh:233
ctype volume() const
Returns volume (area) of quadrilateral.
Definition: boundarygrid.hh:372
int i
Index of the vertex.
Definition: boundarygrid.hh:70
void add(const Quad &quad)
Add a quad to the grid.
Definition: boundarygrid.cpp:85
friend std::ostream & operator<<(std::ostream &os, const Quad &q)
Print to a stream.
Definition: boundarygrid.hh:113
GlobalCoordinate global(const LocalCoordinate &local) const
Map from local coordinates to global coordinates.
Definition: boundarygrid.hh:394
static BoundaryGrid uniform(const FaceCoord &min, const FaceCoord &max, int k1, int k2, bool dc=false)
Establish an uniform quad grid.
Definition: boundarygrid.cpp:24
int corners() const
Returns number of corners.
Definition: boundarygrid.hh:366
friend std::ostream & operator<<(std::ostream &os, const BoundaryGrid &g)
Print to a stream.
Definition: boundarygrid.hh:247
const Dune::FieldMatrix< ctype, coorddimension, mydimension > jacobianInverseTransposed(const LocalCoordinate &local) const
Returns the inverse, transposed Jacobian.
Definition: boundarygrid.hh:475
std::vector< bool > fixNodes
Whether or not a given node is marked as fixed.
Definition: boundarygrid.hh:241
Dune::FieldVector< ctype, mydimension > LocalCoordinate
Domain type.
Definition: boundarygrid.hh:315
GlobalCoordinate center() const
Returns center of quadrilateral.
Definition: boundarygrid.hh:378
size_t size() const
Obtain the number of quads in the grid.
Definition: boundarygrid.hh:154
size_t nodes
Total number of nodes on grid.
Definition: boundarygrid.hh:244
bool bilinearSolve(double epsilon, double order, const double *A, const double *B, std::vector< double > &X, std::vector< double > &Y) const
Solves a bi-linear set of equations in x and y.
Definition: boundarygrid.cpp:213
Predicate for sorting vertices.
Definition: boundarygrid.hh:199
Quad & operator[](int index)
Obtain a reference to a quad.
Definition: boundarygrid.hh:135
Predicate for checking if a vertex falls within a quads bounding box.
Definition: boundarygrid.hh:217
Vertex()
Default constructor.
Definition: boundarygrid.hh:67
bool operator()(const Quad &q)
The comparison operator.
Definition: boundarygrid.hh:223
Geometry class for general hexagons.
Definition: boundarygrid.hh:290
BoundaryGrid()
Default (empty) constructor.
Definition: boundarygrid.hh:54
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition: boundarygrid.hh:40
double ctype
Coordinate element type.
Definition: boundarygrid.hh:312
Vertex v[4]
Vertices.
Definition: boundarygrid.hh:108
Quad()
Default constructor.
Definition: boundarygrid.hh:91
int Q4inv(FaceCoord &res, const Quad &q, const FaceCoord &coord, double epsZero, double epsOut) const
Find the local coordinates of a given point within a given quad.
Definition: boundarygrid.cpp:150
HexGeometry(const BoundaryGrid::Quad &q)
Construct integration element.
Definition: boundarygrid.hh:350
bool operator()(const Vertex &q1, const Vertex &q2)
The comparison operator.
Definition: boundarygrid.hh:205
static void extract(FaceCoord &res, const GlobalCoordinate &coord, int dir)
Helper function for extracting given 2D coordinates from a 3D coordinate.
Definition: boundarygrid.cpp:70