00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
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
00082 static const scalar edgeTol = 1E-3;
00083
00084
00085
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
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
00180
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
00193
00194
00195
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
00209 otherCellI = cellI;
00210
00211 break;
00212 }
00213 }
00214
00215 if (otherCellI == -1)
00216 {
00217
00218
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
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
00256
00257
00258
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
00286 return -1;
00287 }
00288 }
00289
00290
00291
00292
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
00307
00308 return -1;
00309 }
00310 }
00311 }
00312
00313 return axisIndex;
00314 }
00315
00316
00317
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
00335
00336
00337 bool readDict = args.optionFound("dict");
00338 bool overwrite = args.optionFound("overwrite");
00339
00340
00341 labelList refCells;
00342
00343
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
00378 refCells.setSize(mesh.nCells());
00379
00380 forAll(mesh.cells(), cellI)
00381 {
00382 refCells[cellI] = cellI;
00383 }
00384
00385
00386
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
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
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
00452 multiDirRefinement multiRef(mesh, refCells, refineDict);
00453
00454
00455
00456 if (overwrite)
00457 {
00458 mesh.setInstance(oldInstance);
00459 }
00460 mesh.write();
00461
00462
00463
00464
00465 const labelListList& oldToNew = multiRef.addedCells();
00466
00467
00468
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
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
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
00539
00540 printEdgeStats(mesh);
00541
00542 Info<< "End\n" << endl;
00543
00544 return 0;
00545 }
00546
00547
00548