My Project
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 #include <dune/grid/common/mcmgmapper.hh>
22 #include <dune/grid/common/defaultgridview.hh>
23 
24 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
25 
26 
27 #include <vector>
28 
29 namespace Opm {
30 namespace Elasticity {
31 
33 class BoundaryGrid {
34  public:
36  typedef Dune::FieldVector<double,3> GlobalCoordinate;
38  typedef Dune::FieldVector<double,2> FaceCoord;
39 
47  // in natural order along the first direction
48  static BoundaryGrid uniform(const FaceCoord& min, const FaceCoord& max,
49  int k1, int k2, bool dc=false);
50 
53 
55  virtual ~BoundaryGrid() {}
56 
58  // on a boundary
59  class Quad;
60 
62  class Vertex {
63  public:
65  Vertex() : i(-1), c(0), q(0), fixed(false) {}
66 
68  int i;
69 
72 
74  Quad* q;
75 
77  bool fixed;
78 
80  bool operator==(const Vertex& v2)
81  {
82  return hypot(v2.c[0]-c[0],v2.c[1]-c[1]) < 1.e-8;
83  }
84  };
86  class Quad {
87  public:
89  Quad()
90  {
91  v[0].i = v[1].i = v[2].i = v[3].i = 0;
92  v[0].c = v[1].c = v[2].c = v[3].c = 0.f;
93  }
98  FaceCoord pos(double xi, double eta) const;
99 
103  std::vector<double> evalBasis(double xi, double eta) const;
104 
106  Vertex v[4];
109  protected:
111  friend std::ostream& operator <<(std::ostream& os, const Quad& q)
112  {
113  os << "Nodes: " << q.v[0].i << "," << q.v[1].i << "," << q.v[2].i << "," << q.v[3].i << std::endl;
114  os << "(" << q.v[0].c << ")(" << q.v[1].c << ")(" << q.v[2].c << ")(" << q.v[3].c << ")";
115  return os;
116  }
117  };
118 
121  void add(const Quad& quad);
122 
123  void addToColumn(size_t col, const Quad& quad)
124  {
125  if (col >= colGrids.size())
126  colGrids.resize(col+1);
127  colGrids[col].push_back(quad);
128  }
129 
133  Quad& operator[](int index)
134  {
135  return grid[index];
136  }
137 
141  const Quad& operator[](int index) const
142  {
143  return grid[index];
144  }
145 
146  const Quad& getQuad(int col, int index) const
147  {
148  return colGrids[col][index];
149  }
150 
152  size_t size() const
153  {
154  return grid.size();
155  }
156 
157  size_t colSize(int i) const
158  {
159  return colGrids[i].size();
160  }
161 
164  size_t totalNodes() const
165  {
166  return nodes;
167  }
168 
172  bool isFixed(int node) const
173  {
174  return fixNodes[node];
175  }
176 
180  bool find(Vertex& res, const Vertex& coord) const;
181 
186  static void extract(FaceCoord& res,
187  const GlobalCoordinate& coord, int dir);
188 
193  static void extract(Vertex& res,
194  const Dune::FieldVector<double,3>& coord, int dir);
195 
197  struct VertexLess {
200  VertexLess(int comp) : dir(comp) {}
201 
203  bool operator()(const Vertex& q1, const Vertex& q2)
204  {
205  if (dir >= 0)
206  return q1.c[dir] < q2.c[dir];
207  return q1.i < q2.i;
208  }
209 
211  int dir;
212  };
213 
218  BoundedPredicate(const FaceCoord& coord_) : coord(coord_) {}
219 
221  bool operator()(const Quad& q)
222  {
223  double eps = 1.e-8;
224  return (coord[0] >= q.bb[0][0]-eps &&
225  coord[0] <= q.bb[1][0]+eps &&
226  coord[1] >= q.bb[0][1]-eps &&
227  coord[1] <= q.bb[1][1]+eps);
228  }
229 
232  };
233  protected:
235  std::vector<Quad> grid;
236  std::vector<std::vector<Quad> > colGrids;
237 
239  std::vector<bool> fixNodes;
240 
242  size_t nodes;
243 
245  friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
246  {
247  for (size_t i=0;i<g.size();++i)
248  os << g[i] << std::endl;
249  return os;
250  }
251 
261  bool bilinearSolve(double epsilon, double order,
262  const double* A, const double* B,
263  std::vector<double>& X,
264  std::vector<double>& Y) const;
265 
273  bool cubicSolve(double eps, double A, double B, double C,
274  double D, std::vector<double>& X) const;
275 
282  int Q4inv(FaceCoord& res, const Quad& q, const FaceCoord& coord,
283  double epsZero, double epsOut) const;
284 };
285 
287  template<int mydim, int cdim, class GridImp>
289 {
290 };
291 
293  template<int cdim, class GridImp>
294 class HexGeometry<2, cdim, GridImp>
295 {
296  public:
298  enum { dimension = 2};
299 
301  enum { mydimension = 2};
302 
304  enum { coorddimension = cdim };
305 
307  enum {dimensionworld = 2 };
308 
310  typedef double ctype;
311 
313  typedef Dune::FieldVector<ctype,mydimension> LocalCoordinate;
314 
316  typedef Dune::FieldVector<ctype,coorddimension> GlobalCoordinate;
317 
319  typedef Dune::FieldMatrix<ctype,coorddimension,mydimension> Jacobian;
320 
322  typedef Dune::FieldMatrix<ctype,mydimension,coorddimension> JacobianTransposed;
323 
328  HexGeometry(const BoundaryGrid::Quad& q, const GridImp& gv, int dir)
329  {
330  using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridImp>>;
331  Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
332  typename LeafGridView::template Codim<3>::Iterator start=gv.leafGridView().template begin<3>();
333  const typename LeafGridView::template Codim<3>::Iterator itend = gv.leafGridView().template end<3>();
334  for (int i=0;i<4;++i) {
335  typename LeafGridView::template Codim<3>::Iterator it=start;
336  for (; it != itend; ++it) {
337  if (mapper.index(*it) == q.v[i].i)
338  break;
339  }
340  BoundaryGrid::extract(c[i],it->geometry().corner(0),dir);
341  }
342 
343  m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
344  }
345 
349  {
350  for (int i=0;i<4;++i)
351  c[i] = q.v[i].c;
352  m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
353  }
354 
356  Dune::GeometryType type() const
357  {
358  Dune::GeometryType t;
359  t = Dune::GeometryTypes::cube(mydimension);
360  return t;
361  }
362 
364  int corners() const
365  {
366  return 4;
367  }
368 
370  ctype volume() const
371  {
372  return m_volume;
373  }
374 
377  {
378  LocalCoordinate local;
379  local = .5f;
380  return Global(local);
381  }
382 
385  GlobalCoordinate corner(int cor) const
386  {
387  return c[cor];
388  }
389 
393  {
394  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
395  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
396  uvw[0] -= local;
397  // Access pattern for uvw matching ordering of corners.
398  const int pat[4][2] = {{ 0, 0 },
399  { 1, 0 },
400  { 0, 1 },
401  { 1, 1 }};
402  GlobalCoordinate xyz(0.0);
403  for (int i = 0; i < 4; ++i) {
404  GlobalCoordinate corner_contrib = corner(i);
405  double factor = 1.0;
406  for (int j = 0; j < 2; ++j) {
407  factor *= uvw[pat[i][j]][j];
408  }
409  corner_contrib *= factor;
410  xyz += corner_contrib;
411  }
412  return xyz;
413  }
414 
418  {
419  const ctype epsilon = 1e-10;
420  auto refElement = Dune::ReferenceElements<ctype,2>::cube();
421  LocalCoordinate x = refElement.position(0,0);
422  LocalCoordinate dx;
423  do {
424  // DF^n dx^n = F^n, x^{n+1} -= dx^n
425  JacobianTransposed JT = jacobianTransposed(x);
426  GlobalCoordinate z = global(x);
427  z -= y;
428  Dune::Impl::FieldMatrixHelper<double>::template xTRightInvA<2, 2>(JT, z, dx );
429  x -= dx;
430  } while (dx.two_norm2() > epsilon*epsilon);
431  return x;
432  }
433 
436  const Dune::FieldMatrix<ctype, mydimension, coorddimension>
438  {
439  // uvw = { (1-u, 1-v, (u, v) }
440  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
441  uvw[0] -= local;
442  // Access pattern for uvw matching ordering of corners.
443  const int pat[4][3] = {{ 0, 0 },
444  { 1, 0 },
445  { 0, 1 },
446  { 1, 1 }};
447  Dune::FieldMatrix<ctype, mydimension, coorddimension> Jt(0.0);
448  for (int i = 0; i < 4; ++i) {
449  for (int deriv = 0; deriv < 2; ++deriv) {
450  // This part contributing to dg/du_{deriv}
451  double factor = 1.0;
452  for (int j = 0; j < 2; ++j) {
453  factor *= (j != deriv) ? uvw[pat[i][j]][j]
454  : (pat[i][j] == 0 ? -1.0 : 1.0);
455  }
456  GlobalCoordinate corner_contrib = corner(i);
457  corner_contrib *= factor;
458  Jt[deriv] += corner_contrib; // using FieldMatrix row access.
459  }
460  }
461  return Jt;
462  }
463 
466  const Dune::FieldMatrix<ctype, coorddimension, mydimension>
468  {
469  Dune::FieldMatrix<ctype, coorddimension, mydimension> Jti = jacobianTransposed(local);
470  Jti.invert();
471  return Jti;
472  }
473 
477  {
478  Dune::FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local);
479  return Dune::Impl::FieldMatrixHelper<double>::template sqrtDetAAT<2, 2>(Jt);
480  }
481 
482  private:
484  GlobalCoordinate c[4];
485 
487  ctype m_volume;
488 };
489 
492 BoundaryGrid::Vertex minXminY(std::vector<BoundaryGrid::Vertex>& in);
493 
496 BoundaryGrid::Vertex maxXminY(std::vector<BoundaryGrid::Vertex>& in);
497 
500 BoundaryGrid::Vertex maxXmaxY(std::vector<BoundaryGrid::Vertex>& in);
501 
504 BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in);
505 
506 }
507 }
508 
509 #endif
A class describing a linear, quadrilateral element.
Definition: boundarygrid.hh:86
friend std::ostream & operator<<(std::ostream &os, const Quad &q)
Print to a stream.
Definition: boundarygrid.hh:111
Quad()
Default constructor.
Definition: boundarygrid.hh:89
Vertex v[4]
Vertices.
Definition: boundarygrid.hh:106
FaceCoord bb[2]
Bounding box.
Definition: boundarygrid.hh:108
std::vector< double > evalBasis(double xi, double eta) const
Evaluate the basis functions.
Definition: boundarygrid.cpp:340
FaceCoord pos(double xi, double eta) const
Return the physical coordinates corresponding to the given local coordinates.
Definition: boundarygrid.cpp:329
A class describing a 2D vertex.
Definition: boundarygrid.hh:62
bool operator==(const Vertex &v2)
Comparison operator.
Definition: boundarygrid.hh:80
Quad * q
The quad containing the vertex (if search has been done)
Definition: boundarygrid.hh:74
bool fixed
Whether or not this node is fixed.
Definition: boundarygrid.hh:77
FaceCoord c
Coordinates of the vertex (2D)
Definition: boundarygrid.hh:71
int i
Index of the vertex.
Definition: boundarygrid.hh:68
Vertex()
Default constructor.
Definition: boundarygrid.hh:65
A class describing a quad grid.
Definition: boundarygrid.hh:33
std::vector< Quad > grid
Our quadrilateral elements.
Definition: boundarygrid.hh:235
std::vector< bool > fixNodes
Whether or not a given node is marked as fixed.
Definition: boundarygrid.hh:239
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
Quad & operator[](int index)
Obtain a reference to a quad.
Definition: boundarygrid.hh:133
const Quad & operator[](int index) const
Obtain a const reference to a quad.
Definition: boundarygrid.hh:141
virtual ~BoundaryGrid()
Default (empty) destructor.
Definition: boundarygrid.hh:55
Dune::FieldVector< double, 3 > GlobalCoordinate
A coordinate on the underlying 3D grid.
Definition: boundarygrid.hh:36
size_t totalNodes() const
Return the total number of nodes on the grid when known.
Definition: boundarygrid.hh:164
size_t nodes
Total number of nodes on grid.
Definition: boundarygrid.hh:242
BoundaryGrid()
Default (empty) constructor.
Definition: boundarygrid.hh:52
void add(const Quad &quad)
Add a quad to the grid.
Definition: boundarygrid.cpp:85
bool find(Vertex &res, const Vertex &coord) const
Locate the position of a vertex on the grid.
Definition: boundarygrid.cpp:117
friend std::ostream & operator<<(std::ostream &os, const BoundaryGrid &g)
Print to a stream.
Definition: boundarygrid.hh:245
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
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
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition: boundarygrid.hh:38
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
bool isFixed(int node) const
Check if a given node is marked as fixed.
Definition: boundarygrid.hh:172
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
size_t size() const
Obtain the number of quads in the grid.
Definition: boundarygrid.hh:152
Dune::FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: boundarygrid.hh:319
HexGeometry(const BoundaryGrid::Quad &q, const GridImp &gv, int dir)
Construct integration element extracted from a 3D grid.
Definition: boundarygrid.hh:328
int corners() const
Returns number of corners.
Definition: boundarygrid.hh:364
const Dune::FieldMatrix< ctype, coorddimension, mydimension > jacobianInverseTransposed(const LocalCoordinate &local) const
Returns the inverse, transposed Jacobian.
Definition: boundarygrid.hh:467
GlobalCoordinate center() const
Returns center of quadrilateral.
Definition: boundarygrid.hh:376
LocalCoordinate local(const GlobalCoordinate &y) const
Map from global coordinates to local coordinates.
Definition: boundarygrid.hh:417
Dune::FieldVector< ctype, mydimension > LocalCoordinate
Domain type.
Definition: boundarygrid.hh:313
Dune::GeometryType type() const
Returns entity type (a 2D cube)
Definition: boundarygrid.hh:356
ctype volume() const
Returns volume (area) of quadrilateral.
Definition: boundarygrid.hh:370
double ctype
Coordinate element type.
Definition: boundarygrid.hh:310
Dune::FieldVector< ctype, coorddimension > GlobalCoordinate
Range type.
Definition: boundarygrid.hh:316
Dune::FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: boundarygrid.hh:322
const Dune::FieldMatrix< ctype, mydimension, coorddimension > jacobianTransposed(const LocalCoordinate &local) const
Return the transposed jacobian.
Definition: boundarygrid.hh:437
GlobalCoordinate corner(int cor) const
Returns coordinates to requested corner.
Definition: boundarygrid.hh:385
ctype integrationElement(const LocalCoordinate &local) const
Returns the integration element (|J'*J|)^(1/2)
Definition: boundarygrid.hh:476
HexGeometry(const BoundaryGrid::Quad &q)
Construct integration element.
Definition: boundarygrid.hh:348
GlobalCoordinate global(const LocalCoordinate &local) const
Map from local coordinates to global coordinates.
Definition: boundarygrid.hh:392
Geometry class for general hexagons.
Definition: boundarygrid.hh:289
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
Predicate for checking if a vertex falls within a quads bounding box.
Definition: boundarygrid.hh:215
bool operator()(const Quad &q)
The comparison operator.
Definition: boundarygrid.hh:221
FaceCoord coord
The coordinates to check.
Definition: boundarygrid.hh:231
BoundedPredicate(const FaceCoord &coord_)
Default constructor.
Definition: boundarygrid.hh:218
Predicate for sorting vertices.
Definition: boundarygrid.hh:197
VertexLess(int comp)
Default constructor.
Definition: boundarygrid.hh:200
int dir
Compare using this direction, if -1 compare using index values.
Definition: boundarygrid.hh:211
bool operator()(const Vertex &q1, const Vertex &q2)
The comparison operator.
Definition: boundarygrid.hh:203