My Project
GridAdapter.hpp
1 /*
2  Copyright 2010 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_GRIDADAPTER_HEADER_INCLUDED
21 #define OPM_GRIDADAPTER_HEADER_INCLUDED
22 
23 
25 #include <opm/grid/utility/ErrorMacros.hpp>
26 #include <opm/grid/CpGrid.hpp>
27 #include <stdexcept>
28 #include <vector>
29 
31 {
32 public:
37  template <class Grid>
38  void init(const Grid& grid)
39  {
40  buildTopology(grid);
41  buildGeometry(grid);
42  buildGlobalCell(grid);
43  copyCartDims(grid);
44  }
45 
46  UnstructuredGrid* c_grid()
47  {
48  return &g_;
49  }
51  const UnstructuredGrid* c_grid() const
52  {
53  return &g_;
54  }
55 
56  // ------ Forwarding the same interface that init() expects ------
57  //
58  // This is only done in order to verify that init() works correctly.
59  //
60 
61  enum { dimension = 3 }; // This is actually a hack used for testing (dim is a runtime parameter).
62 
63  struct Vector
64  {
65  explicit Vector(const double* source)
66  {
67  for (int i = 0; i < dimension; ++i) {
68  data[i] = source[i];
69  }
70  }
71  double& operator[] (const int ix)
72  {
73  return data[ix];
74  }
75  double operator[] (const int ix) const
76  {
77  return data[ix];
78  }
79  double data[dimension];
80  };
81 
82  // Topology
83  int numCells() const
84  {
85  return g_.number_of_cells;
86  }
87  int numFaces() const
88  {
89  return g_.number_of_faces;
90  }
91  int numVertices() const
92  {
93  return g_.number_of_nodes;
94  }
95 
96  int numCellFaces(int cell) const
97  {
98  return cell_facepos_[cell + 1] - cell_facepos_[cell];
99  }
100  int cellFace(int cell, int local_index) const
101  {
102  return cell_faces_[cell_facepos_[cell] + local_index];
103  }
104  int faceCell(int face, int local_index) const
105  {
106  return face_cells_[2*face + local_index];
107  }
108  int numFaceVertices(int face) const
109  {
110  return face_nodepos_[face + 1] - face_nodepos_[face];
111  }
112  int faceVertex(int face, int local_index) const
113  {
114  return face_nodes_[face_nodepos_[face] + local_index];
115  }
116 
117  // Geometry
118  Vector vertexPosition(int vertex) const
119  {
120  return Vector(&node_coordinates_[g_.dimensions*vertex]);
121  }
122  double faceArea(int face) const
123  {
124  return face_areas_[face];
125  }
126  Vector faceCentroid(int face) const
127  {
128  return Vector(&face_centroids_[g_.dimensions*face]);
129  }
130  Vector faceNormal(int face) const
131  {
132  Vector fn(&face_normals_[g_.dimensions*face]);
133  // We must renormalize since the stored normals are
134  // 'unit normal * face area'.
135  double invfa = 1.0 / faceArea(face);
136  for (int i = 0; i < dimension; ++i) {
137  fn[i] *= invfa;
138  }
139  return fn;
140  }
141  double cellVolume(int cell) const
142  {
143  return cell_volumes_[cell];
144  }
145  Vector cellCentroid(int cell) const
146  {
147  return Vector(&cell_centroids_[g_.dimensions*cell]);
148  }
149 
150  bool operator==(const GridAdapter& other)
151  {
152  return face_nodes_ == other.face_nodes_
153  && face_nodepos_ == other.face_nodepos_
154  && face_cells_ == other.face_cells_
155  && cell_faces_ == other.cell_faces_
156  && cell_facepos_ == other.cell_facepos_
157  && node_coordinates_ == other.node_coordinates_
158  && face_centroids_ == other.face_centroids_
159  && face_areas_ == other.face_areas_
160  && face_normals_ == other.face_normals_
161  && cell_centroids_ == other.cell_centroids_
162  && cell_volumes_ == other.cell_volumes_;
163  }
164  // make a grid which looks periodic but do not have 2 half faces for each
165  // periodic boundary
166  void makeQPeriodic(const std::vector<int>& hf_ind,const std::vector<int>& periodic_cells){
167  for(int i=0; i<int(hf_ind.size());++i){
168  //std::array<int,2> cells;
169  int& cell0=face_cells_[2*cell_faces_[ hf_ind[i] ]+0];
170  int& cell1=face_cells_[2*cell_faces_[ hf_ind[i] ]+1];
171  assert(periodic_cells[2*i+1]>=0);
172  if(periodic_cells[2*i+0] == cell0){
173  assert(cell1==-1);
174  cell1=periodic_cells[2*i+1];
175  }else{
176  assert(periodic_cells[2*i+0] == cell1);
177  assert(cell0==-1);
178  cell0=periodic_cells[2*i+1];
179  }
180  }
181  }
182 private:
183  UnstructuredGrid g_;
184  // Topology storage.
185  std::vector<int> face_nodes_;
186  std::vector<int> face_nodepos_;
187  std::vector<int> face_cells_;
188  std::vector<int> cell_faces_;
189  std::vector<int> cell_facepos_;
190  // Geometry storage.
191  std::vector<double> node_coordinates_;
192  std::vector<double> face_centroids_;
193  std::vector<double> face_areas_;
194  std::vector<double> face_normals_;
195  std::vector<double> cell_centroids_;
196  std::vector<double> cell_volumes_;
197  // The global cell information
198  std::vector<int> global_cell_;
200  void buildGlobalCell(const Dune::CpGrid& grid)
201  {
202  bool all_active=true;
203  int old_cell=-1;
204  global_cell_.resize(grid.numCells());
205  for(int c=0; c<grid.numCells(); ++c)
206  {
207  int new_cell=global_cell_[c]=grid.globalCell()[c];
208  all_active = all_active && (new_cell==old_cell+1);
209  old_cell=new_cell;
210  }
211  if(all_active){
212  g_.global_cell=0;
213  // really release memory
214  std::vector<int>().swap(global_cell_);
215  }
216  else
217  g_.global_cell=&(global_cell_[0]);
218  }
220  template<class Grid>
221  void buildGlobalCell(const Grid& /*grid*/)
222  {
223  g_.global_cell=0;
224  }
226  void copyCartDims(const Dune::CpGrid& grid)
227  {
228  for(int i=0; i<3; ++i)
229  g_.cartdims[i] = grid.logicalCartesianSize()[i];
230  }
232  template <class Grid>
233  void copyCartDims(const Grid& /*grid*/)
234  {}
236  template <class Grid>
237  void buildTopology(const Grid& grid)
238  {
239  // Face topology.
240  int num_cells = grid.numCells();
241  int num_faces = grid.numFaces();
242  face_nodepos_.resize(num_faces + 1);
243  int facenodecount = 0;
244  for (int f = 0; f < num_faces; ++f) {
245  face_nodepos_[f] = facenodecount;
246  facenodecount += grid.numFaceVertices(f);
247  }
248  face_nodepos_.back() = facenodecount;
249  face_nodes_.resize(facenodecount);
250  for (int f = 0; f < num_faces; ++f) {
251  for (int local = 0; local < grid.numFaceVertices(f); ++local) {
252  face_nodes_[face_nodepos_[f] + local] = grid.faceVertex(f, local);
253  }
254  }
255  face_cells_.resize(2*num_faces);
256  for (int f = 0; f < num_faces; ++f) {
257  face_cells_[2*f] = grid.faceCell(f, 0);
258  face_cells_[2*f + 1] = grid.faceCell(f, 1);
259  }
260 
261  // Cell topology.
262  int cellfacecount = 0;
263  cell_facepos_.resize(num_cells + 1);
264  for (int c = 0; c < num_cells; ++c) {
265  cell_facepos_[c] = cellfacecount;
266  cellfacecount += grid.numCellFaces(c);
267  }
268  cell_facepos_.back() = cellfacecount;
269  cell_faces_.resize(cellfacecount);
270  for (int c = 0; c < num_cells; ++c) {
271  for (int local = 0; local < grid.numCellFaces(c); ++local) {
272  cell_faces_[cell_facepos_[c] + local] = grid.cellFace(c, local);
273  }
274  }
275 
276  // Set C grid members.
277  g_.dimensions = Grid::dimension;
278  g_.number_of_cells = grid.numCells();
279  g_.number_of_faces = grid.numFaces();
280  g_.number_of_nodes = grid.numVertices();
281  g_.face_nodes = &face_nodes_[0];
282  g_.face_nodepos = &face_nodepos_[0];
283  g_.face_cells = &face_cells_[0];
284  g_.cell_faces = &cell_faces_[0];
285  g_.cell_facepos = &cell_facepos_[0];
286  }
287 
288 
291  template <class Grid>
292  void buildGeometry(const Grid& grid)
293  {
294  // Node geometry.
295  int num_cells = grid.numCells();
296  int num_nodes = grid.numVertices();
297  int num_faces = grid.numFaces();
298  int dim = Grid::dimension;
299  node_coordinates_.resize(dim*num_nodes);
300  for (int n = 0; n < num_nodes; ++n) {
301  for (int dd = 0; dd < dim; ++dd) {
302  node_coordinates_[dim*n + dd] = grid.vertexPosition(n)[dd];
303  }
304  }
305 
306  // Face geometry.
307  face_centroids_.resize(dim*num_faces);
308  face_areas_.resize(num_faces);
309  face_normals_.resize(dim*num_faces);
310  for (int f = 0; f < num_faces; ++f) {
311  face_areas_[f] = grid.faceArea(f);
312  for (int dd = 0; dd < dim; ++dd) {
313  face_centroids_[dim*f + dd] = grid.faceCentroid(f)[dd];
314  face_normals_[dim*f + dd] = grid.faceNormal(f)[dd]*face_areas_[f];
315  }
316  }
317 
318  // Cell geometry.
319  cell_centroids_.resize(dim*num_cells);
320  cell_volumes_.resize(num_cells);
321  for (int c = 0; c < num_cells; ++c) {
322  cell_volumes_[c] = grid.cellVolume(c);
323  for (int dd = 0; dd < dim; ++dd) {
324  cell_centroids_[dim*c + dd] = grid.cellCentroid(c)[dd];
325  }
326  }
327 
328  // Set C grid members.
329  g_.node_coordinates = &node_coordinates_[0];
330  g_.face_centroids = &face_centroids_[0];
331  g_.face_areas = &face_areas_[0];
332  g_.face_normals = &face_normals_[0];
333  g_.cell_centroids = &cell_centroids_[0];
334  g_.cell_volumes = &cell_volumes_[0];
335  }
336 };
337 
338 
339 #endif // OPM_GRIDADAPTER_HEADER_INCLUDED
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
[ provides Dune::Grid ]
Definition: CpGrid.hpp:207
Definition: GridAdapter.hpp:31
void init(const Grid &grid)
Initialize the grid.
Definition: GridAdapter.hpp:38
const UnstructuredGrid * c_grid() const
Access the underlying C grid.
Definition: GridAdapter.hpp:51
Definition: GridAdapter.hpp:64
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition: UnstructuredGrid.h:99
int * face_nodes
Contains for each face, the indices of its adjacent nodes.
Definition: UnstructuredGrid.h:121
int * face_nodepos
For a face f, face_nodepos[f] contains the starting index for f's nodes in the face_nodes array.
Definition: UnstructuredGrid.h:127
int number_of_cells
The number of cells in the grid.
Definition: UnstructuredGrid.h:109
double * face_areas
Exact or approximate face areas.
Definition: UnstructuredGrid.h:173
int * cell_faces
Contains for each cell, the indices of its adjacent faces.
Definition: UnstructuredGrid.h:146
int number_of_faces
The number of faces in the grid.
Definition: UnstructuredGrid.h:111
double * cell_centroids
Exact or approximate cell centroids, stored consecutively for each cell.
Definition: UnstructuredGrid.h:192
int * cell_facepos
For a cell c, cell_facepos[c] contains the starting index for c's faces in the cell_faces array.
Definition: UnstructuredGrid.h:152
double * cell_volumes
Exact or approximate cell volumes.
Definition: UnstructuredGrid.h:197
int * face_cells
For a face f, face_cells[2*f] and face_cells[2*f + 1] contain the cell indices of the cells adjacent ...
Definition: UnstructuredGrid.h:138
double * face_centroids
Exact or approximate face centroids, stored consecutively for each face.
Definition: UnstructuredGrid.h:168
int number_of_nodes
The number of nodes in the grid.
Definition: UnstructuredGrid.h:113
int dimensions
The topological and geometrical dimensionality of the grid.
Definition: UnstructuredGrid.h:106
int cartdims[3]
Contains the size of the logical cartesian structure (if any) of the grid.
Definition: UnstructuredGrid.h:227
int * global_cell
If non-null, this array contains the logical cartesian indices (in a lexicographic ordering) of each ...
Definition: UnstructuredGrid.h:214
double * face_normals
Exact or approximate face normals, stored consecutively for each face.
Definition: UnstructuredGrid.h:184
double * node_coordinates
Node coordinates, stored consecutively for each node.
Definition: UnstructuredGrid.h:160