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{
32public:
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 }
182private:
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:210
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