My Project
CellQuadrature.hpp
1/*
2 Copyright 2012 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_CELLQUADRATURE_HEADER_INCLUDED
21#define OPM_CELLQUADRATURE_HEADER_INCLUDED
22
24#include <opm/common/ErrorMacros.hpp>
25#include <algorithm>
26#include <cmath>
27
28namespace Opm
29{
30
31
32 namespace {
35 inline double determinantOf(const double* a0,
36 const double* a1,
37 const double* a2)
38 {
39 return
40 a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
41 a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
42 a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
43 }
44
47 inline double tetVolume(const double* p0,
48 const double* p1,
49 const double* p2,
50 const double* p3)
51 {
52 double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
53 double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
54 double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
55 return std::fabs(determinantOf(a, b, c) / 6.0);
56 }
57
60 inline double triangleArea2d(const double* p0,
61 const double* p1,
62 const double* p2)
63 {
64 double a[2] = { p1[0] - p0[0], p1[1] - p0[1] };
65 double b[2] = { p2[0] - p0[0], p2[1] - p0[1] };
66 double a_cross_b = a[0]*b[1] - a[1]*b[0];
67 return 0.5*std::fabs(a_cross_b);
68 }
69 } // anonymous namespace
70
71
120 {
121 public:
123 const int cell,
124 const int degree)
125 : grid_(grid), cell_(cell), degree_(degree)
126 {
127 if (grid.dimensions > 3) {
128 OPM_THROW(std::runtime_error, "CellQuadrature only implemented for up to 3 dimensions.");
129 }
130 if (degree > 2) {
131 OPM_THROW(std::runtime_error, "CellQuadrature exact for polynomial degrees > 1 not implemented.");
132 }
133 }
134
135 int numQuadPts() const
136 {
137 if (degree_ < 2 || grid_.dimensions == 1) {
138 return 1;
139 }
140 // Degree 2 case.
141 if (grid_.dimensions == 2) {
142 return 3*(grid_.cell_facepos[cell_ + 1] - grid_.cell_facepos[cell_]);
143 }
144 assert(grid_.dimensions == 3);
145 int sumnodes = 0;
146 for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
147 const int face = grid_.cell_faces[hf];
148 sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
149 }
150 return 4*sumnodes;
151 }
152
153 void quadPtCoord(const int index, double* coord) const
154 {
155 const int dim = grid_.dimensions;
156 const double* cc = grid_.cell_centroids + dim*cell_;
157 if (degree_ < 2) {
158 std::copy(cc, cc + dim, coord);
159 return;
160 }
161 // Degree 2 case.
162 if (dim == 2) {
163 if (index % 3 == 0) {
164 // Boundary midpoint. This is the face centroid.
165 const int hface = grid_.cell_facepos[cell_] + index/3;
166 const int face = grid_.cell_faces[hface];
167 const double* fc = grid_.face_centroids + dim*face;
168 std::copy(fc, fc + dim, coord);
169 } else {
170 // Interiour midpoint. This is the average of the
171 // cell centroid and a face node (they should
172 // always have two nodes in 2d).
173 const int hface = grid_.cell_facepos[cell_] + index/3;
174 const int face = grid_.cell_faces[hface];
175 const int nodeoff = (index % 3) - 1; // == 0 or 1
176 const int node = grid_.face_nodes[grid_.face_nodepos[face] + nodeoff];
177 const double* nc = grid_.node_coordinates + dim*node;
178 for (int dd = 0; dd < dim; ++dd) {
179 coord[dd] = 0.5*(nc[dd] + cc[dd]);
180 }
181 }
182 return;
183 }
184 assert(dim == 3);
185 int tetindex = index / 4;
186 const int subindex = index % 4;
187 const double* nc = grid_.node_coordinates;
188 for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
189 const int face = grid_.cell_faces[hf];
190 const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
191 if (nfn <= tetindex) {
192 // Our tet is not associated with this face.
193 tetindex -= nfn;
194 continue;
195 }
196 const double* fc = grid_.face_centroids + dim*face;
197 const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
198 const int node0 = fnodes[tetindex];
199 const int node1 = fnodes[(tetindex + 1) % nfn];
200 const double* n0c = nc + dim*node0;
201 const double* n1c = nc + dim*node1;
202 const double a = 0.138196601125010515179541316563436;
203 // Barycentric coordinates of our point in the tet.
204 double baryc[4] = { a, a, a, a };
205 baryc[subindex] = 1.0 - 3.0*a;
206 for (int dd = 0; dd < dim; ++dd) {
207 coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
208 }
209 return;
210 }
211 OPM_THROW(std::runtime_error, "Should never reach this point.");
212 }
213
214 double quadPtWeight(const int index) const
215 {
216 if (degree_ < 2) {
217 return grid_.cell_volumes[cell_];
218 }
219 // Degree 2 case.
220 const int dim = grid_.dimensions;
221 const double* cc = grid_.cell_centroids + dim*cell_;
222 if (dim == 2) {
223 const int hface = grid_.cell_facepos[cell_] + index/3;
224 const int face = grid_.cell_faces[hface];
225 const int* nptr = grid_.face_nodes + grid_.face_nodepos[face];
226 const double* nc0 = grid_.node_coordinates + dim*nptr[0];
227 const double* nc1 = grid_.node_coordinates + dim*nptr[1];
228 return triangleArea2d(nc0, nc1, cc)/3.0;
229 }
230 assert(dim == 3);
231 int tetindex = index / 4;
232 const double* nc = grid_.node_coordinates;
233 for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
234 const int face = grid_.cell_faces[hf];
235 const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
236 if (nfn <= tetindex) {
237 // Our tet is not associated with this face.
238 tetindex -= nfn;
239 continue;
240 }
241 const double* fc = grid_.face_centroids + dim*face;
242 const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
243 const int node0 = fnodes[tetindex];
244 const int node1 = fnodes[(tetindex + 1) % nfn];
245 const double* n0c = nc + dim*node0;
246 const double* n1c = nc + dim*node1;
247 return 0.25*tetVolume(cc, fc, n0c, n1c);
248 }
249 OPM_THROW(std::runtime_error, "Should never reach this point.");
250 }
251
252 private:
253 const UnstructuredGrid& grid_;
254 const int cell_;
255 const int degree_;
256 };
257
258} // namespace Opm
259
260#endif // OPM_CELLQUADRATURE_HEADER_INCLUDED
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
A class providing numerical quadrature for cells.
Definition: CellQuadrature.hpp:120
Holds the implementation of the CpGrid as a pimple.
Definition: CellQuadrature.hpp:29
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 * cell_faces
Contains for each cell, the indices of its adjacent faces.
Definition: UnstructuredGrid.h:146
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
double * face_centroids
Exact or approximate face centroids, stored consecutively for each face.
Definition: UnstructuredGrid.h:168
int dimensions
The topological and geometrical dimensionality of the grid.
Definition: UnstructuredGrid.h:106
double * node_coordinates
Node coordinates, stored consecutively for each node.
Definition: UnstructuredGrid.h:160