My Project
TransTpfa_impl.hpp
1 #include <opm/common/utility/numeric/blas_lapack.h>
3 #include <opm/grid/GridHelpers.hpp>
4 
5 #include <cmath>
6 
7 namespace Dune
8 {
9 class CpGrid;
10 }
11 
12 namespace Opm
13 {
14 namespace UgGridHelpers
15 {
16 int dimensions(const Dune::CpGrid&);
17 
18 double faceArea(const Dune::CpGrid&, int);
19 }
20 }
21 
22 namespace
23 {
24 inline const double* multiplyFaceNormalWithArea(const Dune::CpGrid& grid, int face_index, const double* in)
25 {
26  int d=Opm::UgGridHelpers::dimensions(grid);
27  double* out=new double[d];
28  double area=Opm::UgGridHelpers::faceArea(grid, face_index);
29 
30  for(int i=0;i<d;++i)
31  out[i]=in[i]*area;
32  return out;
33 }
34 
35 inline void maybeFreeFaceNormal(const Dune::CpGrid&, const double* array)
36 {
37  delete[] array;
38 }
39 
40 inline const double* multiplyFaceNormalWithArea(const UnstructuredGrid&, int, const double* in)
41 {
42  return in;
43 }
44 
45 inline void maybeFreeFaceNormal(const UnstructuredGrid&, const double*)
46 {}
47 }
48 
49 /* ---------------------------------------------------------------------- */
50 /* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
51 /* ---------------------------------------------------------------------- */
52 template<class Grid>
53 void
54 tpfa_htrans_compute(const Grid* G, const double *perm, double *htrans)
55 /* ---------------------------------------------------------------------- */
56 {
57  using namespace Opm::UgGridHelpers;
58  int d, j;
59  double s, dist, denom;
60 
61  double Kn[3];
62  typename CellCentroidTraits<Grid>::IteratorType cc = beginCellCentroids(*G);
63  typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
64  typename FaceCellTraits<Grid>::Type face_cells = faceCells(*G);
65 
66  const double *n;
67  const double *K;
68 
69  MAT_SIZE_T nrows, ncols, ldA, incx, incy;
70  double a1, a2;
71 
72  d = dimensions(*G);
73 
74  nrows = ncols = ldA = d;
75  incx = incy = 1 ;
76  a1 = 1.0; a2 = 0.0 ;
77 
78  for (int c =0, i = 0; c < numCells(*G); c++) {
79  K = perm + (c * d * d);
80 
81  typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
82  FaceRow faces = c2f[c];
83 
84  for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
85  f!=end; ++f, ++i)
86  {
87  s = 2.0*(face_cells(*f, 0) == c) - 1.0;
88  n = faceNormal(*G, *f);
89  const double* nn=multiplyFaceNormalWithArea(*G, *f, n);
90  const double* fc = &(faceCentroid(*G, *f)[0]);
91  dgemv_("No Transpose", &nrows, &ncols,
92  &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy);
93  maybeFreeFaceNormal(*G, nn);
94 
95  htrans[i] = denom = 0.0;
96  for (j = 0; j < d; j++) {
97  dist = fc[j] - getCoordinate(cc, j);
98 
99  htrans[i] += s * dist * Kn[j];
100  denom += dist * dist;
101  }
102 
103  assert (denom > 0);
104  htrans[i] /= denom;
105  htrans[i] = std::abs(htrans[i]);
106  }
107  // Move to next cell centroid.
108  cc = increment(cc, 1, d);
109  }
110 }
111 
112 
113 /* ---------------------------------------------------------------------- */
114 template<class Grid>
115 void
116 tpfa_trans_compute(const Grid* G, const double *htrans, double *trans)
117 /* ---------------------------------------------------------------------- */
118 {
119  using namespace Opm::UgGridHelpers;
120 
121  for (int f = 0; f < numFaces(*G); f++) {
122  trans[f] = 0.0;
123  }
124 
125  typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
126 
127  for (int c = 0, i = 0; c < numCells(*G); c++) {
128  typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
129  FaceRow faces = c2f[c];
130 
131  for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
132  f!=end; ++f, ++i)
133  {
134  trans[*f] += 1.0 / htrans[i];
135  }
136  }
137 
138  for (int f = 0; f < numFaces(*G); f++) {
139  trans[f] = 1.0 / trans[f];
140  }
141 }
142 
143 
144 /* ---------------------------------------------------------------------- */
145  template<class Grid>
146 void
148  const double *totmob,
149  const double *htrans,
150  double *trans)
151 /* ---------------------------------------------------------------------- */
152 {
153  using namespace Opm::UgGridHelpers;
154 
155  for (int f = 0; f < numFaces(*G); f++) {
156  trans[f] = 0.0;
157  }
158 
159  typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
160 
161  for (int c = 0, i = 0; c < numCells(*G); c++) {
162  typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
163  FaceRow faces = c2f[c];
164 
165  for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
166  f!=end; ++f, ++i)
167  {
168  trans[*f] += 1.0 / (totmob[c] * htrans[i]);
169  }
170  }
171 
172  for (int f = 0; f < numFaces(*G); f++) {
173  trans[f] = 1.0 / trans[f];
174  }
175 }
[ provides Dune::Grid ]
Definition: CpGrid.hpp:207
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
T area(const Point< T, 2 > *c)
Computes the area of a 2-dimensional triangle.
Definition: Volumes.hpp:118
Holds the implementation of the CpGrid as a pimple.
Definition: CellQuadrature.hpp:29
Maps the grid type to the associated type of the cell to faces mapping.
Definition: GridHelpers.hpp:277
Traits of the cell centroids of a grid.
Definition: GridHelpers.hpp:126
Traits of the face to attached cell mappping of a grid.
Definition: GridHelpers.hpp:334
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition: UnstructuredGrid.h:99
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,...
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.
Definition: trans_tpfa.c:19
void tpfa_trans_compute(struct UnstructuredGrid *G, const double *htrans, double *trans)
Compute two-point transmissibilities from one-sided transmissibilities.
Definition: trans_tpfa.c:74