polyMeshTetDecomposition.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) 2011-2013 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  (
314  "labelList"
315  "polyMeshTetDecomposition::findFaceBasePts"
316  "("
317  "const polyMesh& mesh, "
318  "scalar tol, "
319  "bool report"
320  ")"
321  ) << "Coupled face base point exchange failure for face "
322  << fI
323  << abort(FatalError);
324  }
325 
326  if (bFTetBasePtI < 1)
327  {
328  // If the base point is -1, it should be left as such to
329  // indicate a problem, if it is 0, then no action is required.
330 
331  continue;
332  }
333 
334  label patchI = mesh.boundaryMesh().patchID()[bFI];
335 
336  if (patches[patchI].coupled())
337  {
338  const coupledPolyPatch& pp =
339  refCast<const coupledPolyPatch>(patches[patchI]);
340 
341  // Calculated base points on coupled faces are those of
342  // the owner patch face. They need to be reindexed to for
343  // the non-owner face, which has the opposite order.
344 
345  // So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
346 
347  // i.e.:
348 
349  // owner coupledPolyPatch face
350  // face (a b c d e f)
351  // fPtI 0 1 2 3 4 5
352  // +
353  // #
354 
355  // neighbour coupledPolyPatch face
356  // face (a f e d c b)
357  // fPtI 0 1 2 3 4 5
358  // +
359  // #
360  // +: 6 - 1 = 5
361  // #: 6 - 2 = 4
362 
363  if (!pp.owner())
364  {
365  bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
366  }
367  }
368  }
369 
370  return tetBasePtIs;
371 }
372 
373 
375 (
376  const polyMesh& mesh,
377  scalar tol,
378  const bool report,
379  labelHashSet* setPtr
380 )
381 {
382  const labelList& own = mesh.faceOwner();
383  const labelList& nei = mesh.faceNeighbour();
384  const polyBoundaryMesh& patches = mesh.boundaryMesh();
385 
386  const vectorField& cc = mesh.cellCentres();
387  const vectorField& fc = mesh.faceCentres();
388 
389  // Calculate coupled cell centre
390  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
391 
392  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
393  {
394  neiCc[faceI - mesh.nInternalFaces()] = cc[own[faceI]];
395  }
396 
397  syncTools::swapBoundaryFacePositions(mesh, neiCc);
398 
399  const faceList& fcs = mesh.faces();
400 
401  const pointField& p = mesh.points();
402 
403  label nErrorTets = 0;
404 
405  forAll(fcs, faceI)
406  {
407  const face& f = fcs[faceI];
408 
409  forAll(f, fPtI)
410  {
411  scalar tetQual = tetPointRef
412  (
413  p[f[fPtI]],
414  p[f.nextLabel(fPtI)],
415  fc[faceI],
416  cc[own[faceI]]
417  ).quality();
418 
419  if (tetQual > -tol)
420  {
421  if (setPtr)
422  {
423  setPtr->insert(faceI);
424  }
425 
426  nErrorTets++;
427  break; // no need to check other tets
428  }
429  }
430 
431  if (mesh.isInternalFace(faceI))
432  {
433  // Create the neighbour tet - it will have positive volume
434  const face& f = fcs[faceI];
435 
436  forAll(f, fPtI)
437  {
438  scalar tetQual = tetPointRef
439  (
440  p[f[fPtI]],
441  p[f.nextLabel(fPtI)],
442  fc[faceI],
443  cc[nei[faceI]]
444  ).quality();
445 
446  if (tetQual < tol)
447  {
448  if (setPtr)
449  {
450  setPtr->insert(faceI);
451  }
452 
453  nErrorTets++;
454  break;
455  }
456  }
457 
458  if (findSharedBasePoint(mesh, faceI, tol, report) == -1)
459  {
460  if (setPtr)
461  {
462  setPtr->insert(faceI);
463  }
464 
465  nErrorTets++;
466  }
467  }
468  else
469  {
470  label patchI = patches.patchID()[faceI - mesh.nInternalFaces()];
471 
472  if (patches[patchI].coupled())
473  {
474  if
475  (
476  findSharedBasePoint
477  (
478  mesh,
479  faceI,
480  neiCc[faceI - mesh.nInternalFaces()],
481  tol,
482  report
483  ) == -1
484  )
485  {
486  if (setPtr)
487  {
488  setPtr->insert(faceI);
489  }
490 
491  nErrorTets++;
492  }
493  }
494  else
495  {
496  if (findBasePoint(mesh, faceI, tol, report) == -1)
497  {
498  if (setPtr)
499  {
500  setPtr->insert(faceI);
501  }
502 
503  nErrorTets++;
504  }
505  }
506  }
507  }
508 
509  reduce(nErrorTets, sumOp<label>());
510 
511  if (nErrorTets > 0)
512  {
513  if (report)
514  {
515  Info<< " ***Error in face tets: "
516  << nErrorTets << " faces with low quality or negative volume "
517  << "decomposition tets." << endl;
518  }
519 
520  return true;
521  }
522  else
523  {
524  if (report)
525  {
526  Info<< " Face tets OK." << endl;
527  }
528 
529  return false;
530  }
531 }
532 
533 
535 (
536  const polyMesh& mesh,
537  label fI,
538  label cI
539 )
540 {
541  static label nWarnings = 0;
542  static const label maxWarnings = 100;
543 
544  const faceList& pFaces = mesh.faces();
545  const labelList& pOwner = mesh.faceOwner();
546 
547  const face& f = pFaces[fI];
548 
549  label nTets = f.size() - 2;
550 
551  List<tetIndices> faceTets(nTets);
552 
553  bool own = (pOwner[fI] == cI);
554 
555  label tetBasePtI = mesh.tetBasePtIs()[fI];
556 
557  if (tetBasePtI == -1)
558  {
559  if (nWarnings < maxWarnings)
560  {
561  WarningIn
562  (
563  "List<tetIndices> "
564  "polyMeshTetDecomposition::faceTetIndices"
565  "("
566  "const polyMesh&, "
567  "label, "
568  "label"
569  ")"
570  ) << "No base point for face " << fI << ", " << f
571  << ", produces a valid tet decomposition."
572  << endl;
573  nWarnings++;
574  }
575  if (nWarnings == maxWarnings)
576  {
577  Warning<< "Suppressing any further warnings." << endl;
578  nWarnings++;
579  }
580 
581  tetBasePtI = 0;
582  }
583 
584  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
585  {
586  tetIndices& faceTetIs = faceTets[tetPtI - 1];
587 
588  label facePtI = (tetPtI + tetBasePtI) % f.size();
589  label otherFacePtI = f.fcIndex(facePtI);
590 
591  faceTetIs.cell() = cI;
592 
593  faceTetIs.face() = fI;
594 
595  faceTetIs.faceBasePt() = tetBasePtI;
596 
597  if (own)
598  {
599  faceTetIs.facePtA() = facePtI;
600  faceTetIs.facePtB() = otherFacePtI;
601  }
602  else
603  {
604  faceTetIs.facePtA() = otherFacePtI;
605  faceTetIs.facePtB() = facePtI;
606  }
607 
608  faceTetIs.tetPt() = tetPtI;
609  }
610 
611  return faceTets;
612 }
613 
614 
616 (
617  const polyMesh& mesh,
618  const label fI,
619  const label cI,
620  const label tetPtI
621 )
622 {
623  static label nWarnings = 0;
624  static const label maxWarnings = 100;
625 
626  const face& f = mesh.faces()[fI];
627  bool own = (mesh.faceOwner()[fI] == cI);
628  label tetBasePtI = mesh.tetBasePtIs()[fI];
629  if (tetBasePtI == -1)
630  {
631  if (nWarnings < maxWarnings)
632  {
633  WarningIn
634  (
635  "tetIndices "
636  "polyMeshTetDecomposition::triangleTetIndices"
637  "("
638  "const polyMesh&, "
639  "label, "
640  "label, "
641  "label"
642  ")"
643  ) << "No base point for face " << fI << ", " << f
644  << ", produces a valid tet decomposition."
645  << endl;
646  nWarnings++;
647  }
648  if (nWarnings == maxWarnings)
649  {
650  Warning<< "Suppressing any further warnings." << endl;
651  nWarnings++;
652  }
653 
654  tetBasePtI = 0;
655  }
656 
657  tetIndices faceTetIs;
658 
659  label facePtI = (tetPtI + tetBasePtI) % f.size();
660  label otherFacePtI = f.fcIndex(facePtI);
661 
662  faceTetIs.cell() = cI;
663 
664  faceTetIs.face() = fI;
665 
666  faceTetIs.faceBasePt() = tetBasePtI;
667 
668  if (own)
669  {
670  faceTetIs.facePtA() = facePtI;
671  faceTetIs.facePtB() = otherFacePtI;
672  }
673  else
674  {
675  faceTetIs.facePtA() = otherFacePtI;
676  faceTetIs.facePtB() = facePtI;
677  }
678 
679  faceTetIs.tetPt() = tetPtI;
680 
681  return faceTetIs;
682 }
683 
684 
686 (
687  const polyMesh& mesh,
688  label cI
689 )
690 {
691  const faceList& pFaces = mesh.faces();
692  const cellList& pCells = mesh.cells();
693 
694  const cell& thisCell = pCells[cI];
695 
696  label nTets = 0;
697 
698  forAll(thisCell, cFI)
699  {
700  nTets += pFaces[thisCell[cFI]].size() - 2;
701  }
702 
703  DynamicList<tetIndices> cellTets(nTets);
704 
705  forAll(thisCell, cFI)
706  {
707  label fI = thisCell[cFI];
708 
709  cellTets.append(faceTetIndices(mesh, fI, cI));
710  }
711 
712  return cellTets;
713 }
714 
715 
717 (
718  const polyMesh& mesh,
719  label cI,
720  const point& pt
721 )
722 {
723  const faceList& pFaces = mesh.faces();
724  const cellList& pCells = mesh.cells();
725 
726  const cell& thisCell = pCells[cI];
727 
728  tetIndices tetContainingPt;
729 
730 
731  forAll(thisCell, cFI)
732  {
733  label fI = thisCell[cFI];
734  const face& f = pFaces[fI];
735 
736  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
737  {
738  // Get tetIndices of face triangle
739  tetIndices faceTetIs
740  (
741  triangleTetIndices
742  (
743  mesh,
744  fI,
745  cI,
746  tetPtI
747  )
748  );
749 
750  // Check if inside
751  if (faceTetIs.tet(mesh).inside(pt))
752  {
753  tetContainingPt = faceTetIs;
754  break;
755  }
756  }
757 
758  if (tetContainingPt.cell() != -1)
759  {
760  break;
761  }
762  }
763 
764  return tetContainingPt;
765 }
766 
767 
768 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label cell() const
Return the cell.
Definition: tetIndicesI.H:30
static tetIndices triangleTetIndices(const polyMesh &mesh, label fI, label cI, const label tetPtI)
Return the tet decomposition of the given triangle of the given face.
label nFaces() const
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
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.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
labelList f(nPoints)
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
Find a suitable base point for each face for decomposition.
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:864
label faceBasePt() const
Return the face base point.
Definition: tetIndicesI.H:42
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
static const scalar minTetQuality
Minimum tetrahedron quality.
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 bool checkFaceTets(const polyMesh &mesh, scalar tol=minTetQuality, const bool report=false, labelHashSet *setPtr=NULL)
Check face-decomposition tet volume.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const cellList & cells() const
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
Definition: tetrahedronI.H:456
const vectorField & cellCentres() const
messageStream Info
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
tetrahedron< point, const point & > tetPointRef
Definition: tetrahedron.H:78
patches[0]
messageStream Warning
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:231
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
Definition: tetIndicesI.H:66
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
Foam::polyBoundaryMesh.
volScalarField & p
Definition: createFields.H:51
#define forAll(list, i)
Definition: UList.H:421
A tetrahedron primitive.
Definition: tetrahedron.H:62
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
label nInternalFaces() const
const labelList & patchID() const
Per boundary face label the patch index.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
error FatalError
label facePtA() const
Return face point A.
Definition: tetIndicesI.H:48
virtual bool owner() const =0
Does this side own the patch ?
A List obtained as a section of another List.
Definition: SubList.H:53
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
label tetPt() const
Return the characterising tetPtI.
Definition: tetIndicesI.H:60
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
label facePtB() const
Return face point B.
Definition: tetIndicesI.H:54
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const vectorField & faceCentres() const
label face() const
Return the face.
Definition: tetIndicesI.H:36
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116