polyMeshTetDecomposition.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) 2011-2018 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 
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 // Note: the use of this tolerance is ad-hoc, there may be extreme
31 // cases where the resulting tetrahedra still have particle tracking
32 // problems, or tets with lower quality may track OK.
33 const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = sqr(small);
34 
35 
36 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
37 
39 (
40  const polyMesh& mesh,
41  label fI,
42  const point& nCc,
43  scalar tol,
44  bool report
45 )
46 {
47  const faceList& pFaces = mesh.faces();
48  const pointField& pPts = mesh.points();
49  const vectorField& pC = mesh.cellCentres();
50  const labelList& pOwner = mesh.faceOwner();
51 
52  const face& f = pFaces[fI];
53 
54  label oCI = pOwner[fI];
55 
56  const point& oCc = pC[oCI];
57 
58  List<scalar> tetQualities(2, 0.0);
59 
60  forAll(f, faceBasePtI)
61  {
62  scalar thisBaseMinTetQuality = vGreat;
63 
64  const point& tetBasePt = pPts[f[faceBasePtI]];
65 
66  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
67  {
68  label facePtI = (tetPtI + faceBasePtI) % f.size();
69  label otherFacePtI = f.fcIndex(facePtI);
70 
71  {
72  // owner cell tet
73  label ptAI = f[facePtI];
74  label ptBI = f[otherFacePtI];
75 
76  const point& pA = pPts[ptAI];
77  const point& pB = pPts[ptBI];
78 
79  tetPointRef tet(oCc, tetBasePt, pA, pB);
80 
81  tetQualities[0] = tet.quality();
82  }
83 
84  {
85  // neighbour cell tet
86  label ptAI = f[otherFacePtI];
87  label ptBI = f[facePtI];
88 
89  const point& pA = pPts[ptAI];
90  const point& pB = pPts[ptBI];
91 
92  tetPointRef tet(nCc, tetBasePt, pA, pB);
93 
94  tetQualities[1] = tet.quality();
95  }
96 
97  if (min(tetQualities) < thisBaseMinTetQuality)
98  {
99  thisBaseMinTetQuality = min(tetQualities);
100  }
101  }
102 
103  if (thisBaseMinTetQuality > tol)
104  {
105  return faceBasePtI;
106  }
107  }
108 
109  // If a base point hasn't triggered a return by now, then there is
110  // none that can produce a good decomposition
111  return -1;
112 }
113 
114 
116 (
117  const polyMesh& mesh,
118  label fI,
119  scalar tol,
120  bool report
121 )
122 {
123  return findSharedBasePoint
124  (
125  mesh,
126  fI,
127  mesh.cellCentres()[mesh.faceNeighbour()[fI]],
128  tol,
129  report
130  );
131 }
132 
133 
135 (
136  const polyMesh& mesh,
137  label fI,
138  scalar tol,
139  bool report
140 )
141 {
142  const faceList& pFaces = mesh.faces();
143  const pointField& pPts = mesh.points();
144  const vectorField& pC = mesh.cellCentres();
145  const labelList& pOwner = mesh.faceOwner();
146 
147  const face& f = pFaces[fI];
148 
149  label cI = pOwner[fI];
150 
151  bool own = (pOwner[fI] == cI);
152 
153  const point& cC = pC[cI];
154 
155  forAll(f, faceBasePtI)
156  {
157  scalar thisBaseMinTetQuality = vGreat;
158 
159  const point& tetBasePt = pPts[f[faceBasePtI]];
160 
161  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
162  {
163  label facePtI = (tetPtI + faceBasePtI) % f.size();
164  label otherFacePtI = f.fcIndex(facePtI);
165 
166  label ptAI = -1;
167  label ptBI = -1;
168 
169  if (own)
170  {
171  ptAI = f[facePtI];
172  ptBI = f[otherFacePtI];
173  }
174  else
175  {
176  ptAI = f[otherFacePtI];
177  ptBI = f[facePtI];
178  }
179 
180  const point& pA = pPts[ptAI];
181  const point& pB = pPts[ptBI];
182 
183  tetPointRef tet(cC, tetBasePt, pA, pB);
184 
185  scalar tetQuality = tet.quality();
186 
187  if (tetQuality < thisBaseMinTetQuality)
188  {
189  thisBaseMinTetQuality = tetQuality;
190  }
191  }
192 
193  if (thisBaseMinTetQuality > tol)
194  {
195  return faceBasePtI;
196  }
197  }
198 
199  // If a base point hasn't triggered a return by now, then there is
200  // none that can produce a good decomposition
201  return -1;
202 }
203 
204 
206 (
207  const polyMesh& mesh,
208  scalar tol,
209  bool report
210 )
211 {
212  const labelList& pOwner = mesh.faceOwner();
213  const vectorField& pC = mesh.cellCentres();
214 
215  // Find a suitable base point for each face, considering both
216  // cells for interface faces or those on coupled patches
217 
218  labelList tetBasePtIs(mesh.nFaces(), -1);
219 
220  label nInternalFaces = mesh.nInternalFaces();
221 
222  for (label fI = 0; fI < nInternalFaces; fI++)
223  {
224  tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
225  }
226 
227  pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
228 
229  for(label facei = nInternalFaces; facei < mesh.nFaces(); facei++)
230  {
231  neighbourCellCentres[facei - nInternalFaces] =
232  pC[pOwner[facei]];
233  }
234 
235  syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
236 
237  const polyBoundaryMesh& patches = mesh.boundaryMesh();
238 
239  SubList<label> boundaryFaceTetBasePtIs
240  (
241  tetBasePtIs,
242  mesh.nFaces() - nInternalFaces,
243  nInternalFaces
244  );
245 
246  for
247  (
248  label fI = nInternalFaces, bFI = 0;
249  fI < mesh.nFaces();
250  fI++, bFI++
251  )
252  {
253  label patchi =
254  mesh.boundaryMesh().patchID()[bFI];
255 
256  if (patches[patchi].coupled())
257  {
258  const coupledPolyPatch& pp =
259  refCast<const coupledPolyPatch>(patches[patchi]);
260 
261  if (pp.owner())
262  {
263  boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
264  (
265  mesh,
266  fI,
267  neighbourCellCentres[bFI],
268  tol,
269  report
270  );
271  }
272  else
273  {
274  // Assign -2, to distinguish from a failed base point
275  // find, which returns -1.
276  boundaryFaceTetBasePtIs[bFI] = -2;
277  }
278  }
279  else
280  {
281  boundaryFaceTetBasePtIs[bFI] = findBasePoint
282  (
283  mesh,
284  fI,
285  tol,
286  report
287  );
288  }
289  }
290 
291  // maxEqOp will replace the -2 values on the neighbour patches
292  // with the result from the owner base point find.
293 
294  syncTools::syncBoundaryFaceList
295  (
296  mesh,
297  boundaryFaceTetBasePtIs,
299  );
300 
301  for
302  (
303  label fI = nInternalFaces, bFI = 0;
304  fI < mesh.nFaces();
305  fI++, bFI++
306  )
307  {
308  label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
309 
310  if (bFTetBasePtI == -2)
311  {
313  << "Coupled face base point exchange failure for face "
314  << fI
315  << abort(FatalError);
316  }
317 
318  if (bFTetBasePtI < 1)
319  {
320  // If the base point is -1, it should be left as such to
321  // indicate a problem, if it is 0, then no action is required.
322 
323  continue;
324  }
325 
326  label patchi = mesh.boundaryMesh().patchID()[bFI];
327 
328  if (patches[patchi].coupled())
329  {
330  const coupledPolyPatch& pp =
331  refCast<const coupledPolyPatch>(patches[patchi]);
332 
333  // Calculated base points on coupled faces are those of
334  // the owner patch face. They need to be reindexed to for
335  // the non-owner face, which has the opposite order.
336 
337  // So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
338 
339  // i.e.:
340 
341  // owner coupledPolyPatch face
342  // face (a b c d e f)
343  // fPtI 0 1 2 3 4 5
344  // +
345  // #
346 
347  // neighbour coupledPolyPatch face
348  // face (a f e d c b)
349  // fPtI 0 1 2 3 4 5
350  // +
351  // #
352  // +: 6 - 1 = 5
353  // #: 6 - 2 = 4
354 
355  if (!pp.owner())
356  {
357  bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
358  }
359  }
360  }
361 
362  return tetBasePtIs;
363 }
364 
365 
367 (
368  const polyMesh& mesh,
369  scalar tol,
370  const bool report,
371  labelHashSet* setPtr
372 )
373 {
374  const labelList& own = mesh.faceOwner();
375  const labelList& nei = mesh.faceNeighbour();
376  const polyBoundaryMesh& patches = mesh.boundaryMesh();
377 
378  const vectorField& cc = mesh.cellCentres();
379  const vectorField& fc = mesh.faceCentres();
380 
381  // Calculate coupled cell centre
382  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
383 
384  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
385  {
386  neiCc[facei - mesh.nInternalFaces()] = cc[own[facei]];
387  }
388 
389  syncTools::swapBoundaryFacePositions(mesh, neiCc);
390 
391  const faceList& fcs = mesh.faces();
392 
393  const pointField& p = mesh.points();
394 
395  label nErrorTets = 0;
396 
397  forAll(fcs, facei)
398  {
399  const face& f = fcs[facei];
400 
401  forAll(f, fPtI)
402  {
403  scalar tetQual = tetPointRef
404  (
405  p[f[fPtI]],
406  p[f.nextLabel(fPtI)],
407  fc[facei],
408  cc[own[facei]]
409  ).quality();
410 
411  if (tetQual > -tol)
412  {
413  if (setPtr)
414  {
415  setPtr->insert(facei);
416  }
417 
418  nErrorTets++;
419  break; // no need to check other tets
420  }
421  }
422 
423  if (mesh.isInternalFace(facei))
424  {
425  // Create the neighbour tet - it will have positive volume
426  const face& f = fcs[facei];
427 
428  forAll(f, fPtI)
429  {
430  scalar tetQual = tetPointRef
431  (
432  p[f[fPtI]],
433  p[f.nextLabel(fPtI)],
434  fc[facei],
435  cc[nei[facei]]
436  ).quality();
437 
438  if (tetQual < tol)
439  {
440  if (setPtr)
441  {
442  setPtr->insert(facei);
443  }
444 
445  nErrorTets++;
446  break;
447  }
448  }
449 
450  if (findSharedBasePoint(mesh, facei, tol, report) == -1)
451  {
452  if (setPtr)
453  {
454  setPtr->insert(facei);
455  }
456 
457  nErrorTets++;
458  }
459  }
460  else
461  {
462  label patchi = patches.patchID()[facei - mesh.nInternalFaces()];
463 
464  if (patches[patchi].coupled())
465  {
466  if
467  (
468  findSharedBasePoint
469  (
470  mesh,
471  facei,
472  neiCc[facei - mesh.nInternalFaces()],
473  tol,
474  report
475  ) == -1
476  )
477  {
478  if (setPtr)
479  {
480  setPtr->insert(facei);
481  }
482 
483  nErrorTets++;
484  }
485  }
486  else
487  {
488  if (findBasePoint(mesh, facei, tol, report) == -1)
489  {
490  if (setPtr)
491  {
492  setPtr->insert(facei);
493  }
494 
495  nErrorTets++;
496  }
497  }
498  }
499  }
500 
501  reduce(nErrorTets, sumOp<label>());
502 
503  if (nErrorTets > 0)
504  {
505  if (report)
506  {
507  Info<< " ***Error in face tets: "
508  << nErrorTets << " faces with low quality or negative volume "
509  << "decomposition tets." << endl;
510  }
511 
512  return true;
513  }
514  else
515  {
516  if (report)
517  {
518  Info<< " Face tets OK." << endl;
519  }
520 
521  return false;
522  }
523 }
524 
525 
527 (
528  const polyMesh& mesh,
529  label fI,
530  label cI
531 )
532 {
533  const faceList& pFaces = mesh.faces();
534 
535  const face& f = pFaces[fI];
536 
537  label nTets = f.size() - 2;
538 
539  List<tetIndices> faceTets(nTets);
540 
541  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++)
542  {
543  faceTets[tetPtI - 1] = tetIndices(cI, fI, tetPtI);
544  }
545 
546  return faceTets;
547 }
548 
549 
551 (
552  const polyMesh& mesh,
553  label cI
554 )
555 {
556  const faceList& pFaces = mesh.faces();
557  const cellList& pCells = mesh.cells();
558 
559  const cell& thisCell = pCells[cI];
560 
561  label nTets = 0;
562 
563  forAll(thisCell, cFI)
564  {
565  nTets += pFaces[thisCell[cFI]].size() - 2;
566  }
567 
568  DynamicList<tetIndices> cellTets(nTets);
569 
570  forAll(thisCell, cFI)
571  {
572  label fI = thisCell[cFI];
573 
574  cellTets.append(faceTetIndices(mesh, fI, cI));
575  }
576 
577  return cellTets;
578 }
579 
580 
582 (
583  const polyMesh& mesh,
584  label cI,
585  const point& pt
586 )
587 {
588  const faceList& pFaces = mesh.faces();
589  const cellList& pCells = mesh.cells();
590 
591  const cell& thisCell = pCells[cI];
592 
593  tetIndices tetContainingPt;
594 
595 
596  forAll(thisCell, cFI)
597  {
598  label fI = thisCell[cFI];
599  const face& f = pFaces[fI];
600 
601  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
602  {
603  // Get tetIndices of face triangle
604  tetIndices faceTetIs(cI, fI, tetPtI);
605 
606  // Check if inside
607  if (faceTetIs.tet(mesh).inside(pt))
608  {
609  tetContainingPt = faceTetIs;
610  break;
611  }
612  }
613 
614  if (tetContainingPt.cell() != -1)
615  {
616  break;
617  }
618  }
619 
620  return tetContainingPt;
621 }
622 
623 
624 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
A tetrahedron primitive.
Definition: tetrahedron.H:61
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
Find a suitable base point for each face for decomposition.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1047
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const labelList & patchID() const
Per boundary face label the patch index.
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
const cellList & cells() const
patches[0]
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
Definition: tetIndicesI.H:100
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1003
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
static bool checkFaceTets(const polyMesh &mesh, scalar tol=minTetQuality, const bool report=false, labelHashSet *setPtr=nullptr)
Check face-decomposition tet volume.
A List obtained as a section of another List.
Definition: SubList.H:53
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
label cell() const
Return the cell.
Definition: tetIndicesI.H:28
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:230
static label findSharedBasePoint(const polyMesh &mesh, label fI, const point &nCc, scalar tol, bool report=false)
Find the first suitable base point to use for a minimum.
virtual bool owner() const =0
Does this side own the patch ?
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1041
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1028
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
Definition: tetrahedronI.H:413
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const vectorField & faceCentres() const
label patchi
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
messageStream Info
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
volScalarField & p
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
static const scalar minTetQuality
Minimum tetrahedron quality.
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.