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
00063
00064
00065
00066
00067
00068
00069
00070
00071
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
00134
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,
00143 const label f1,
00144
00145 const vector& n0,
00146 const vector& n1
00147 )
00148 {
00149 const labelList& own = mesh.faceOwner();
00150
00151 bool sameFaceOrder = !((own[f0] == cellI) ^ (own[f1] == cellI));
00152
00153
00154
00155 scalar normalCosAngle = n0 & n1;
00156
00157 if (sameFaceOrder)
00158 {
00159 normalCosAngle = -normalCosAngle;
00160 }
00161
00162
00163
00164
00165
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
00179 if (fcCosAngle <= 0)
00180 {
00181
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
00199 if (fcCosAngle > 0)
00200 {
00201
00202 return true;
00203 }
00204 else
00205 {
00206
00207 if (normalCosAngle > cosAngle)
00208 {
00209 return false;
00210 }
00211 else
00212 {
00213 return true;
00214 }
00215 }
00216 }
00217 }
00218
00219
00220
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
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
00241
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
00264 rightI = faceI;
00265 rightFp = fp1;
00266
00267 if (leftI != -1)
00268 {
00269
00270 break;
00271 }
00272 }
00273 }
00274 else
00275 {
00276 if (fp1 != -1)
00277 {
00278
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
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
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
00356
00357
00358 planeN += 0.01*halfNorm;
00359
00360 planeN /= mag(planeN);
00361
00362
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
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
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
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
00418 const cellModel& hex = *(cellModeller::lookup("hex"));
00419
00420
00421 edgeVertex ev(mesh);
00422
00423
00424
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
00492 halfNorm = 0.5*(n0 - n1);
00493 }
00494 else
00495 {
00496
00497
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
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
00544 break;
00545 }
00546 }
00547 }
00548 }
00549
00550 cutCells.shrink();
00551 cellLoops.shrink();
00552 cellEdgeWeights.shrink();
00553 }
00554
00555
00556
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
00609 geomCellLooper cellCutter(mesh);
00610
00611 geomCellLooper::setSnapTol(edgeTol);
00612
00613
00614 cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
00615
00616 if (readSet)
00617 {
00618
00619 cellSet cells(mesh, args.option("set"));
00620
00621 cellsToCut = cells;
00622 }
00623
00624 while (true)
00625 {
00626 if (!readSet)
00627 {
00628
00629 for (label cellI = 0; cellI < mesh.nCells(); cellI++)
00630 {
00631 cellsToCut.insert(cellI);
00632 }
00633 }
00634
00635
00636
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
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
00668 cellCuts cuts
00669 (
00670 mesh,
00671 cutCells,
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
00684 forAll(cuts.cellLoops(), cellI)
00685 {
00686 if (cuts.cellLoops()[cellI].size())
00687 {
00688
00689
00690 cellsToCut.erase(cellI);
00691 }
00692 }
00693
00694
00695 polyTopoChange meshMod(mesh);
00696
00697
00698 meshCutter cutter(mesh);
00699
00700
00701 cutter.setRefinement(cuts, meshMod);
00702
00703
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
00719 cutter.updateMesh(morphMap());
00720
00721
00722 cellsToCut.updateMesh(morphMap());
00723
00724 Info<< "Remaining:" << cellsToCut.size() << endl;
00725
00726
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