1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "surfaceSets.H"
28 #include "triSurface.H"
29 #include "triSurfaceSearch.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"
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 //Foam::scalar Foam::surfaceSets::minEdgeLen
47 // const primitiveMesh& mesh,
51 // const edgeList& edges = mesh.edges();
53 // const pointField& points = mesh.points();
55 // const labelList& pEdges = mesh.pointEdges()[pointI];
57 // scalar minLen = GREAT;
61 // minLen = min(minLen, edges[pEdges[i]].mag(points));
67 //// Returns true if cell uses at least one selected point
68 //bool Foam::surfaceSets::usesPoint
70 // const primitiveMesh& mesh,
71 // const boolList& selectedPoint,
75 // const labelList& cFaces = mesh.cells()[cellI];
77 // forAll(cFaces, cFaceI)
79 // label faceI = cFaces[cFaceI];
81 // const face& f = mesh.faces()[faceI];
85 // if (selectedPoint[f[fp]])
96 //// Remove cells in allCells which are connected to other cells in allCells
97 //// by outside vertices only. Since these outside vertices will be moved onto
98 //// a surface they might result in flat cells.
99 //Foam::label Foam::surfaceSets::removeHangingCells
101 // const primitiveMesh& mesh,
102 // const triSurfaceSearch& querySurf,
103 // labelHashSet& internalCells
106 // const pointField& points = mesh.points();
107 // const cellList& cells = mesh.cells();
108 // const faceList& faces = mesh.faces();
110 // // Determine cells that have all points on the boundary.
111 // labelHashSet flatCandidates(getHangingCells(mesh, internalCells));
113 // // All boundary points will become visible after subsetting and will be
115 // // to surface. Calculate new volume for cells using only these points and
116 // // check if it does not become too small.
118 // // Get points used by flatCandidates.
119 // labelHashSet outsidePoints(flatCandidates.size());
121 // // Snap outside points to surface
122 // pointField newPoints(points);
124 // forAllConstIter(labelHashSet, flatCandidates, iter)
126 // const cell& cFaces = cells[iter.key()];
128 // forAll(cFaces, cFaceI)
130 // const face& f = faces[cFaces[cFaceI]];
134 // label pointI = f[fp];
136 // if (outsidePoints.insert(pointI))
138 // // Calculate new position for this outside point
139 // tmp<pointField> tnearest =
140 // querySurf.calcNearest(pointField(1, points[pointI]));
141 // newPoints[pointI] = tnearest()[0];
148 // // Calculate new volume for mixed cells
149 // label nRemoved = 0;
150 // forAllConstIter(labelHashSet, flatCandidates, iter)
152 // label cellI = iter.key();
154 // const cell& cll = cells[cellI];
156 // scalar newVol = cll.mag(newPoints, faces);
157 // scalar oldVol = mesh.cellVolumes()[cellI];
159 // if (newVol < 0.1 * oldVol)
161 // internalCells.erase(cellI);
170 //// Select all points out of pointSet where the distance to the surface
171 //// is less than a factor times a local length scale (minimum length of
172 //// connected edges)
173 //void Foam::surfaceSets::getNearPoints
175 // const primitiveMesh& mesh,
176 // const triSurface&,
177 // const triSurfaceSearch& querySurf,
178 // const scalar edgeFactor,
179 // const pointSet& candidateSet,
180 // pointSet& nearPointSet
183 // if (edgeFactor <= 0)
188 // labelList candidates(candidateSet.toc());
190 // pointField candidatePoints(candidates.size());
191 // forAll(candidates, i)
193 // candidatePoints[i] = mesh.points()[candidates[i]];
196 // tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
197 // const pointField& nearest = tnearest();
199 // const pointField& points = mesh.points();
201 // forAll(candidates, i)
203 // label pointI = candidates[i];
205 // scalar minLen = minEdgeLen(mesh, pointI);
207 // scalar dist = mag(nearest[i] - points[pointI]);
209 // if (dist < edgeFactor * minLen)
211 // nearPointSet.insert(pointI);
217 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
220 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
222 void Foam::surfaceSets::getSurfaceSets
224 const polyMesh& mesh,
227 const triSurfaceSearch& querySurf,
228 const pointField& outsidePts,
230 const label nCutLayers,
232 labelHashSet& inside,
233 labelHashSet& outside,
237 // Construct search engine on mesh
238 meshSearch queryMesh(mesh, true);
240 // Cut faces with surface and classify cells
241 cellClassification cellType
251 // Trim cutCells so they are max nCutLayers away (as seen in point-cell
252 // walk) from outside cells.
253 cellType.trimCutCells
256 cellClassification::OUTSIDE,
257 cellClassification::INSIDE
261 forAll(cellType, cellI)
263 label cType = cellType[cellI];
265 if (cType == cellClassification::CUT)
269 else if (cType == cellClassification::INSIDE)
271 inside.insert(cellI);
273 else if (cType == cellClassification::OUTSIDE)
275 outside.insert(cellI);
281 Foam::labelHashSet Foam::surfaceSets::getHangingCells
283 const primitiveMesh& mesh,
284 const labelHashSet& internalCells
287 const cellList& cells = mesh.cells();
288 const faceList& faces = mesh.faces();
291 // Divide points into
292 // -referenced by internal only
293 // -referenced by outside only
296 List<pointStatus> pointSide(mesh.nPoints(), NOTSET);
298 for (label cellI = 0; cellI < mesh.nCells(); cellI++)
300 if (internalCells.found(cellI))
302 // Inside cell. Mark all vertices seen from this cell.
303 const labelList& cFaces = cells[cellI];
305 forAll(cFaces, cFaceI)
307 const face& f = faces[cFaces[cFaceI]];
311 label pointI = f[fp];
313 if (pointSide[pointI] == NOTSET)
315 pointSide[pointI] = INSIDE;
317 else if (pointSide[pointI] == OUTSIDE)
319 pointSide[pointI] = MIXED;
323 // mixed or inside stay same
331 const labelList& cFaces = cells[cellI];
333 forAll(cFaces, cFaceI)
335 const face& f = faces[cFaces[cFaceI]];
339 label pointI = f[fp];
341 if (pointSide[pointI] == NOTSET)
343 pointSide[pointI] = OUTSIDE;
345 else if (pointSide[pointI] == INSIDE)
347 pointSide[pointI] = MIXED;
351 // mixed or outside stay same
359 //OFstream mixedStr("mixed.obj");
361 //forAll(pointSide, pointI)
363 // if (pointSide[pointI] == MIXED)
365 // const point& pt = points[pointI];
367 // mixedStr << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
373 // Determine cells using mixed points only
375 labelHashSet mixedOnlyCells(internalCells.size());
377 forAllConstIter(labelHashSet, internalCells, iter)
379 const label cellI = iter.key();
380 const cell& cFaces = cells[cellI];
382 label usesMixedOnly = true;
386 const face& f = faces[cFaces[i]];
390 if (pointSide[f[fp]] != MIXED)
392 usesMixedOnly = false;
404 mixedOnlyCells.insert(cellI);
408 return mixedOnlyCells;
412 //void Foam::surfaceSets::writeSurfaceSets
414 // const polyMesh& mesh,
415 // const fileName& surfName,
416 // const triSurface& surf,
417 // const triSurfaceSearch& querySurf,
418 // const pointField& outsidePts,
419 // const scalar edgeFactor
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);
427 // // Get inside/outside/cut cells
442 // Pout<< "rawInside:" << rawInside.size() << endl;
447 // nRemoved = removeHangingCells(mesh, querySurf, rawInside);
450 // << "Removed " << nRemoved
451 // << " rawInside cells that have all their points on the outside"
454 // while (nRemoved != 0);
456 // Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
457 // << rawInside.instance()/rawInside.local()/rawInside.name()
459 // rawInside.write();
463 // // Select outside cells
464 // surfaceToCell outsideSource
471 // false, // includeCut
472 // false, // includeInside
473 // true, // includeOutside
474 // -GREAT, // nearDist
475 // -GREAT // curvature
478 // outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
480 // Pout<< "rawOutside:" << rawOutside.size() << endl;
484 // nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
487 // << "Removed " << nRemoved
488 // << " rawOutside cells that have all their points on the outside"
491 // while (nRemoved != 0);
493 // Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
494 // << rawOutside.instance()/rawOutside.local()/rawOutside.name()
496 // rawOutside.write();
499 // // Select cut cells by negating inside and outside set.
500 // cutCells.invert(mesh.nCells());
502 // cellToCell deleteInsideSource(mesh, rawInside.name());
504 // deleteInsideSource.applyToSet(topoSetSource::DELETE, cutCells);
505 // Pout<< "Writing cut cells (" << cutCells.size() << ") to cellSet "
506 // << cutCells.instance()/cutCells.local()/cutCells.name()
512 // // Remove cells with points too close to surface.
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);
521 // // Get all points that are too close to surface.
522 // pointSet nearPoints(mesh, "nearPoints", cutPoints.size());
535 // << "Selected " << nearPoints.size()
536 // << " points that are closer than " << edgeFactor
537 // << " times the local minimum lengthscale to the surface"
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);
548 // // Start off from copy of rawInside, rawOutside
549 // cellSet inside(mesh, "inside", rawInside);
550 // cellSet outside(mesh, "outside", rawOutside);
552 // pToCell.applyToSet(topoSetSource::DELETE, inside);
553 // pToCell.applyToSet(topoSetSource::DELETE, outside);
556 // << "Removed " << rawInside.size() - inside.size()
557 // << " inside cells that are too close to the surface" << endl;
560 // << "Removed " << rawOutside.size() - outside.size()
561 // << " inside cells that are too close to the surface" << nl << endl;
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.
574 // nRemoved = removeHangingCells(mesh, querySurf, inside);
577 // << "Removed " << nRemoved
578 // << " inside cells that have all their points on the outside"
581 // while (nRemoved != 0);
584 // nRemoved = removeHangingCells(mesh, querySurf, outside);
587 // << "Removed " << nRemoved
588 // << " outside cells that have all their points on the inside"
591 // while (nRemoved != 0);
599 // Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
600 // << inside.instance()/inside.local()/inside.name()
605 // Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
606 // << outside.instance()/outside.local()/outside.name()
613 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
616 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
619 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
622 // ************************************************************************* //