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-2019 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  const point& cC,
42  const label fI,
43  const bool isOwner,
44  const label faceBasePtI
45 )
46 {
47  // Does fan decomposition of face (starting at faceBasePti) and determines
48  // min quality over all resulting tets.
49 
50  const pointField& pPts = mesh.points();
51  const face& f = mesh.faces()[fI];
52  const point& tetBasePt = pPts[f[faceBasePtI]];
53 
54  scalar thisBaseMinTetQuality = vGreat;
55 
56  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
57  {
58  label facePtI = (tetPtI + faceBasePtI) % f.size();
59  label otherFacePtI = f.fcIndex(facePtI);
60 
61  label ptAI = -1;
62  label ptBI = -1;
63 
64  if (isOwner)
65  {
66  ptAI = f[facePtI];
67  ptBI = f[otherFacePtI];
68  }
69  else
70  {
71  ptAI = f[otherFacePtI];
72  ptBI = f[facePtI];
73  }
74 
75  const point& pA = pPts[ptAI];
76  const point& pB = pPts[ptBI];
77 
78  tetPointRef tet(cC, tetBasePt, pA, pB);
79 
80  scalar tetQuality = tet.quality();
81 
82  if (tetQuality < thisBaseMinTetQuality)
83  {
84  thisBaseMinTetQuality = tetQuality;
85  }
86  }
87  return thisBaseMinTetQuality;
88 }
89 
90 
92 (
93  const polyMesh& mesh,
94  label fI,
95  const point& nCc,
96  scalar tol,
97  bool report
98 )
99 {
100  const faceList& pFaces = mesh.faces();
101  const vectorField& pC = mesh.cellCentres();
102  const labelList& pOwner = mesh.faceOwner();
103 
104  const face& f = pFaces[fI];
105 
106  label oCI = pOwner[fI];
107 
108  const point& oCc = pC[oCI];
109 
110  forAll(f, faceBasePtI)
111  {
112  scalar minQ = minQuality(mesh, oCc, fI, true, faceBasePtI);
113  minQ = min(minQ, minQuality(mesh, nCc, fI, false, faceBasePtI));
114 
115  if (minQ > tol)
116  {
117  return faceBasePtI;
118  }
119  }
120 
121  // If a base point hasn't triggered a return by now, then there is
122  // none that can produce a good decomposition
123  return -1;
124 }
125 
126 
128 (
129  const polyMesh& mesh,
130  label fI,
131  scalar tol,
132  bool report
133 )
134 {
135  return findSharedBasePoint
136  (
137  mesh,
138  fI,
139  mesh.cellCentres()[mesh.faceNeighbour()[fI]],
140  tol,
141  report
142  );
143 }
144 
145 
147 (
148  const polyMesh& mesh,
149  label fI,
150  scalar tol,
151  bool report
152 )
153 {
154  const faceList& pFaces = mesh.faces();
155  const vectorField& pC = mesh.cellCentres();
156  const labelList& pOwner = mesh.faceOwner();
157 
158  const face& f = pFaces[fI];
159 
160  label cI = pOwner[fI];
161 
162  bool own = (pOwner[fI] == cI);
163 
164  const point& cC = pC[cI];
165 
166  forAll(f, faceBasePtI)
167  {
168  scalar quality = minQuality(mesh, cC, fI, own, faceBasePtI);
169 
170  if (quality > tol)
171  {
172  return faceBasePtI;
173  }
174  }
175 
176  // If a base point hasn't triggered a return by now, then there is
177  // none that can produce a good decomposition
178  return -1;
179 }
180 
181 
183 (
184  const polyMesh& mesh,
185  scalar tol,
186  bool report
187 )
188 {
189  const labelList& pOwner = mesh.faceOwner();
190  const vectorField& pC = mesh.cellCentres();
191 
192  // Find a suitable base point for each face, considering both
193  // cells for interface faces or those on coupled patches
194 
195  labelList tetBasePtIs(mesh.nFaces(), -1);
196 
197  label nInternalFaces = mesh.nInternalFaces();
198 
199  for (label fI = 0; fI < nInternalFaces; fI++)
200  {
201  tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
202  }
203 
204  pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
205 
206  for(label facei = nInternalFaces; facei < mesh.nFaces(); facei++)
207  {
208  neighbourCellCentres[facei - nInternalFaces] =
209  pC[pOwner[facei]];
210  }
211 
212  syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
213 
214  const polyBoundaryMesh& patches = mesh.boundaryMesh();
215 
216  SubList<label> boundaryFaceTetBasePtIs
217  (
218  tetBasePtIs,
219  mesh.nFaces() - nInternalFaces,
220  nInternalFaces
221  );
222 
223  for
224  (
225  label fI = nInternalFaces, bFI = 0;
226  fI < mesh.nFaces();
227  fI++, bFI++
228  )
229  {
230  label patchi =
231  mesh.boundaryMesh().patchID()[bFI];
232 
233  if (patches[patchi].coupled())
234  {
235  const coupledPolyPatch& pp =
236  refCast<const coupledPolyPatch>(patches[patchi]);
237 
238  if (pp.owner())
239  {
240  boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
241  (
242  mesh,
243  fI,
244  neighbourCellCentres[bFI],
245  tol,
246  report
247  );
248  }
249  else
250  {
251  // Assign -2, to distinguish from a failed base point
252  // find, which returns -1.
253  boundaryFaceTetBasePtIs[bFI] = -2;
254  }
255  }
256  else
257  {
258  boundaryFaceTetBasePtIs[bFI] = findBasePoint
259  (
260  mesh,
261  fI,
262  tol,
263  report
264  );
265  }
266  }
267 
268  // maxEqOp will replace the -2 values on the neighbour patches
269  // with the result from the owner base point find.
270 
271  syncTools::syncBoundaryFaceList
272  (
273  mesh,
274  boundaryFaceTetBasePtIs,
276  );
277 
278  for
279  (
280  label fI = nInternalFaces, bFI = 0;
281  fI < mesh.nFaces();
282  fI++, bFI++
283  )
284  {
285  label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
286 
287  if (bFTetBasePtI == -2)
288  {
290  << "Coupled face base point exchange failure for face "
291  << fI
292  << abort(FatalError);
293  }
294 
295  if (bFTetBasePtI < 1)
296  {
297  // If the base point is -1, it should be left as such to
298  // indicate a problem, if it is 0, then no action is required.
299 
300  continue;
301  }
302 
303  label patchi = mesh.boundaryMesh().patchID()[bFI];
304 
305  if (patches[patchi].coupled())
306  {
307  const coupledPolyPatch& pp =
308  refCast<const coupledPolyPatch>(patches[patchi]);
309 
310  // Calculated base points on coupled faces are those of
311  // the owner patch face. They need to be reindexed to for
312  // the non-owner face, which has the opposite order.
313 
314  // So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
315 
316  // i.e.:
317 
318  // owner coupledPolyPatch face
319  // face (a b c d e f)
320  // fPtI 0 1 2 3 4 5
321  // +
322  // #
323 
324  // neighbour coupledPolyPatch face
325  // face (a f e d c b)
326  // fPtI 0 1 2 3 4 5
327  // +
328  // #
329  // +: 6 - 1 = 5
330  // #: 6 - 2 = 4
331 
332  if (!pp.owner())
333  {
334  bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
335  }
336  }
337  }
338 
339  return tetBasePtIs;
340 }
341 
342 
344 (
345  const polyMesh& mesh,
346  scalar tol,
347  const bool report,
348  labelHashSet* setPtr
349 )
350 {
351  const labelList& own = mesh.faceOwner();
352  const labelList& nei = mesh.faceNeighbour();
353  const polyBoundaryMesh& patches = mesh.boundaryMesh();
354 
355  const vectorField& cc = mesh.cellCentres();
356  const vectorField& fc = mesh.faceCentres();
357 
358  // Calculate coupled cell centre
359  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
360 
361  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
362  {
363  neiCc[facei - mesh.nInternalFaces()] = cc[own[facei]];
364  }
365 
366  syncTools::swapBoundaryFacePositions(mesh, neiCc);
367 
368  const faceList& fcs = mesh.faces();
369 
370  const pointField& p = mesh.points();
371 
372  label nErrorTets = 0;
373 
374  forAll(fcs, facei)
375  {
376  const face& f = fcs[facei];
377 
378  forAll(f, fPtI)
379  {
380  scalar tetQual = tetPointRef
381  (
382  p[f[fPtI]],
383  p[f.nextLabel(fPtI)],
384  fc[facei],
385  cc[own[facei]]
386  ).quality();
387 
388  if (tetQual > -tol)
389  {
390  if (setPtr)
391  {
392  setPtr->insert(facei);
393  }
394 
395  nErrorTets++;
396  break; // no need to check other tets
397  }
398  }
399 
400  if (mesh.isInternalFace(facei))
401  {
402  // Create the neighbour tet - it will have positive volume
403  const face& f = fcs[facei];
404 
405  forAll(f, fPtI)
406  {
407  scalar tetQual = tetPointRef
408  (
409  p[f[fPtI]],
410  p[f.nextLabel(fPtI)],
411  fc[facei],
412  cc[nei[facei]]
413  ).quality();
414 
415  if (tetQual < tol)
416  {
417  if (setPtr)
418  {
419  setPtr->insert(facei);
420  }
421 
422  nErrorTets++;
423  break;
424  }
425  }
426 
427  if (findSharedBasePoint(mesh, facei, tol, report) == -1)
428  {
429  if (setPtr)
430  {
431  setPtr->insert(facei);
432  }
433 
434  nErrorTets++;
435  }
436  }
437  else
438  {
439  label patchi = patches.patchID()[facei - mesh.nInternalFaces()];
440 
441  if (patches[patchi].coupled())
442  {
443  if
444  (
445  findSharedBasePoint
446  (
447  mesh,
448  facei,
449  neiCc[facei - mesh.nInternalFaces()],
450  tol,
451  report
452  ) == -1
453  )
454  {
455  if (setPtr)
456  {
457  setPtr->insert(facei);
458  }
459 
460  nErrorTets++;
461  }
462  }
463  else
464  {
465  if (findBasePoint(mesh, facei, tol, report) == -1)
466  {
467  if (setPtr)
468  {
469  setPtr->insert(facei);
470  }
471 
472  nErrorTets++;
473  }
474  }
475  }
476  }
477 
478  reduce(nErrorTets, sumOp<label>());
479 
480  if (nErrorTets > 0)
481  {
482  if (report)
483  {
484  Info<< " ***Error in face tets: "
485  << nErrorTets << " faces with low quality or negative volume "
486  << "decomposition tets." << endl;
487  }
488 
489  return true;
490  }
491  else
492  {
493  if (report)
494  {
495  Info<< " Face tets OK." << endl;
496  }
497 
498  return false;
499  }
500 }
501 
502 
504 (
505  const polyMesh& mesh,
506  label fI,
507  label cI
508 )
509 {
510  const faceList& pFaces = mesh.faces();
511 
512  const face& f = pFaces[fI];
513 
514  label nTets = f.size() - 2;
515 
516  List<tetIndices> faceTets(nTets);
517 
518  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++)
519  {
520  faceTets[tetPtI - 1] = tetIndices(cI, fI, tetPtI);
521  }
522 
523  return faceTets;
524 }
525 
526 
528 (
529  const polyMesh& mesh,
530  label cI
531 )
532 {
533  const faceList& pFaces = mesh.faces();
534  const cellList& pCells = mesh.cells();
535 
536  const cell& thisCell = pCells[cI];
537 
538  label nTets = 0;
539 
540  forAll(thisCell, cFI)
541  {
542  nTets += pFaces[thisCell[cFI]].size() - 2;
543  }
544 
545  DynamicList<tetIndices> cellTets(nTets);
546 
547  forAll(thisCell, cFI)
548  {
549  label fI = thisCell[cFI];
550 
551  cellTets.append(faceTetIndices(mesh, fI, cI));
552  }
553 
554  return move(cellTets);
555 }
556 
557 
559 (
560  const polyMesh& mesh,
561  label cI,
562  const point& pt
563 )
564 {
565  const faceList& pFaces = mesh.faces();
566  const cellList& pCells = mesh.cells();
567 
568  const cell& thisCell = pCells[cI];
569 
570  tetIndices tetContainingPt;
571 
572 
573  forAll(thisCell, cFI)
574  {
575  label fI = thisCell[cFI];
576  const face& f = pFaces[fI];
577 
578  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
579  {
580  // Get tetIndices of face triangle
581  tetIndices faceTetIs(cI, fI, tetPtI);
582 
583  // Check if inside
584  if (faceTetIs.tet(mesh).inside(pt))
585  {
586  tetContainingPt = faceTetIs;
587  break;
588  }
589  }
590 
591  if (tetContainingPt.cell() != -1)
592  {
593  break;
594  }
595  }
596 
597  return tetContainingPt;
598 }
599 
600 
601 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
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:434
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:323
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static scalar minQuality(const polyMesh &mesh, const point &cC, const label fI, const bool isOwner, const label faceBasePtI)
Given a face and cc and starting index for triangulation determine.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const labelList & patchID() const
Per boundary face label the patch index.
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:100
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:111
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
Definition: tetIndicesI.H:98
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:1131
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:1169
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
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:1156
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.
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, &mergedCyclicPolyPatch::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
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
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.