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