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 #include "boundaryCutter.H"
00029 #include <OpenFOAM/polyMesh.H>
00030 #include <dynamicMesh/polyTopoChange.H>
00031 #include <dynamicMesh/polyAddCell.H>
00032 #include <dynamicMesh/polyAddFace.H>
00033 #include <dynamicMesh/polyAddPoint.H>
00034 #include <dynamicMesh/polyModifyFace.H>
00035 #include <dynamicMesh/polyModifyPoint.H>
00036 #include <OpenFOAM/mapPolyMesh.H>
00037 #include <meshTools/meshTools.H>
00038
00039
00040
00041 namespace Foam
00042 {
00043
00044 defineTypeNameAndDebug(boundaryCutter, 0);
00045
00046 }
00047
00048
00049
00050
00051
00052
00053
00054 void Foam::boundaryCutter::getFaceInfo
00055 (
00056 const label faceI,
00057 label& patchID,
00058 label& zoneID,
00059 label& zoneFlip
00060 ) const
00061 {
00062 patchID = -1;
00063
00064 if (!mesh_.isInternalFace(faceI))
00065 {
00066 patchID = mesh_.boundaryMesh().whichPatch(faceI);
00067 }
00068
00069 zoneID = mesh_.faceZones().whichZone(faceI);
00070
00071 zoneFlip = false;
00072
00073 if (zoneID >= 0)
00074 {
00075 const faceZone& fZone = mesh_.faceZones()[zoneID];
00076
00077 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
00078 }
00079 }
00080
00081
00082
00083
00084 Foam::face Foam::boundaryCutter::addEdgeCutsToFace
00085 (
00086 const label faceI,
00087 const Map<labelList>& edgeToAddedPoints
00088 ) const
00089 {
00090 const edgeList& edges = mesh_.edges();
00091 const face& f = mesh_.faces()[faceI];
00092 const labelList& fEdges = mesh_.faceEdges()[faceI];
00093
00094
00095 DynamicList<label> newFace(2 * f.size());
00096
00097 forAll(f, fp)
00098 {
00099
00100 newFace.append(f[fp]);
00101
00102
00103 label v1 = f.nextLabel(fp);
00104
00105 label edgeI = meshTools::findEdge(edges, fEdges, f[fp], v1);
00106
00107 Map<labelList>::const_iterator fnd = edgeToAddedPoints.find(edgeI);
00108
00109 if (fnd != edgeToAddedPoints.end())
00110 {
00111
00112 const labelList& addedPoints = fnd();
00113
00114 if (edges[edgeI].start() == f[fp])
00115 {
00116
00117 forAll(addedPoints, i)
00118 {
00119 newFace.append(addedPoints[i]);
00120 }
00121 }
00122 else
00123 {
00124
00125 forAllReverse(addedPoints, i)
00126 {
00127 newFace.append(addedPoints[i]);
00128 }
00129 }
00130 }
00131 }
00132
00133 face returnFace;
00134 returnFace.transfer(newFace);
00135
00136 if (debug)
00137 {
00138 Pout<< "addEdgeCutsToFace:" << nl
00139 << " from : " << f << nl
00140 << " to : " << returnFace << endl;
00141 }
00142
00143 return returnFace;
00144 }
00145
00146
00147 void Foam::boundaryCutter::addFace
00148 (
00149 const label faceI,
00150 const face& newFace,
00151
00152 bool& modifiedFace,
00153 polyTopoChange& meshMod
00154 ) const
00155 {
00156
00157 label patchID, zoneID, zoneFlip;
00158 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
00159 label own = mesh_.faceOwner()[faceI];
00160 label masterPoint = mesh_.faces()[faceI][0];
00161
00162 if (!modifiedFace)
00163 {
00164 meshMod.setAction
00165 (
00166 polyModifyFace
00167 (
00168 newFace,
00169 faceI,
00170 own,
00171 -1,
00172 false,
00173 patchID,
00174 false,
00175 zoneID,
00176 zoneFlip
00177 )
00178 );
00179
00180 modifiedFace = true;
00181 }
00182 else
00183 {
00184 meshMod.setAction
00185 (
00186 polyAddFace
00187 (
00188 newFace,
00189 own,
00190 -1,
00191 masterPoint,
00192 -1,
00193 -1,
00194 false,
00195 patchID,
00196 zoneID,
00197 zoneFlip
00198 )
00199 );
00200 }
00201 }
00202
00203
00204
00205
00206 bool Foam::boundaryCutter::splitFace
00207 (
00208 const label faceI,
00209 const Map<point>& pointToPos,
00210 const Map<labelList>& edgeToAddedPoints,
00211 polyTopoChange& meshMod
00212 ) const
00213 {
00214 const edgeList& edges = mesh_.edges();
00215 const face& f = mesh_.faces()[faceI];
00216 const labelList& fEdges = mesh_.faceEdges()[faceI];
00217
00218
00219 label nSplitEdges = 0;
00220 label nModPoints = 0;
00221 label nTotalSplits = 0;
00222
00223 forAll(f, fp)
00224 {
00225 if (pointToPos.found(f[fp]))
00226 {
00227 nModPoints++;
00228 nTotalSplits++;
00229 }
00230
00231
00232 label nextV = f.nextLabel(fp);
00233
00234 label edgeI = meshTools::findEdge(edges, fEdges, f[fp], nextV);
00235
00236 Map<labelList>::const_iterator fnd = edgeToAddedPoints.find(edgeI);
00237
00238 if (fnd != edgeToAddedPoints.end())
00239 {
00240 nSplitEdges++;
00241 nTotalSplits += fnd().size();
00242 }
00243 }
00244
00245 if (debug)
00246 {
00247 Pout<< "Face:" << faceI
00248 << " nModPoints:" << nModPoints
00249 << " nSplitEdges:" << nSplitEdges
00250 << " nTotalSplits:" << nTotalSplits << endl;
00251 }
00252
00253 if (nSplitEdges == 0 && nModPoints == 0)
00254 {
00255 FatalErrorIn("boundaryCutter::splitFace") << "Problem : face:" << faceI
00256 << " nSplitEdges:" << nSplitEdges
00257 << " nTotalSplits:" << nTotalSplits
00258 << abort(FatalError);
00259 return false;
00260 }
00261 else if (nSplitEdges + nModPoints == 1)
00262 {
00263
00264
00265 Warning << "Face " << faceI << " has only one edge cut " << endl;
00266 return false;
00267 }
00268 else
00269 {
00270
00271
00272
00273
00274
00275
00276 label patchID, zoneID, zoneFlip;
00277 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
00278
00279
00280 face extendedFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
00281
00282
00283 label startFp = -1;
00284
00285 forAll(extendedFace, fp)
00286 {
00287 if (extendedFace[fp] >= mesh_.nPoints())
00288 {
00289 startFp = fp;
00290 break;
00291 }
00292 }
00293
00294 if (startFp == -1)
00295 {
00296
00297 forAll(extendedFace, fp)
00298 {
00299 if (pointToPos.found(extendedFace[fp]))
00300 {
00301 startFp = fp;
00302 break;
00303 }
00304 }
00305 }
00306
00307 if (startFp == -1)
00308 {
00309 FatalErrorIn("boundaryCutter::splitFace")
00310 << "Problem" << abort(FatalError);
00311 }
00312
00313
00314
00315 bool modifiedFace = false;
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 DynamicList<label> newFace(extendedFace.size());
00340
00341 label fp = startFp;
00342
00343 forAll(extendedFace, i)
00344 {
00345 label pointI = extendedFace[fp];
00346
00347 newFace.append(pointI);
00348
00349 if
00350 (
00351 newFace.size() > 2
00352 && (
00353 pointI >= mesh_.nPoints()
00354 || pointToPos.found(pointI)
00355 )
00356 )
00357 {
00358
00359 face tmpFace;
00360 tmpFace.transfer(newFace);
00361
00362
00363 addFace(faceI, tmpFace, modifiedFace, meshMod);
00364
00365
00366 newFace.append(extendedFace[startFp]);
00367 newFace.append(extendedFace[fp]);
00368 }
00369
00370 fp = (fp+1) % extendedFace.size();
00371 }
00372
00373
00374 if (newFace.size() > 2)
00375 {
00376
00377 face tmpFace;
00378 tmpFace.transfer(newFace);
00379
00380
00381 addFace(faceI, tmpFace, modifiedFace, meshMod);
00382 }
00383
00384
00385 return true;
00386 }
00387 }
00388
00389
00390
00391
00392
00393 Foam::boundaryCutter::boundaryCutter(const polyMesh& mesh)
00394 :
00395 mesh_(mesh),
00396 edgeAddedPoints_(),
00397 faceAddedPoint_()
00398 {}
00399
00400
00401
00402
00403 Foam::boundaryCutter::~boundaryCutter()
00404 {}
00405
00406
00407
00408
00409 void Foam::boundaryCutter::setRefinement
00410 (
00411 const Map<point>& pointToPos,
00412 const Map<List<point> >& edgeToCuts,
00413 const Map<labelPair>& faceToSplit,
00414 const Map<point>& faceToFeaturePoint,
00415 polyTopoChange& meshMod
00416 )
00417 {
00418
00419 edgeAddedPoints_.clear();
00420
00421 faceAddedPoint_.clear();
00422 faceAddedPoint_.resize(faceToFeaturePoint.size());
00423
00424
00425
00426
00427
00428
00429
00430 forAllConstIter(Map<point>, pointToPos, iter)
00431 {
00432 meshMod.setAction
00433 (
00434 polyModifyPoint
00435 (
00436 iter.key(),
00437 iter(),
00438 false,
00439 -1,
00440 true
00441 )
00442 );
00443 }
00444
00445
00446
00447
00448
00449
00450
00451 Map<labelList> edgeToAddedPoints(edgeToCuts.size());
00452
00453 forAllConstIter(Map<List<point> >, edgeToCuts, iter)
00454 {
00455 label edgeI = iter.key();
00456
00457 const edge& e = mesh_.edges()[edgeI];
00458
00459
00460 const List<point>& cuts = iter();
00461
00462 forAll(cuts, cutI)
00463 {
00464
00465 const point& featurePoint = cuts[cutI];
00466
00467 label addedPointI =
00468 meshMod.setAction
00469 (
00470 polyAddPoint
00471 (
00472 featurePoint,
00473 e.start(),
00474 -1,
00475 true
00476 )
00477 );
00478
00479 Map<labelList>::iterator fnd = edgeToAddedPoints.find(edgeI);
00480
00481 if (fnd != edgeToAddedPoints.end())
00482 {
00483 labelList& addedPoints = fnd();
00484
00485 label sz = addedPoints.size();
00486 addedPoints.setSize(sz+1);
00487 addedPoints[sz] = addedPointI;
00488 }
00489 else
00490 {
00491 edgeToAddedPoints.insert(edgeI, labelList(1, addedPointI));
00492 }
00493
00494 if (debug)
00495 {
00496 Pout<< "Added point " << addedPointI << " for edge " << edgeI
00497 << " with cuts:" << edgeToAddedPoints[edgeI] << endl;
00498 }
00499 }
00500 }
00501
00502
00503
00504
00505
00506
00507 forAllConstIter(Map<point>, faceToFeaturePoint, iter)
00508 {
00509 label faceI = iter.key();
00510
00511 const face& f = mesh_.faces()[faceI];
00512
00513 if (faceToSplit.found(faceI))
00514 {
00515 FatalErrorIn("boundaryCutter::setRefinement")
00516 << "Face " << faceI << " vertices " << f
00517 << " is both marked for face-centre decomposition and"
00518 << " diagonal splitting."
00519 << abort(FatalError);
00520 }
00521
00522 if (mesh_.isInternalFace(faceI))
00523 {
00524 FatalErrorIn("boundaryCutter::setRefinement")
00525 << "Face " << faceI << " vertices " << f
00526 << " is not an external face. Cannot split it"
00527 << abort(FatalError);
00528 }
00529
00530 label addedPointI =
00531 meshMod.setAction
00532 (
00533 polyAddPoint
00534 (
00535 iter(),
00536 f[0],
00537 -1,
00538 true
00539 )
00540 );
00541 faceAddedPoint_.insert(faceI, addedPointI);
00542
00543 if (debug)
00544 {
00545 Pout<< "Added point " << addedPointI << " for feature point "
00546 << iter() << " on face " << faceI << " with centre "
00547 << mesh_.faceCentres()[faceI] << endl;
00548 }
00549 }
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 boolList faceUptodate(mesh_.nFaces(), false);
00560
00561
00562
00563 forAllConstIter(Map<label>, faceAddedPoint_, iter)
00564 {
00565 label faceI = iter.key();
00566
00567
00568 face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
00569
00570 label addedPointI = iter();
00571
00572
00573 label patchID, zoneID, zoneFlip;
00574 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
00575 label own = mesh_.faceOwner()[faceI];
00576 label masterPoint = mesh_.faces()[faceI][0];
00577
00578
00579
00580 face tri(3);
00581
00582 forAll(newFace, fp)
00583 {
00584 label nextV = newFace.nextLabel(fp);
00585
00586 tri[0] = newFace[fp];
00587 tri[1] = nextV;
00588 tri[2] = addedPointI;
00589
00590 if (fp == 0)
00591 {
00592
00593 meshMod.setAction
00594 (
00595 polyModifyFace
00596 (
00597 tri,
00598 faceI,
00599 own,
00600 -1,
00601 false,
00602 patchID,
00603 false,
00604 zoneID,
00605 zoneFlip
00606 )
00607 );
00608 }
00609 else
00610 {
00611
00612 meshMod.setAction
00613 (
00614 polyAddFace
00615 (
00616 tri,
00617 own,
00618 -1,
00619 masterPoint,
00620 -1,
00621 -1,
00622 false,
00623 patchID,
00624 zoneID,
00625 zoneFlip
00626 )
00627 );
00628 }
00629 }
00630
00631 faceUptodate[faceI] = true;
00632 }
00633
00634
00635
00636 forAllConstIter(Map<labelPair>, faceToSplit, iter)
00637 {
00638 label faceI = iter.key();
00639
00640 const face& f = mesh_.faces()[faceI];
00641
00642 if (faceAddedPoint_.found(faceI))
00643 {
00644 FatalErrorIn("boundaryCutter::setRefinement")
00645 << "Face " << faceI << " vertices " << f
00646 << " is both marked for face-centre decomposition and"
00647 << " diagonal splitting."
00648 << abort(FatalError);
00649 }
00650
00651
00652
00653 face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
00654
00655
00656 label patchID, zoneID, zoneFlip;
00657 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
00658 label own = mesh_.faceOwner()[faceI];
00659 label masterPoint = mesh_.faces()[faceI][0];
00660
00661
00662 const labelPair& diag = iter();
00663
00664 label fp0 = findIndex(newFace, f[diag[0]]);
00665 label fp1 = findIndex(newFace, f[diag[1]]);
00666
00667 if (fp0 == -1 || fp1 == -1 || fp0 == fp1)
00668 {
00669 FatalErrorIn("boundaryCutter::setRefinement")
00670 << "Problem : Face " << faceI << " vertices " << f
00671 << " newFace:" << newFace << " diagonal:" << f[diag[0]]
00672 << ' ' << f[diag[1]]
00673 << abort(FatalError);
00674 }
00675
00676
00677
00678
00679 DynamicList<label> newVerts(newFace.size());
00680
00681
00682 label fp = fp0;
00683
00684 do
00685 {
00686 newVerts.append(newFace[fp]);
00687
00688 fp = (fp == newFace.size()-1 ? 0 : fp+1);
00689
00690 } while (fp != fp1);
00691
00692 newVerts.append(newFace[fp1]);
00693
00694
00695
00696 meshMod.setAction
00697 (
00698 polyModifyFace
00699 (
00700 face(newVerts.shrink()),
00701 faceI,
00702 own,
00703 -1,
00704 false,
00705 patchID,
00706 false,
00707 zoneID,
00708 zoneFlip
00709 )
00710 );
00711
00712
00713 newVerts.clear();
00714
00715
00716
00717 do
00718 {
00719 newVerts.append(newFace[fp]);
00720
00721 fp = (fp == newFace.size()-1 ? 0 : fp+1);
00722
00723 } while (fp != fp0);
00724
00725 newVerts.append(newFace[fp0]);
00726
00727
00728 meshMod.setAction
00729 (
00730 polyAddFace
00731 (
00732 face(newVerts.shrink()),
00733 own,
00734 -1,
00735 masterPoint,
00736 -1,
00737 -1,
00738 false,
00739 patchID,
00740 zoneID,
00741 zoneFlip
00742 )
00743 );
00744
00745 faceUptodate[faceI] = true;
00746 }
00747
00748
00749
00750
00751 forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
00752 {
00753 label edgeI = iter.key();
00754
00755 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
00756
00757 forAll(eFaces, i)
00758 {
00759 label faceI = eFaces[i];
00760
00761 if (!faceUptodate[faceI] && !mesh_.isInternalFace(faceI))
00762 {
00763
00764 if (splitFace(faceI, pointToPos, edgeToAddedPoints, meshMod))
00765 {
00766
00767 faceUptodate[faceI] = true;
00768 }
00769 }
00770 }
00771 }
00772
00773
00774
00775
00776
00777 forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
00778 {
00779 label edgeI = iter.key();
00780
00781 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
00782
00783 forAll(eFaces, i)
00784 {
00785 label faceI = eFaces[i];
00786
00787 if (!faceUptodate[faceI])
00788 {
00789
00790 face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
00791
00792 label own = mesh_.faceOwner()[faceI];
00793
00794 label nei = -1;
00795
00796 if (mesh_.isInternalFace(faceI))
00797 {
00798 nei = mesh_.faceNeighbour()[faceI];
00799 }
00800
00801 label patchID, zoneID, zoneFlip;
00802 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
00803
00804 meshMod.setAction
00805 (
00806 polyModifyFace
00807 (
00808 newFace,
00809 faceI,
00810 own,
00811 nei,
00812 false,
00813 patchID,
00814 false,
00815 zoneID,
00816 zoneFlip
00817 )
00818 );
00819
00820 faceUptodate[faceI] = true;
00821 }
00822 }
00823 }
00824
00825
00826
00827 edgeAddedPoints_.resize(edgeToCuts.size());
00828
00829 forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
00830 {
00831 edgeAddedPoints_.insert(mesh_.edges()[iter.key()], iter());
00832 }
00833 }
00834
00835
00836 void Foam::boundaryCutter::updateMesh(const mapPolyMesh& morphMap)
00837 {
00838
00839
00840
00841
00842
00843
00844 {
00845
00846 Map<label> newAddedPoints(faceAddedPoint_.size());
00847
00848 forAllConstIter(Map<label>, faceAddedPoint_, iter)
00849 {
00850 label oldFaceI = iter.key();
00851
00852 label newFaceI = morphMap.reverseFaceMap()[oldFaceI];
00853
00854 label oldPointI = iter();
00855
00856 label newPointI = morphMap.reversePointMap()[oldPointI];
00857
00858 if (newFaceI >= 0 && newPointI >= 0)
00859 {
00860 newAddedPoints.insert(newFaceI, newPointI);
00861 }
00862 }
00863
00864
00865 faceAddedPoint_.transfer(newAddedPoints);
00866 }
00867
00868
00869
00870
00871
00872
00873
00874 {
00875
00876 HashTable<labelList, edge, Hash<edge> >
00877 newEdgeAddedPoints(edgeAddedPoints_.size());
00878
00879 for
00880 (
00881 HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
00882 edgeAddedPoints_.begin();
00883 iter != edgeAddedPoints_.end();
00884 ++iter
00885 )
00886 {
00887 const edge& e = iter.key();
00888
00889 label newStart = morphMap.reversePointMap()[e.start()];
00890
00891 label newEnd = morphMap.reversePointMap()[e.end()];
00892
00893 if (newStart >= 0 && newEnd >= 0)
00894 {
00895 const labelList& addedPoints = iter();
00896
00897 labelList newAddedPoints(addedPoints.size());
00898 label newI = 0;
00899
00900 forAll(addedPoints, i)
00901 {
00902 label newAddedPointI =
00903 morphMap.reversePointMap()[addedPoints[i]];
00904
00905 if (newAddedPointI >= 0)
00906 {
00907 newAddedPoints[newI++] = newAddedPointI;
00908 }
00909 }
00910 if (newI > 0)
00911 {
00912 newAddedPoints.setSize(newI);
00913
00914 edge newE = edge(newStart, newEnd);
00915
00916 newEdgeAddedPoints.insert(newE, newAddedPoints);
00917 }
00918 }
00919 }
00920
00921
00922 edgeAddedPoints_.transfer(newEdgeAddedPoints);
00923 }
00924 }
00925
00926
00927