surfaceSets.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-2018 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 "surfaceSets.H"
27 #include "polyMesh.H"
28 #include "triSurface.H"
29 #include "triSurfaceSearch.H"
30 #include "pointSet.H"
31 #include "cellSet.H"
32 #include "surfaceToCell.H"
33 #include "cellToPoint.H"
34 #include "cellToCell.H"
35 #include "pointToCell.H"
36 #include "meshSearch.H"
37 #include "cellClassification.H"
38 
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 //Foam::scalar Foam::surfaceSets::minEdgeLen
46 //(
47 // const primitiveMesh& mesh,
48 // const label pointi
49 //)
50 //{
51 // const edgeList& edges = mesh.edges();
52 //
53 // const pointField& points = mesh.points();
54 //
55 // const labelList& pEdges = mesh.pointEdges()[pointi];
56 //
57 // scalar minLen = great;
58 //
59 // forAll(pEdges, i)
60 // {
61 // minLen = min(minLen, edges[pEdges[i]].mag(points));
62 // }
63 // return minLen;
64 //}
65 //
66 //
68 //bool Foam::surfaceSets::usesPoint
69 //(
70 // const primitiveMesh& mesh,
71 // const boolList& selectedPoint,
72 // const label celli
73 //)
74 //{
75 // const labelList& cFaces = mesh.cells()[celli];
76 //
77 // forAll(cFaces, cFacei)
78 // {
79 // label facei = cFaces[cFacei];
80 //
81 // const face& f = mesh.faces()[facei];
82 //
83 // forAll(f, fp)
84 // {
85 // if (selectedPoint[f[fp]])
86 // {
87 // return true;
88 // }
89 // }
90 // }
91 // return false;
92 //}
93 
94 
95 
99 //Foam::label Foam::surfaceSets::removeHangingCells
100 //(
101 // const primitiveMesh& mesh,
102 // const triSurfaceSearch& querySurf,
103 // labelHashSet& internalCells
104 //)
105 //{
106 // const pointField& points = mesh.points();
107 // const cellList& cells = mesh.cells();
108 // const faceList& faces = mesh.faces();
109 //
110 // // Determine cells that have all points on the boundary.
111 // labelHashSet flatCandidates(getHangingCells(mesh, internalCells));
112 //
113 // // All boundary points will become visible after subsetting and will be
114 // // snapped
115 // // to surface. Calculate new volume for cells using only these points and
116 // // check if it does not become too small.
117 //
118 // // Get points used by flatCandidates.
119 // labelHashSet outsidePoints(flatCandidates.size());
120 //
121 // // Snap outside points to surface
122 // pointField newPoints(points);
123 //
124 // forAllConstIter(labelHashSet, flatCandidates, iter)
125 // {
126 // const cell& cFaces = cells[iter.key()];
127 //
128 // forAll(cFaces, cFacei)
129 // {
130 // const face& f = faces[cFaces[cFacei]];
131 //
132 // forAll(f, fp)
133 // {
134 // label pointi = f[fp];
135 //
136 // if (outsidePoints.insert(pointi))
137 // {
138 // // Calculate new position for this outside point
139 // tmp<pointField> tnearest =
140 // querySurf.calcNearest(pointField(1, points[pointi]));
141 // newPoints[pointi] = tnearest()[0];
142 // }
143 // }
144 // }
145 // }
146 //
147 //
148 // // Calculate new volume for mixed cells
149 // label nRemoved = 0;
150 // forAllConstIter(labelHashSet, flatCandidates, iter)
151 // {
152 // label celli = iter.key();
153 //
154 // const cell& cll = cells[celli];
155 //
156 // scalar newVol = cll.mag(newPoints, faces);
157 // scalar oldVol = mesh.cellVolumes()[celli];
158 //
159 // if (newVol < 0.1 * oldVol)
160 // {
161 // internalCells.erase(celli);
162 // nRemoved++;
163 // }
164 // }
165 //
166 // return nRemoved;
167 //}
168 
169 
173 //void Foam::surfaceSets::getNearPoints
174 //(
175 // const primitiveMesh& mesh,
176 // const triSurface&,
177 // const triSurfaceSearch& querySurf,
178 // const scalar edgeFactor,
179 // const pointSet& candidateSet,
180 // pointSet& nearPointSet
181 //)
182 //{
183 // if (edgeFactor <= 0)
184 // {
185 // return;
186 // }
187 //
188 // labelList candidates(candidateSet.toc());
189 //
190 // pointField candidatePoints(candidates.size());
191 // forAll(candidates, i)
192 // {
193 // candidatePoints[i] = mesh.points()[candidates[i]];
194 // }
195 //
196 // tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
197 // const pointField& nearest = tnearest();
198 //
199 // const pointField& points = mesh.points();
200 //
201 // forAll(candidates, i)
202 // {
203 // label pointi = candidates[i];
204 //
205 // scalar minLen = minEdgeLen(mesh, pointi);
206 //
207 // scalar dist = mag(nearest[i] - points[pointi]);
208 //
209 // if (dist < edgeFactor * minLen)
210 // {
211 // nearPointSet.insert(pointi);
212 // }
213 // }
214 //}
215 
216 
217 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
218 
219 
220 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
221 
223 (
224  const polyMesh& mesh,
225  const fileName&,
226  const triSurface&,
227  const triSurfaceSearch& querySurf,
228  const pointField& outsidePts,
229 
230  const label nCutLayers,
231 
232  labelHashSet& inside,
233  labelHashSet& outside,
234  labelHashSet& cut
235 )
236 {
237  // Construct search engine on mesh
238  meshSearch queryMesh(mesh);
239 
240  // Cut faces with surface and classify cells
241  cellClassification cellType
242  (
243  mesh,
244  queryMesh,
245  querySurf,
246  outsidePts
247  );
248 
249  if (nCutLayers > 0)
250  {
251  // Trim cutCells so they are max nCutLayers away (as seen in point-cell
252  // walk) from outside cells.
253  cellType.trimCutCells
254  (
255  nCutLayers,
256  cellClassification::OUTSIDE,
257  cellClassification::INSIDE
258  );
259  }
260 
261  forAll(cellType, celli)
262  {
263  label cType = cellType[celli];
264 
265  if (cType == cellClassification::CUT)
266  {
267  cut.insert(celli);
268  }
269  else if (cType == cellClassification::INSIDE)
270  {
271  inside.insert(celli);
272  }
273  else if (cType == cellClassification::OUTSIDE)
274  {
275  outside.insert(celli);
276  }
277  }
278 }
279 
280 
282 (
283  const primitiveMesh& mesh,
284  const labelHashSet& internalCells
285 )
286 {
287  const cellList& cells = mesh.cells();
288  const faceList& faces = mesh.faces();
289 
290 
291  // Divide points into
292  // -referenced by internal only
293  // -referenced by outside only
294  // -mixed
295 
296  List<pointStatus> pointSide(mesh.nPoints(), NOTSET);
297 
298  for (label celli = 0; celli < mesh.nCells(); celli++)
299  {
300  if (internalCells.found(celli))
301  {
302  // Inside cell. Mark all vertices seen from this cell.
303  const labelList& cFaces = cells[celli];
304 
305  forAll(cFaces, cFacei)
306  {
307  const face& f = faces[cFaces[cFacei]];
308 
309  forAll(f, fp)
310  {
311  label pointi = f[fp];
312 
313  if (pointSide[pointi] == NOTSET)
314  {
315  pointSide[pointi] = INSIDE;
316  }
317  else if (pointSide[pointi] == OUTSIDE)
318  {
319  pointSide[pointi] = MIXED;
320  }
321  else
322  {
323  // mixed or inside stay same
324  }
325  }
326  }
327  }
328  else
329  {
330  // Outside cell
331  const labelList& cFaces = cells[celli];
332 
333  forAll(cFaces, cFacei)
334  {
335  const face& f = faces[cFaces[cFacei]];
336 
337  forAll(f, fp)
338  {
339  label pointi = f[fp];
340 
341  if (pointSide[pointi] == NOTSET)
342  {
343  pointSide[pointi] = OUTSIDE;
344  }
345  else if (pointSide[pointi] == INSIDE)
346  {
347  pointSide[pointi] = MIXED;
348  }
349  else
350  {
351  // mixed or outside stay same
352  }
353  }
354  }
355  }
356  }
357 
358 
359  // OFstream mixedStr("mixed.obj");
360  //
361  // forAll(pointSide, pointi)
362  //{
363  // if (pointSide[pointi] == MIXED)
364  // {
365  // const point& pt = points[pointi];
366  //
367  // mixedStr << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
368  // << endl;
369  // }
370  //}
371 
372 
373  // Determine cells using mixed points only
374 
375  labelHashSet mixedOnlyCells(internalCells.size());
376 
377  forAllConstIter(labelHashSet, internalCells, iter)
378  {
379  const label celli = iter.key();
380  const cell& cFaces = cells[celli];
381 
382  label usesMixedOnly = true;
383 
384  forAll(cFaces, i)
385  {
386  const face& f = faces[cFaces[i]];
387 
388  forAll(f, fp)
389  {
390  if (pointSide[f[fp]] != MIXED)
391  {
392  usesMixedOnly = false;
393  break;
394  }
395  }
396 
397  if (!usesMixedOnly)
398  {
399  break;
400  }
401  }
402  if (usesMixedOnly)
403  {
404  mixedOnlyCells.insert(celli);
405  }
406  }
407 
408  return mixedOnlyCells;
409 }
410 
411 
412 //void Foam::surfaceSets::writeSurfaceSets
413 //(
414 // const polyMesh& mesh,
415 // const fileName& surfName,
416 // const triSurface& surf,
417 // const triSurfaceSearch& querySurf,
418 // const pointField& outsidePts,
419 // const scalar edgeFactor
420 //)
421 //{
422 // // Cellsets for inside/outside determination
423 // cellSet rawInside(mesh, "rawInside", mesh.nCells()/10);
424 // cellSet rawOutside(mesh, "rawOutside", mesh.nCells()/10);
425 // cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
426 //
427 // // Get inside/outside/cut cells
428 // getSurfaceSets
429 // (
430 // mesh,
431 // surfName,
432 // surf,
433 // querySurf,
434 // outsidePts,
435 //
436 // rawInside,
437 // rawOutside,
438 // cutCells
439 // );
440 //
441 //
442 // Pout<< "rawInside:" << rawInside.size() << endl;
443 //
444 // label nRemoved;
445 // do
446 // {
447 // nRemoved = removeHangingCells(mesh, querySurf, rawInside);
448 //
449 // Pout<< nl
450 // << "Removed " << nRemoved
451 // << " rawInside cells that have all their points on the outside"
452 // << endl;
453 // }
454 // while (nRemoved != 0);
455 //
456 // Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
457 // << rawInside.instance()/rawInside.local()/rawInside.name()
458 // << endl << endl;
459 // rawInside.write();
460 //
461 //
462 //
463 // // Select outside cells
464 // surfaceToCell outsideSource
465 // (
466 // mesh,
467 // surfName,
468 // surf,
469 // querySurf,
470 // outsidePts,
471 // false, // includeCut
472 // false, // includeInside
473 // true, // includeOutside
474 // -great, // nearDist
475 // -great // curvature
476 // );
477 //
478 // outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
479 //
480 // Pout<< "rawOutside:" << rawOutside.size() << endl;
481 //
482 // do
483 // {
484 // nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
485 //
486 // Pout<< nl
487 // << "Removed " << nRemoved
488 // << " rawOutside cells that have all their points on the outside"
489 // << endl;
490 // }
491 // while (nRemoved != 0);
492 //
493 // Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
494 // << rawOutside.instance()/rawOutside.local()/rawOutside.name()
495 // << endl << endl;
496 // rawOutside.write();
497 //
498 //
499 // // Select cut cells by negating inside and outside set.
500 // cutCells.invert(mesh.nCells());
501 //
502 // cellToCell deleteInsideSource(mesh, rawInside.name());
503 //
504 // deleteInsideSource.applyToSet(topoSetSource::DELETE, cutCells);
505 // Pout<< "Writing cut cells (" << cutCells.size() << ") to cellSet "
506 // << cutCells.instance()/cutCells.local()/cutCells.name()
507 // << endl << endl;
508 // cutCells.write();
509 //
510 //
511 // //
512 // // Remove cells with points too close to surface.
513 // //
514 //
515 //
516 // // Get all points in cutCells.
517 // pointSet cutPoints(mesh, "cutPoints", 4*cutCells.size());
518 // cellToPoint cutSource(mesh, "cutCells", cellToPoint::ALL);
519 // cutSource.applyToSet(topoSetSource::NEW, cutPoints);
520 //
521 // // Get all points that are too close to surface.
522 // pointSet nearPoints(mesh, "nearPoints", cutPoints.size());
523 //
524 // getNearPoints
525 // (
526 // mesh,
527 // surf,
528 // querySurf,
529 // edgeFactor,
530 // cutPoints,
531 // nearPoints
532 // );
533 //
534 // Pout<< nl
535 // << "Selected " << nearPoints.size()
536 // << " points that are closer than " << edgeFactor
537 // << " times the local minimum lengthscale to the surface"
538 // << nl << endl;
539 //
540 //
541 // // Remove cells that use any of the points in nearPoints
542 // // from the inside and outsideCells.
543 // nearPoints.write();
544 // pointToCell pToCell(mesh, nearPoints.name(), pointToCell::ANY);
545 //
546 //
547 //
548 // // Start off from copy of rawInside, rawOutside
549 // cellSet inside(mesh, "inside", rawInside);
550 // cellSet outside(mesh, "outside", rawOutside);
551 //
552 // pToCell.applyToSet(topoSetSource::DELETE, inside);
553 // pToCell.applyToSet(topoSetSource::DELETE, outside);
554 //
555 // Pout<< nl
556 // << "Removed " << rawInside.size() - inside.size()
557 // << " inside cells that are too close to the surface" << endl;
558 //
559 // Pout<< nl
560 // << "Removed " << rawOutside.size() - outside.size()
561 // << " inside cells that are too close to the surface" << nl << endl;
562 //
563 //
564 //
565 // //
566 // // Remove cells with one or no faces on rest of cellSet. Note: Problem is
567 // // not these cells an sich but rather that all of their vertices will be
568 // // outside vertices and thus projected onto surface. Which might (if they
569 // // project onto same surface) result in flattened cells.
570 // //
571 //
572 // do
573 // {
574 // nRemoved = removeHangingCells(mesh, querySurf, inside);
575 //
576 // Pout<< nl
577 // << "Removed " << nRemoved
578 // << " inside cells that have all their points on the outside"
579 // << endl;
580 // }
581 // while (nRemoved != 0);
582 // do
583 // {
584 // nRemoved = removeHangingCells(mesh, querySurf, outside);
585 //
586 // Pout<< nl
587 // << "Removed " << nRemoved
588 // << " outside cells that have all their points on the inside"
589 // << endl;
590 // }
591 // while (nRemoved != 0);
592 //
593 //
594 // //
595 // // Write
596 // //
597 //
598 //
599 // Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
600 // << inside.instance()/inside.local()/inside.name()
601 // << endl << endl;
602 //
603 // inside.write();
604 //
605 // Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
606 // << outside.instance()/outside.local()/outside.name()
607 // << endl << endl;
608 //
609 // outside.write();
610 //}
611 
612 
613 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
614 
615 
616 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
617 
618 
619 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
620 
621 
622 // ************************************************************************* //
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
A class for handling file names.
Definition: fileName.H:79
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
label nCells() const
const cellList & cells() const
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
static labelHashSet getHangingCells(const primitiveMesh &mesh, const labelHashSet &internalCells)
Get cells using points on &#39;outside&#39; only.
Definition: surfaceSets.C:282
Helper class to search on triSurface.
const cellShapeList & cells
&#39;Cuts&#39; a mesh with a surface.
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
static void getSurfaceSets(const polyMesh &mesh, const fileName &surfName, const triSurface &surf, const triSurfaceSearch &querySurf, const pointField &outsidePts, const label nCutLayers, labelHashSet &inside, labelHashSet &outside, labelHashSet &cut)
Divide cells into cut,inside and outside.
Definition: surfaceSets.C:223
labelList f(nPoints)
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
virtual const faceList & faces() const =0
Return faces.
label nPoints() const
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
Triangulated surface description with patch information.
Definition: triSurface.H:66