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: ************************ //