primitiveMeshTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "primitiveMeshTools.H"
27 #include "syncTools.H"
28 #include "pyramidPointFaceRef.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
33 (
34  const primitiveMesh& mesh,
35  const pointField& p,
36  const vectorField& fCtrs,
37  const vectorField& fAreas,
38 
39  const label facei,
40  const point& ownCc,
41  const point& neiCc
42 )
43 {
44  vector Cpf = fCtrs[facei] - ownCc;
45  vector d = neiCc - ownCc;
46 
47  // Skewness vector
48  vector sv =
49  Cpf
50  - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
51  vector svHat = sv/(mag(sv) + rootVSmall);
52 
53  // Normalisation distance calculated as the approximate distance
54  // from the face centre to the edge of the face in the direction
55  // of the skewness
56  scalar fd = 0.2*mag(d) + rootVSmall;
57  const face& f = mesh.faces()[facei];
58  forAll(f, pi)
59  {
60  fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
61  }
62 
63  // Normalised skewness
64  return mag(sv)/fd;
65 }
66 
67 
69 (
70  const primitiveMesh& mesh,
71  const pointField& p,
72  const vectorField& fCtrs,
73  const vectorField& fAreas,
74 
75  const label facei,
76  const point& ownCc
77 )
78 {
79  vector Cpf = fCtrs[facei] - ownCc;
80 
81  vector normal = fAreas[facei];
82  normal /= mag(normal) + rootVSmall;
83  vector d = normal*(normal & Cpf);
84 
85 
86  // Skewness vector
87  vector sv =
88  Cpf
89  - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
90  vector svHat = sv/(mag(sv) + rootVSmall);
91 
92  // Normalisation distance calculated as the approximate distance
93  // from the face centre to the edge of the face in the direction
94  // of the skewness
95  scalar fd = 0.4*mag(d) + rootVSmall;
96  const face& f = mesh.faces()[facei];
97  forAll(f, pi)
98  {
99  fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
100  }
101 
102  // Normalised skewness
103  return mag(sv)/fd;
104 }
105 
106 
108 (
109  const point& ownCc,
110  const point& neiCc,
111  const vector& s
112 )
113 {
114  vector d = neiCc - ownCc;
115 
116  return (d & s)/(mag(d)*mag(s) + rootVSmall);
117 }
118 
119 
120 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
121 
123 (
124  const primitiveMesh& mesh,
125  const vectorField& areas,
126  const vectorField& cc
127 )
128 {
129  const labelList& own = mesh.faceOwner();
130  const labelList& nei = mesh.faceNeighbour();
131 
132  tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
133  scalarField& ortho = tortho.ref();
134 
135  // Internal faces
136  forAll(nei, facei)
137  {
138  ortho[facei] = faceOrthogonality
139  (
140  cc[own[facei]],
141  cc[nei[facei]],
142  areas[facei]
143  );
144  }
145 
146  return tortho;
147 }
148 
149 
151 (
152  const primitiveMesh& mesh,
153  const pointField& p,
154  const vectorField& fCtrs,
155  const vectorField& fAreas,
156  const vectorField& cellCtrs
157 )
158 {
159  const labelList& own = mesh.faceOwner();
160  const labelList& nei = mesh.faceNeighbour();
161 
162  tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
163  scalarField& skew = tskew.ref();
164 
165  forAll(nei, facei)
166  {
167  skew[facei] = faceSkewness
168  (
169  mesh,
170  p,
171  fCtrs,
172  fAreas,
173 
174  facei,
175  cellCtrs[own[facei]],
176  cellCtrs[nei[facei]]
177  );
178  }
179 
180 
181  // Boundary faces: consider them to have only skewness error.
182  // (i.e. treat as if mirror cell on other side)
183 
184  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
185  {
186  skew[facei] = boundaryFaceSkewness
187  (
188  mesh,
189  p,
190  fCtrs,
191  fAreas,
192  facei,
193  cellCtrs[own[facei]]
194  );
195  }
196 
197  return tskew;
198 }
199 
200 
202 (
203  const primitiveMesh& mesh,
204  const pointField& points,
205  const vectorField& ctrs,
206 
207  scalarField& ownPyrVol,
208  scalarField& neiPyrVol
209 )
210 {
211  const labelList& own = mesh.faceOwner();
212  const labelList& nei = mesh.faceNeighbour();
213  const faceList& f = mesh.faces();
214 
215  ownPyrVol.setSize(mesh.nFaces());
216  neiPyrVol.setSize(mesh.nInternalFaces());
217 
218  forAll(f, facei)
219  {
220  // Create the owner pyramid
221  ownPyrVol[facei] = -pyramidPointFaceRef
222  (
223  f[facei],
224  ctrs[own[facei]]
225  ).mag(points);
226 
227  if (mesh.isInternalFace(facei))
228  {
229  // Create the neighbour pyramid - it will have positive volume
230  neiPyrVol[facei] = pyramidPointFaceRef
231  (
232  f[facei],
233  ctrs[nei[facei]]
234  ).mag(points);
235  }
236  }
237 }
238 
239 
241 (
242  const primitiveMesh& mesh,
243  const Vector<label>& meshD,
244  const vectorField& areas,
245  const scalarField& vols,
246 
247  scalarField& openness,
248  scalarField& aratio
249 )
250 {
251  const labelList& own = mesh.faceOwner();
252  const labelList& nei = mesh.faceNeighbour();
253 
254  // Loop through cell faces and sum up the face area vectors for each cell.
255  // This should be zero in all vector components
256 
257  vectorField sumClosed(mesh.nCells(), Zero);
258  vectorField sumMagClosed(mesh.nCells(), Zero);
259 
260  forAll(own, facei)
261  {
262  // Add to owner
263  sumClosed[own[facei]] += areas[facei];
264  sumMagClosed[own[facei]] += cmptMag(areas[facei]);
265  }
266 
267  forAll(nei, facei)
268  {
269  // Subtract from neighbour
270  sumClosed[nei[facei]] -= areas[facei];
271  sumMagClosed[nei[facei]] += cmptMag(areas[facei]);
272  }
273 
274 
275  label nDims = 0;
276  for (direction dir = 0; dir < vector::nComponents; dir++)
277  {
278  if (meshD[dir] == 1)
279  {
280  nDims++;
281  }
282  }
283 
284 
285  // Check the sums
286  openness.setSize(mesh.nCells());
287  aratio.setSize(mesh.nCells());
288 
289  forAll(sumClosed, celli)
290  {
291  scalar maxOpenness = 0;
292 
293  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
294  {
295  maxOpenness = max
296  (
297  maxOpenness,
298  mag(sumClosed[celli][cmpt])
299  /(sumMagClosed[celli][cmpt] + rootVSmall)
300  );
301  }
302  openness[celli] = maxOpenness;
303 
304  // Calculate the aspect ration as the maximum of Cartesian component
305  // aspect ratio to the total area hydraulic area aspect ratio
306  scalar minCmpt = vGreat;
307  scalar maxCmpt = -vGreat;
308  for (direction dir = 0; dir < vector::nComponents; dir++)
309  {
310  if (meshD[dir] == 1)
311  {
312  minCmpt = min(minCmpt, sumMagClosed[celli][dir]);
313  maxCmpt = max(maxCmpt, sumMagClosed[celli][dir]);
314  }
315  }
316 
317  scalar aspectRatio = maxCmpt/(minCmpt + rootVSmall);
318  if (nDims == 3)
319  {
320  scalar v = max(rootVSmall, vols[celli]);
321 
322  aspectRatio = max
323  (
324  aspectRatio,
325  1.0/6.0*cmptSum(sumMagClosed[celli])/pow(v, 2.0/3.0)
326  );
327  }
328 
329  aratio[celli] = aspectRatio;
330  }
331 }
332 
333 
335 (
336  const scalar maxSin,
337  const primitiveMesh& mesh,
338  const pointField& p,
339  const vectorField& faceAreas
340 )
341 {
342  const faceList& fcs = mesh.faces();
343 
344  vectorField faceNormals(faceAreas);
345  faceNormals /= mag(faceNormals) + rootVSmall;
346 
347  tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
348  scalarField& faceAngles = tfaceAngles.ref();
349 
350 
351  forAll(fcs, facei)
352  {
353  const face& f = fcs[facei];
354 
355  // Get edge from f[0] to f[size-1];
356  vector ePrev(p[f.first()] - p[f.last()]);
357  scalar magEPrev = mag(ePrev);
358  ePrev /= magEPrev + rootVSmall;
359 
360  scalar maxEdgeSin = 0.0;
361 
362  forAll(f, fp0)
363  {
364  // Get vertex after fp
365  label fp1 = f.fcIndex(fp0);
366 
367  // Normalised vector between two consecutive points
368  vector e10(p[f[fp1]] - p[f[fp0]]);
369  scalar magE10 = mag(e10);
370  e10 /= magE10 + rootVSmall;
371 
372  if (magEPrev > small && magE10 > small)
373  {
374  vector edgeNormal = ePrev ^ e10;
375  scalar magEdgeNormal = mag(edgeNormal);
376 
377  if (magEdgeNormal < maxSin)
378  {
379  // Edges (almost) aligned -> face is ok.
380  }
381  else
382  {
383  // Check normal
384  edgeNormal /= magEdgeNormal;
385 
386  if ((edgeNormal & faceNormals[facei]) < small)
387  {
388  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
389  }
390  }
391  }
392 
393  ePrev = e10;
394  magEPrev = magE10;
395  }
396 
397  faceAngles[facei] = maxEdgeSin;
398  }
399 
400  return tfaceAngles;
401 }
402 
403 
405 (
406  const primitiveMesh& mesh,
407  const pointField& p,
408  const vectorField& fCtrs,
409  const vectorField& faceAreas
410 )
411 {
412  const faceList& fcs = mesh.faces();
413 
414  // Areas are calculated as the sum of areas. (see
415  // primitiveMeshFaceCentresAndAreas.C)
416  scalarField magAreas(mag(faceAreas));
417 
418  tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
419  scalarField& faceFlatness = tfaceFlatness.ref();
420 
421 
422  forAll(fcs, facei)
423  {
424  const face& f = fcs[facei];
425 
426  if (f.size() > 3 && magAreas[facei] > rootVSmall)
427  {
428  const point& fc = fCtrs[facei];
429 
430  // Calculate the sum of magnitude of areas and compare to magnitude
431  // of sum of areas.
432 
433  scalar sumA = 0.0;
434 
435  forAll(f, fp)
436  {
437  const point& thisPoint = p[f[fp]];
438  const point& nextPoint = p[f.nextLabel(fp)];
439 
440  // Triangle around fc.
441  vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
442  sumA += mag(n);
443  }
444 
445  faceFlatness[facei] = magAreas[facei]/(sumA + rootVSmall);
446  }
447  }
448 
449  return tfaceFlatness;
450 }
451 
452 
454 (
455  const primitiveMesh& mesh,
456  const Vector<label>& meshD,
457  const vectorField& faceAreas,
458  const PackedBoolList& internalOrCoupledFace
459 )
460 {
461  // Determine number of dimensions and (for 2D) missing dimension
462  label nDims = 0;
463  label twoD = -1;
464  for (direction dir = 0; dir < vector::nComponents; dir++)
465  {
466  if (meshD[dir] == 1)
467  {
468  nDims++;
469  }
470  else
471  {
472  twoD = dir;
473  }
474  }
475 
476  tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
477  scalarField& cellDeterminant = tcellDeterminant.ref();
478 
479  const cellList& c = mesh.cells();
480 
481  if (nDims == 1)
482  {
483  cellDeterminant = 1.0;
484  }
485  else
486  {
487  forAll(c, celli)
488  {
489  const labelList& curFaces = c[celli];
490 
491  // Calculate local normalisation factor
492  scalar avgArea = 0;
493 
494  label nInternalFaces = 0;
495 
496  forAll(curFaces, i)
497  {
498  if (internalOrCoupledFace[curFaces[i]])
499  {
500  avgArea += mag(faceAreas[curFaces[i]]);
501 
502  nInternalFaces++;
503  }
504  }
505 
506  if (nInternalFaces == 0)
507  {
508  cellDeterminant[celli] = 0;
509  }
510  else
511  {
512  avgArea /= nInternalFaces;
513 
514  symmTensor areaTensor(Zero);
515 
516  forAll(curFaces, i)
517  {
518  if (internalOrCoupledFace[curFaces[i]])
519  {
520  areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
521  }
522  }
523 
524  if (nDims == 2)
525  {
526  // Add the missing eigenvector (such that it does not
527  // affect the determinant)
528  if (twoD == 0)
529  {
530  areaTensor.xx() = 1;
531  }
532  else if (twoD == 1)
533  {
534  areaTensor.yy() = 1;
535  }
536  else
537  {
538  areaTensor.zz() = 1;
539  }
540  }
541 
542  cellDeterminant[celli] = mag(det(areaTensor));
543  }
544  }
545  }
546 
547  return tcellDeterminant;
548 }
549 
550 
551 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
label nInternalFaces() const
dimensionedTensor skew(const dimensionedTensor &dt)
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
label nFaces() const
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
const Cmpt & xx() const
Definition: SymmTensorI.H:87
dimensionedSymmTensor sqr(const dimensionedVector &dv)
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label nCells() const
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:500
const cellList & cells() const
dimensionedScalar det(const dimensionedSphericalTensor &dt)
static tmp< scalarField > faceConcavity(const scalar maxSin, const primitiveMesh &mesh, const pointField &p, const vectorField &faceAreas)
Generate face concavity field. Returns per face the (sin of the)
T & first()
Return the first element of the list.
Definition: UListI.H:114
const Cmpt & yy() const
Definition: SymmTensorI.H:105
const dimensionedScalar c
Speed of light in a vacuum.
static tmp< scalarField > faceSkewness(const primitiveMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
static scalar boundaryFaceSkewness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &fAreas, const label facei, const point &ownCc)
Skewness of single boundary face.
const Cmpt & zz() const
Definition: SymmTensorI.H:117
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:275
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const zero Zero
Definition: zero.H:97
pyramid< point, const point &, const face & > pyramidPointFaceRef
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
volScalarField scalarField(fieldObject, mesh)
static void facePyramidVolume(const primitiveMesh &mesh, const pointField &points, const vectorField &cellCtrs, scalarField &ownPyrVol, scalarField &neiPyrVol)
Generate face pyramid volume fields.
labelList f(nPoints)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void setSize(const label)
Reset size of List.
Definition: List.C:281
A bit-packed bool list.
virtual const faceList & faces() const =0
Return faces.
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
static tmp< scalarField > faceOrthogonality(const primitiveMesh &mesh, const vectorField &fAreas, const vectorField &cellCtrs)
Generate non-orthogonality field (internal faces only)
A class for managing temporary objects.
Definition: PtrList.H:53
static void cellClosedness(const primitiveMesh &mesh, const Vector< label > &meshD, const vectorField &areas, const scalarField &vols, scalarField &openness, scalarField &aratio)
Generate cell openness and cell ascpect ratio field.
T & last()
Return the last element of the list.
Definition: UListI.H:128
static tmp< scalarField > faceFlatness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &faceAreas)
Generate face flatness field. Compares the individual triangles&#39;.
static tmp< scalarField > cellDeterminant(const primitiveMesh &mesh, const Vector< label > &directions, const vectorField &faceAreas, const PackedBoolList &internalOrCoupledFace)
Generate cell determinant field.