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