ENH: autoLayerDriver: better layering information message
[OpenFOAM-2.0.x.git] / src / meshTools / surfaceSets / surfaceSets.C
blob47f553517a4bdc72e6182445a922ae616ac11aeb
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
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
19     for more details.
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"
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"
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
45 //Foam::scalar Foam::surfaceSets::minEdgeLen
46 //(
47 //    const primitiveMesh& mesh,
48 //    const label pointI
49 //)
50 //{
51 //    const edgeList& edges = mesh.edges();
53 //    const pointField& points = mesh.points();
55 //    const labelList& pEdges = mesh.pointEdges()[pointI];
57 //    scalar minLen = GREAT;
59 //    forAll(pEdges, i)
60 //    {
61 //        minLen = min(minLen, edges[pEdges[i]].mag(points));
62 //    }
63 //    return minLen;
64 //}
67 //// Returns true if cell uses at least one selected point
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];
77 //    forAll(cFaces, cFaceI)
78 //    {
79 //        label faceI = cFaces[cFaceI];
81 //        const face& f = mesh.faces()[faceI];
83 //        forAll(f, fp)
84 //        {
85 //            if (selectedPoint[f[fp]])
86 //            {
87 //                return true;
88 //            }
89 //        }
90 //    }
91 //    return false;
92 //}
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
114 //    // snapped
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)
125 //    {
126 //        const cell& cFaces = cells[iter.key()];
128 //        forAll(cFaces, cFaceI)
129 //        {
130 //            const face& f = faces[cFaces[cFaceI]];
132 //            forAll(f, fp)
133 //            {
134 //                label pointI = f[fp];
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 //    }
148 //    // Calculate new volume for mixed cells
149 //    label nRemoved = 0;
150 //    forAllConstIter(labelHashSet, flatCandidates, iter)
151 //    {
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)
160 //        {
161 //            internalCells.erase(cellI);
162 //            nRemoved++;
163 //        }
164 //    }
166 //    return nRemoved;
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)
184 //    {
185 //        return;
186 //    }
188 //    labelList candidates(candidateSet.toc());
190 //    pointField candidatePoints(candidates.size());
191 //    forAll(candidates, i)
192 //    {
193 //        candidatePoints[i] = mesh.points()[candidates[i]];
194 //    }
196 //    tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
197 //    const pointField& nearest = tnearest();
199 //    const pointField& points = mesh.points();
201 //    forAll(candidates, i)
202 //    {
203 //        label pointI = candidates[i];
205 //        scalar minLen = minEdgeLen(mesh, pointI);
207 //        scalar dist = mag(nearest[i] - points[pointI]);
209 //        if (dist < edgeFactor * minLen)
210 //        {
211 //            nearPointSet.insert(pointI);
212 //        }
213 //    }
217 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
220 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
222 void Foam::surfaceSets::getSurfaceSets
224     const polyMesh& mesh,
225     const fileName&,
226     const triSurface&,
227     const triSurfaceSearch& querySurf,
228     const pointField& outsidePts,
230     const label nCutLayers,
232     labelHashSet& inside,
233     labelHashSet& outside,
234     labelHashSet& cut
237     // Construct search engine on mesh
238     meshSearch queryMesh(mesh, true);
240     // Cut faces with surface and classify cells
241     cellClassification cellType
242     (
243         mesh,
244         queryMesh,
245         querySurf,
246         outsidePts
247     );
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     }
261     forAll(cellType, cellI)
262     {
263         label cType = cellType[cellI];
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     }
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
294     // -mixed
296     List<pointStatus> pointSide(mesh.nPoints(), NOTSET);
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];
305             forAll(cFaces, cFaceI)
306             {
307                 const face& f = faces[cFaces[cFaceI]];
309                 forAll(f, fp)
310                 {
311                     label pointI = f[fp];
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];
333             forAll(cFaces, cFaceI)
334             {
335                 const face& f = faces[cFaces[cFaceI]];
337                 forAll(f, fp)
338                 {
339                     label pointI = f[fp];
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     }
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     //}
373     // Determine cells using mixed points only
375     labelHashSet mixedOnlyCells(internalCells.size());
377     forAllConstIter(labelHashSet, internalCells, iter)
378     {
379         const label cellI = iter.key();
380         const cell& cFaces = cells[cellI];
382         label usesMixedOnly = true;
384         forAll(cFaces, i)
385         {
386             const face& f = faces[cFaces[i]];
388             forAll(f, fp)
389             {
390                 if (pointSide[f[fp]] != MIXED)
391                 {
392                     usesMixedOnly = false;
393                     break;
394                 }
395             }
397             if (!usesMixedOnly)
398             {
399                 break;
400             }
401         }
402         if (usesMixedOnly)
403         {
404             mixedOnlyCells.insert(cellI);
405         }
406     }
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
428 //    getSurfaceSets
429 //    (
430 //        mesh,
431 //        surfName,
432 //        surf,
433 //        querySurf,
434 //        outsidePts,
436 //        rawInside,
437 //        rawOutside,
438 //        cutCells
439 //    );
442 //    Pout<< "rawInside:" << rawInside.size() << endl;
444 //    label nRemoved;
445 //    do
446 //    {
447 //        nRemoved = removeHangingCells(mesh, querySurf, rawInside);
449 //        Pout<< nl
450 //            << "Removed " << nRemoved
451 //            << " rawInside cells that have all their points on the outside"
452 //            << endl;
453 //    }
454 //    while (nRemoved != 0);
456 //    Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
457 //        << rawInside.instance()/rawInside.local()/rawInside.name()
458 //        << endl << endl;
459 //    rawInside.write();
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 //    );
478 //    outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
480 //    Pout<< "rawOutside:" << rawOutside.size() << endl;
482 //    do
483 //    {
484 //        nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
486 //        Pout<< nl
487 //            << "Removed " << nRemoved
488 //            << " rawOutside cells that have all their points on the outside"
489 //            << endl;
490 //    }
491 //    while (nRemoved != 0);
493 //    Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
494 //        << rawOutside.instance()/rawOutside.local()/rawOutside.name()
495 //        << endl << endl;
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()
507 //        << endl << endl;
508 //    cutCells.write();
511 //    //
512 //    // Remove cells with points too close to surface.
513 //    //
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());
524 //    getNearPoints
525 //    (
526 //        mesh,
527 //        surf,
528 //        querySurf,
529 //        edgeFactor,
530 //        cutPoints,
531 //        nearPoints
532 //    );
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;
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);
555 //    Pout<< nl
556 //        << "Removed " << rawInside.size() - inside.size()
557 //        << " inside cells that are too close to the surface" << endl;
559 //    Pout<< nl
560 //        << "Removed " << rawOutside.size() - outside.size()
561 //        << " inside cells that are too close to the surface" << nl << endl;
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 //    //
572 //    do
573 //    {
574 //        nRemoved = removeHangingCells(mesh, querySurf, inside);
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);
586 //        Pout<< nl
587 //            << "Removed " << nRemoved
588 //            << " outside cells that have all their points on the inside"
589 //            << endl;
590 //    }
591 //    while (nRemoved != 0);
594 //    //
595 //    // Write
596 //    //
599 //    Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
600 //        << inside.instance()/inside.local()/inside.name()
601 //        << endl << endl;
603 //    inside.write();
605 //    Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
606 //        << outside.instance()/outside.local()/outside.name()
607 //        << endl << endl;
609 //    outside.write();
613 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
616 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
619 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
622 // ************************************************************************* //