cellClassification.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 
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  << "outsidePoint " << outsidePts[outsidePtI]
292  << " is not inside any cell"
293  << nl << "It might be on a face or outside the geometry"
294  << exit(FatalError);
295  }
296 
297  if (celli >= 0)
298  {
299  cellInfoList[celli] = cellInfo(cellClassification::OUTSIDE);
300 
301  // Mark faces of celli
302  const labelList& myFaces = mesh_.cells()[celli];
303  forAll(myFaces, myFacei)
304  {
305  outsideFacesMap.insert(myFaces[myFacei]);
306  }
307  }
308  }
309 
310 
311  //
312  // Mark faces to start wave from
313  //
314 
315  labelList changedFaces(outsideFacesMap.toc());
316 
317  List<cellInfo> changedFacesInfo
318  (
319  changedFaces.size(),
321  );
322 
323  MeshWave<cellInfo> cellInfoCalc
324  (
325  mesh_,
326  changedFaces, // Labels of changed faces
327  changedFacesInfo, // Information on changed faces
328  cellInfoList, // Information on all cells
329  mesh_.globalData().nTotalCells()+1 // max iterations
330  );
331 
332  // Get information out of cellInfoList
333  const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
334 
335  forAll(allInfo, celli)
336  {
337  label t = allInfo[celli].type();
338 
340  {
342  }
343  operator[](celli) = t;
344  }
345 }
346 
347 
348 void Foam::cellClassification::classifyPoints
349 (
350  const label meshType,
351  const labelList& cellType,
352  List<pointStatus>& pointSide
353 ) const
354 {
355  pointSide.setSize(mesh_.nPoints());
356 
357  forAll(mesh_.pointCells(), pointi)
358  {
359  const labelList& pCells = mesh_.pointCells()[pointi];
360 
361  pointSide[pointi] = UNSET;
362 
363  forAll(pCells, i)
364  {
365  label type = cellType[pCells[i]];
366 
367  if (type == meshType)
368  {
369  if (pointSide[pointi] == UNSET)
370  {
371  pointSide[pointi] = MESH;
372  }
373  else if (pointSide[pointi] == NONMESH)
374  {
375  pointSide[pointi] = MIXED;
376 
377  break;
378  }
379  }
380  else
381  {
382  if (pointSide[pointi] == UNSET)
383  {
384  pointSide[pointi] = NONMESH;
385  }
386  else if (pointSide[pointi] == MESH)
387  {
388  pointSide[pointi] = MIXED;
389 
390  break;
391  }
392  }
393  }
394  }
395 }
396 
397 
398 bool Foam::cellClassification::usesMixedPointsOnly
399 (
400  const List<pointStatus>& pointSide,
401  const label celli
402 ) const
403 {
404  const faceList& faces = mesh_.faces();
405 
406  const cell& cFaces = mesh_.cells()[celli];
407 
408  forAll(cFaces, cFacei)
409  {
410  const face& f = faces[cFaces[cFacei]];
411 
412  forAll(f, fp)
413  {
414  if (pointSide[f[fp]] != MIXED)
415  {
416  return false;
417  }
418  }
419  }
420 
421  // All points are mixed.
422  return true;
423 }
424 
425 
426 void Foam::cellClassification::getMeshOutside
427 (
428  const label meshType,
429  faceList& outsideFaces,
430  labelList& outsideOwner
431 ) const
432 {
433  const labelList& own = mesh_.faceOwner();
434  const labelList& nbr = mesh_.faceNeighbour();
435  const faceList& faces = mesh_.faces();
436 
437  outsideFaces.setSize(mesh_.nFaces());
438  outsideOwner.setSize(mesh_.nFaces());
439  label outsideI = 0;
440 
441  // Get faces on interface between meshType and non-meshType
442 
443  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
444  {
445  label ownType = operator[](own[facei]);
446  label nbrType = operator[](nbr[facei]);
447 
448  if (ownType == meshType && nbrType != meshType)
449  {
450  outsideFaces[outsideI] = faces[facei];
451  outsideOwner[outsideI] = own[facei]; // meshType cell
452  outsideI++;
453  }
454  else if (ownType != meshType && nbrType == meshType)
455  {
456  outsideFaces[outsideI] = faces[facei];
457  outsideOwner[outsideI] = nbr[facei]; // meshType cell
458  outsideI++;
459  }
460  }
461 
462  // Get faces on outside of real mesh with cells of meshType.
463 
464  for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
465  {
466  if (operator[](own[facei]) == meshType)
467  {
468  outsideFaces[outsideI] = faces[facei];
469  outsideOwner[outsideI] = own[facei]; // meshType cell
470  outsideI++;
471  }
472  }
473  outsideFaces.setSize(outsideI);
474  outsideOwner.setSize(outsideI);
475 }
476 
477 
478 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
479 
480 // Construct from mesh and surface and point(s) on outside
482 (
483  const polyMesh& mesh,
484  const meshSearch& meshQuery,
485  const triSurfaceSearch& surfQuery,
486  const pointField& outsidePoints
487 )
488 :
490  mesh_(mesh)
491 {
492  markCells
493  (
494  meshQuery,
495  markFaces(surfQuery),
496  outsidePoints
497  );
498 }
499 
500 
501 // Construct from mesh and meshType.
503 (
504  const polyMesh& mesh,
505  const labelList& cellType
506 )
507 :
508  labelList(cellType),
509  mesh_(mesh)
510 {
511  if (mesh_.nCells() != size())
512  {
514  << "Number of elements of cellType argument is not equal to the"
515  << " number of cells" << abort(FatalError);
516  }
517 }
518 
519 
520 // Copy constructor
522 :
523  labelList(cType),
524  mesh_(cType.mesh())
525 {}
526 
527 
528 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
529 
530 
531 // Makes cutCells further than nLayers away from meshType to
532 // fillType. Returns number of cells changed.
534 (
535  const label nLayers,
536  const label meshType,
537  const label fillType
538 )
539 {
540  // Temporary cell type for growing.
541  labelList newCellType(*this);
542 
543 // // Split types into outside and rest
544 // forAll(*this, celli)
545 // {
546 // label type = operator[](celli);
547 //
548 // if (type == meshType)
549 // {
550 // newCellType[celli] = type;
551 // }
552 // else
553 // {
554 // newCellType[celli] = fillType;
555 // }
556 // }
557 
558  newCellType = *this;
559 
560 
561  // Do point-cell-point walk on newCellType out from meshType cells
562  for (label iter = 0; iter < nLayers; iter++)
563  {
564  // Get status of points: visible from meshType (pointSide == MESH)
565  // or fillType/cutCells cells (pointSide == NONMESH) or
566  // both (pointSide == MIXED)
567  List<pointStatus> pointSide(mesh_.nPoints());
568  classifyPoints(meshType, newCellType, pointSide);
569 
570  // Grow layer of outside cells
571  forAll(pointSide, pointi)
572  {
573  if (pointSide[pointi] == MIXED)
574  {
575  // Make cut
576  const labelList& pCells = mesh_.pointCells()[pointi];
577 
578  forAll(pCells, i)
579  {
580  label type = newCellType[pCells[i]];
581 
582  if (type == cellClassification::CUT)
583  {
584  // Found cut cell sharing point with
585  // mesh type cell. Grow.
586  newCellType[pCells[i]] = meshType;
587  }
588  }
589  }
590  }
591  }
592 
593  // Merge newCellType into *this:
594  // - leave meshType cells intact
595  // - leave nonMesh cells intact
596  // - make cutcells fillType if they were not marked in newCellType
597 
598  label nChanged = 0;
599 
600  forAll(newCellType, celli)
601  {
602  if (operator[](celli) == cellClassification::CUT)
603  {
604  if (newCellType[celli] != meshType)
605  {
606  // Cell was cutCell but further than nLayers away from
607  // meshType. Convert to fillType.
608  operator[](celli) = fillType;
609  nChanged++;
610  }
611  }
612  }
613 
614  return nChanged;
615 }
616 
617 
618 // Grow surface by pointNeighbours of type meshType
620 (
621  const label meshType,
622  const label fillType
623 )
624 {
625  boolList hasMeshType(mesh_.nPoints(), false);
626 
627  // Mark points used by meshType cells
628 
629  forAll(mesh_.pointCells(), pointi)
630  {
631  const labelList& myCells = mesh_.pointCells()[pointi];
632 
633  // Check if one of cells has meshType
634  forAll(myCells, myCelli)
635  {
636  label type = operator[](myCells[myCelli]);
637 
638  if (type == meshType)
639  {
640  hasMeshType[pointi] = true;
641 
642  break;
643  }
644  }
645  }
646 
647  // Change neighbours of marked points
648 
649  label nChanged = 0;
650 
651  forAll(hasMeshType, pointi)
652  {
653  if (hasMeshType[pointi])
654  {
655  const labelList& myCells = mesh_.pointCells()[pointi];
656 
657  forAll(myCells, myCelli)
658  {
659  if (operator[](myCells[myCelli]) != meshType)
660  {
661  operator[](myCells[myCelli]) = fillType;
662 
663  nChanged++;
664  }
665  }
666  }
667  }
668  return nChanged;
669 }
670 
671 
672 // Check all points used by cells of meshType for use of at least one point
673 // which is not on the outside. If all points are on the outside of the mesh
674 // this would probably result in a flattened cell when projecting the cell
675 // onto the surface.
677 (
678  const label meshType,
679  const label fillType,
680  const label maxIter
681 )
682 {
683  label nTotChanged = 0;
684 
685  for (label iter = 0; iter < maxIter; iter++)
686  {
687  label nChanged = 0;
688 
689  // Get status of points: visible from meshType or non-meshType cells
690  List<pointStatus> pointSide(mesh_.nPoints());
691  classifyPoints(meshType, *this, pointSide);
692 
693  // Check all cells using mixed point type for whether they use mixed
694  // points only. Note: could probably speed this up by counting number
695  // of mixed verts per face and mixed faces per cell or something?
696  forAll(pointSide, pointi)
697  {
698  if (pointSide[pointi] == MIXED)
699  {
700  const labelList& pCells = mesh_.pointCells()[pointi];
701 
702  forAll(pCells, i)
703  {
704  label celli = pCells[i];
705 
706  if (operator[](celli) == meshType)
707  {
708  if (usesMixedPointsOnly(pointSide, celli))
709  {
710  operator[](celli) = fillType;
711 
712  nChanged++;
713  }
714  }
715  }
716  }
717  }
718  nTotChanged += nChanged;
719 
720  Pout<< "removeHangingCells : changed " << nChanged
721  << " hanging cells" << endl;
722 
723  if (nChanged == 0)
724  {
725  break;
726  }
727  }
728 
729  return nTotChanged;
730 }
731 
732 
734 (
735  const label meshType,
736  const label fillType,
737  const label maxIter
738 )
739 {
740  label nTotChanged = 0;
741 
742  for (label iter = 0; iter < maxIter; iter++)
743  {
744  // Get interface between meshType cells and non-meshType cells as a list
745  // of faces and for each face the cell which is the meshType.
746  faceList outsideFaces;
747  labelList outsideOwner;
748 
749  getMeshOutside(meshType, outsideFaces, outsideOwner);
750 
751  // Build primitivePatch out of it and check it for problems.
752  primitiveFacePatch fp(outsideFaces, mesh_.points());
753 
754  const labelListList& edgeFaces = fp.edgeFaces();
755 
756  label nChanged = 0;
757 
758  // Check all edgeFaces for non-manifoldness
759 
760  forAll(edgeFaces, edgeI)
761  {
762  const labelList& eFaces = edgeFaces[edgeI];
763 
764  if (eFaces.size() > 2)
765  {
766  // patch connected through pinched edge. Remove first face using
767  // edge (and not yet changed)
768 
769  forAll(eFaces, i)
770  {
771  label patchFacei = eFaces[i];
772 
773  label ownerCell = outsideOwner[patchFacei];
774 
775  if (operator[](ownerCell) == meshType)
776  {
777  operator[](ownerCell) = fillType;
778 
779  nChanged++;
780 
781  break;
782  }
783  }
784  }
785  }
786 
787  nTotChanged += nChanged;
788 
789  Pout<< "fillRegionEdges : changed " << nChanged
790  << " cells using multiply connected edges" << endl;
791 
792  if (nChanged == 0)
793  {
794  break;
795  }
796  }
797 
798  return nTotChanged;
799 }
800 
801 
803 (
804  const label meshType,
805  const label fillType,
806  const label maxIter
807 )
808 {
809  label nTotChanged = 0;
810 
811  for (label iter = 0; iter < maxIter; iter++)
812  {
813  // Get interface between meshType cells and non-meshType cells as a list
814  // of faces and for each face the cell which is the meshType.
815  faceList outsideFaces;
816  labelList outsideOwner;
817 
818  getMeshOutside(meshType, outsideFaces, outsideOwner);
819 
820  // Build primitivePatch out of it and check it for problems.
821  primitiveFacePatch fp(outsideFaces, mesh_.points());
822 
823  labelHashSet nonManifoldPoints;
824 
825  // Check for non-manifold points.
826  fp.checkPointManifold(false, &nonManifoldPoints);
827 
828  const Map<label>& meshPointMap = fp.meshPointMap();
829 
830  label nChanged = 0;
831 
832  forAllConstIter(labelHashSet, nonManifoldPoints, iter)
833  {
834  // Find a face on fp using point and remove it.
835  const label patchPointi = meshPointMap[iter.key()];
836 
837  const labelList& pFaces = fp.pointFaces()[patchPointi];
838 
839  // Remove any face using conflicting point. Does first face which
840  // has not yet been done. Could be more intelligent and decide which
841  // one would be best to remove.
842  forAll(pFaces, i)
843  {
844  const label patchFacei = pFaces[i];
845  const label ownerCell = outsideOwner[patchFacei];
846 
847  if (operator[](ownerCell) == meshType)
848  {
849  operator[](ownerCell) = fillType;
850 
851  nChanged++;
852  break;
853  }
854  }
855  }
856 
857  nTotChanged += nChanged;
858 
859  Pout<< "fillRegionPoints : changed " << nChanged
860  << " cells using multiply connected points" << endl;
861 
862  if (nChanged == 0)
863  {
864  break;
865  }
866  }
867 
868  return nTotChanged;
869 }
870 
871 
873 {
874  os << "Cells:" << size() << endl
875  << " notset : " << count(*this, NOTSET) << endl
876  << " cut : " << count(*this, CUT) << endl
877  << " inside : " << count(*this, INSIDE) << endl
878  << " outside : " << count(*this, OUTSIDE) << endl;
879 }
880 
881 
882 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
883 
885 {
887 }
888 
889 
890 // ************************************************************************* //
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
Definition: meshSearch.H:57
label trimCutCells(const label nLayers, const label meshType, const label fillType)
#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 const Form max
Definition: VectorSpace.H:115
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
T & operator[](const label)
Return element of UList.
Definition: UListI.H:167
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
List< face > faceList
Definition: faceListFwd.H:43
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void operator=(const cellClassification &)
static const Form min
Definition: VectorSpace.H:116
Various functions to operate on Lists.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:208
cellClassification(const polyMesh &mesh, const meshSearch &meshQuery, const triSurfaceSearch &surfQuery, const pointField &outsidePoints)
Construct from mesh and surface and point(s) on outside.
Helper class to search on triSurface.
A list of faces which address into the list of points.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
&#39;Cuts&#39; a mesh with a surface.
List< edge > edgeList
Definition: edgeList.H:38
label fillRegionPoints(const label meshType, const label fillType, const label maxIter)
Find regionPoints and fill all neighbours. Iterate until nothing.
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
const labelListList & pointCells() const
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const char nl
Definition: Ostream.H:265
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
void operator=(const UList< label > &)
Assignment to UList operator. Takes linear time.
Definition: List.C:376
void setSize(const label)
Reset size of List.
Definition: List.C:281
vector point
Point is a vector.
Definition: point.H:41
label growSurface(const label meshType, const label fillType)
Sets vertex neighbours of meshType cells to fillType.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
label fillRegionEdges(const label meshType, const label fillType, const label maxIter)
Find regionEdges and fill one neighbour. Iterate until nothing.
void writeStats(Ostream &os) const
Write statistics on cell types to Ostream.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
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 fillHangingCells(const label meshType, const label fillType, const label maxIter)
Find hanging cells (cells with all points on outside) and set their.
label size() const
Return the number of elements in the UList.
Definition: ListI.H:170
Namespace for OpenFOAM.