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 #include "topoCellLooper.H"
00027 #include <meshTools/cellFeatures.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mathematicalConstants.H>
00030 #include <OpenFOAM/DynamicList.H>
00031 #include <OpenFOAM/ListOps.H>
00032 #include <meshTools/meshTools.H>
00033 #include <OpenFOAM/hexMatcher.H>
00034
00035 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00036
00037
00038
00039
00040 namespace Foam
00041 {
00042 defineTypeNameAndDebug(topoCellLooper, 0);
00043 addToRunTimeSelectionTable(cellLooper, topoCellLooper, word);
00044 }
00045
00046
00047 const Foam::scalar Foam::topoCellLooper::featureCos =
00048 Foam::cos(10.0 * mathematicalConstant::pi/180.0);
00049
00050
00051
00052
00053
00054 template <class T>
00055 void Foam::topoCellLooper::subsetList
00056 (
00057 const label startI,
00058 const label freeI,
00059 DynamicList<T>& lst
00060 )
00061 {
00062 if (startI == 0)
00063 {
00064
00065
00066 if (freeI < 0)
00067 {
00068 FatalErrorIn("topoCellLooper::subsetList")
00069 << "startI:" << startI << " freeI:" << freeI
00070 << " lst:" << lst << abort(FatalError);
00071 }
00072 lst.setCapacity(freeI);
00073 }
00074 else
00075 {
00076
00077 label newI = 0;
00078 for (label elemI = startI; elemI < freeI; elemI++)
00079 {
00080 lst[newI++] = lst[elemI];
00081 }
00082
00083 if ((freeI - startI) < 0)
00084 {
00085 FatalErrorIn("topoCellLooper::subsetList")
00086 << "startI:" << startI << " freeI:" << freeI
00087 << " lst:" << lst << abort(FatalError);
00088 }
00089
00090 lst.setCapacity(freeI - startI);
00091 }
00092 }
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 void Foam::topoCellLooper::walkFace
00137 (
00138 const cellFeatures& features,
00139 const label faceI,
00140 const label startEdgeI,
00141 const label startVertI,
00142 const label nFeaturePts,
00143
00144 label& edgeI,
00145 label& vertI
00146 ) const
00147 {
00148 const labelList& fEdges = mesh().faceEdges()[faceI];
00149
00150 edgeI = startEdgeI;
00151
00152 vertI = startVertI;
00153
00154
00155 label nVisited = 0;
00156
00157 if (vertI == -1)
00158 {
00159
00160 vertI = mesh().edges()[edgeI].start();
00161
00162 if (features.isFeatureVertex(faceI, vertI))
00163 {
00164 nVisited++;
00165 }
00166 }
00167
00168 if ((edgeI == -1) || !meshTools::edgeOnFace(mesh(), faceI, edgeI))
00169 {
00170
00171
00172 edgeI = getFirstVertEdge(faceI, vertI);
00173 }
00174
00175
00176
00177 do
00178 {
00179
00180 edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
00181
00182 if (nVisited == nFeaturePts)
00183 {
00184 break;
00185 }
00186
00187 vertI = mesh().edges()[edgeI].otherVertex(vertI);
00188
00189 if (features.isFeatureVertex(faceI, vertI))
00190 {
00191 nVisited++;
00192 }
00193 }
00194 while (true);
00195 }
00196
00197
00198
00199
00200
00201 Foam::labelList Foam::topoCellLooper::getSuperEdge
00202 (
00203 const cellFeatures& features,
00204 const label faceI,
00205 const label startEdgeI,
00206 const label startVertI
00207 ) const
00208 {
00209 const labelList& fEdges = mesh().faceEdges()[faceI];
00210
00211 labelList superVerts(fEdges.size());
00212 label superVertI = 0;
00213
00214
00215 label edgeI = startEdgeI;
00216
00217 label vertI = startVertI;
00218
00219 superVerts[superVertI++] = vertI;
00220
00221 label prevEdgeI = -1;
00222
00223 do
00224 {
00225 vertI = mesh().edges()[edgeI].otherVertex(vertI);
00226
00227 superVerts[superVertI++] = vertI;
00228
00229 prevEdgeI = edgeI;
00230
00231 edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
00232 }
00233 while (!features.isFeaturePoint(prevEdgeI, edgeI));
00234
00235 superVerts.setSize(superVertI);
00236
00237 return superVerts;
00238 }
00239
00240
00241
00242
00243 Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge
00244 (
00245 const vector& refDir,
00246 const label cellI,
00247 const cellFeatures& features
00248 ) const
00249 {
00250 const labelList& cEdges = mesh().cellEdges()[cellI];
00251
00252 const point& ctr = mesh().cellCentres()[cellI];
00253
00254 label cutEdgeI = -1;
00255 scalar maxCos = -GREAT;
00256
00257 forAll(cEdges, cEdgeI)
00258 {
00259 label edgeI = cEdges[cEdgeI];
00260
00261 if (!features.isFeatureEdge(edgeI))
00262 {
00263 const edge& e = mesh().edges()[edgeI];
00264
00265
00266 vector e0 = mesh().points()[e.start()] - ctr;
00267 vector e1 = mesh().points()[e.end()] - ctr;
00268
00269 vector n = e0 ^ e1;
00270 n /= mag(n);
00271
00272 scalar cosAngle = mag(refDir & n);
00273
00274 if (cosAngle > maxCos)
00275 {
00276 maxCos = cosAngle;
00277
00278 cutEdgeI = edgeI;
00279 }
00280 }
00281 }
00282
00283 return cutEdgeI;
00284 }
00285
00286
00287
00288
00289
00290 void Foam::topoCellLooper::walkAcrossFace
00291 (
00292 const cellFeatures& features,
00293 const label faceI,
00294 const label startEdgeI,
00295 const label startVertI,
00296 const label nFeats,
00297
00298 label& edgeI,
00299 label& vertI
00300 ) const
00301 {
00302 label oppositeVertI = -1;
00303 label oppositeEdgeI = -1;
00304
00305
00306 walkFace
00307 (
00308 features,
00309 faceI,
00310 startEdgeI,
00311 startVertI,
00312 nFeats,
00313
00314 oppositeEdgeI,
00315 oppositeVertI
00316 );
00317
00318
00319
00320 labelList superEdge =
00321 getSuperEdge
00322 (
00323 features,
00324 faceI,
00325 oppositeEdgeI,
00326 oppositeVertI
00327 );
00328
00329 label sz = superEdge.size();
00330
00331 if (sz == 2)
00332 {
00333
00334
00335
00336 vertI = -1;
00337 edgeI = oppositeEdgeI;
00338 }
00339 else if (sz == 3)
00340 {
00341 vertI = superEdge[1];
00342 edgeI = -1;
00343 }
00344 else
00345 {
00346
00347 label index = sz/2;
00348
00349 if (debug)
00350 {
00351 Pout<< " Don't know what to do. Stepped to non-feature point "
00352 << "at index " << index << " in superEdge:" << superEdge
00353 << endl;
00354 }
00355
00356 vertI = superEdge[index];
00357 edgeI = -1;
00358 }
00359 }
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 void Foam::topoCellLooper::walkSplitHex
00377 (
00378 const label cellI,
00379 const cellFeatures& features,
00380 const label fromFaceI,
00381 const label fromEdgeI,
00382 const label fromVertI,
00383
00384 DynamicList<label>& loop,
00385 DynamicList<scalar>& loopWeights
00386 ) const
00387 {
00388
00389 label faceI = fromFaceI;
00390 label edgeI = fromEdgeI;
00391 label vertI = fromVertI;
00392
00393 do
00394 {
00395 if (debug)
00396 {
00397 Pout<< "Entering walk with : cell:" << cellI << " face:" << faceI;
00398 if (faceI != -1)
00399 {
00400 Pout<< " verts:" << mesh().faces()[faceI];
00401 }
00402 Pout<< " edge:" << edgeI;
00403 if (edgeI != -1)
00404 {
00405 Pout<< " verts:" << mesh().edges()[edgeI];
00406 }
00407 Pout<< " vert:" << vertI << endl;
00408 }
00409
00410 label startLoop = -1;
00411
00412 if
00413 (
00414 (vertI != -1)
00415 && (
00416 (startLoop =
00417 findIndex
00418 (
00419 loop,
00420 vertToEVert(vertI)
00421 )
00422 )
00423 != -1
00424 )
00425 )
00426 {
00427
00428 label firstFree = loop.size();
00429
00430 subsetList(startLoop, firstFree, loop);
00431 subsetList(startLoop, firstFree, loopWeights);
00432
00433 break;
00434 }
00435 if
00436 (
00437 (edgeI != -1)
00438 && (
00439 (startLoop =
00440 findIndex
00441 (
00442 loop,
00443 edgeToEVert(edgeI)
00444 )
00445 )
00446 != -1
00447 )
00448 )
00449 {
00450
00451 label firstFree = loop.size();
00452
00453 subsetList(startLoop, firstFree, loop);
00454 subsetList(startLoop, firstFree, loopWeights);
00455
00456 break;
00457 }
00458
00459
00460 if (vertI == -1)
00461 {
00462
00463 if (edgeI == -1)
00464 {
00465 FatalErrorIn("walkSplitHex") << "Illegal edge and vert"
00466 << abort(FatalError);
00467 }
00468
00469 loop.append(edgeToEVert(edgeI));
00470 loopWeights.append(0.5);
00471
00472
00473 faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
00474
00475 if (debug)
00476 {
00477 Pout<< " stepped across edge " << mesh().edges()[edgeI]
00478 << " to face " << faceI << " verts:"
00479 << mesh().faces()[faceI] << endl;
00480 }
00481
00482 label nextEdgeI = -1;
00483 label nextVertI = -1;
00484
00485
00486 walkAcrossFace
00487 (
00488 features,
00489 faceI,
00490 edgeI,
00491 vertI,
00492 2,
00493
00494 nextEdgeI,
00495 nextVertI
00496 );
00497
00498 edgeI = nextEdgeI;
00499 vertI = nextVertI;
00500 }
00501 else
00502 {
00503
00504
00505 loop.append(vertToEVert(vertI));
00506 loopWeights.append(-GREAT);
00507
00508 if (edgeI == -1)
00509 {
00510
00511
00512 labelList nextEdges = getVertEdgesNonFace
00513 (
00514 cellI,
00515 faceI,
00516 vertI
00517 );
00518
00519 if (nextEdges.empty())
00520 {
00521
00522 const labelList& pFaces = mesh().pointFaces()[vertI];
00523
00524 forAll(pFaces, pFaceI)
00525 {
00526 label thisFaceI = pFaces[pFaceI];
00527
00528 if
00529 (
00530 (thisFaceI != faceI)
00531 && meshTools::faceOnCell(mesh(), cellI, thisFaceI)
00532 )
00533 {
00534 faceI = thisFaceI;
00535 break;
00536 }
00537 }
00538
00539 if (debug)
00540 {
00541 Pout<< " stepped from non-edge vertex " << vertI
00542 << " to face " << faceI << " verts:"
00543 << mesh().faces()[faceI]
00544 << " since candidate edges:" << nextEdges << endl;
00545 }
00546
00547 label nextEdgeI = -1;
00548 label nextVertI = -1;
00549
00550 walkAcrossFace
00551 (
00552 features,
00553 faceI,
00554 edgeI,
00555 vertI,
00556 2,
00557
00558 nextEdgeI,
00559 nextVertI
00560 );
00561
00562 edgeI = nextEdgeI;
00563 vertI = nextVertI;
00564 }
00565 else if (nextEdges.size() == 1)
00566 {
00567
00568 edgeI = nextEdges[0];
00569
00570 if (debug)
00571 {
00572 Pout<< " stepped from non-edge vertex " << vertI
00573 << " along edge " << edgeI << " verts:"
00574 << mesh().edges()[edgeI]
00575 << " out of candidate edges:"
00576 << nextEdges << endl;
00577 }
00578
00579 vertI = mesh().edges()[edgeI].otherVertex(vertI);
00580
00581 faceI = -1;
00582 }
00583 else
00584 {
00585
00586
00587
00588 label index = nextEdges.size()/2;
00589
00590 edgeI = nextEdges[index];
00591
00592 if (debug)
00593 {
00594 Pout<< " stepped from non-edge vertex " << vertI
00595 << " along edge " << edgeI << " verts:"
00596 << mesh().edges()[edgeI]
00597 << " out of candidate edges:" << nextEdges << endl;
00598 }
00599
00600 vertI = mesh().edges()[edgeI].otherVertex(vertI);
00601
00602 faceI = -1;
00603 }
00604 }
00605 else
00606 {
00607
00608 labelList nextFaces =
00609 getVertFacesNonEdge
00610 (
00611 cellI,
00612 edgeI,
00613 vertI
00614 );
00615
00616 if (nextFaces.size() == 1)
00617 {
00618
00619 faceI = nextFaces[0];
00620
00621 label nextEdgeI = -1;
00622 label nextVertI = -1;
00623
00624 walkAcrossFace
00625 (
00626 features,
00627 faceI,
00628 edgeI,
00629 vertI,
00630 2,
00631
00632 nextEdgeI,
00633 nextVertI
00634 );
00635
00636 edgeI = nextEdgeI;
00637 vertI = nextVertI;
00638 }
00639 else if (nextFaces.size() == 2)
00640 {
00641
00642 faceI = -1;
00643
00644 edgeI =
00645 meshTools::getSharedEdge
00646 (
00647 mesh(),
00648 nextFaces[0],
00649 nextFaces[1]
00650 );
00651
00652 vertI = mesh().edges()[edgeI].otherVertex(vertI);
00653 }
00654 else
00655 {
00656 FatalErrorIn("walkFromVert") << "Not yet implemented"
00657 << "Choosing from more than "
00658 << "two candidates:" << nextFaces
00659 << " when coming from vertex " << vertI << " on cell "
00660 << cellI << abort(FatalError);
00661 }
00662 }
00663 }
00664
00665 if (debug)
00666 {
00667 Pout<< "Walked to : face:" << faceI;
00668 if (faceI != -1)
00669 {
00670 Pout<< " verts:" << mesh().faces()[faceI];
00671 }
00672 Pout<< " edge:" << edgeI;
00673 if (edgeI != -1)
00674 {
00675 Pout<< " verts:" << mesh().edges()[edgeI];
00676 }
00677 Pout<< " vert:" << vertI << endl;
00678 }
00679 }
00680 while (true);
00681 }
00682
00683
00684
00685
00686
00687 Foam::topoCellLooper::topoCellLooper(const polyMesh& mesh)
00688 :
00689 hexCellLooper(mesh)
00690 {}
00691
00692
00693
00694
00695 Foam::topoCellLooper::~topoCellLooper()
00696 {}
00697
00698
00699
00700
00701 bool Foam::topoCellLooper::cut
00702 (
00703 const vector& refDir,
00704 const label cellI,
00705 const boolList& vertIsCut,
00706 const boolList& edgeIsCut,
00707 const scalarField& edgeWeight,
00708
00709 labelList& loop,
00710 scalarField& loopWeights
00711 ) const
00712 {
00713 if (mesh().cellShapes()[cellI].model() == hex_)
00714 {
00715
00716 return
00717 hexCellLooper::cut
00718 (
00719 refDir,
00720 cellI,
00721 vertIsCut,
00722 edgeIsCut,
00723 edgeWeight,
00724 loop,
00725 loopWeights
00726 );
00727 }
00728 else
00729 {
00730 cellFeatures superCell(mesh(), featureCos, cellI);
00731
00732 if (hexMatcher().isA(superCell.faces()))
00733 {
00734 label edgeI =
00735 getAlignedNonFeatureEdge
00736 (
00737 refDir,
00738 cellI,
00739 superCell
00740 );
00741
00742 label vertI = -1;
00743
00744 label faceI = -1;
00745
00746 if (edgeI != -1)
00747 {
00748
00749 vertI = mesh().edges()[edgeI].start();
00750 }
00751 else
00752 {
00753
00754
00755 edgeI = getMisAlignedEdge(refDir, cellI);
00756
00757
00758 label face0;
00759 label face1;
00760 meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
00761
00762 faceI = face0;
00763 }
00764
00765 label nEstCuts = 2*mesh().cells()[cellI].size();
00766
00767 DynamicList<label> localLoop(nEstCuts);
00768 DynamicList<scalar> localLoopWeights(nEstCuts);
00769
00770 walkSplitHex
00771 (
00772 cellI,
00773 superCell,
00774 faceI,
00775 edgeI,
00776 vertI,
00777
00778 localLoop,
00779 localLoopWeights
00780 );
00781
00782 if (localLoop.size() <=2)
00783 {
00784 return false;
00785 }
00786 else
00787 {
00788 loop.transfer(localLoop);
00789 loopWeights.transfer(localLoopWeights);
00790
00791 return true;
00792 }
00793 }
00794 else
00795 {
00796
00797 return hexCellLooper::cut
00798 (
00799 refDir,
00800 cellI,
00801 vertIsCut,
00802 edgeIsCut,
00803 edgeWeight,
00804 loop,
00805 loopWeights
00806 );
00807 }
00808 }
00809 }
00810
00811
00812 bool Foam::topoCellLooper::cut
00813 (
00814 const plane& cutPlane,
00815 const label cellI,
00816 const boolList& vertIsCut,
00817 const boolList& edgeIsCut,
00818 const scalarField& edgeWeight,
00819
00820 labelList& loop,
00821 scalarField& loopWeights
00822 ) const
00823 {
00824
00825 return
00826 hexCellLooper::cut
00827 (
00828 cutPlane,
00829 cellI,
00830 vertIsCut,
00831 edgeIsCut,
00832 edgeWeight,
00833
00834 loop,
00835 loopWeights
00836 );
00837 }
00838
00839
00840