cellClassification.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 
26 #include "cellClassification.H"
27 #include "triSurfaceSearch.H"
28 #include "indexedOctree.H"
29 #include "treeDataFace.H"
30 #include "meshSearch.H"
31 #include "cellInfo.H"
32 #include "polyMesh.H"
33 #include "MeshWave.H"
34 #include "ListOps.H"
35 #include "meshTools.H"
36 #include "cpuTime.H"
37 #include "triSurface.H"
38 #include "globalMeshData.H"
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 namespace Foam
43 {
44 defineTypeNameAndDebug(cellClassification, 0);
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 Foam::label Foam::cellClassification::count
51 (
52  const labelList& elems,
53  const label elem
54 )
55 {
56  label cnt = 0;
57 
58  forAll(elems, elemI)
59  {
60  if (elems[elemI] == elem)
61  {
62  cnt++;
63  }
64  }
65  return cnt;
66 
67 }
68 
69 
70 // Mark all faces that are cut by the surface. Two pass:
71 // Pass1: mark all mesh edges that intersect surface (quick since triangle
72 // pierce test).
73 // Pass2: Check for all point neighbours of these faces whether any of their
74 // faces are pierced.
75 Foam::boolList Foam::cellClassification::markFaces
76 (
77  const triSurfaceSearch& search
78 ) const
79 {
80  cpuTime timer;
81 
82  boolList cutFace(mesh_.nFaces(), false);
83 
84  label nCutFaces = 0;
85 
86  // Intersect mesh edges with surface (is fast) and mark all faces that
87  // use edge.
88 
89  forAll(mesh_.edges(), edgeI)
90  {
91  if (debug && (edgeI % 10000 == 0))
92  {
93  Pout<< "Intersecting mesh edge " << edgeI << " with surface"
94  << endl;
95  }
96 
97  const edge& e = mesh_.edges()[edgeI];
98 
99  const point& p0 = mesh_.points()[e.start()];
100  const point& p1 = mesh_.points()[e.end()];
101 
102  pointIndexHit pHit(search.tree().findLineAny(p0, p1));
103 
104  if (pHit.hit())
105  {
106  const labelList& myFaces = mesh_.edgeFaces()[edgeI];
107 
108  forAll(myFaces, myFaceI)
109  {
110  label faceI = myFaces[myFaceI];
111 
112  if (!cutFace[faceI])
113  {
114  cutFace[faceI] = true;
115 
116  nCutFaces++;
117  }
118  }
119  }
120  }
121 
122  if (debug)
123  {
124  Pout<< "Intersected edges of mesh with surface in = "
125  << timer.cpuTimeIncrement() << " s\n" << endl << endl;
126  }
127 
128  //
129  // Construct octree on faces that have not yet been marked as cut
130  //
131 
132  labelList allFaces(mesh_.nFaces() - nCutFaces);
133 
134  label allFaceI = 0;
135 
136  forAll(cutFace, faceI)
137  {
138  if (!cutFace[faceI])
139  {
140  allFaces[allFaceI++] = faceI;
141  }
142  }
143 
144  if (debug)
145  {
146  Pout<< "Testing " << allFaceI << " faces for piercing by surface"
147  << endl;
148  }
149 
150  treeBoundBox allBb(mesh_.points());
151 
152  // Extend domain slightly (also makes it 3D if was 2D)
153  scalar tol = 1e-6 * allBb.avgDim();
154 
155  point& bbMin = allBb.min();
156  bbMin.x() -= tol;
157  bbMin.y() -= tol;
158  bbMin.z() -= tol;
159 
160  point& bbMax = allBb.max();
161  bbMax.x() += 2*tol;
162  bbMax.y() += 2*tol;
163  bbMax.z() += 2*tol;
164 
165  indexedOctree<treeDataFace> faceTree
166  (
167  treeDataFace(false, mesh_, allFaces),
168  allBb, // overall search domain
169  8, // maxLevel
170  10, // leafsize
171  3.0 // duplicity
172  );
173 
174  const triSurface& surf = search.surface();
175  const edgeList& edges = surf.edges();
176  const pointField& localPoints = surf.localPoints();
177 
178  label nAddFaces = 0;
179 
180  forAll(edges, edgeI)
181  {
182  if (debug && (edgeI % 10000 == 0))
183  {
184  Pout<< "Intersecting surface edge " << edgeI
185  << " with mesh faces" << endl;
186  }
187  const edge& e = edges[edgeI];
188 
189  const point& start = localPoints[e.start()];
190  const point& end = localPoints[e.end()];
191 
192  vector edgeNormal(end - start);
193  const scalar edgeMag = mag(edgeNormal);
194  const vector smallVec = 1e-9*edgeNormal;
195 
196  edgeNormal /= edgeMag+VSMALL;
197 
198  // Current start of pierce test
199  point pt = start;
200 
201  while (true)
202  {
203  pointIndexHit pHit(faceTree.findLine(pt, end));
204 
205  if (!pHit.hit())
206  {
207  break;
208  }
209  else
210  {
211  label faceI = faceTree.shapes().faceLabels()[pHit.index()];
212 
213  if (!cutFace[faceI])
214  {
215  cutFace[faceI] = true;
216 
217  nAddFaces++;
218  }
219 
220  // Restart from previous endpoint
221  pt = pHit.hitPoint() + smallVec;
222 
223  if (((pt-start) & edgeNormal) >= edgeMag)
224  {
225  break;
226  }
227  }
228  }
229  }
230 
231  if (debug)
232  {
233  Pout<< "Detected an additional " << nAddFaces << " faces cut"
234  << endl;
235 
236  Pout<< "Intersected edges of surface with mesh faces in = "
237  << timer.cpuTimeIncrement() << " s\n" << endl << endl;
238  }
239 
240  return cutFace;
241 }
242 
243 
244 // Determine faces cut by surface and use to divide cells into types. See
245 // cellInfo. All cells reachable from outsidePts are considered to be of type
246 // 'outside'
247 void Foam::cellClassification::markCells
248 (
249  const meshSearch& queryMesh,
250  const boolList& piercedFace,
251  const pointField& outsidePts
252 )
253 {
254  // Use meshwave to partition mesh, starting from outsidePt
255 
256 
257  // Construct null; sets type to NOTSET
258  List<cellInfo> cellInfoList(mesh_.nCells());
259 
260  // Mark cut cells first
261  forAll(piercedFace, faceI)
262  {
263  if (piercedFace[faceI])
264  {
265  cellInfoList[mesh_.faceOwner()[faceI]] =
266  cellInfo(cellClassification::CUT);
267 
268  if (mesh_.isInternalFace(faceI))
269  {
270  cellInfoList[mesh_.faceNeighbour()[faceI]] =
271  cellInfo(cellClassification::CUT);
272  }
273  }
274  }
275 
276  //
277  // Mark cells containing outside points as being outside
278  //
279 
280  // Coarse guess number of faces
281  labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
282 
283  forAll(outsidePts, outsidePtI)
284  {
285  // Use linear search for points.
286  label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
287 
288  if (returnReduce(cellI, maxOp<label>()) == -1)
289  {
291  (
292  "List<cellClassification::cType> markCells"
293  "(const meshSearch&, const boolList&, const pointField&)"
294  ) << "outsidePoint " << outsidePts[outsidePtI]
295  << " is not inside any cell"
296  << nl << "It might be on a face or outside the geometry"
297  << exit(FatalError);
298  }
299 
300  if (cellI >= 0)
301  {
302  cellInfoList[cellI] = cellInfo(cellClassification::OUTSIDE);
303 
304  // Mark faces of cellI
305  const labelList& myFaces = mesh_.cells()[cellI];
306  forAll(myFaces, myFaceI)
307  {
308  outsideFacesMap.insert(myFaces[myFaceI]);
309  }
310  }
311  }
312 
313 
314  //
315  // Mark faces to start wave from
316  //
317 
318  labelList changedFaces(outsideFacesMap.toc());
319 
320  List<cellInfo> changedFacesInfo
321  (
322  changedFaces.size(),
324  );
325 
326  MeshWave<cellInfo> cellInfoCalc
327  (
328  mesh_,
329  changedFaces, // Labels of changed faces
330  changedFacesInfo, // Information on changed faces
331  cellInfoList, // Information on all cells
332  mesh_.globalData().nTotalCells()+1 // max iterations
333  );
334 
335  // Get information out of cellInfoList
336  const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
337 
338  forAll(allInfo, cellI)
339  {
340  label t = allInfo[cellI].type();
341 
343  {
345  }
346  operator[](cellI) = t;
347  }
348 }
349 
350 
351 void Foam::cellClassification::classifyPoints
352 (
353  const label meshType,
354  const labelList& cellType,
355  List<pointStatus>& pointSide
356 ) const
357 {
358  pointSide.setSize(mesh_.nPoints());
359 
360  forAll(mesh_.pointCells(), pointI)
361  {
362  const labelList& pCells = mesh_.pointCells()[pointI];
363 
364  pointSide[pointI] = UNSET;
365 
366  forAll(pCells, i)
367  {
368  label type = cellType[pCells[i]];
369 
370  if (type == meshType)
371  {
372  if (pointSide[pointI] == UNSET)
373  {
374  pointSide[pointI] = MESH;
375  }
376  else if (pointSide[pointI] == NONMESH)
377  {
378  pointSide[pointI] = MIXED;
379 
380  break;
381  }
382  }
383  else
384  {
385  if (pointSide[pointI] == UNSET)
386  {
387  pointSide[pointI] = NONMESH;
388  }
389  else if (pointSide[pointI] == MESH)
390  {
391  pointSide[pointI] = MIXED;
392 
393  break;
394  }
395  }
396  }
397  }
398 }
399 
400 
401 bool Foam::cellClassification::usesMixedPointsOnly
402 (
403  const List<pointStatus>& pointSide,
404  const label cellI
405 ) const
406 {
407  const faceList& faces = mesh_.faces();
408 
409  const cell& cFaces = mesh_.cells()[cellI];
410 
411  forAll(cFaces, cFaceI)
412  {
413  const face& f = faces[cFaces[cFaceI]];
414 
415  forAll(f, fp)
416  {
417  if (pointSide[f[fp]] != MIXED)
418  {
419  return false;
420  }
421  }
422  }
423 
424  // All points are mixed.
425  return true;
426 }
427 
428 
429 void Foam::cellClassification::getMeshOutside
430 (
431  const label meshType,
432  faceList& outsideFaces,
433  labelList& outsideOwner
434 ) const
435 {
436  const labelList& own = mesh_.faceOwner();
437  const labelList& nbr = mesh_.faceNeighbour();
438  const faceList& faces = mesh_.faces();
439 
440  outsideFaces.setSize(mesh_.nFaces());
441  outsideOwner.setSize(mesh_.nFaces());
442  label outsideI = 0;
443 
444  // Get faces on interface between meshType and non-meshType
445 
446  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
447  {
448  label ownType = operator[](own[faceI]);
449  label nbrType = operator[](nbr[faceI]);
450 
451  if (ownType == meshType && nbrType != meshType)
452  {
453  outsideFaces[outsideI] = faces[faceI];
454  outsideOwner[outsideI] = own[faceI]; // meshType cell
455  outsideI++;
456  }
457  else if (ownType != meshType && nbrType == meshType)
458  {
459  outsideFaces[outsideI] = faces[faceI];
460  outsideOwner[outsideI] = nbr[faceI]; // meshType cell
461  outsideI++;
462  }
463  }
464 
465  // Get faces on outside of real mesh with cells of meshType.
466 
467  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
468  {
469  if (operator[](own[faceI]) == meshType)
470  {
471  outsideFaces[outsideI] = faces[faceI];
472  outsideOwner[outsideI] = own[faceI]; // meshType cell
473  outsideI++;
474  }
475  }
476  outsideFaces.setSize(outsideI);
477  outsideOwner.setSize(outsideI);
478 }
479 
480 
481 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
482 
483 // Construct from mesh and surface and point(s) on outside
485 (
486  const polyMesh& mesh,
487  const meshSearch& meshQuery,
488  const triSurfaceSearch& surfQuery,
489  const pointField& outsidePoints
490 )
491 :
493  mesh_(mesh)
494 {
495  markCells
496  (
497  meshQuery,
498  markFaces(surfQuery),
499  outsidePoints
500  );
501 }
502 
503 
504 // Construct from mesh and meshType.
506 (
507  const polyMesh& mesh,
508  const labelList& cellType
509 )
510 :
511  labelList(cellType),
512  mesh_(mesh)
513 {
514  if (mesh_.nCells() != size())
515  {
517  (
518  "cellClassification::cellClassification"
519  "(const polyMesh&, const labelList&)"
520  ) << "Number of elements of cellType argument is not equal to the"
521  << " number of cells" << abort(FatalError);
522  }
523 }
524 
525 
526 // Construct as copy
528 :
529  labelList(cType),
530  mesh_(cType.mesh())
531 {}
532 
533 
534 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
535 
536 
537 // Makes cutCells further than nLayers away from meshType to
538 // fillType. Returns number of cells changed.
540 (
541  const label nLayers,
542  const label meshType,
543  const label fillType
544 )
545 {
546  // Temporary cell type for growing.
547  labelList newCellType(*this);
548 
549 // // Split types into outside and rest
550 // forAll(*this, cellI)
551 // {
552 // label type = operator[](cellI);
553 //
554 // if (type == meshType)
555 // {
556 // newCellType[cellI] = type;
557 // }
558 // else
559 // {
560 // newCellType[cellI] = fillType;
561 // }
562 // }
563 
564  newCellType = *this;
565 
566 
567  // Do point-cell-point walk on newCellType out from meshType cells
568  for (label iter = 0; iter < nLayers; iter++)
569  {
570  // Get status of points: visible from meshType (pointSide == MESH)
571  // or fillType/cutCells cells (pointSide == NONMESH) or
572  // both (pointSide == MIXED)
573  List<pointStatus> pointSide(mesh_.nPoints());
574  classifyPoints(meshType, newCellType, pointSide);
575 
576  // Grow layer of outside cells
577  forAll(pointSide, pointI)
578  {
579  if (pointSide[pointI] == MIXED)
580  {
581  // Make cut
582  const labelList& pCells = mesh_.pointCells()[pointI];
583 
584  forAll(pCells, i)
585  {
586  label type = newCellType[pCells[i]];
587 
588  if (type == cellClassification::CUT)
589  {
590  // Found cut cell sharing point with
591  // mesh type cell. Grow.
592  newCellType[pCells[i]] = meshType;
593  }
594  }
595  }
596  }
597  }
598 
599  // Merge newCellType into *this:
600  // - leave meshType cells intact
601  // - leave nonMesh cells intact
602  // - make cutcells fillType if they were not marked in newCellType
603 
604  label nChanged = 0;
605 
606  forAll(newCellType, cellI)
607  {
608  if (operator[](cellI) == cellClassification::CUT)
609  {
610  if (newCellType[cellI] != meshType)
611  {
612  // Cell was cutCell but further than nLayers away from
613  // meshType. Convert to fillType.
614  operator[](cellI) = fillType;
615  nChanged++;
616  }
617  }
618  }
619 
620  return nChanged;
621 }
622 
623 
624 // Grow surface by pointNeighbours of type meshType
626 (
627  const label meshType,
628  const label fillType
629 )
630 {
631  boolList hasMeshType(mesh_.nPoints(), false);
632 
633  // Mark points used by meshType cells
634 
635  forAll(mesh_.pointCells(), pointI)
636  {
637  const labelList& myCells = mesh_.pointCells()[pointI];
638 
639  // Check if one of cells has meshType
640  forAll(myCells, myCellI)
641  {
642  label type = operator[](myCells[myCellI]);
643 
644  if (type == meshType)
645  {
646  hasMeshType[pointI] = true;
647 
648  break;
649  }
650  }
651  }
652 
653  // Change neighbours of marked points
654 
655  label nChanged = 0;
656 
657  forAll(hasMeshType, pointI)
658  {
659  if (hasMeshType[pointI])
660  {
661  const labelList& myCells = mesh_.pointCells()[pointI];
662 
663  forAll(myCells, myCellI)
664  {
665  if (operator[](myCells[myCellI]) != meshType)
666  {
667  operator[](myCells[myCellI]) = fillType;
668 
669  nChanged++;
670  }
671  }
672  }
673  }
674  return nChanged;
675 }
676 
677 
678 // Check all points used by cells of meshType for use of at least one point
679 // which is not on the outside. If all points are on the outside of the mesh
680 // this would probably result in a flattened cell when projecting the cell
681 // onto the surface.
683 (
684  const label meshType,
685  const label fillType,
686  const label maxIter
687 )
688 {
689  label nTotChanged = 0;
690 
691  for (label iter = 0; iter < maxIter; iter++)
692  {
693  label nChanged = 0;
694 
695  // Get status of points: visible from meshType or non-meshType cells
696  List<pointStatus> pointSide(mesh_.nPoints());
697  classifyPoints(meshType, *this, pointSide);
698 
699  // Check all cells using mixed point type for whether they use mixed
700  // points only. Note: could probably speed this up by counting number
701  // of mixed verts per face and mixed faces per cell or something?
702  forAll(pointSide, pointI)
703  {
704  if (pointSide[pointI] == MIXED)
705  {
706  const labelList& pCells = mesh_.pointCells()[pointI];
707 
708  forAll(pCells, i)
709  {
710  label cellI = pCells[i];
711 
712  if (operator[](cellI) == meshType)
713  {
714  if (usesMixedPointsOnly(pointSide, cellI))
715  {
716  operator[](cellI) = fillType;
717 
718  nChanged++;
719  }
720  }
721  }
722  }
723  }
724  nTotChanged += nChanged;
725 
726  Pout<< "removeHangingCells : changed " << nChanged
727  << " hanging cells" << endl;
728 
729  if (nChanged == 0)
730  {
731  break;
732  }
733  }
734 
735  return nTotChanged;
736 }
737 
738 
740 (
741  const label meshType,
742  const label fillType,
743  const label maxIter
744 )
745 {
746  label nTotChanged = 0;
747 
748  for (label iter = 0; iter < maxIter; iter++)
749  {
750  // Get interface between meshType cells and non-meshType cells as a list
751  // of faces and for each face the cell which is the meshType.
752  faceList outsideFaces;
753  labelList outsideOwner;
754 
755  getMeshOutside(meshType, outsideFaces, outsideOwner);
756 
757  // Build primitivePatch out of it and check it for problems.
758  primitiveFacePatch fp(outsideFaces, mesh_.points());
759 
760  const labelListList& edgeFaces = fp.edgeFaces();
761 
762  label nChanged = 0;
763 
764  // Check all edgeFaces for non-manifoldness
765 
766  forAll(edgeFaces, edgeI)
767  {
768  const labelList& eFaces = edgeFaces[edgeI];
769 
770  if (eFaces.size() > 2)
771  {
772  // patch connected through pinched edge. Remove first face using
773  // edge (and not yet changed)
774 
775  forAll(eFaces, i)
776  {
777  label patchFaceI = eFaces[i];
778 
779  label ownerCell = outsideOwner[patchFaceI];
780 
781  if (operator[](ownerCell) == meshType)
782  {
783  operator[](ownerCell) = fillType;
784 
785  nChanged++;
786 
787  break;
788  }
789  }
790  }
791  }
792 
793  nTotChanged += nChanged;
794 
795  Pout<< "fillRegionEdges : changed " << nChanged
796  << " cells using multiply connected edges" << endl;
797 
798  if (nChanged == 0)
799  {
800  break;
801  }
802  }
803 
804  return nTotChanged;
805 }
806 
807 
809 (
810  const label meshType,
811  const label fillType,
812  const label maxIter
813 )
814 {
815  label nTotChanged = 0;
816 
817  for (label iter = 0; iter < maxIter; iter++)
818  {
819  // Get interface between meshType cells and non-meshType cells as a list
820  // of faces and for each face the cell which is the meshType.
821  faceList outsideFaces;
822  labelList outsideOwner;
823 
824  getMeshOutside(meshType, outsideFaces, outsideOwner);
825 
826  // Build primitivePatch out of it and check it for problems.
827  primitiveFacePatch fp(outsideFaces, mesh_.points());
828 
829  labelHashSet nonManifoldPoints;
830 
831  // Check for non-manifold points.
832  fp.checkPointManifold(false, &nonManifoldPoints);
833 
834  const Map<label>& meshPointMap = fp.meshPointMap();
835 
836  label nChanged = 0;
837 
838  forAllConstIter(labelHashSet, nonManifoldPoints, iter)
839  {
840  // Find a face on fp using point and remove it.
841  const label patchPointI = meshPointMap[iter.key()];
842 
843  const labelList& pFaces = fp.pointFaces()[patchPointI];
844 
845  // Remove any face using conflicting point. Does first face which
846  // has not yet been done. Could be more intelligent and decide which
847  // one would be best to remove.
848  forAll(pFaces, i)
849  {
850  const label patchFaceI = pFaces[i];
851  const label ownerCell = outsideOwner[patchFaceI];
852 
853  if (operator[](ownerCell) == meshType)
854  {
855  operator[](ownerCell) = fillType;
856 
857  nChanged++;
858  break;
859  }
860  }
861  }
862 
863  nTotChanged += nChanged;
864 
865  Pout<< "fillRegionPoints : changed " << nChanged
866  << " cells using multiply connected points" << endl;
867 
868  if (nChanged == 0)
869  {
870  break;
871  }
872  }
873 
874  return nTotChanged;
875 }
876 
877 
879 {
880  os << "Cells:" << size() << endl
881  << " notset : " << count(*this, NOTSET) << endl
882  << " cut : " << count(*this, CUT) << endl
883  << " inside : " << count(*this, INSIDE) << endl
884  << " outside : " << count(*this, OUTSIDE) << endl;
885 }
886 
887 
888 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
889 
891 {
893 }
894 
895 
896 // ************************************************************************* //
cellClassification(const polyMesh &mesh, const meshSearch &meshQuery, const triSurfaceSearch &surfQuery, const pointField &outsidePoints)
Construct from mesh and surface and point(s) on outside.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
vector point
Point is a vector.
Definition: point.H:41
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
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 growSurface(const label meshType, const label fillType)
Sets vertex neighbours of meshType cells to fillType.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
label nCells() const
Various functions to operate on Lists.
const labelListList & pointCells() const
static const Vector max
Definition: Vector.H:82
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
label fillRegionPoints(const label meshType, const label fillType, const label maxIter)
Find regionPoints and fill all neighbours. Iterate until nothing.
List< face > faceList
Definition: faceListFwd.H:43
void operator=(const cellClassification &)
void writeStats(Ostream &os) const
Write statistics on cell types to Ostream.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
T & operator[](const label)
Return element of UList.
Definition: UListI.H:163
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
#define forAll(list, i)
Definition: UList.H:421
const Cmpt & x() const
Definition: VectorI.H:65
&#39;Cuts&#39; a mesh with a surface.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
Definition: meshSearch.H:57
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label fillRegionEdges(const label meshType, const label fillType, const label maxIter)
Find regionEdges and fill one neighbour. Iterate until nothing.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
A list of faces which address into the list of points.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static const Vector min
Definition: Vector.H:83
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
label size() const
Return the number of elements in the UList.
Definition: ListI.H:83
label trimCutCells(const label nLayers, const label meshType, const label fillType)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label fillHangingCells(const label meshType, const label fillType, const label maxIter)
Find hanging cells (cells with all points on outside) and set their.
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
label nPoints() const
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Helper class to search on triSurface.
void operator=(const UList< label > &)
Assignment from UList operator. Takes linear time.