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 #include "polyDualMesh.H"
00030 #include <meshTools/meshTools.H>
00031 #include <OpenFOAM/OFstream.H>
00032 #include <OpenFOAM/Time.H>
00033 #include <OpenFOAM/SortableList.H>
00034 #include <meshTools/pointSet.H>
00035
00036
00037
00038 namespace Foam
00039 {
00040 defineTypeNameAndDebug(polyDualMesh, 0);
00041 }
00042
00043
00044
00045
00046
00047
00048 Foam::labelList Foam::polyDualMesh::getFaceOrder
00049 (
00050 const labelList& faceOwner,
00051 const labelList& faceNeighbour,
00052 const cellList& cells,
00053 label& nInternalFaces
00054 )
00055 {
00056 labelList oldToNew(faceOwner.size(), -1);
00057
00058
00059 label newFaceI = 0;
00060
00061 forAll(cells, cellI)
00062 {
00063 const labelList& cFaces = cells[cellI];
00064
00065 SortableList<label> nbr(cFaces.size());
00066
00067 forAll(cFaces, i)
00068 {
00069 label faceI = cFaces[i];
00070
00071 label nbrCellI = faceNeighbour[faceI];
00072
00073 if (nbrCellI != -1)
00074 {
00075
00076 if (nbrCellI == cellI)
00077 {
00078 nbrCellI = faceOwner[faceI];
00079 }
00080
00081 if (cellI < nbrCellI)
00082 {
00083
00084 nbr[i] = nbrCellI;
00085 }
00086 else
00087 {
00088
00089 nbr[i] = -1;
00090 }
00091 }
00092 else
00093 {
00094
00095 nbr[i] = -1;
00096 }
00097 }
00098
00099 nbr.sort();
00100
00101 forAll(nbr, i)
00102 {
00103 if (nbr[i] != -1)
00104 {
00105 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
00106 }
00107 }
00108 }
00109
00110 nInternalFaces = newFaceI;
00111
00112 Pout<< "nInternalFaces:" << nInternalFaces << endl;
00113 Pout<< "nFaces:" << faceOwner.size() << endl;
00114 Pout<< "nCells:" << cells.size() << endl;
00115
00116
00117 for (label faceI = newFaceI; faceI < faceOwner.size(); faceI++)
00118 {
00119 oldToNew[faceI] = faceI;
00120 }
00121
00122
00123
00124 forAll(oldToNew, faceI)
00125 {
00126 if (oldToNew[faceI] == -1)
00127 {
00128 FatalErrorIn
00129 (
00130 "polyDualMesh::getFaceOrder"
00131 "(const labelList&, const labelList&, const label) const"
00132 ) << "Did not determine new position"
00133 << " for face " << faceI
00134 << abort(FatalError);
00135 }
00136 }
00137
00138 return oldToNew;
00139 }
00140
00141
00142
00143
00144
00145
00146
00147 void Foam::polyDualMesh::getPointEdges
00148 (
00149 const primitivePatch& patch,
00150 const label faceI,
00151 const label pointI,
00152 label& e0,
00153 label& e1
00154 )
00155 {
00156 const labelList& fEdges = patch.faceEdges()[faceI];
00157 const face& f = patch.localFaces()[faceI];
00158
00159 e0 = -1;
00160 e1 = -1;
00161
00162 forAll(fEdges, i)
00163 {
00164 label edgeI = fEdges[i];
00165
00166 const edge& e = patch.edges()[edgeI];
00167
00168 if (e[0] == pointI)
00169 {
00170
00171 label index = findIndex(f, pointI);
00172
00173 if (f.nextLabel(index) == e[1])
00174 {
00175 e1 = edgeI;
00176 }
00177 else
00178 {
00179 e0 = edgeI;
00180 }
00181
00182 if (e0 != -1 && e1 != -1)
00183 {
00184 return;
00185 }
00186 }
00187 else if (e[1] == pointI)
00188 {
00189
00190 label index = findIndex(f, pointI);
00191
00192 if (f.nextLabel(index) == e[0])
00193 {
00194 e1 = edgeI;
00195 }
00196 else
00197 {
00198 e0 = edgeI;
00199 }
00200
00201 if (e0 != -1 && e1 != -1)
00202 {
00203 return;
00204 }
00205 }
00206 }
00207
00208 FatalErrorIn("getPointEdges") << "Cannot find two edges on face:" << faceI
00209 << " vertices:" << patch.localFaces()[faceI]
00210 << " that uses point:" << pointI
00211 << abort(FatalError);
00212 }
00213
00214
00215
00216 Foam::labelList Foam::polyDualMesh::collectPatchSideFace
00217 (
00218 const polyPatch& patch,
00219 const label patchToDualOffset,
00220 const labelList& edgeToDualPoint,
00221 const labelList& pointToDualPoint,
00222 const label pointI,
00223
00224 label& edgeI
00225 )
00226 {
00227
00228
00229
00230 label meshPointI = patch.meshPoints()[pointI];
00231 const labelList& pFaces = patch.pointFaces()[pointI];
00232
00233 DynamicList<label> dualFace;
00234
00235 if (pointToDualPoint[meshPointI] >= 0)
00236 {
00237
00238 dualFace.setCapacity(pFaces.size()+2+1);
00239
00240 dualFace.append(pointToDualPoint[meshPointI]);
00241 }
00242 else
00243 {
00244 dualFace.setCapacity(pFaces.size()+2);
00245 }
00246
00247
00248 if (edgeToDualPoint[patch.meshEdges()[edgeI]] < 0)
00249 {
00250 FatalErrorIn("polyDualMesh::collectPatchSideFace") << edgeI
00251 << abort(FatalError);
00252 }
00253
00254 dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
00255
00256 label faceI = patch.edgeFaces()[edgeI][0];
00257
00258
00259 bool reverseFace;
00260
00261 label e0, e1;
00262 getPointEdges(patch, faceI, pointI, e0, e1);
00263
00264 if (e0 == edgeI)
00265 {
00266 reverseFace = true;
00267 }
00268 else
00269 {
00270 reverseFace = false;
00271 }
00272
00273 while (true)
00274 {
00275
00276 dualFace.append(faceI + patchToDualOffset);
00277
00278
00279 label e0, e1;
00280 getPointEdges(patch, faceI, pointI, e0, e1);
00281
00282 if (e0 == edgeI)
00283 {
00284 edgeI = e1;
00285 }
00286 else
00287 {
00288 edgeI = e0;
00289 }
00290
00291 if (edgeToDualPoint[patch.meshEdges()[edgeI]] >= 0)
00292 {
00293
00294 dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
00295 }
00296
00297 const labelList& eFaces = patch.edgeFaces()[edgeI];
00298
00299 if (eFaces.size() == 1)
00300 {
00301
00302 break;
00303 }
00304
00305
00306 if (eFaces[0] == faceI)
00307 {
00308 faceI = eFaces[1];
00309 }
00310 else
00311 {
00312 faceI = eFaces[0];
00313 }
00314 }
00315
00316 dualFace.shrink();
00317
00318 if (reverseFace)
00319 {
00320 reverse(dualFace);
00321 }
00322
00323 return dualFace;
00324 }
00325
00326
00327
00328
00329
00330 void Foam::polyDualMesh::collectPatchInternalFace
00331 (
00332 const polyPatch& patch,
00333 const label patchToDualOffset,
00334 const labelList& edgeToDualPoint,
00335
00336 const label pointI,
00337 const label startEdgeI,
00338
00339 labelList& dualFace2,
00340 labelList& featEdgeIndices2
00341 )
00342 {
00343
00344 const labelList& meshEdges = patch.meshEdges();
00345 const labelList& pFaces = patch.pointFaces()[pointI];
00346
00347
00348 DynamicList<label> dualFace(pFaces.size());
00349
00350 DynamicList<label> featEdgeIndices(dualFace.size());
00351
00352
00353 label edgeI = startEdgeI;
00354 label faceI = patch.edgeFaces()[edgeI][0];
00355
00356
00357 bool reverseFace;
00358
00359 label e0, e1;
00360 getPointEdges(patch, faceI, pointI, e0, e1);
00361
00362 if (e0 == edgeI)
00363 {
00364 reverseFace = true;
00365 }
00366 else
00367 {
00368 reverseFace = false;
00369 }
00370
00371 while (true)
00372 {
00373
00374 dualFace.append(faceI + patchToDualOffset);
00375
00376
00377 label e0, e1;
00378 getPointEdges(patch, faceI, pointI, e0, e1);
00379
00380 if (e0 == edgeI)
00381 {
00382 edgeI = e1;
00383 }
00384 else
00385 {
00386 edgeI = e0;
00387 }
00388
00389 if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
00390 {
00391
00392 dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
00393 featEdgeIndices.append(dualFace.size()-1);
00394 }
00395
00396 if (edgeI == startEdgeI)
00397 {
00398 break;
00399 }
00400
00401
00402 const labelList& eFaces = patch.edgeFaces()[edgeI];
00403
00404 if (eFaces[0] == faceI)
00405 {
00406 faceI = eFaces[1];
00407 }
00408 else
00409 {
00410 faceI = eFaces[0];
00411 }
00412 }
00413
00414 dualFace2.transfer(dualFace);
00415
00416 featEdgeIndices2.transfer(featEdgeIndices);
00417
00418 if (reverseFace)
00419 {
00420 reverse(dualFace2);
00421
00422
00423 forAll(featEdgeIndices2, i)
00424 {
00425 featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
00426 }
00427
00428 reverse(featEdgeIndices2);
00429 }
00430 }
00431
00432
00433 void Foam::polyDualMesh::splitFace
00434 (
00435 const polyPatch& patch,
00436 const labelList& pointToDualPoint,
00437
00438 const label pointI,
00439 const labelList& dualFace,
00440 const labelList& featEdgeIndices,
00441
00442 DynamicList<face>& dualFaces,
00443 DynamicList<label>& dualOwner,
00444 DynamicList<label>& dualNeighbour,
00445 DynamicList<label>& dualRegion
00446 )
00447 {
00448
00449
00450
00451 label meshPointI = patch.meshPoints()[pointI];
00452
00453 if (pointToDualPoint[meshPointI] >= 0)
00454 {
00455
00456
00457 if (featEdgeIndices.size() < 2)
00458 {
00459
00460 dualFaces.append(face(dualFace));
00461 dualOwner.append(meshPointI);
00462 dualNeighbour.append(-1);
00463 dualRegion.append(patch.index());
00464 }
00465 else
00466 {
00467
00468
00469
00470 forAll(featEdgeIndices, i)
00471 {
00472 label startFp = featEdgeIndices[i];
00473
00474 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
00475
00476
00477 label sz = 0;
00478
00479 if (endFp > startFp)
00480 {
00481 sz = endFp - startFp + 2;
00482 }
00483 else
00484 {
00485 sz = endFp + dualFace.size() - startFp + 2;
00486 }
00487 face subFace(sz);
00488
00489
00490 subFace[0] = pointToDualPoint[patch.meshPoints()[pointI]];
00491
00492
00493 for (label subFp = 1; subFp < subFace.size(); subFp++)
00494 {
00495 subFace[subFp] = dualFace[startFp];
00496
00497 startFp = (startFp+1) % dualFace.size();
00498 }
00499
00500 dualFaces.append(face(subFace));
00501 dualOwner.append(meshPointI);
00502 dualNeighbour.append(-1);
00503 dualRegion.append(patch.index());
00504 }
00505 }
00506 }
00507 else
00508 {
00509
00510 if (featEdgeIndices.size() < 2)
00511 {
00512
00513 dualFaces.append(face(dualFace));
00514 dualOwner.append(meshPointI);
00515 dualNeighbour.append(-1);
00516 dualRegion.append(patch.index());
00517 }
00518 else
00519 {
00520
00521
00522
00523
00524
00525 DynamicList<label> subFace(dualFace.size());
00526
00527 forAll(featEdgeIndices, featI)
00528 {
00529 label startFp = featEdgeIndices[featI];
00530 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
00531
00532 label fp = startFp;
00533
00534 while (true)
00535 {
00536 label vertI = dualFace[fp];
00537
00538 subFace.append(vertI);
00539
00540 if (fp == endFp)
00541 {
00542 break;
00543 }
00544
00545 fp = dualFace.fcIndex(fp);
00546 }
00547
00548 if (subFace.size() > 2)
00549 {
00550
00551 subFace.shrink();
00552
00553 dualFaces.append(face(subFace));
00554 dualOwner.append(meshPointI);
00555 dualNeighbour.append(-1);
00556 dualRegion.append(patch.index());
00557
00558 subFace.clear();
00559 }
00560 }
00561
00562 if (subFace.size() > 2)
00563 {
00564
00565 subFace.shrink();
00566
00567 dualFaces.append(face(subFace));
00568 dualOwner.append(meshPointI);
00569 dualNeighbour.append(-1);
00570 dualRegion.append(patch.index());
00571
00572 subFace.clear();
00573 }
00574 }
00575 }
00576 }
00577
00578
00579
00580 void Foam::polyDualMesh::dualPatch
00581 (
00582 const polyPatch& patch,
00583 const label patchToDualOffset,
00584 const labelList& edgeToDualPoint,
00585 const labelList& pointToDualPoint,
00586
00587 const pointField& dualPoints,
00588
00589 DynamicList<face>& dualFaces,
00590 DynamicList<label>& dualOwner,
00591 DynamicList<label>& dualNeighbour,
00592 DynamicList<label>& dualRegion
00593 )
00594 {
00595 const labelList& meshEdges = patch.meshEdges();
00596
00597
00598
00599
00600
00601
00602 labelList doneEdgeSide(meshEdges.size(), 0);
00603
00604 boolList donePoint(patch.nPoints(), false);
00605
00606
00607
00608
00609
00610 forAll(doneEdgeSide, patchEdgeI)
00611 {
00612 const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
00613
00614 if (eFaces.size() == 1)
00615 {
00616 const edge& e = patch.edges()[patchEdgeI];
00617
00618 forAll(e, eI)
00619 {
00620 label bitMask = 1<<eI;
00621
00622 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
00623 {
00624
00625
00626
00627 label pointI = e[eI];
00628
00629 label edgeI = patchEdgeI;
00630 labelList dualFace
00631 (
00632 collectPatchSideFace
00633 (
00634 patch,
00635 patchToDualOffset,
00636 edgeToDualPoint,
00637 pointToDualPoint,
00638
00639 pointI,
00640 edgeI
00641 )
00642 );
00643
00644
00645 if (patch.edges()[edgeI][0] == pointI)
00646 {
00647 doneEdgeSide[edgeI] |= 1;
00648 }
00649 else
00650 {
00651 doneEdgeSide[edgeI] |= 2;
00652 }
00653
00654 dualFaces.append(face(dualFace));
00655 dualOwner.append(patch.meshPoints()[pointI]);
00656 dualNeighbour.append(-1);
00657 dualRegion.append(patch.index());
00658
00659 doneEdgeSide[patchEdgeI] |= bitMask;
00660 donePoint[pointI] = true;
00661 }
00662 }
00663 }
00664 }
00665
00666
00667
00668
00669
00670
00671 forAll(donePoint, pointI)
00672 {
00673 if (!donePoint[pointI])
00674 {
00675 labelList dualFace, featEdgeIndices;
00676
00677 collectPatchInternalFace
00678 (
00679 patch,
00680 patchToDualOffset,
00681 edgeToDualPoint,
00682 pointI,
00683 patch.pointEdges()[pointI][0],
00684
00685 dualFace,
00686 featEdgeIndices
00687 );
00688
00689
00690
00692
00693
00694
00695
00696
00697 splitFace
00698 (
00699 patch,
00700 pointToDualPoint,
00701 pointI,
00702 dualFace,
00703 featEdgeIndices,
00704
00705 dualFaces,
00706 dualOwner,
00707 dualNeighbour,
00708 dualRegion
00709 );
00710
00711 donePoint[pointI] = true;
00712 }
00713 }
00714 }
00715
00716
00717 void Foam::polyDualMesh::calcDual
00718 (
00719 const polyMesh& mesh,
00720 const labelList& featureEdges,
00721 const labelList& featurePoints
00722 )
00723 {
00724 const polyBoundaryMesh& patches = mesh.boundaryMesh();
00725 const label nIntFaces = mesh.nInternalFaces();
00726
00727
00728
00729 primitivePatch allBoundary
00730 (
00731 SubList<face>
00732 (
00733 mesh.faces(),
00734 mesh.nFaces() - nIntFaces,
00735 nIntFaces
00736 ),
00737 mesh.points()
00738 );
00739
00740
00741 labelList meshEdges
00742 (
00743 allBoundary.meshEdges
00744 (
00745 mesh.edges(),
00746 mesh.pointEdges()
00747 )
00748 );
00749
00750 {
00751 pointSet nonManifoldPoints
00752 (
00753 mesh,
00754 "nonManifoldPoints",
00755 meshEdges.size() / 100
00756 );
00757
00758 allBoundary.checkPointManifold(true, &nonManifoldPoints);
00759
00760 if (nonManifoldPoints.size())
00761 {
00762 nonManifoldPoints.write();
00763
00764 FatalErrorIn
00765 (
00766 "polyDualMesh::calcDual(const polyMesh&, const labelList&"
00767 ", const labelList&)"
00768 ) << "There are " << nonManifoldPoints.size() << " points where"
00769 << " the outside of the mesh is non-manifold." << nl
00770 << "Such a mesh cannot be converted to a dual." << nl
00771 << "Writing points at non-manifold sites to pointSet "
00772 << nonManifoldPoints.name()
00773 << exit(FatalError);
00774 }
00775 }
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788 pointField dualPoints
00789 (
00790 mesh.nCells()
00791 + mesh.nFaces() - nIntFaces
00792 + featureEdges.size()
00793 + featurePoints.size()
00794 );
00795
00796 label dualPointI = 0;
00797
00798
00799
00800 const pointField& cellCentres = mesh.cellCentres();
00801
00802 cellPoint_.setSize(cellCentres.size());
00803
00804 forAll(cellCentres, cellI)
00805 {
00806 cellPoint_[cellI] = dualPointI;
00807 dualPoints[dualPointI++] = cellCentres[cellI];
00808 }
00809
00810
00811 const pointField& faceCentres = mesh.faceCentres();
00812
00813 boundaryFacePoint_.setSize(mesh.nFaces() - nIntFaces);
00814
00815 for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
00816 {
00817 boundaryFacePoint_[faceI - nIntFaces] = dualPointI;
00818 dualPoints[dualPointI++] = faceCentres[faceI];
00819 }
00820
00821
00822
00823
00824
00825 labelList edgeToDualPoint(mesh.nEdges(), -2);
00826
00827 forAll(meshEdges, patchEdgeI)
00828 {
00829 label edgeI = meshEdges[patchEdgeI];
00830
00831 edgeToDualPoint[edgeI] = -1;
00832 }
00833
00834 forAll(featureEdges, i)
00835 {
00836 label edgeI = featureEdges[i];
00837
00838 const edge& e = mesh.edges()[edgeI];
00839
00840 edgeToDualPoint[edgeI] = dualPointI;
00841 dualPoints[dualPointI++] = e.centre(mesh.points());
00842 }
00843
00844
00845
00846
00847
00848
00849
00850
00851 labelList pointToDualPoint(mesh.nPoints(), -3);
00852
00853 forAll(patches, patchI)
00854 {
00855 const labelList& meshPoints = patches[patchI].meshPoints();
00856
00857 forAll(meshPoints, i)
00858 {
00859 pointToDualPoint[meshPoints[i]] = -2;
00860 }
00861
00862 const labelListList& loops = patches[patchI].edgeLoops();
00863
00864 forAll(loops, i)
00865 {
00866 const labelList& loop = loops[i];
00867
00868 forAll(loop, j)
00869 {
00870 pointToDualPoint[meshPoints[loop[j]]] = -1;
00871 }
00872 }
00873 }
00874
00875 forAll(featurePoints, i)
00876 {
00877 label pointI = featurePoints[i];
00878
00879 pointToDualPoint[pointI] = dualPointI;
00880 dualPoints[dualPointI++] = mesh.points()[pointI];
00881 }
00882
00883
00884
00885
00886
00887 DynamicList<face> dynDualFaces(mesh.nEdges());
00888 DynamicList<label> dynDualOwner(mesh.nEdges());
00889 DynamicList<label> dynDualNeighbour(mesh.nEdges());
00890 DynamicList<label> dynDualRegion(mesh.nEdges());
00891
00892
00893
00894
00895
00896 forAll(meshEdges, patchEdgeI)
00897 {
00898 label edgeI = meshEdges[patchEdgeI];
00899
00900 const edge& e = mesh.edges()[edgeI];
00901
00902 label owner = -1;
00903 label neighbour = -1;
00904
00905 if (e[0] < e[1])
00906 {
00907 owner = e[0];
00908 neighbour = e[1];
00909 }
00910 else
00911 {
00912 owner = e[1];
00913 neighbour = e[0];
00914 }
00915
00916
00917 const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
00918
00919 if (patchFaces.size() != 2)
00920 {
00921 FatalErrorIn("polyDualMesh::calcDual")
00922 << "Cannot handle edges with " << patchFaces.size()
00923 << " connected boundary faces."
00924 << abort(FatalError);
00925 }
00926
00927 label face0 = patchFaces[0] + nIntFaces;
00928 const face& f0 = mesh.faces()[face0];
00929
00930 label face1 = patchFaces[1] + nIntFaces;
00931
00932
00933
00934
00935 label startFaceI = -1;
00936 label endFaceI = -1;
00937
00938 label index = findIndex(f0, neighbour);
00939
00940 if (f0.nextLabel(index) == owner)
00941 {
00942 startFaceI = face0;
00943 endFaceI = face1;
00944 }
00945 else
00946 {
00947 startFaceI = face1;
00948 endFaceI = face0;
00949 }
00950
00951
00952
00953 DynamicList<label> dualFace;
00954
00955 if (edgeToDualPoint[edgeI] >= 0)
00956 {
00957
00958 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
00959
00960 dualFace.append(edgeToDualPoint[edgeI]);
00961 }
00962 else
00963 {
00964 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
00965 }
00966
00967
00968 dualFace.append(mesh.nCells() + startFaceI - nIntFaces);
00969
00970 label cellI = mesh.faceOwner()[startFaceI];
00971 label faceI = startFaceI;
00972
00973 while (true)
00974 {
00975
00976 dualFace.append(cellI);
00977
00978
00979 label f0, f1;
00980 meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
00981
00982 if (f0 == faceI)
00983 {
00984 faceI = f1;
00985 }
00986 else
00987 {
00988 faceI = f0;
00989 }
00990
00991
00992 if (faceI == endFaceI)
00993 {
00994 break;
00995 }
00996
00997 if (mesh.faceOwner()[faceI] == cellI)
00998 {
00999 cellI = mesh.faceNeighbour()[faceI];
01000 }
01001 else
01002 {
01003 cellI = mesh.faceOwner()[faceI];
01004 }
01005 }
01006
01007
01008 dualFace.append(mesh.nCells() + endFaceI - nIntFaces);
01009
01010 dynDualFaces.append(face(dualFace.shrink()));
01011 dynDualOwner.append(owner);
01012 dynDualNeighbour.append(neighbour);
01013 dynDualRegion.append(-1);
01014
01015 {
01016
01017 const face& f = dynDualFaces[dynDualFaces.size()-1];
01018 vector n = f.normal(dualPoints);
01019 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
01020 {
01021 WarningIn("calcDual") << "Incorrect orientation"
01022 << " on boundary edge:" << edgeI
01023 << mesh.points()[mesh.edges()[edgeI][0]]
01024 << mesh.points()[mesh.edges()[edgeI][1]]
01025 << endl;
01026 }
01027 }
01028 }
01029
01030
01031
01032
01033
01034 forAll(edgeToDualPoint, edgeI)
01035 {
01036 if (edgeToDualPoint[edgeI] == -2)
01037 {
01038
01039
01040
01041
01042 const edge& e = mesh.edges()[edgeI];
01043
01044 label owner = -1;
01045 label neighbour = -1;
01046
01047 if (e[0] < e[1])
01048 {
01049 owner = e[0];
01050 neighbour = e[1];
01051 }
01052 else
01053 {
01054 owner = e[1];
01055 neighbour = e[0];
01056 }
01057
01058
01059 const labelList& eCells = mesh.edgeCells()[edgeI];
01060
01061 label cellI = eCells[0];
01062
01063
01064 label face0, face1;
01065 meshTools::getEdgeFaces(mesh, cellI, edgeI, face0, face1);
01066
01067
01068
01069 const face& f0 = mesh.faces()[face0];
01070
01071 label index = findIndex(f0, neighbour);
01072
01073 bool f0OrderOk = (f0.nextLabel(index) == owner);
01074
01075 label startFaceI = -1;
01076
01077 if (f0OrderOk == (mesh.faceOwner()[face0] == cellI))
01078 {
01079 startFaceI = face0;
01080 }
01081 else
01082 {
01083 startFaceI = face1;
01084 }
01085
01086
01087 DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
01088
01089 label faceI = startFaceI;
01090
01091 while (true)
01092 {
01093
01094 dualFace.append(cellI);
01095
01096
01097 label f0, f1;
01098 meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
01099
01100 if (f0 == faceI)
01101 {
01102 faceI = f1;
01103 }
01104 else
01105 {
01106 faceI = f0;
01107 }
01108
01109
01110 if (faceI == startFaceI)
01111 {
01112 break;
01113 }
01114
01115 if (mesh.faceOwner()[faceI] == cellI)
01116 {
01117 cellI = mesh.faceNeighbour()[faceI];
01118 }
01119 else
01120 {
01121 cellI = mesh.faceOwner()[faceI];
01122 }
01123 }
01124
01125 dynDualFaces.append(face(dualFace.shrink()));
01126 dynDualOwner.append(owner);
01127 dynDualNeighbour.append(neighbour);
01128 dynDualRegion.append(-1);
01129
01130 {
01131
01132 const face& f = dynDualFaces[dynDualFaces.size()-1];
01133 vector n = f.normal(dualPoints);
01134 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
01135 {
01136 WarningIn("calcDual") << "Incorrect orientation"
01137 << " on internal edge:" << edgeI
01138 << mesh.points()[mesh.edges()[edgeI][0]]
01139 << mesh.points()[mesh.edges()[edgeI][1]]
01140 << endl;
01141 }
01142 }
01143 }
01144 }
01145
01146
01147 if (debug)
01148 {
01149 dynDualFaces.shrink();
01150 dynDualOwner.shrink();
01151 dynDualNeighbour.shrink();
01152 dynDualRegion.shrink();
01153
01154 OFstream str("dualInternalFaces.obj");
01155 Pout<< "polyDualMesh::calcDual : dumping internal faces to "
01156 << str.name() << endl;
01157
01158 forAll(dualPoints, dualPointI)
01159 {
01160 meshTools::writeOBJ(str, dualPoints[dualPointI]);
01161 }
01162 forAll(dynDualFaces, dualFaceI)
01163 {
01164 const face& f = dynDualFaces[dualFaceI];
01165
01166 str<< 'f';
01167 forAll(f, fp)
01168 {
01169 str<< ' ' << f[fp]+1;
01170 }
01171 str<< nl;
01172 }
01173 }
01174
01175 const label nInternalFaces = dynDualFaces.size();
01176
01177
01178
01179
01180 forAll(patches, patchI)
01181 {
01182 const polyPatch& pp = patches[patchI];
01183
01184 dualPatch
01185 (
01186 pp,
01187 mesh.nCells() + pp.start() - nIntFaces,
01188 edgeToDualPoint,
01189 pointToDualPoint,
01190
01191 dualPoints,
01192
01193 dynDualFaces,
01194 dynDualOwner,
01195 dynDualNeighbour,
01196 dynDualRegion
01197 );
01198 }
01199
01200
01201
01202
01203 faceList dualFaces(dynDualFaces.shrink(), true);
01204 dynDualFaces.clear();
01205
01206 labelList dualOwner(dynDualOwner.shrink(), true);
01207 dynDualOwner.clear();
01208
01209 labelList dualNeighbour(dynDualNeighbour.shrink(), true);
01210 dynDualNeighbour.clear();
01211
01212 labelList dualRegion(dynDualRegion.shrink(), true);
01213 dynDualRegion.clear();
01214
01215
01216
01217
01218 if (debug)
01219 {
01220 OFstream str("dualFaces.obj");
01221 Pout<< "polyDualMesh::calcDual : dumping all faces to "
01222 << str.name() << endl;
01223
01224 forAll(dualPoints, dualPointI)
01225 {
01226 meshTools::writeOBJ(str, dualPoints[dualPointI]);
01227 }
01228 forAll(dualFaces, dualFaceI)
01229 {
01230 const face& f = dualFaces[dualFaceI];
01231
01232 str<< 'f';
01233 forAll(f, fp)
01234 {
01235 str<< ' ' << f[fp]+1;
01236 }
01237 str<< nl;
01238 }
01239 }
01240
01241
01242 cellList dualCells(mesh.nPoints());
01243
01244 forAll(dualCells, cellI)
01245 {
01246 dualCells[cellI].setSize(0);
01247 }
01248
01249 forAll(dualOwner, faceI)
01250 {
01251 label cellI = dualOwner[faceI];
01252
01253 labelList& cFaces = dualCells[cellI];
01254
01255 label sz = cFaces.size();
01256 cFaces.setSize(sz+1);
01257 cFaces[sz] = faceI;
01258 }
01259 forAll(dualNeighbour, faceI)
01260 {
01261 label cellI = dualNeighbour[faceI];
01262
01263 if (cellI != -1)
01264 {
01265 labelList& cFaces = dualCells[cellI];
01266
01267 label sz = cFaces.size();
01268 cFaces.setSize(sz+1);
01269 cFaces[sz] = faceI;
01270 }
01271 }
01272
01273
01274
01275
01276 label dummy;
01277
01278 labelList oldToNew
01279 (
01280 getFaceOrder
01281 (
01282 dualOwner,
01283 dualNeighbour,
01284 dualCells,
01285 dummy
01286 )
01287 );
01288
01289
01290 inplaceReorder(oldToNew, dualFaces);
01291 inplaceReorder(oldToNew, dualOwner);
01292 inplaceReorder(oldToNew, dualNeighbour);
01293 inplaceReorder(oldToNew, dualRegion);
01294 forAll(dualCells, cellI)
01295 {
01296 inplaceRenumber(oldToNew, dualCells[cellI]);
01297 }
01298
01299
01300
01301 labelList patchSizes(patches.size(), 0);
01302
01303 forAll(dualRegion, faceI)
01304 {
01305 if (dualRegion[faceI] >= 0)
01306 {
01307 patchSizes[dualRegion[faceI]]++;
01308 }
01309 }
01310
01311 labelList patchStarts(patches.size(), 0);
01312
01313 label faceI = nInternalFaces;
01314
01315 forAll(patches, patchI)
01316 {
01317 patchStarts[patchI] = faceI;
01318
01319 faceI += patchSizes[patchI];
01320 }
01321
01322
01323 Pout<< "nFaces:" << dualFaces.size()
01324 << " patchSizes:" << patchSizes
01325 << " patchStarts:" << patchStarts
01326 << endl;
01327
01328
01329
01330 List<polyPatch*> dualPatches(patches.size());
01331
01332 forAll(patches, patchI)
01333 {
01334 const polyPatch& pp = patches[patchI];
01335
01336 dualPatches[patchI] = pp.clone
01337 (
01338 boundaryMesh(),
01339 patchI,
01340 0,
01341 0
01342 ).ptr();
01343 }
01344 addPatches(dualPatches);
01345
01346
01347 resetPrimitives
01348 (
01349 xferMove(dualPoints),
01350 xferMove(dualFaces),
01351 xferMove(dualOwner),
01352 xferMove(dualNeighbour),
01353 patchSizes,
01354 patchStarts
01355 );
01356 }
01357
01358
01359
01360
01361
01362 Foam::polyDualMesh::polyDualMesh(const IOobject& io)
01363 :
01364 polyMesh(io),
01365 cellPoint_
01366 (
01367 IOobject
01368 (
01369 "cellPoint",
01370 time().findInstance(meshDir(), "cellPoint"),
01371 meshSubDir,
01372 *this,
01373 IOobject::MUST_READ,
01374 IOobject::NO_WRITE
01375 )
01376 ),
01377 boundaryFacePoint_
01378 (
01379 IOobject
01380 (
01381 "boundaryFacePoint",
01382 time().findInstance(meshDir(), "boundaryFacePoint"),
01383 meshSubDir,
01384 *this,
01385 IOobject::MUST_READ,
01386 IOobject::NO_WRITE
01387 )
01388 )
01389 {}
01390
01391
01392
01393 Foam::polyDualMesh::polyDualMesh
01394 (
01395 const polyMesh& mesh,
01396 const labelList& featureEdges,
01397 const labelList& featurePoints
01398 )
01399 :
01400 polyMesh
01401 (
01402 mesh,
01403 xferCopy(pointField()),
01404 xferCopy(faceList()),
01405 xferCopy(cellList())
01406 ),
01407 cellPoint_
01408 (
01409 IOobject
01410 (
01411 "cellPoint",
01412 time().findInstance(meshDir(), "faces"),
01413 meshSubDir,
01414 *this,
01415 IOobject::NO_READ,
01416 IOobject::AUTO_WRITE
01417 ),
01418 labelList(mesh.nCells(), -1)
01419 ),
01420 boundaryFacePoint_
01421 (
01422 IOobject
01423 (
01424 "boundaryFacePoint",
01425 time().findInstance(meshDir(), "faces"),
01426 meshSubDir,
01427 *this,
01428 IOobject::NO_READ,
01429 IOobject::AUTO_WRITE
01430 ),
01431 labelList(mesh.nFaces() - mesh.nInternalFaces())
01432 )
01433 {
01434 calcDual(mesh, featureEdges, featurePoints);
01435 }
01436
01437
01438
01439 Foam::polyDualMesh::polyDualMesh
01440 (
01441 const polyMesh& mesh,
01442 const scalar featureCos
01443 )
01444 :
01445 polyMesh
01446 (
01447 mesh,
01448 xferCopy(pointField()),
01449 xferCopy(faceList()),
01450 xferCopy(cellList())
01451 ),
01452 cellPoint_
01453 (
01454 IOobject
01455 (
01456 "cellPoint",
01457 time().findInstance(meshDir(), "faces"),
01458 meshSubDir,
01459 *this,
01460 IOobject::NO_READ,
01461 IOobject::AUTO_WRITE
01462 ),
01463 labelList(mesh.nCells(), -1)
01464 ),
01465 boundaryFacePoint_
01466 (
01467 IOobject
01468 (
01469 "boundaryFacePoint",
01470 time().findInstance(meshDir(), "faces"),
01471 meshSubDir,
01472 *this,
01473 IOobject::NO_READ,
01474 IOobject::AUTO_WRITE
01475 ),
01476 labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
01477 )
01478 {
01479 labelList featureEdges, featurePoints;
01480
01481 calcFeatures(mesh, featureCos, featureEdges, featurePoints);
01482 calcDual(mesh, featureEdges, featurePoints);
01483 }
01484
01485
01486 void Foam::polyDualMesh::calcFeatures
01487 (
01488 const polyMesh& mesh,
01489 const scalar featureCos,
01490 labelList& featureEdges,
01491 labelList& featurePoints
01492 )
01493 {
01494
01495 primitivePatch allBoundary
01496 (
01497 SubList<face>
01498 (
01499 mesh.faces(),
01500 mesh.nFaces() - mesh.nInternalFaces(),
01501 mesh.nInternalFaces()
01502 ),
01503 mesh.points()
01504 );
01505
01506
01507 labelList allRegion(allBoundary.size());
01508
01509 const polyBoundaryMesh& patches = mesh.boundaryMesh();
01510
01511 forAll(patches, patchI)
01512 {
01513 const polyPatch& pp = patches[patchI];
01514
01515 forAll(pp, i)
01516 {
01517 allRegion[i + pp.start() - mesh.nInternalFaces()] = patchI;
01518 }
01519 }
01520
01521
01522
01523
01524
01525 const labelListList& edgeFaces = allBoundary.edgeFaces();
01526 const vectorField& faceNormals = allBoundary.faceNormals();
01527 const labelList& meshPoints = allBoundary.meshPoints();
01528
01529 boolList isFeatureEdge(edgeFaces.size(), false);
01530
01531 forAll(edgeFaces, edgeI)
01532 {
01533 const labelList& eFaces = edgeFaces[edgeI];
01534
01535 if (eFaces.size() != 2)
01536 {
01537
01538 const edge& e = allBoundary.edges()[edgeI];
01539
01540 WarningIn("polyDualMesh::calcFeatures") << "Edge "
01541 << meshPoints[e[0]] << ' ' << meshPoints[e[1]]
01542 << " coords:" << mesh.points()[meshPoints[e[0]]]
01543 << mesh.points()[meshPoints[e[1]]]
01544 << " has more than 2 faces connected to it:"
01545 << eFaces.size() << endl;
01546
01547 isFeatureEdge[edgeI] = true;
01548 }
01549 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
01550 {
01551 isFeatureEdge[edgeI] = true;
01552 }
01553 else if
01554 (
01555 (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
01556 < featureCos
01557 )
01558 {
01559 isFeatureEdge[edgeI] = true;
01560 }
01561 }
01562
01563
01564
01565
01566
01567 const labelListList& pointEdges = allBoundary.pointEdges();
01568
01569 DynamicList<label> allFeaturePoints(pointEdges.size());
01570
01571 forAll(pointEdges, pointI)
01572 {
01573 const labelList& pEdges = pointEdges[pointI];
01574
01575 label nFeatEdges = 0;
01576
01577 forAll(pEdges, i)
01578 {
01579 if (isFeatureEdge[pEdges[i]])
01580 {
01581 nFeatEdges++;
01582 }
01583 }
01584 if (nFeatEdges > 2)
01585 {
01586
01587 allFeaturePoints.append(allBoundary.meshPoints()[pointI]);
01588 }
01589 }
01590 featurePoints.transfer(allFeaturePoints);
01591
01592 if (debug)
01593 {
01594 OFstream str("featurePoints.obj");
01595
01596 Pout<< "polyDualMesh::calcFeatures : dumping feature points to "
01597 << str.name() << endl;
01598
01599 forAll(featurePoints, i)
01600 {
01601 label pointI = featurePoints[i];
01602 meshTools::writeOBJ(str, mesh.points()[pointI]);
01603 }
01604 }
01605
01606
01607
01608 labelList meshEdges
01609 (
01610 allBoundary.meshEdges
01611 (
01612 mesh.edges(),
01613 mesh.cellEdges(),
01614 SubList<label>
01615 (
01616 mesh.faceOwner(),
01617 allBoundary.size(),
01618 mesh.nInternalFaces()
01619 )
01620 )
01621 );
01622
01623 DynamicList<label> allFeatureEdges(isFeatureEdge.size());
01624 forAll(isFeatureEdge, edgeI)
01625 {
01626 if (isFeatureEdge[edgeI])
01627 {
01628
01629 allFeatureEdges.append(meshEdges[edgeI]);
01630 }
01631 }
01632 featureEdges.transfer(allFeatureEdges);
01633 }
01634
01635
01636
01637
01638 Foam::polyDualMesh::~polyDualMesh()
01639 {}
01640
01641
01642