FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

surfaceSets.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025 
00026 \*---------------------------------------------------------------------------*/
00027 
00028 #include "surfaceSets.H"
00029 #include <OpenFOAM/polyMesh.H>
00030 #include <triSurface/triSurface.H>
00031 #include <meshTools/triSurfaceSearch.H>
00032 #include <meshTools/pointSet.H>
00033 #include <meshTools/cellSet.H>
00034 #include <meshTools/surfaceToCell.H>
00035 #include <meshTools/cellToPoint.H>
00036 #include <meshTools/cellToCell.H>
00037 #include <meshTools/pointToCell.H>
00038 #include <meshTools/meshSearch.H>
00039 #include <meshTools/cellClassification.H>
00040 
00041 
00042 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00043 
00044 
00045 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00046 
00047 //Foam::scalar Foam::surfaceSets::minEdgeLen
00048 //(
00049 //    const primitiveMesh& mesh,
00050 //    const label pointI
00051 //)
00052 //{
00053 //    const edgeList& edges = mesh.edges();
00054 //
00055 //    const pointField& points = mesh.points();
00056 //
00057 //    const labelList& pEdges = mesh.pointEdges()[pointI];
00058 //
00059 //    scalar minLen = GREAT;
00060 //
00061 //    forAll(pEdges, i)
00062 //    {
00063 //        minLen = min(minLen, edges[pEdges[i]].mag(points));
00064 //    }
00065 //    return minLen;
00066 //}
00067 //
00068 //
00070 //bool Foam::surfaceSets::usesPoint
00071 //(
00072 //    const primitiveMesh& mesh,
00073 //    const boolList& selectedPoint,
00074 //    const label cellI
00075 //)
00076 //{
00077 //    const labelList& cFaces = mesh.cells()[cellI];
00078 //
00079 //    forAll(cFaces, cFaceI)
00080 //    {
00081 //        label faceI = cFaces[cFaceI];
00082 //
00083 //        const face& f = mesh.faces()[faceI];
00084 //
00085 //        forAll(f, fp)
00086 //        {
00087 //            if (selectedPoint[f[fp]])
00088 //            {
00089 //                return true;
00090 //            }
00091 //        }
00092 //    }
00093 //    return false;
00094 //}
00095 
00096 
00097 
00101 //Foam::label Foam::surfaceSets::removeHangingCells
00102 //(
00103 //    const primitiveMesh& mesh,
00104 //    const triSurfaceSearch& querySurf,
00105 //    labelHashSet& internalCells
00106 //)
00107 //{
00108 //    const pointField& points = mesh.points();
00109 //    const cellList& cells = mesh.cells();
00110 //    const faceList& faces = mesh.faces();
00111 //
00112 //    // Determine cells that have all points on the boundary.
00113 //    labelHashSet flatCandidates(getHangingCells(mesh, internalCells));
00114 //
00115 //    // All boundary points will become visible after subsetting and will be
00116 //    // snapped
00117 //    // to surface. Calculate new volume for cells using only these points and
00118 //    // check if it does not become too small.
00119 //
00120 //    // Get points used by flatCandidates.
00121 //    labelHashSet outsidePoints(flatCandidates.size());
00122 //
00123 //    // Snap outside points to surface
00124 //    pointField newPoints(points);
00125 //
00126 //    for
00127 //    (
00128 //        labelHashSet::const_iterator iter = flatCandidates.begin();
00129 //        iter != flatCandidates.end();
00130 //        ++iter
00131 //    )
00132 //    {
00133 //        const cell& cFaces = cells[iter.key()];
00134 //
00135 //        forAll(cFaces, cFaceI)
00136 //        {
00137 //            const face& f = faces[cFaces[cFaceI]];
00138 //
00139 //            forAll(f, fp)
00140 //            {
00141 //                label pointI = f[fp];
00142 //
00143 //                if (outsidePoints.insert(pointI))
00144 //                {
00145 //                    // Calculate new position for this outside point
00146 //                    tmp<pointField> tnearest =
00147 //                        querySurf.calcNearest(pointField(1, points[pointI]));
00148 //                    newPoints[pointI] = tnearest()[0];
00149 //                }
00150 //            }
00151 //        }               
00152 //    }
00153 //
00154 //
00155 //    // Calculate new volume for mixed cells
00156 //    label nRemoved = 0;
00157 //    for
00158 //    (
00159 //        labelHashSet::const_iterator iter = flatCandidates.begin();
00160 //        iter != flatCandidates.end();
00161 //        ++iter
00162 //    )
00163 //    {
00164 //        label cellI = iter.key();
00165 //    
00166 //        const cell& cll = cells[cellI];
00167 //
00168 //        scalar newVol = cll.mag(newPoints, faces);
00169 //        scalar oldVol = mesh.cellVolumes()[cellI];
00170 //
00171 //        if (newVol < 0.1 * oldVol)
00172 //        {
00173 //            internalCells.erase(cellI);
00174 //            nRemoved++;
00175 //        }
00176 //    }
00177 //
00178 //    return nRemoved;
00179 //}
00180 
00181 
00184 //void Foam::surfaceSets::getNearPoints
00185 //(
00186 //    const primitiveMesh& mesh,
00187 //    const triSurface&,
00188 //    const triSurfaceSearch& querySurf,
00189 //    const scalar edgeFactor,
00190 //    const pointSet& candidateSet,
00191 //    pointSet& nearPointSet
00192 //)
00193 //{
00194 //    if (edgeFactor <= 0)
00195 //    {
00196 //        return;
00197 //    }
00198 //
00199 //    labelList candidates(candidateSet.toc());
00200 //
00201 //    pointField candidatePoints(candidates.size());
00202 //    forAll(candidates, i)
00203 //    {
00204 //        candidatePoints[i] = mesh.points()[candidates[i]];
00205 //    }
00206 //
00207 //    tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
00208 //    const pointField& nearest = tnearest();
00209 //
00210 //    const pointField& points = mesh.points();
00211 //
00212 //    forAll(candidates, i)
00213 //    {
00214 //        label pointI = candidates[i];
00215 //
00216 //        scalar minLen = minEdgeLen(mesh, pointI);
00217 //
00218 //        scalar dist = mag(nearest[i] - points[pointI]);
00219 //
00220 //        if (dist < edgeFactor * minLen)
00221 //        {
00222 //            nearPointSet.insert(pointI);
00223 //        }
00224 //    }
00225 //}
00226 
00227 
00228 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00229 
00230 
00231 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00232 
00233 void Foam::surfaceSets::getSurfaceSets
00234 (
00235     const polyMesh& mesh,
00236     const fileName&,
00237     const triSurface&,
00238     const triSurfaceSearch& querySurf,
00239     const pointField& outsidePts,
00240 
00241     const label nCutLayers,
00242 
00243     labelHashSet& inside,
00244     labelHashSet& outside,
00245     labelHashSet& cut
00246 )
00247 {
00248     // Construct search engine on mesh
00249     meshSearch queryMesh(mesh, true);
00250 
00251     // Cut faces with surface and classify cells
00252     cellClassification cellType
00253     (
00254         mesh,
00255         queryMesh,
00256         querySurf,
00257         outsidePts
00258     );
00259 
00260     if (nCutLayers > 0)
00261     {
00262         // Trim cutCells so they are max nCutLayers away (as seen in point-cell
00263         // walk) from outside cells.
00264         cellType.trimCutCells
00265         (
00266             nCutLayers,
00267             cellClassification::OUTSIDE, 
00268             cellClassification::INSIDE
00269         );
00270     }
00271 
00272     forAll(cellType, cellI)
00273     {
00274         label cType = cellType[cellI];
00275 
00276         if (cType == cellClassification::CUT)
00277         {
00278             cut.insert(cellI);
00279         }
00280         else if (cType == cellClassification::INSIDE)
00281         {
00282             inside.insert(cellI);
00283         }
00284         else if (cType == cellClassification::OUTSIDE)
00285         {
00286             outside.insert(cellI);
00287         }
00288     }
00289 }
00290 
00291 
00292 Foam::labelHashSet Foam::surfaceSets::getHangingCells
00293 (
00294     const primitiveMesh& mesh,
00295     const labelHashSet& internalCells
00296 )
00297 {
00298     const cellList& cells = mesh.cells();
00299     const faceList& faces = mesh.faces();
00300 
00301 
00302     // Divide points into
00303     // -referenced by internal only
00304     // -referenced by outside only
00305     // -mixed
00306 
00307     List<pointStatus> pointSide(mesh.nPoints(), NOTSET);
00308 
00309     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
00310     {
00311         if (internalCells.found(cellI))
00312         {
00313             // Inside cell. Mark all vertices seen from this cell.
00314             const labelList& cFaces = cells[cellI];
00315 
00316             forAll(cFaces, cFaceI)
00317             {
00318                 const face& f = faces[cFaces[cFaceI]];
00319 
00320                 forAll(f, fp)
00321                 {
00322                     label pointI = f[fp];
00323 
00324                     if (pointSide[pointI] == NOTSET)
00325                     {
00326                         pointSide[pointI] = INSIDE;
00327                     }
00328                     else if (pointSide[pointI] == OUTSIDE)
00329                     {
00330                         pointSide[pointI] = MIXED;
00331                     }
00332                     else
00333                     {
00334                         // mixed or inside stay same
00335                     }
00336                 }
00337             }
00338         }
00339         else
00340         {
00341             // Outside cell
00342             const labelList& cFaces = cells[cellI];
00343 
00344             forAll(cFaces, cFaceI)
00345             {
00346                 const face& f = faces[cFaces[cFaceI]];
00347 
00348                 forAll(f, fp)
00349                 {
00350                     label pointI = f[fp];
00351 
00352                     if (pointSide[pointI] == NOTSET)
00353                     {
00354                         pointSide[pointI] = OUTSIDE;
00355                     }
00356                     else if (pointSide[pointI] == INSIDE)
00357                     {
00358                         pointSide[pointI] = MIXED;
00359                     }
00360                     else
00361                     {
00362                         // mixed or outside stay same
00363                     }
00364                 }
00365             }
00366         }
00367     }
00368 
00369 
00370     //OFstream mixedStr("mixed.obj");
00371     //
00372     //forAll(pointSide, pointI)
00373     //{
00374     //    if (pointSide[pointI] == MIXED)
00375     //    {
00376     //        const point& pt = points[pointI];
00377     //
00378     //        mixedStr << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
00379     //            << endl;
00380     //    }
00381     //}
00382 
00383 
00384     // Determine cells using mixed points only
00385 
00386     labelHashSet mixedOnlyCells(internalCells.size());
00387 
00388     for
00389     (
00390         labelHashSet::const_iterator iter = internalCells.begin();
00391         iter != internalCells.end();
00392         ++iter
00393     )
00394     {
00395         label cellI = iter.key();
00396 
00397         const cell& cFaces = cells[cellI];
00398 
00399         label usesMixedOnly = true;
00400 
00401         forAll(cFaces, i)
00402         {
00403             const face& f = faces[cFaces[i]];
00404 
00405             forAll(f, fp)
00406             {
00407                 if (pointSide[f[fp]] != MIXED)
00408                 {
00409                     usesMixedOnly = false;
00410                     break;
00411                 }
00412             }
00413 
00414             if (!usesMixedOnly)
00415             {
00416                 break;
00417             }
00418         }
00419         if (usesMixedOnly)
00420         {
00421             mixedOnlyCells.insert(cellI);
00422         }
00423     }
00424 
00425     return mixedOnlyCells;
00426 }
00427 
00428 
00429 //void Foam::surfaceSets::writeSurfaceSets
00430 //(
00431 //    const polyMesh& mesh,
00432 //    const fileName& surfName,
00433 //    const triSurface& surf,
00434 //    const triSurfaceSearch& querySurf,
00435 //    const pointField& outsidePts,
00436 //    const scalar edgeFactor
00437 //)
00438 //{
00439 //    // Cellsets for inside/outside determination
00440 //    cellSet rawInside(mesh, "rawInside", mesh.nCells()/10);
00441 //    cellSet rawOutside(mesh, "rawOutside", mesh.nCells()/10);
00442 //    cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
00443 //
00444 //    // Get inside/outside/cut cells
00445 //    getSurfaceSets
00446 //    (
00447 //        mesh,
00448 //        surfName,
00449 //        surf,
00450 //        querySurf,
00451 //        outsidePts,
00452 //
00453 //        rawInside,
00454 //        rawOutside,
00455 //        cutCells
00456 //    );
00457 //
00458 //
00459 //    Pout<< "rawInside:" << rawInside.size() << endl;
00460 //
00461 //    label nRemoved;
00462 //    do
00463 //    {
00464 //        nRemoved = removeHangingCells(mesh, querySurf, rawInside);
00465 //
00466 //        Pout<< nl
00467 //            << "Removed " << nRemoved
00468 //            << " rawInside cells that have all their points on the outside"
00469 //            << endl;
00470 //    }
00471 //    while (nRemoved != 0);
00472 //
00473 //    Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
00474 //        << rawInside.instance()/rawInside.local()/rawInside.name()
00475 //        << endl << endl;
00476 //    rawInside.write();
00477 //
00478 //
00479 //
00480 //    // Select outside cells
00481 //    surfaceToCell outsideSource
00482 //    (
00483 //        mesh,
00484 //        surfName,
00485 //        surf,
00486 //        querySurf, 
00487 //        outsidePts,
00488 //        false,          // includeCut
00489 //        false,          // includeInside
00490 //        true,           // includeOutside
00491 //        -GREAT,         // nearDist
00492 //        -GREAT          // curvature
00493 //    );
00494 //
00495 //    outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
00496 //
00497 //    Pout<< "rawOutside:" << rawOutside.size() << endl;
00498 //
00499 //    do
00500 //    {
00501 //        nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
00502 //
00503 //        Pout<< nl
00504 //            << "Removed " << nRemoved
00505 //            << " rawOutside cells that have all their points on the outside"
00506 //            << endl;
00507 //    }
00508 //    while (nRemoved != 0);
00509 //
00510 //    Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
00511 //        << rawOutside.instance()/rawOutside.local()/rawOutside.name()
00512 //        << endl << endl;
00513 //    rawOutside.write();
00514 //
00515 //
00516 //    // Select cut cells by negating inside and outside set.
00517 //    cutCells.invert(mesh.nCells());
00518 //
00519 //    cellToCell deleteInsideSource(mesh, rawInside.name());
00520 //
00521 //    deleteInsideSource.applyToSet(topoSetSource::DELETE, cutCells);
00522 //    Pout<< "Writing cut cells (" << cutCells.size() << ") to cellSet "
00523 //        << cutCells.instance()/cutCells.local()/cutCells.name()
00524 //        << endl << endl;
00525 //    cutCells.write();
00526 //
00527 //
00528 //    //
00529 //    // Remove cells with points too close to surface.
00530 //    //
00531 //
00532 //
00533 //    // Get all points in cutCells.
00534 //    pointSet cutPoints(mesh, "cutPoints", 4*cutCells.size());
00535 //    cellToPoint cutSource(mesh, "cutCells", cellToPoint::ALL);
00536 //    cutSource.applyToSet(topoSetSource::NEW, cutPoints);
00537 //
00538 //    // Get all points that are too close to surface.
00539 //    pointSet nearPoints(mesh, "nearPoints", cutPoints.size());
00540 //
00541 //    getNearPoints
00542 //    (
00543 //        mesh,
00544 //        surf,
00545 //        querySurf,
00546 //        edgeFactor,
00547 //        cutPoints,
00548 //        nearPoints
00549 //    );
00550 //
00551 //    Pout<< nl
00552 //        << "Selected " << nearPoints.size()
00553 //        << " points that are closer than " << edgeFactor
00554 //        << " times the local minimum lengthscale to the surface"
00555 //        << nl << endl;
00556 //
00557 //
00558 //    // Remove cells that use any of the points in nearPoints
00559 //    // from the inside and outsideCells.
00560 //    nearPoints.write();
00561 //    pointToCell pToCell(mesh, nearPoints.name(), pointToCell::ANY);
00562 //
00563 //
00564 //
00565 //    // Start off from copy of rawInside, rawOutside
00566 //    cellSet inside(mesh, "inside", rawInside);
00567 //    cellSet outside(mesh, "outside", rawOutside);
00568 //
00569 //    pToCell.applyToSet(topoSetSource::DELETE, inside);
00570 //    pToCell.applyToSet(topoSetSource::DELETE, outside);
00571 //
00572 //    Pout<< nl
00573 //        << "Removed " << rawInside.size() - inside.size()
00574 //        << " inside cells that are too close to the surface" << endl;
00575 //
00576 //    Pout<< nl
00577 //        << "Removed " << rawOutside.size() - outside.size()
00578 //        << " inside cells that are too close to the surface" << nl << endl;
00579 //
00580 //
00581 //
00582 //    //
00583 //    // Remove cells with one or no faces on rest of cellSet. Note: Problem is
00584 //    // not these cells an sich but rather that all of their vertices will be
00585 //    // outside vertices and thus projected onto surface. Which might (if they
00586 //    // project onto same surface) result in flattened cells.
00587 //    //
00588 //
00589 //    do
00590 //    {
00591 //        nRemoved = removeHangingCells(mesh, querySurf, inside);
00592 //
00593 //        Pout<< nl
00594 //            << "Removed " << nRemoved
00595 //            << " inside cells that have all their points on the outside"
00596 //            << endl;
00597 //    }
00598 //    while (nRemoved != 0);
00599 //    do
00600 //    {
00601 //        nRemoved = removeHangingCells(mesh, querySurf, outside);
00602 //
00603 //        Pout<< nl
00604 //            << "Removed " << nRemoved
00605 //            << " outside cells that have all their points on the inside"
00606 //            << endl;
00607 //    }
00608 //    while (nRemoved != 0);
00609 //    
00610 //
00611 //    //
00612 //    // Write
00613 //    //
00614 //
00615 //
00616 //    Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
00617 //        << inside.instance()/inside.local()/inside.name()
00618 //        << endl << endl;
00619 //
00620 //    inside.write();
00621 //
00622 //    Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
00623 //        << outside.instance()/outside.local()/outside.name()
00624 //        << endl << endl;
00625 //    
00626 //    outside.write();
00627 //}
00628 
00629 
00630 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
00631 
00632 
00633 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
00634 
00635 
00636 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
00637 
00638 
00639 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines