20 #ifndef OPM_CELLQUADRATURE_HEADER_INCLUDED
21 #define OPM_CELLQUADRATURE_HEADER_INCLUDED
24 #include <opm/grid/utility/ErrorMacros.hpp>
35 inline double determinantOf(
const double* a0,
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]);
47 inline double tetVolume(
const double* p0,
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);
60 inline double triangleArea2d(
const double* p0,
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);
125 : grid_(grid), cell_(cell), degree_(degree)
128 OPM_THROW(std::runtime_error,
"CellQuadrature only implemented for up to 3 dimensions.");
131 OPM_THROW(std::runtime_error,
"CellQuadrature exact for polynomial degrees > 1 not implemented.");
135 int numQuadPts()
const
153 void quadPtCoord(
const int index,
double* coord)
const
158 std::copy(cc, cc + dim, coord);
163 if (index % 3 == 0) {
168 std::copy(fc, fc + dim, coord);
175 const int nodeoff = (index % 3) - 1;
178 for (
int dd = 0; dd < dim; ++dd) {
179 coord[dd] = 0.5*(nc[dd] + cc[dd]);
185 int tetindex = index / 4;
186 const int subindex = index % 4;
191 if (nfn <= tetindex) {
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;
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];
211 OPM_THROW(std::runtime_error,
"Should never reach this point.");
214 double quadPtWeight(
const int index)
const
228 return triangleArea2d(nc0, nc1, cc)/3.0;
231 int tetindex = index / 4;
236 if (nfn <= tetindex) {
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);
249 OPM_THROW(std::runtime_error,
"Should never reach this point.");
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