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 "meshTools.H"
00027 #include <OpenFOAM/polyMesh.H>
00028 #include <OpenFOAM/hexMatcher.H>
00029
00030
00031
00032
00033 bool Foam::meshTools::visNormal
00034 (
00035 const vector& n,
00036 const vectorField& faceNormals,
00037 const labelList& faceLabels
00038 )
00039 {
00040 forAll(faceLabels, i)
00041 {
00042 if ((faceNormals[faceLabels[i]] & n) < SMALL)
00043 {
00044
00045 return false;
00046 }
00047 }
00048
00049 return true;
00050 }
00051
00052
00053 Foam::vectorField Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)
00054 {
00055 vectorField octantNormal(8);
00056 octantNormal[mXmYmZ] = vector(-1, -1, -1);
00057 octantNormal[pXmYmZ] = vector(1, -1, -1);
00058 octantNormal[mXpYmZ] = vector(-1, 1, -1);
00059 octantNormal[pXpYmZ] = vector(1, 1, -1);
00060
00061 octantNormal[mXmYpZ] = vector(-1, -1, 1);
00062 octantNormal[pXmYpZ] = vector(1, -1, 1);
00063 octantNormal[mXpYpZ] = vector(-1, 1, 1);
00064 octantNormal[pXpYpZ] = vector(1, 1, 1);
00065
00066 octantNormal /= mag(octantNormal);
00067
00068
00069 vectorField pn(pp.nPoints());
00070
00071 const vectorField& faceNormals = pp.faceNormals();
00072 const vectorField& pointNormals = pp.pointNormals();
00073 const labelListList& pointFaces = pp.pointFaces();
00074
00075 forAll(pointFaces, pointI)
00076 {
00077 const labelList& pFaces = pointFaces[pointI];
00078
00079 if (visNormal(pointNormals[pointI], faceNormals, pFaces))
00080 {
00081 pn[pointI] = pointNormals[pointI];
00082 }
00083 else
00084 {
00085 WarningIn
00086 (
00087 "Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)"
00088 ) << "Average point normal not visible for point:"
00089 << pp.meshPoints()[pointI] << endl;
00090
00091 label visOctant =
00092 mXmYmZMask
00093 | pXmYmZMask
00094 | mXpYmZMask
00095 | pXpYmZMask
00096 | mXmYpZMask
00097 | pXmYpZMask
00098 | mXpYpZMask
00099 | pXpYpZMask;
00100
00101 forAll(pFaces, i)
00102 {
00103 const vector& n = faceNormals[pFaces[i]];
00104
00105 if (n.x() > SMALL)
00106 {
00107
00108 visOctant &= ~mXmYmZMask;
00109 visOctant &= ~mXmYpZMask;
00110 visOctant &= ~mXpYmZMask;
00111 visOctant &= ~mXpYpZMask;
00112 }
00113 else if (n.x() < -SMALL)
00114 {
00115
00116 visOctant &= ~pXmYmZMask;
00117 visOctant &= ~pXmYpZMask;
00118 visOctant &= ~pXpYmZMask;
00119 visOctant &= ~pXpYpZMask;
00120 }
00121
00122 if (n.y() > SMALL)
00123 {
00124 visOctant &= ~mXmYmZMask;
00125 visOctant &= ~mXmYpZMask;
00126 visOctant &= ~pXmYmZMask;
00127 visOctant &= ~pXmYpZMask;
00128 }
00129 else if (n.y() < -SMALL)
00130 {
00131 visOctant &= ~mXpYmZMask;
00132 visOctant &= ~mXpYpZMask;
00133 visOctant &= ~pXpYmZMask;
00134 visOctant &= ~pXpYpZMask;
00135 }
00136
00137 if (n.z() > SMALL)
00138 {
00139 visOctant &= ~mXmYmZMask;
00140 visOctant &= ~mXpYmZMask;
00141 visOctant &= ~pXmYmZMask;
00142 visOctant &= ~pXpYmZMask;
00143 }
00144 else if (n.z() < -SMALL)
00145 {
00146 visOctant &= ~mXmYpZMask;
00147 visOctant &= ~mXpYpZMask;
00148 visOctant &= ~pXmYpZMask;
00149 visOctant &= ~pXpYpZMask;
00150 }
00151 }
00152
00153 label visI = -1;
00154
00155 label mask = 1;
00156
00157 for (label octant = 0; octant < 8; octant++)
00158 {
00159 if (visOctant & mask)
00160 {
00161
00162
00163 visI = octant;
00164
00165 break;
00166 }
00167 mask <<= 1;
00168 }
00169
00170 if (visI != -1)
00171 {
00172
00173 pn[pointI] = octantNormal[visI];
00174 }
00175 else
00176 {
00177 pn[pointI] = vector::zero;
00178
00179 WarningIn
00180 (
00181 "Foam::meshTools::calcBoxPointNormals"
00182 "(const primitivePatch& pp)"
00183 ) << "No visible octant for point:" << pp.meshPoints()[pointI]
00184 << " cooord:" << pp.points()[pp.meshPoints()[pointI]] << nl
00185 << "Normal set to " << pn[pointI] << endl;
00186 }
00187 }
00188 }
00189
00190 return pn;
00191 }
00192
00193
00194 Foam::vector Foam::meshTools::normEdgeVec
00195 (
00196 const primitiveMesh& mesh,
00197 const label edgeI
00198 )
00199 {
00200 vector eVec = mesh.edges()[edgeI].vec(mesh.points());
00201
00202 eVec /= mag(eVec);
00203
00204 return eVec;
00205 }
00206
00207
00208 void Foam::meshTools::writeOBJ
00209 (
00210 Ostream& os,
00211 const point& pt
00212 )
00213 {
00214 os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
00215 }
00216
00217
00218 void Foam::meshTools::writeOBJ
00219 (
00220 Ostream& os,
00221 const faceList& faces,
00222 const pointField& points,
00223 const labelList& faceLabels
00224 )
00225 {
00226 Map<label> foamToObj(4*faceLabels.size());
00227
00228 label vertI = 0;
00229
00230 forAll(faceLabels, i)
00231 {
00232 const face& f = faces[faceLabels[i]];
00233
00234 forAll(f, fp)
00235 {
00236 if (foamToObj.insert(f[fp], vertI))
00237 {
00238 writeOBJ(os, points[f[fp]]);
00239 vertI++;
00240 }
00241 }
00242
00243 os << 'l';
00244 forAll(f, fp)
00245 {
00246 os << ' ' << foamToObj[f[fp]]+1;
00247 }
00248 os << ' ' << foamToObj[f[0]]+1 << endl;
00249 }
00250 }
00251
00252
00253 void Foam::meshTools::writeOBJ
00254 (
00255 Ostream& os,
00256 const faceList& faces,
00257 const pointField& points
00258 )
00259 {
00260 labelList allFaces(faces.size());
00261 forAll(allFaces, i)
00262 {
00263 allFaces[i] = i;
00264 }
00265 writeOBJ(os, faces, points, allFaces);
00266 }
00267
00268
00269 void Foam::meshTools::writeOBJ
00270 (
00271 Ostream& os,
00272 const cellList& cells,
00273 const faceList& faces,
00274 const pointField& points,
00275 const labelList& cellLabels
00276 )
00277 {
00278 labelHashSet usedFaces(4*cellLabels.size());
00279
00280 forAll(cellLabels, i)
00281 {
00282 const cell& cFaces = cells[cellLabels[i]];
00283
00284 forAll(cFaces, j)
00285 {
00286 usedFaces.insert(cFaces[j]);
00287 }
00288 }
00289
00290 writeOBJ(os, faces, points, usedFaces.toc());
00291 }
00292
00293
00294 bool Foam::meshTools::edgeOnCell
00295 (
00296 const primitiveMesh& mesh,
00297 const label cellI,
00298 const label edgeI
00299 )
00300 {
00301 return findIndex(mesh.edgeCells(edgeI), cellI) != -1;
00302 }
00303
00304
00305 bool Foam::meshTools::edgeOnFace
00306 (
00307 const primitiveMesh& mesh,
00308 const label faceI,
00309 const label edgeI
00310 )
00311 {
00312 return findIndex(mesh.faceEdges(faceI), edgeI) != -1;
00313 }
00314
00315
00316
00317 bool Foam::meshTools::faceOnCell
00318 (
00319 const primitiveMesh& mesh,
00320 const label cellI,
00321 const label faceI
00322 )
00323 {
00324 if (mesh.isInternalFace(faceI))
00325 {
00326 if
00327 (
00328 (mesh.faceOwner()[faceI] == cellI)
00329 || (mesh.faceNeighbour()[faceI] == cellI)
00330 )
00331 {
00332 return true;
00333 }
00334 }
00335 else
00336 {
00337 if (mesh.faceOwner()[faceI] == cellI)
00338 {
00339 return true;
00340 }
00341 }
00342 return false;
00343 }
00344
00345
00346 Foam::label Foam::meshTools::findEdge
00347 (
00348 const edgeList& edges,
00349 const labelList& candidates,
00350 const label v0,
00351 const label v1
00352 )
00353 {
00354 forAll(candidates, i)
00355 {
00356 label edgeI = candidates[i];
00357
00358 const edge& e = edges[edgeI];
00359
00360 if ((e[0] == v0 && e[1] == v1) || (e[0] == v1 && e[1] == v0))
00361 {
00362 return edgeI;
00363 }
00364 }
00365 return -1;
00366 }
00367
00368
00369 Foam::label Foam::meshTools::findEdge
00370 (
00371 const primitiveMesh& mesh,
00372 const label v0,
00373 const label v1
00374 )
00375 {
00376 const edgeList& edges = mesh.edges();
00377
00378 const labelList& v0Edges = mesh.pointEdges()[v0];
00379
00380 forAll(v0Edges, i)
00381 {
00382 label edgeI = v0Edges[i];
00383
00384 const edge& e = edges[edgeI];
00385
00386 if ((e.start() == v1) || (e.end() == v1))
00387 {
00388 return edgeI;
00389 }
00390 }
00391 return -1;
00392 }
00393
00394
00395 Foam::label Foam::meshTools::getSharedEdge
00396 (
00397 const primitiveMesh& mesh,
00398 const label f0,
00399 const label f1
00400 )
00401 {
00402 const labelList& f0Edges = mesh.faceEdges()[f0];
00403 const labelList& f1Edges = mesh.faceEdges()[f1];
00404
00405 forAll(f0Edges, f0EdgeI)
00406 {
00407 label edge0 = f0Edges[f0EdgeI];
00408
00409 forAll(f1Edges, f1EdgeI)
00410 {
00411 label edge1 = f1Edges[f1EdgeI];
00412
00413 if (edge0 == edge1)
00414 {
00415 return edge0;
00416 }
00417 }
00418 }
00419 FatalErrorIn
00420 (
00421 "meshTools::getSharedEdge(const primitiveMesh&, const label"
00422 ", const label)"
00423 ) << "Faces " << f0 << " and " << f1 << " do not share an edge"
00424 << abort(FatalError);
00425
00426 return -1;
00427
00428 }
00429
00430
00431 Foam::label Foam::meshTools::getSharedFace
00432 (
00433 const primitiveMesh& mesh,
00434 const label cell0I,
00435 const label cell1I
00436 )
00437 {
00438 const cell& cFaces = mesh.cells()[cell0I];
00439
00440 forAll(cFaces, cFaceI)
00441 {
00442 label faceI = cFaces[cFaceI];
00443
00444 if
00445 (
00446 mesh.isInternalFace(faceI)
00447 && (
00448 mesh.faceOwner()[faceI] == cell1I
00449 || mesh.faceNeighbour()[faceI] == cell1I
00450 )
00451 )
00452 {
00453 return faceI;
00454 }
00455 }
00456
00457
00458 FatalErrorIn
00459 (
00460 "meshTools::getSharedFace(const primitiveMesh&, const label"
00461 ", const label)"
00462 ) << "No common face for"
00463 << " cell0I:" << cell0I << " faces:" << cFaces
00464 << " cell1I:" << cell1I << " faces:"
00465 << mesh.cells()[cell1I]
00466 << abort(FatalError);
00467
00468 return -1;
00469 }
00470
00471
00472
00473 void Foam::meshTools::getEdgeFaces
00474 (
00475 const primitiveMesh& mesh,
00476 const label cellI,
00477 const label edgeI,
00478 label& face0,
00479 label& face1
00480 )
00481 {
00482 const labelList& eFaces = mesh.edgeFaces(edgeI);
00483
00484 face0 = -1;
00485 face1 = -1;
00486
00487 forAll(eFaces, eFaceI)
00488 {
00489 label faceI = eFaces[eFaceI];
00490
00491 if (faceOnCell(mesh, cellI, faceI))
00492 {
00493 if (face0 == -1)
00494 {
00495 face0 = faceI;
00496 }
00497 else
00498 {
00499 face1 = faceI;
00500
00501 return;
00502 }
00503 }
00504 }
00505
00506 if ((face0 == -1) || (face1 == -1))
00507 {
00508 FatalErrorIn
00509 (
00510 "meshTools::getEdgeFaces(const primitiveMesh&, const label"
00511 ", const label, label&, label&"
00512 ) << "Can not find faces using edge " << mesh.edges()[edgeI]
00513 << " on cell " << cellI << abort(FatalError);
00514 }
00515 }
00516
00517
00518
00519 Foam::label Foam::meshTools::otherEdge
00520 (
00521 const primitiveMesh& mesh,
00522 const labelList& edgeLabels,
00523 const label thisEdgeI,
00524 const label thisVertI
00525 )
00526 {
00527 forAll(edgeLabels, edgeLabelI)
00528 {
00529 label edgeI = edgeLabels[edgeLabelI];
00530
00531 if (edgeI != thisEdgeI)
00532 {
00533 const edge& e = mesh.edges()[edgeI];
00534
00535 if ((e.start() == thisVertI) || (e.end() == thisVertI))
00536 {
00537 return edgeI;
00538 }
00539 }
00540 }
00541
00542 FatalErrorIn
00543 (
00544 "meshTools::otherEdge(const primitiveMesh&, const labelList&"
00545 ", const label, const label)"
00546 ) << "Can not find edge in "
00547 << UIndirectList<edge>(mesh.edges(), edgeLabels)()
00548 << " connected to edge "
00549 << thisEdgeI << " with vertices " << mesh.edges()[thisEdgeI]
00550 << " on side " << thisVertI << abort(FatalError);
00551
00552 return -1;
00553 }
00554
00555
00556
00557 Foam::label Foam::meshTools::otherFace
00558 (
00559 const primitiveMesh& mesh,
00560 const label cellI,
00561 const label faceI,
00562 const label edgeI
00563 )
00564 {
00565 label face0;
00566 label face1;
00567
00568 getEdgeFaces(mesh, cellI, edgeI, face0, face1);
00569
00570 if (face0 == faceI)
00571 {
00572 return face1;
00573 }
00574 else
00575 {
00576 return face0;
00577 }
00578 }
00579
00580
00581
00582 Foam::label Foam::meshTools::otherCell
00583 (
00584 const primitiveMesh& mesh,
00585 const label otherCellI,
00586 const label faceI
00587 )
00588 {
00589 if (!mesh.isInternalFace(faceI))
00590 {
00591 FatalErrorIn
00592 (
00593 "meshTools::otherCell(const primitiveMesh&, const label"
00594 ", const label)"
00595 ) << "Face " << faceI << " is not internal"
00596 << abort(FatalError);
00597 }
00598
00599 label newCellI = mesh.faceOwner()[faceI];
00600
00601 if (newCellI == otherCellI)
00602 {
00603 newCellI = mesh.faceNeighbour()[faceI];
00604 }
00605 return newCellI;
00606 }
00607
00608
00609
00610
00611 Foam::label Foam::meshTools::walkFace
00612 (
00613 const primitiveMesh& mesh,
00614 const label faceI,
00615 const label startEdgeI,
00616 const label startVertI,
00617 const label nEdges
00618 )
00619 {
00620 const labelList& fEdges = mesh.faceEdges(faceI);
00621
00622 label edgeI = startEdgeI;
00623
00624 label vertI = startVertI;
00625
00626 for (label iter = 0; iter < nEdges; iter++)
00627 {
00628 edgeI = otherEdge(mesh, fEdges, edgeI, vertI);
00629
00630 vertI = mesh.edges()[edgeI].otherVertex(vertI);
00631 }
00632
00633 return edgeI;
00634 }
00635
00636
00637 void Foam::meshTools::constrainToMeshCentre
00638 (
00639 const polyMesh& mesh,
00640 point& pt
00641 )
00642 {
00643 const Vector<label>& dirs = mesh.geometricD();
00644
00645 const point& min = mesh.bounds().min();
00646 const point& max = mesh.bounds().max();
00647
00648 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00649 {
00650 if (dirs[cmpt] == -1)
00651 {
00652 pt[cmpt] = 0.5*(min[cmpt] + max[cmpt]);
00653 }
00654 }
00655 }
00656
00657
00658 void Foam::meshTools::constrainToMeshCentre
00659 (
00660 const polyMesh& mesh,
00661 pointField& pts
00662 )
00663 {
00664 const Vector<label>& dirs = mesh.geometricD();
00665
00666 const point& min = mesh.bounds().min();
00667 const point& max = mesh.bounds().max();
00668
00669 bool isConstrained = false;
00670 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00671 {
00672 if (dirs[cmpt] == -1)
00673 {
00674 isConstrained = true;
00675 break;
00676 }
00677 }
00678
00679 if (isConstrained)
00680 {
00681 forAll(pts, i)
00682 {
00683 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00684 {
00685 if (dirs[cmpt] == -1)
00686 {
00687 pts[i][cmpt] = 0.5*(min[cmpt] + max[cmpt]);
00688 }
00689 }
00690 }
00691 }
00692 }
00693
00694
00695
00696 void Foam::meshTools::constrainDirection
00697 (
00698 const polyMesh& mesh,
00699 const Vector<label>& dirs,
00700 vector& d
00701 )
00702 {
00703 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00704 {
00705 if (dirs[cmpt] == -1)
00706 {
00707 d[cmpt] = 0.0;
00708 }
00709 }
00710 }
00711
00712
00713 void Foam::meshTools::constrainDirection
00714 (
00715 const polyMesh& mesh,
00716 const Vector<label>& dirs,
00717 vectorField& d
00718 )
00719 {
00720 bool isConstrained = false;
00721 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00722 {
00723 if (dirs[cmpt] == -1)
00724 {
00725 isConstrained = true;
00726 break;
00727 }
00728 }
00729
00730 if (isConstrained)
00731 {
00732 forAll(d, i)
00733 {
00734 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00735 {
00736 if (dirs[cmpt] == -1)
00737 {
00738 d[i][cmpt] = 0.0;
00739 }
00740 }
00741 }
00742 }
00743 }
00744
00745
00746 void Foam::meshTools::getParallelEdges
00747 (
00748 const primitiveMesh& mesh,
00749 const label cellI,
00750 const label e0,
00751 label& e1,
00752 label& e2,
00753 label& e3
00754 )
00755 {
00756
00757 label faceI = meshTools::otherFace(mesh, cellI, -1, e0);
00758
00759
00760 e1 = meshTools::walkFace(mesh, faceI, e0, mesh.edges()[e0].end(), 2);
00761
00762 faceI = meshTools::otherFace(mesh, cellI, faceI, e1);
00763
00764 e2 = meshTools::walkFace(mesh, faceI, e1, mesh.edges()[e1].end(), 2);
00765
00766 faceI = meshTools::otherFace(mesh, cellI, faceI, e2);
00767
00768 e3 = meshTools::walkFace(mesh, faceI, e2, mesh.edges()[e2].end(), 2);
00769 }
00770
00771
00772 Foam::vector Foam::meshTools::edgeToCutDir
00773 (
00774 const primitiveMesh& mesh,
00775 const label cellI,
00776 const label startEdgeI
00777 )
00778 {
00779 if (!hexMatcher().isA(mesh, cellI))
00780 {
00781 FatalErrorIn
00782 (
00783 "Foam::meshTools::getCutDir(const label, const label)"
00784 ) << "Not a hex : cell:" << cellI << abort(FatalError);
00785 }
00786
00787
00788 vector avgVec(normEdgeVec(mesh, startEdgeI));
00789
00790 label edgeI = startEdgeI;
00791
00792 label faceI = -1;
00793
00794 for (label i = 0; i < 3; i++)
00795 {
00796
00797 faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
00798
00799 vector eVec(normEdgeVec(mesh, edgeI));
00800
00801 if ((eVec & avgVec) > 0)
00802 {
00803 avgVec += eVec;
00804 }
00805 else
00806 {
00807 avgVec -= eVec;
00808 }
00809
00810 label vertI = mesh.edges()[edgeI].end();
00811
00812 edgeI = meshTools::walkFace(mesh, faceI, edgeI, vertI, 2);
00813 }
00814
00815 avgVec /= mag(avgVec) + VSMALL;
00816
00817 return avgVec;
00818 }
00819
00820
00821
00822 Foam::label Foam::meshTools::cutDirToEdge
00823 (
00824 const primitiveMesh& mesh,
00825 const label cellI,
00826 const vector& cutDir
00827 )
00828 {
00829 if (!hexMatcher().isA(mesh, cellI))
00830 {
00831 FatalErrorIn
00832 (
00833 "Foam::meshTools::getCutDir(const label, const vector&)"
00834 ) << "Not a hex : cell:" << cellI << abort(FatalError);
00835 }
00836
00837 const labelList& cEdges = mesh.cellEdges()[cellI];
00838
00839 labelHashSet doneEdges(2*cEdges.size());
00840
00841 scalar maxCos = -GREAT;
00842 label maxEdgeI = -1;
00843
00844 for (label i = 0; i < 4; i++)
00845 {
00846 forAll(cEdges, cEdgeI)
00847 {
00848 label e0 = cEdges[cEdgeI];
00849
00850 if (!doneEdges.found(e0))
00851 {
00852 vector avgDir(edgeToCutDir(mesh, cellI, e0));
00853
00854 scalar cosAngle = mag(avgDir & cutDir);
00855
00856 if (cosAngle > maxCos)
00857 {
00858 maxCos = cosAngle;
00859 maxEdgeI = e0;
00860 }
00861
00862
00863 label e1, e2, e3;
00864 getParallelEdges(mesh, cellI, e0, e1, e2, e3);
00865
00866 doneEdges.insert(e0);
00867 doneEdges.insert(e1);
00868 doneEdges.insert(e2);
00869 doneEdges.insert(e3);
00870 }
00871 }
00872 }
00873
00874 forAll(cEdges, cEdgeI)
00875 {
00876 if (!doneEdges.found(cEdges[cEdgeI]))
00877 {
00878 FatalErrorIn
00879 (
00880 "meshTools::cutDirToEdge(const label, const vector&)"
00881 ) << "Cell:" << cellI << " edges:" << cEdges << endl
00882 << "Edge:" << cEdges[cEdgeI] << " not yet handled"
00883 << abort(FatalError);
00884 }
00885 }
00886
00887 if (maxEdgeI == -1)
00888 {
00889 FatalErrorIn
00890 (
00891 "meshTools::cutDirToEdge(const label, const vector&)"
00892 ) << "Problem : did not find edge aligned with " << cutDir
00893 << " on cell " << cellI << abort(FatalError);
00894 }
00895
00896 return maxEdgeI;
00897 }
00898
00899
00900