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

refineMesh.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     refineMesh
00026 
00027 Description
00028     Utility to refine cells in multiple directions.
00029 
00030     Either supply -all option to refine all cells (3D refinement for 3D
00031     cases; 2D for 2D cases) or reads a refineMeshDict with
00032     - cellSet to refine
00033     - directions to refine
00034 
00035 Usage
00036 
00037     - refineMesh [OPTIONS]
00038 
00039     @param -dict <specify refineMeshDict>\n
00040     Refine according to refineMeshDict.
00041 
00042     @param -overwrite \n
00043     Overwrite existing data.
00044 
00045     @param -case <dir>\n
00046     Case directory.
00047 
00048     @param -parallel \n
00049     Run in parallel.
00050 
00051     @param -help \n
00052     Display help message.
00053 
00054     @param -doc \n
00055     Display Doxygen API documentation page for this application.
00056 
00057     @param -srcDoc \n
00058     Display Doxygen source documentation page for this application.
00059 
00060 \*---------------------------------------------------------------------------*/
00061 
00062 #include <OpenFOAM/argList.H>
00063 #include <OpenFOAM/polyMesh.H>
00064 #include <OpenFOAM/Time.H>
00065 #include <dynamicMesh/undoableMeshCutter.H>
00066 #include <dynamicMesh/hexCellLooper.H>
00067 #include <meshTools/cellSet.H>
00068 #include <meshTools/twoDPointCorrector.H>
00069 #include <dynamicMesh/directions.H>
00070 #include <OpenFOAM/OFstream.H>
00071 #include <dynamicMesh/multiDirRefinement.H>
00072 #include <OpenFOAM/labelIOList.H>
00073 #include <OpenFOAM/wedgePolyPatch.H>
00074 #include <OpenFOAM/plane.H>
00075 
00076 using namespace Foam;
00077 
00078 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00079 
00080 
00081 // Max cos angle for edges to be considered aligned with axis.
00082 static const scalar edgeTol = 1E-3;
00083 
00084 
00085 // Calculate some edge statistics on mesh.
00086 void printEdgeStats(const primitiveMesh& mesh)
00087 {
00088     label nX = 0;
00089     label nY = 0;
00090     label nZ = 0;
00091 
00092     scalar minX = GREAT;
00093     scalar maxX = -GREAT;
00094     vector x(1, 0, 0);
00095 
00096     scalar minY = GREAT;
00097     scalar maxY = -GREAT;
00098     vector y(0, 1, 0);
00099 
00100     scalar minZ = GREAT;
00101     scalar maxZ = -GREAT;
00102     vector z(0, 0, 1);
00103 
00104     scalar minOther = GREAT;
00105     scalar maxOther = -GREAT;
00106 
00107     const edgeList& edges = mesh.edges();
00108 
00109     forAll(edges, edgeI)
00110     {
00111         const edge& e = edges[edgeI];
00112 
00113         vector eVec(e.vec(mesh.points()));
00114 
00115         scalar eMag = mag(eVec);
00116 
00117         eVec /= eMag;
00118 
00119         if (mag(eVec & x) > 1-edgeTol)
00120         {
00121             minX = min(minX, eMag);
00122             maxX = max(maxX, eMag);
00123             nX++;
00124         }
00125         else if (mag(eVec & y) > 1-edgeTol)
00126         {
00127             minY = min(minY, eMag);
00128             maxY = max(maxY, eMag);
00129             nY++;
00130         }
00131         else if (mag(eVec & z) > 1-edgeTol)
00132         {
00133             minZ = min(minZ, eMag);
00134             maxZ = max(maxZ, eMag);
00135             nZ++;
00136         }
00137         else
00138         {
00139             minOther = min(minOther, eMag);
00140             maxOther = max(maxOther, eMag);
00141         }
00142     }
00143 
00144     Pout<< "Mesh edge statistics:" << endl
00145         << "    x aligned :  number:" << nX << "\tminLen:" << minX
00146         << "\tmaxLen:" << maxX << endl
00147         << "    y aligned :  number:" << nY << "\tminLen:" << minY
00148         << "\tmaxLen:" << maxY << endl
00149         << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
00150         << "\tmaxLen:" << maxZ << endl
00151         << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
00152         << "\tminLen:" << minOther
00153         << "\tmaxLen:" << maxOther << endl << endl;
00154 }
00155 
00156 
00157 // Return index of coordinate axis.
00158 label axis(const vector& normal)
00159 {
00160     label axisIndex = -1;
00161 
00162     if (mag(normal & point(1, 0, 0)) > (1-edgeTol))
00163     {
00164         axisIndex = 0;
00165     }
00166     else if (mag(normal & point(0, 1, 0)) > (1-edgeTol))
00167     {
00168         axisIndex = 1;
00169     }
00170     else if (mag(normal & point(0, 0, 1)) > (1-edgeTol))
00171     {
00172         axisIndex = 2;
00173     }
00174 
00175     return axisIndex;
00176 }
00177 
00178 
00179 //- Returns -1 or cartesian coordinate component (0=x, 1=y, 2=z) of normal
00180 //  in case of 2D mesh
00181 label twoDNess(const polyMesh& mesh)
00182 {
00183     const pointField& ctrs = mesh.cellCentres();
00184 
00185     if (ctrs.size() < 2)
00186     {
00187         return -1;
00188     }
00189 
00190 
00191     //
00192     // 1. All cell centres on single plane aligned with x, y or z
00193     //
00194 
00195     // Determine 3 points to base plane on.
00196     vector vec10 = ctrs[1] - ctrs[0];
00197     vec10 /= mag(vec10);
00198 
00199     label otherCellI = -1;
00200 
00201     for (label cellI = 2; cellI < ctrs.size(); cellI++)
00202     {
00203         vector vec(ctrs[cellI] - ctrs[0]);
00204         vec /= mag(vec);
00205 
00206         if (mag(vec & vec10) < 0.9)
00207         {
00208             // ctrs[cellI] not in line with n
00209             otherCellI = cellI;
00210 
00211             break;
00212         }
00213     }
00214 
00215     if (otherCellI == -1)
00216     {
00217         // Cannot find cell to make decent angle with cell0-cell1 vector.
00218         // Note: what to do here? All cells (almost) in one line. Maybe 1D case?
00219         return -1;
00220     }
00221 
00222     plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);
00223 
00224 
00225     forAll(ctrs, cellI)
00226     {
00227         const labelList& cEdges = mesh.cellEdges()[cellI];
00228 
00229         scalar minLen = GREAT;
00230 
00231         forAll(cEdges, i)
00232         {
00233             minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
00234         }
00235 
00236         if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
00237         {
00238             // Centres not in plane
00239             return  -1;
00240         }
00241     }
00242 
00243     label axisIndex = axis(cellPlane.normal());
00244 
00245     if (axisIndex == -1)
00246     {
00247         return axisIndex;
00248     }
00249 
00250 
00251     const polyBoundaryMesh& patches = mesh.boundaryMesh();
00252 
00253 
00254     //
00255     // 2. No edges without points on boundary
00256     //
00257 
00258     // Mark boundary points
00259     boolList boundaryPoint(mesh.points().size(), false);
00260 
00261     forAll(patches, patchI)
00262     {
00263         const polyPatch& patch = patches[patchI];
00264 
00265         forAll(patch, patchFaceI)
00266         {
00267             const face& f = patch[patchFaceI];
00268 
00269             forAll(f, fp)
00270             {
00271                 boundaryPoint[f[fp]] = true;
00272             }
00273         }
00274     }
00275 
00276 
00277     const edgeList& edges = mesh.edges();
00278 
00279     forAll(edges, edgeI)
00280     {
00281         const edge& e = edges[edgeI];
00282 
00283         if (!boundaryPoint[e.start()] && !boundaryPoint[e.end()])
00284         {
00285             // Edge has no point on boundary.
00286             return -1;
00287         }
00288     }
00289 
00290 
00291     // 3. For all non-wedge patches: all faces either perp or aligned with
00292     //    cell-plane normal. (wedge patches already checked upon construction)
00293 
00294     forAll(patches, patchI)
00295     {
00296         const polyPatch& patch = patches[patchI];
00297 
00298         if (!isA<wedgePolyPatch>(patch))
00299         {
00300             const vectorField& n = patch.faceAreas();
00301 
00302             scalarField cosAngle = mag(n/mag(n) & cellPlane.normal());
00303 
00304             if (mag(min(cosAngle) - max(cosAngle)) > 1E-6)
00305             {
00306                 // cosAngle should be either ~1 over all faces (2D front and
00307                 // back) or ~0 (all other patches perp to 2D)
00308                 return -1;
00309             }
00310         }
00311     }
00312 
00313     return axisIndex;
00314 }
00315 
00316 
00317 // Main program:
00318 
00319 int main(int argc, char *argv[])
00320 {
00321     Foam::argList::validOptions.insert("dict", "");
00322     Foam::argList::validOptions.insert("overwrite", "");
00323 
00324 #   include <OpenFOAM/setRootCase.H>
00325 #   include <OpenFOAM/createTime.H>
00326     runTime.functionObjects().off();
00327 #   include <OpenFOAM/createPolyMesh.H>
00328     const word oldInstance = mesh.pointsInstance();
00329 
00330     printEdgeStats(mesh);
00331 
00332 
00333     //
00334     // Read/construct control dictionary
00335     //
00336 
00337     bool readDict = args.optionFound("dict");
00338     bool overwrite = args.optionFound("overwrite");
00339 
00340     // List of cells to refine
00341     labelList refCells;
00342 
00343     // Dictionary to control refinement
00344     dictionary refineDict;
00345 
00346     if (readDict)
00347     {
00348         Info<< "Refining according to refineMeshDict" << nl << endl;
00349 
00350         refineDict =
00351             IOdictionary
00352             (
00353                 IOobject
00354                 (
00355                     "refineMeshDict",
00356                     runTime.system(),
00357                     mesh,
00358                     IOobject::MUST_READ,
00359                     IOobject::NO_WRITE
00360                 )
00361             );
00362 
00363         word setName(refineDict.lookup("set"));
00364 
00365         cellSet cells(mesh, setName);
00366 
00367         Pout<< "Read " << cells.size() << " cells from cellSet "
00368             << cells.instance()/cells.local()/cells.name()
00369             << endl << endl;
00370 
00371         refCells = cells.toc();
00372     }
00373     else
00374     {
00375         Info<< "Refining all cells" << nl << endl;
00376 
00377         // Select all cells
00378         refCells.setSize(mesh.nCells());
00379 
00380         forAll(mesh.cells(), cellI)
00381         {
00382             refCells[cellI] = cellI;
00383         }
00384 
00385 
00386         // Set refinement directions based on 2D/3D
00387         label axisIndex = twoDNess(mesh);
00388 
00389         if (axisIndex == -1)
00390         {
00391             Info<< "3D case; refining all directions" << nl << endl;
00392 
00393             wordList directions(3);
00394             directions[0] = "tan1";
00395             directions[1] = "tan2";
00396             directions[2] = "normal";
00397             refineDict.add("directions", directions);
00398 
00399             // Use hex cutter
00400             refineDict.add("useHexTopology", "true");
00401         }
00402         else
00403         {
00404             wordList directions(2);
00405 
00406             if (axisIndex == 0)
00407             {
00408                 Info<< "2D case; refining in directions y,z\n" << endl;
00409                 directions[0] = "tan2";
00410                 directions[1] = "normal";
00411             }
00412             else if (axisIndex == 1)
00413             {
00414                 Info<< "2D case; refining in directions x,z\n" << endl;
00415                 directions[0] = "tan1";
00416                 directions[1] = "normal";
00417             }
00418             else
00419             {
00420                 Info<< "2D case; refining in directions x,y\n" << endl;
00421                 directions[0] = "tan1";
00422                 directions[1] = "tan2";
00423             }
00424 
00425             refineDict.add("directions", directions);
00426 
00427             // Use standard cutter
00428             refineDict.add("useHexTopology", "false");
00429         }
00430 
00431         refineDict.add("coordinateSystem", "global");
00432 
00433         dictionary coeffsDict;
00434         coeffsDict.add("tan1", vector(1, 0, 0));
00435         coeffsDict.add("tan2", vector(0, 1, 0));
00436         refineDict.add("globalCoeffs", coeffsDict);
00437 
00438         refineDict.add("geometricCut", "false");
00439         refineDict.add("writeMesh", "false");
00440     }
00441 
00442 
00443     string oldTimeName(runTime.timeName());
00444 
00445     if (!overwrite)
00446     {
00447         runTime++;
00448     }
00449 
00450 
00451     // Multi-directional refinement (does multiple iterations)
00452     multiDirRefinement multiRef(mesh, refCells, refineDict);
00453 
00454 
00455     // Write resulting mesh
00456     if (overwrite)
00457     {
00458         mesh.setInstance(oldInstance);
00459     }
00460     mesh.write();
00461 
00462 
00463     // Get list of cell splits.
00464     // (is for every cell in old mesh the cells they have been split into)
00465     const labelListList& oldToNew = multiRef.addedCells();
00466 
00467 
00468     // Create cellSet with added cells for easy inspection
00469     cellSet newCells(mesh, "refinedCells", refCells.size());
00470 
00471     forAll(oldToNew, oldCellI)
00472     {
00473         const labelList& added = oldToNew[oldCellI];
00474 
00475         forAll(added, i)
00476         {
00477             newCells.insert(added[i]);
00478         }
00479     }
00480 
00481     Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
00482         << newCells.instance()/newCells.local()/newCells.name()
00483         << endl << endl;
00484 
00485     newCells.write();
00486 
00487 
00488 
00489 
00490     //
00491     // Invert cell split to construct map from new to old
00492     //
00493 
00494     labelIOList newToOld
00495     (
00496         IOobject
00497         (
00498             "cellMap",
00499             runTime.timeName(),
00500             polyMesh::meshSubDir,
00501             mesh,
00502             IOobject::NO_READ,
00503             IOobject::AUTO_WRITE
00504         ),
00505         mesh.nCells()
00506     );
00507     newToOld.note() =
00508         "From cells in mesh at "
00509       + runTime.timeName()
00510       + " to cells in mesh at "
00511       + oldTimeName;
00512 
00513 
00514     forAll(oldToNew, oldCellI)
00515     {
00516         const labelList& added = oldToNew[oldCellI];
00517 
00518         if (added.size())
00519         {
00520             forAll(added, i)
00521             {
00522                 newToOld[added[i]] = oldCellI;
00523             }
00524         }
00525         else
00526         {
00527             // Unrefined cell
00528             newToOld[oldCellI] = oldCellI;
00529         }
00530     }
00531 
00532     Info<< "Writing map from new to old cell to "
00533         << newToOld.objectPath() << nl << endl;
00534 
00535     newToOld.write();
00536 
00537 
00538     // Some statistics.
00539 
00540     printEdgeStats(mesh);
00541 
00542     Info<< "End\n" << endl;
00543 
00544     return 0;
00545 }
00546 
00547 
00548 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines