My Project
FaceQuadrature.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_FACEQUADRATURE_HEADER_INCLUDED
21#define OPM_FACEQUADRATURE_HEADER_INCLUDED
22
23#include <opm/common/ErrorMacros.hpp>
25#include <cmath>
26
27namespace Opm
28{
29
30
31 namespace {
35 inline void cross(const double* a, const double* b, double* res)
36 {
37 res[0] = a[1]*b[2] - a[2]*b[1];
38 res[1] = a[2]*b[0] - a[0]*b[2];
39 res[2] = a[0]*b[1] - a[1]*b[0];
40 }
41
44 inline double triangleArea3d(const double* p0,
45 const double* p1,
46 const double* p2)
47 {
48 double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
49 double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
50 double cr[3];
51 cross(a, b, cr);
52 return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
53 }
54 } // anonymous namespace
55
56
83 {
84 public:
86 const int face,
87 const int degree)
88 : grid_(grid), face_(face), degree_(degree)
89 {
90 if (grid_.dimensions > 3) {
91 OPM_THROW(std::runtime_error, "FaceQuadrature only implemented for up to 3 dimensions.");
92 }
93 if (degree_ > 2) {
94 OPM_THROW(std::runtime_error, "FaceQuadrature exact for polynomial degrees > 2 not implemented.");
95 }
96 }
97
98 int numQuadPts() const
99 {
100 if (degree_ < 2 || grid_.dimensions < 2) {
101 return 1;
102 }
103 // Degree 2 case.
104 if (grid_.dimensions == 2) {
105 return 3;
106 } else {
107 return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
108 }
109 }
110
111 void quadPtCoord(const int index, double* coord) const
112 {
113 const int dim = grid_.dimensions;
114 const double* fc = grid_.face_centroids + dim*face_;
115 if (degree_ < 2 || dim < 2) {
116 std::copy(fc, fc + dim, coord);
117 return;
118 }
119 // Degree 2 case.
120 const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
121 const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
122 const double* nc = grid_.node_coordinates;
123 if (dim == 2) {
124 assert(nn == 2);
125 const double* pa[3] = { nc + dim*fnodes[0], fc, nc + dim*fnodes[1] };
126 std::copy(pa[index], pa[index] + dim, coord);
127 return;
128 }
129 assert(dim == 3);
130 if (index < nn) {
131 // Boundary edge midpoint.
132 const int node0 = fnodes[index];
133 const int node1 = fnodes[(index + 1)%nn];
134 for (int dd = 0; dd < dim; ++dd) {
135 coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
136 }
137 } else {
138 // Interiour edge midpoint.
139 // Recall that index is now in [nn, 2*nn).
140 const int node = fnodes[index - nn];
141 for (int dd = 0; dd < dim; ++dd) {
142 coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
143 }
144 }
145 }
146
147 double quadPtWeight(const int index) const
148 {
149 if (degree_ < 2) {
150 return grid_.face_areas[face_];
151 }
152 // Degree 2 case.
153 const int dim = grid_.dimensions;
154 if (dim == 2) {
155 const double simpsonw[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 };
156 return grid_.face_areas[face_]*simpsonw[index];
157 }
158 assert(dim == 3);
159 const double* fc = grid_.face_centroids + dim*face_;
160 const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
161 const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
162 const double* nc = grid_.node_coordinates;
163 if (index < nn) {
164 // Boundary edge midpoint.
165 const int node0 = fnodes[index];
166 const int node1 = fnodes[(index + 1)%nn];
167 const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
168 return area / 3.0;
169 } else {
170 // Interiour edge midpoint.
171 // Recall that index is now in [nn, 2*nn).
172 const int node0 = fnodes[(index - 1) % nn];
173 const int node1 = fnodes[index - nn];
174 const int node2 = fnodes[(index + 1) % nn];
175 const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
176 const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
177 return (area0 + area1) / 3.0;
178 }
179 }
180
181 private:
182 const UnstructuredGrid& grid_;
183 const int face_;
184 const int degree_;
185 };
186
187} // namespace Opm
188
189#endif // OPM_FACEQUADRATURE_HEADER_INCLUDED
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
A class providing numerical quadrature for faces.
Definition: FaceQuadrature.hpp:83
FieldVector< T, 3 > cross(const FieldVector< T, 3 > &a, const FieldVector< T, 3 > &b)
Definition: Volumes.hpp:58
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
double * face_areas
Exact or approximate face areas.
Definition: UnstructuredGrid.h:173
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