primitiveMeshTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2012-2015 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();
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();
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(), vector::zero);
258  vectorField sumMagClosed(mesh.nCells(), vector::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();
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  // Normalized 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();
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();
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 normalization 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(symmTensor::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 // ************************************************************************* //
unsigned char direction
Definition: direction.H:43
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)
label nFaces() const
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedTensor skew(const dimensionedTensor &dt)
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
Templated 3D symmetric tensor derived from VectorSpace adding construction from 6 components...
Definition: SymmTensor.H:53
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;.
T & last()
Return the last element of the list.
Definition: UListI.H:131
A bit-packed bool list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label nCells() const
const cellList & cells() const
virtual const faceList & faces() const =0
Return faces.
T & first()
Return the first element of the list.
Definition: UListI.H:117
const Cmpt & xx() const
Definition: SymmTensorI.H:85
label n
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
void setSize(const label)
Reset size of List.
Definition: List.C:318
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
const Cmpt & yy() const
Definition: SymmTensorI.H:103
const Cmpt & zz() const
Definition: SymmTensorI.H:115
static tmp< scalarField > faceOrthogonality(const primitiveMesh &mesh, const vectorField &fAreas, const vectorField &cellCtrs)
Generate non-orthogonality field (internal faces only)
label nInternalFaces() const
pyramid< point, const point &, const face & > pyramidPointFaceRef
Cmpt cmptSum(const VectorSpace< Form, Cmpt, nCmpt > &vs)
Definition: VectorSpaceI.H:414
static tmp< scalarField > faceSkewness(const primitiveMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
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.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
A normal distribution model.
const dimensionedScalar c
Speed of light in a vacuum.
static tmp< scalarField > cellDeterminant(const primitiveMesh &mesh, const Vector< label > &directions, const vectorField &faceAreas, const PackedBoolList &internalOrCoupledFace)
Generate cell determinant field.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
volScalarField scalarField(fieldObject, mesh)
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.
virtual const labelList & faceOwner() const =0
Face face-owner addresing.
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
A class for managing temporary objects.
Definition: PtrList.H:118
static void facePyramidVolume(const primitiveMesh &mesh, const pointField &points, const vectorField &cellCtrs, scalarField &ownPyrVol, scalarField &neiPyrVol)
Generate face pyramid volume fields.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.