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

splitCells.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-2010 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 Application
00025     splitCells
00026 
00027 Description
00028     Utility to split cells with flat faces.
00029 
00030     Uses a geometric cut with a plane dividing the edge angle into two so
00031     might produce funny cells. For hexes it will use by default a cut from
00032     edge onto opposite edge (i.e. purely topological).
00033 
00034     Options:
00035     - split cells from cellSet only
00036     - use geometric cut for hexes as well
00037 
00038     The angle is the angle between two faces sharing an edge as seen from
00039     inside each cell. So a cube will have all angles 90. If you want
00040     to split cells with cell centre outside use e.g. angle 200
00041 
00042 Usage
00043 
00044     - splitCells [OPTIONS] <edge angle [0..360]>
00045 
00046     @param <edge angle [0..360]> \n
00047     @todo Detailed description of argument.
00048 
00049     @param -geometry \n
00050     Geometry based splitting for hex cells.
00051 
00052     @param -tol <tolerance>\n
00053     Edge snap tolerance.
00054 
00055     @param -set <cellSet name>\n
00056     Cell set to split.
00057 
00058     @param -overwrite \n
00059     Overwrite existing data.
00060 
00061     @param -case <dir>\n
00062     Case directory.
00063 
00064     @param -help \n
00065     Display help message.
00066 
00067     @param -doc \n
00068     Display Doxygen API documentation page for this application.
00069 
00070     @param -srcDoc \n
00071     Display Doxygen source documentation page for this application.
00072 
00073 \*---------------------------------------------------------------------------*/
00074 
00075 #include <OpenFOAM/argList.H>
00076 #include <OpenFOAM/Time.H>
00077 #include <dynamicMesh/polyTopoChange.H>
00078 #include <dynamicMesh/polyTopoChanger.H>
00079 #include <OpenFOAM/mapPolyMesh.H>
00080 #include <OpenFOAM/polyMesh.H>
00081 #include <dynamicMesh/cellCuts.H>
00082 #include <meshTools/cellSet.H>
00083 #include <OpenFOAM/cellModeller.H>
00084 #include <dynamicMesh/meshCutter.H>
00085 #include <OpenFOAM/mathematicalConstants.H>
00086 #include <dynamicMesh/geomCellLooper.H>
00087 #include <OpenFOAM/plane.H>
00088 #include <dynamicMesh/edgeVertex.H>
00089 #include <meshTools/meshTools.H>
00090 #include <OpenFOAM/ListOps.H>
00091 
00092 using namespace Foam;
00093 
00094 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00095 
00096 
00097 labelList pack(const boolList& lst)
00098 {
00099     labelList packedLst(lst.size());
00100     label packedI = 0;
00101 
00102     forAll(lst, i)
00103     {
00104         if (lst[i])
00105         {
00106             packedLst[packedI++] = i;
00107         }
00108     }
00109     packedLst.setSize(packedI);
00110 
00111     return packedLst;
00112 }
00113 
00114 
00115 scalarField pack(const boolList& lst, const scalarField& elems)
00116 {
00117     scalarField packedElems(lst.size());
00118     label packedI = 0;
00119 
00120     forAll(lst, i)
00121     {
00122         if (lst[i])
00123         {
00124             packedElems[packedI++] = elems[i];
00125         }
00126     }
00127     packedElems.setSize(packedI);
00128 
00129     return packedElems;
00130 }
00131 
00132 
00133 // Given sin and cos of max angle between normals calculate whether f0 and f1
00134 // on cellI make larger angle. Uses sinAngle only for quadrant detection.
00135 bool largerAngle
00136 (
00137     const primitiveMesh& mesh,
00138     const scalar cosAngle,
00139     const scalar sinAngle,
00140 
00141     const label cellI,
00142     const label f0,             // face label
00143     const label f1,
00144 
00145     const vector& n0,           // normal at f0
00146     const vector& n1
00147 )
00148 {
00149     const labelList& own = mesh.faceOwner();
00150 
00151     bool sameFaceOrder = !((own[f0] == cellI) ^ (own[f1] == cellI));
00152 
00153     // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
00154     // gives -1.
00155     scalar normalCosAngle = n0 & n1;
00156 
00157     if (sameFaceOrder)
00158     {
00159         normalCosAngle = -normalCosAngle;
00160     }
00161 
00162 
00163     // Get cos between faceCentre and normal vector to determine in
00164     // which quadrant angle is. (Is correct for unwarped faces only!)
00165     // Correct for non-outwards pointing normal.
00166     vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
00167     c1c0 /= mag(c1c0) + VSMALL;
00168 
00169     scalar fcCosAngle = n0 & c1c0;
00170 
00171     if (own[f0] != cellI)
00172     {
00173         fcCosAngle = -fcCosAngle;
00174     }
00175 
00176     if (sinAngle < 0.0)
00177     {
00178         // Looking for concave angles (quadrant 3 or 4)
00179         if (fcCosAngle <= 0)
00180         {
00181             // Angle is convex so smaller.
00182             return false;
00183         }
00184         else
00185         {
00186             if (normalCosAngle < cosAngle)
00187             {
00188                 return false;
00189             }
00190             else
00191             {
00192                 return true;
00193             }
00194         }
00195     }
00196     else
00197     {
00198         // Looking for convex angles (quadrant 1 or 2)
00199         if (fcCosAngle > 0)
00200         {
00201             // Concave angle
00202             return true;
00203         }
00204         else
00205         {
00206             // Convex. Check cos of normal vectors.
00207             if (normalCosAngle > cosAngle)
00208             {
00209                 return false;
00210             }
00211             else
00212             {
00213                 return true;
00214             }
00215         }
00216     }
00217 }
00218 
00219 
00220 // Split hex (and hex only) along edgeI creating two prisms
00221 bool splitHex
00222 (
00223     const polyMesh& mesh,
00224     const label cellI,
00225     const label edgeI,
00226 
00227     DynamicList<label>& cutCells,
00228     DynamicList<labelList>& cellLoops,
00229     DynamicList<scalarField>& cellEdgeWeights
00230 )
00231 {
00232     // cut handling functions
00233     edgeVertex ev(mesh);
00234 
00235     const edgeList& edges = mesh.edges();
00236     const faceList& faces = mesh.faces();
00237 
00238     const edge& e = edges[edgeI];
00239 
00240     // Get faces on the side, i.e. faces not using edge but still using one of
00241     // the edge endpoints.
00242 
00243     label leftI = -1;
00244     label rightI = -1;
00245     label leftFp = -1;
00246     label rightFp = -1;
00247 
00248     const cell& cFaces = mesh.cells()[cellI];
00249 
00250     forAll(cFaces, i)
00251     {
00252         label faceI = cFaces[i];
00253 
00254         const face& f = faces[faceI];
00255 
00256         label fp0 = findIndex(f, e[0]);
00257         label fp1 = findIndex(f, e[1]);
00258 
00259         if (fp0 == -1)
00260         {
00261             if (fp1 != -1)
00262             {
00263                 // Face uses e[1] but not e[0]
00264                 rightI = faceI;
00265                 rightFp = fp1;
00266 
00267                 if (leftI != -1)
00268                 {
00269                     // Have both faces so exit
00270                     break;
00271                 }
00272             }
00273         }
00274         else
00275         {
00276             if (fp1 != -1)
00277             {
00278                 // Face uses both e[1] and e[0]
00279             }
00280             else
00281             {
00282                 leftI = faceI;
00283                 leftFp = fp0;
00284 
00285                 if (rightI != -1)
00286                 {
00287                     break;
00288                 }
00289             }
00290         }
00291     }
00292 
00293     if (leftI == -1 || rightI == -1)
00294     {
00295         FatalErrorIn("splitHex") << "Problem : leftI:" << leftI
00296             << " rightI:" << rightI << abort(FatalError);
00297     }
00298 
00299     // Walk two vertices further on faces.
00300 
00301     const face& leftF = faces[leftI];
00302 
00303     label leftV = leftF[(leftFp + 2) % leftF.size()];
00304 
00305     const face& rightF = faces[rightI];
00306 
00307     label rightV = rightF[(rightFp + 2) % rightF.size()];
00308 
00309     labelList loop(4);
00310     loop[0] = ev.vertToEVert(e[0]);
00311     loop[1] = ev.vertToEVert(leftV);
00312     loop[2] = ev.vertToEVert(rightV);
00313     loop[3] = ev.vertToEVert(e[1]);
00314 
00315     scalarField loopWeights(4);
00316     loopWeights[0] = -GREAT;
00317     loopWeights[1] = -GREAT;
00318     loopWeights[2] = -GREAT;
00319     loopWeights[3] = -GREAT;
00320 
00321     cutCells.append(cellI);
00322     cellLoops.append(loop);
00323     cellEdgeWeights.append(loopWeights);
00324 
00325     return true;
00326 }
00327 
00328 
00329 // Split cellI along edgeI with a plane along halfNorm direction.
00330 bool splitCell
00331 (
00332     const polyMesh& mesh,
00333     const geomCellLooper& cellCutter,
00334 
00335     const label cellI,
00336     const label edgeI,
00337     const vector& halfNorm,
00338 
00339     const boolList& vertIsCut,
00340     const boolList& edgeIsCut,
00341     const scalarField& edgeWeight,
00342 
00343     DynamicList<label>& cutCells,
00344     DynamicList<labelList>& cellLoops,
00345     DynamicList<scalarField>& cellEdgeWeights
00346 )
00347 {
00348     const edge& e = mesh.edges()[edgeI];
00349 
00350     vector eVec = e.vec(mesh.points());
00351     eVec /= mag(eVec);
00352 
00353     vector planeN = eVec ^ halfNorm;
00354 
00355     // Slightly tilt plane to make it not cut edges exactly
00356     // halfway on fully regular meshes (since we want cuts
00357     // to be snapped to vertices)
00358     planeN += 0.01*halfNorm;
00359 
00360     planeN /= mag(planeN);
00361 
00362     // Define plane through edge
00363     plane cutPlane(mesh.points()[e.start()], planeN);
00364 
00365     labelList loop;
00366     scalarField loopWeights;
00367 
00368     if
00369     (
00370         cellCutter.cut
00371         (
00372             cutPlane,
00373             cellI,
00374             vertIsCut,
00375             edgeIsCut,
00376             edgeWeight,
00377             loop,
00378             loopWeights
00379         )
00380     )
00381     {
00382         // Did manage to cut cell. Copy into overall list.
00383         cutCells.append(cellI);
00384         cellLoops.append(loop);
00385         cellEdgeWeights.append(loopWeights);
00386 
00387         return true;
00388     }
00389     else
00390     {
00391         return false;
00392     }
00393 }
00394 
00395 
00396 // Collects cuts for all cells in cellSet
00397 void collectCuts
00398 (
00399     const polyMesh& mesh,
00400     const geomCellLooper& cellCutter,
00401     const bool geometry,
00402     const scalar minCos,
00403     const scalar minSin,
00404     const cellSet& cellsToCut,
00405 
00406     DynamicList<label>& cutCells,
00407     DynamicList<labelList>& cellLoops,
00408     DynamicList<scalarField>& cellEdgeWeights
00409 )
00410 {
00411     // Get data from mesh
00412     const cellShapeList& cellShapes = mesh.cellShapes();
00413     const labelList& own = mesh.faceOwner();
00414     const labelListList& cellEdges = mesh.cellEdges();
00415     const vectorField& faceAreas = mesh.faceAreas();
00416 
00417     // Hex shape
00418     const cellModel& hex = *(cellModeller::lookup("hex"));
00419 
00420     // cut handling functions
00421     edgeVertex ev(mesh);
00422 
00423 
00424     // Cut information per mesh entity
00425     boolList vertIsCut(mesh.nPoints(), false);
00426     boolList edgeIsCut(mesh.nEdges(), false);
00427     scalarField edgeWeight(mesh.nEdges(), -GREAT);
00428 
00429     for
00430     (
00431         cellSet::const_iterator iter = cellsToCut.begin();
00432         iter != cellsToCut.end();
00433         ++iter
00434     )
00435     {
00436         label cellI = iter.key();
00437 
00438         const labelList& cEdges = cellEdges[cellI];
00439 
00440         forAll(cEdges, i)
00441         {
00442             label edgeI = cEdges[i];
00443 
00444             label f0, f1;
00445             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
00446 
00447             vector n0 = faceAreas[f0];
00448             n0 /= mag(n0);
00449 
00450             vector n1 = faceAreas[f1];
00451             n1 /= mag(n1);
00452 
00453             if
00454             (
00455                 largerAngle
00456                 (
00457                     mesh,
00458                     minCos,
00459                     minSin,
00460 
00461                     cellI,
00462                     f0,
00463                     f1,
00464                     n0,
00465                     n1
00466                 )
00467             )
00468             {
00469                 bool splitOk = false;
00470 
00471                 if (!geometry && cellShapes[cellI].model() == hex)
00472                 {
00473                     splitOk =
00474                         splitHex
00475                         (
00476                             mesh,
00477                             cellI,
00478                             edgeI,
00479 
00480                             cutCells,
00481                             cellLoops,
00482                             cellEdgeWeights
00483                         );
00484                 }
00485                 else
00486                 {
00487                     vector halfNorm;
00488 
00489                     if ((own[f0] == cellI) ^ (own[f1] == cellI))
00490                     {
00491                         // Opposite owner orientation
00492                         halfNorm = 0.5*(n0 - n1);
00493                     }
00494                     else
00495                     {
00496                         // Faces have same owner or same neighbour so
00497                         // normals point in same direction
00498                         halfNorm = 0.5*(n0 + n1);
00499                     }
00500 
00501                     splitOk =
00502                         splitCell
00503                         (
00504                             mesh,
00505                             cellCutter,
00506                             cellI,
00507                             edgeI,
00508                             halfNorm,
00509 
00510                             vertIsCut,
00511                             edgeIsCut,
00512                             edgeWeight,
00513 
00514                             cutCells,
00515                             cellLoops,
00516                             cellEdgeWeights
00517                         );
00518                 }
00519 
00520 
00521                 if (splitOk)
00522                 {
00523                     // Update cell/edge/vertex wise info.
00524                     label index = cellLoops.size() - 1;
00525                     const labelList& loop = cellLoops[index];
00526                     const scalarField& loopWeights = cellEdgeWeights[index];
00527 
00528                     forAll(loop, i)
00529                     {
00530                         label cut = loop[i];
00531 
00532                         if (ev.isEdge(cut))
00533                         {
00534                             edgeIsCut[ev.getEdge(cut)] = true;
00535                             edgeWeight[ev.getEdge(cut)] = loopWeights[i];
00536                         }
00537                         else
00538                         {
00539                             vertIsCut[ev.getVertex(cut)] = true;
00540                         }
00541                     }
00542 
00543                     // Stop checking edges for this cell.
00544                     break;
00545                 }
00546             }
00547         }
00548     }
00549 
00550     cutCells.shrink();
00551     cellLoops.shrink();
00552     cellEdgeWeights.shrink();
00553 }
00554 
00555 
00556 // Main program:
00557 
00558 int main(int argc, char *argv[])
00559 {
00560     argList::noParallel();
00561     argList::validOptions.insert("set", "cellSet name");
00562     argList::validOptions.insert("geometry", "");
00563     argList::validOptions.insert("tol", "edge snap tolerance");
00564     argList::validOptions.insert("overwrite", "");
00565     argList::validArgs.append("edge angle [0..360]");
00566 
00567 #   include <OpenFOAM/setRootCase.H>
00568 #   include <OpenFOAM/createTime.H>
00569     runTime.functionObjects().off();
00570 #   include <OpenFOAM/createPolyMesh.H>
00571     const word oldInstance = mesh.pointsInstance();
00572 
00573     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
00574 
00575     scalar radAngle = featureAngle * mathematicalConstant::pi/180.0;
00576     scalar minCos = Foam::cos(radAngle);
00577     scalar minSin = Foam::sin(radAngle);
00578 
00579     bool readSet   = args.optionFound("set");
00580     bool geometry  = args.optionFound("geometry");
00581     bool overwrite = args.optionFound("overwrite");
00582 
00583     scalar edgeTol = 0.2;
00584     args.optionReadIfPresent("tol", edgeTol);
00585 
00586     Info<< "Trying to split cells with internal angles > feature angle\n" << nl
00587         << "featureAngle      : " << featureAngle << nl
00588         << "edge snapping tol : " << edgeTol << nl;
00589     if (readSet)
00590     {
00591         Info<< "candidate cells   : cellSet " << args.option("set") << nl;
00592     }
00593     else
00594     {
00595         Info<< "candidate cells   : all cells" << nl;
00596     }
00597     if (geometry)
00598     {
00599         Info<< "hex cuts          : geometric; using edge tolerance" << nl;
00600     }
00601     else
00602     {
00603         Info<< "hex cuts          : topological; cut to opposite edge" << nl;
00604     }
00605     Info<< endl;
00606 
00607 
00608     // Cell circumference cutter
00609     geomCellLooper cellCutter(mesh);
00610     // Snap all edge cuts close to endpoints to vertices.
00611     geomCellLooper::setSnapTol(edgeTol);
00612 
00613     // Candidate cells to cut
00614     cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
00615 
00616     if (readSet)
00617     {
00618         // Read cells to cut from cellSet
00619         cellSet cells(mesh, args.option("set"));
00620 
00621         cellsToCut = cells;
00622     }
00623 
00624     while (true)
00625     {
00626         if (!readSet)
00627         {
00628             // Try all cells for cutting
00629             for (label cellI = 0; cellI < mesh.nCells(); cellI++)
00630             {
00631                 cellsToCut.insert(cellI);
00632             }
00633         }
00634 
00635 
00636         // Cut information per cut cell
00637         DynamicList<label> cutCells(mesh.nCells()/10 + 10);
00638         DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
00639         DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
00640 
00641         collectCuts
00642         (
00643             mesh,
00644             cellCutter,
00645             geometry,
00646             minCos,
00647             minSin,
00648             cellsToCut,
00649 
00650             cutCells,
00651             cellLoops,
00652             cellEdgeWeights
00653         );
00654 
00655         cellSet cutSet(mesh, "cutSet", cutCells.size());
00656         forAll(cutCells, i)
00657         {
00658             cutSet.insert(cutCells[i]);
00659         }
00660 
00661         // Gets cuts across cells from cuts through edges.
00662         Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
00663             << cutSet.instance()/cutSet.local()/cutSet.name()
00664             << endl << endl;
00665         cutSet.write();
00666 
00667         // Analyze cuts for clashes.
00668         cellCuts cuts
00669         (
00670             mesh,
00671             cutCells,       // cells candidate for cutting
00672             cellLoops,
00673             cellEdgeWeights
00674         );
00675 
00676         Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
00677 
00678         if (cuts.nLoops() == 0)
00679         {
00680             break;
00681         }
00682 
00683         // Remove cut cells from cellsToCut  (Note:only relevant if -readSet)
00684         forAll(cuts.cellLoops(), cellI)
00685         {
00686             if (cuts.cellLoops()[cellI].size())
00687             {
00688                 //Info<< "Removing cut cell " << cellI << " from wishlist"
00689                 //    << endl;
00690                 cellsToCut.erase(cellI);
00691             }
00692         }
00693 
00694         // At least some cells are cut.
00695         polyTopoChange meshMod(mesh);
00696 
00697         // Cutting engine
00698         meshCutter cutter(mesh);
00699 
00700         // Insert mesh refinement into polyTopoChange.
00701         cutter.setRefinement(cuts, meshMod);
00702 
00703         // Do all changes
00704         Info<< "Morphing ..." << endl;
00705 
00706         if (!overwrite)
00707         {
00708             runTime++;
00709         }
00710 
00711         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
00712 
00713         if (morphMap().hasMotionPoints())
00714         {
00715             mesh.movePoints(morphMap().preMotionPoints());
00716         }
00717 
00718         // Update stored labels on meshCutter
00719         cutter.updateMesh(morphMap());
00720 
00721         // Update cellSet
00722         cellsToCut.updateMesh(morphMap());
00723 
00724         Info<< "Remaining:" << cellsToCut.size() << endl;
00725 
00726         // Write resulting mesh
00727         if (overwrite)
00728         {
00729             mesh.setInstance(oldInstance);
00730         }
00731 
00732         Info<< "Writing refined morphMesh to time " << runTime.timeName()
00733             << endl;
00734 
00735         mesh.write();
00736     }
00737 
00738     Info << "End\n" << endl;
00739 
00740     return 0;
00741 }
00742 
00743 
00744 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines