1 #include <opm/core/linalg/blas_lapack.h>
3 #include <opm/core/grid/GridHelpers.hpp>
14 namespace UgGridHelpers
16 int dimensions(
const Dune::CpGrid&);
18 double faceArea(
const Dune::CpGrid&,
int);
25 inline const double* multiplyFaceNormalWithArea(
const Dune::CpGrid& grid,
int face_index,
const double* in)
27 int d=Opm::UgGridHelpers::dimensions(grid);
28 double* out=
new double[d];
29 double area=Opm::UgGridHelpers::faceArea(grid, face_index);
36 inline void maybeFreeFaceNormal(
const Dune::CpGrid&,
const double* array)
40 #endif // HAVE_OPM_GRID
42 inline const double* multiplyFaceNormalWithArea(
const UnstructuredGrid&,
int,
const double* in)
47 inline void maybeFreeFaceNormal(
const UnstructuredGrid&,
const double*)
59 using namespace Opm::UgGridHelpers;
61 double s, dist, denom;
64 typename CellCentroidTraits<Grid>::IteratorType cc = beginCellCentroids(*G);
65 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
66 typename FaceCellTraits<Grid>::Type face_cells = faceCells(*G);
71 MAT_SIZE_T nrows, ncols, ldA, incx, incy;
76 nrows = ncols = ldA = d;
80 for (
int c =0, i = 0; c < numCells(*G); c++) {
81 K = perm + (c * d * d);
83 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
84 FaceRow faces = c2f[c];
86 for(
typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
89 s = 2.0*(face_cells(*f, 0) == c) - 1.0;
90 n = faceNormal(*G, *f);
91 const double* nn=multiplyFaceNormalWithArea(*G, *f, n);
92 const double* fc = &(faceCentroid(*G, *f)[0]);
93 dgemv_(
"No Transpose", &nrows, &ncols,
94 &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy);
95 maybeFreeFaceNormal(*G, nn);
97 htrans[i] = denom = 0.0;
98 for (j = 0; j < d; j++) {
99 dist = fc[j] - getCoordinate(cc, j);
101 htrans[i] += s * dist * Kn[j];
102 denom += dist * dist;
107 htrans[i] = std::abs(htrans[i]);
110 cc = increment(cc, 1, d);
121 using namespace Opm::UgGridHelpers;
123 for (
int f = 0; f < numFaces(*G); f++) {
127 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
129 for (
int c = 0, i = 0; c < numCells(*G); c++) {
130 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
131 FaceRow faces = c2f[c];
133 for(
typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
136 trans[*f] += 1.0 / htrans[i];
140 for (
int f = 0; f < numFaces(*G); f++) {
141 trans[f] = 1.0 / trans[f];
150 const double *totmob,
151 const double *htrans,
155 using namespace Opm::UgGridHelpers;
157 for (
int f = 0; f < numFaces(*G); f++) {
161 typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
163 for (
int c = 0, i = 0; c < numCells(*G); c++) {
164 typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
165 FaceRow faces = c2f[c];
167 for(
typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
170 trans[*f] += 1.0 / (totmob[c] * htrans[i]);
174 for (
int f = 0; f < numFaces(*G); f++) {
175 trans[f] = 1.0 / trans[f];
void tpfa_trans_compute(struct UnstructuredGrid *G, const double *htrans, double *trans)
Compute two-point transmissibilities from one-sided transmissibilities.
Routines to assist in the calculation of two-point transmissibilities.
void tpfa_eff_trans_compute(struct UnstructuredGrid *G, const double *totmob, const double *htrans, double *trans)
Calculate effective two-point transmissibilities from one-sided, total mobility weighted, transmissibilities.
Definition: trans_tpfa.c:103
void tpfa_htrans_compute(struct UnstructuredGrid *G, const double *perm, double *htrans)
Calculate static, one-sided transmissibilities for use in the two-point flux approximation method...