FreeFOAM The Cross-Platform CFD Toolkit
meshTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "meshTools.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/hexMatcher.H>
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
32 // Check if n is in same direction as normals of all faceLabels
34 (
35  const vector& n,
36  const vectorField& faceNormals,
37  const labelList& faceLabels
38 )
39 {
40  forAll(faceLabels, i)
41  {
42  if ((faceNormals[faceLabels[i]] & n) < SMALL)
43  {
44  // Found normal in different direction from n.
45  return false;
46  }
47  }
48 
49  return true;
50 }
51 
52 
54 {
55  vectorField octantNormal(8);
56  octantNormal[mXmYmZ] = vector(-1, -1, -1);
57  octantNormal[pXmYmZ] = vector(1, -1, -1);
58  octantNormal[mXpYmZ] = vector(-1, 1, -1);
59  octantNormal[pXpYmZ] = vector(1, 1, -1);
60 
61  octantNormal[mXmYpZ] = vector(-1, -1, 1);
62  octantNormal[pXmYpZ] = vector(1, -1, 1);
63  octantNormal[mXpYpZ] = vector(-1, 1, 1);
64  octantNormal[pXpYpZ] = vector(1, 1, 1);
65 
66  octantNormal /= mag(octantNormal);
67 
68 
69  vectorField pn(pp.nPoints());
70 
71  const vectorField& faceNormals = pp.faceNormals();
72  const vectorField& pointNormals = pp.pointNormals();
73  const labelListList& pointFaces = pp.pointFaces();
74 
75  forAll(pointFaces, pointI)
76  {
77  const labelList& pFaces = pointFaces[pointI];
78 
79  if (visNormal(pointNormals[pointI], faceNormals, pFaces))
80  {
81  pn[pointI] = pointNormals[pointI];
82  }
83  else
84  {
85  WarningIn
86  (
87  "Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)"
88  ) << "Average point normal not visible for point:"
89  << pp.meshPoints()[pointI] << endl;
90 
91  label visOctant =
93  | pXmYmZMask
94  | mXpYmZMask
95  | pXpYmZMask
96  | mXmYpZMask
97  | pXmYpZMask
98  | mXpYpZMask
99  | pXpYpZMask;
100 
101  forAll(pFaces, i)
102  {
103  const vector& n = faceNormals[pFaces[i]];
104 
105  if (n.x() > SMALL)
106  {
107  // All -x octants become invisible
108  visOctant &= ~mXmYmZMask;
109  visOctant &= ~mXmYpZMask;
110  visOctant &= ~mXpYmZMask;
111  visOctant &= ~mXpYpZMask;
112  }
113  else if (n.x() < -SMALL)
114  {
115  // All +x octants become invisible
116  visOctant &= ~pXmYmZMask;
117  visOctant &= ~pXmYpZMask;
118  visOctant &= ~pXpYmZMask;
119  visOctant &= ~pXpYpZMask;
120  }
121 
122  if (n.y() > SMALL)
123  {
124  visOctant &= ~mXmYmZMask;
125  visOctant &= ~mXmYpZMask;
126  visOctant &= ~pXmYmZMask;
127  visOctant &= ~pXmYpZMask;
128  }
129  else if (n.y() < -SMALL)
130  {
131  visOctant &= ~mXpYmZMask;
132  visOctant &= ~mXpYpZMask;
133  visOctant &= ~pXpYmZMask;
134  visOctant &= ~pXpYpZMask;
135  }
136 
137  if (n.z() > SMALL)
138  {
139  visOctant &= ~mXmYmZMask;
140  visOctant &= ~mXpYmZMask;
141  visOctant &= ~pXmYmZMask;
142  visOctant &= ~pXpYmZMask;
143  }
144  else if (n.z() < -SMALL)
145  {
146  visOctant &= ~mXmYpZMask;
147  visOctant &= ~mXpYpZMask;
148  visOctant &= ~pXmYpZMask;
149  visOctant &= ~pXpYpZMask;
150  }
151  }
152 
153  label visI = -1;
154 
155  label mask = 1;
156 
157  for (label octant = 0; octant < 8; octant++)
158  {
159  if (visOctant & mask)
160  {
161  // Take first visible octant found. Note: could use
162  // combination of visible octants here.
163  visI = octant;
164 
165  break;
166  }
167  mask <<= 1;
168  }
169 
170  if (visI != -1)
171  {
172  // Take a vector in this octant.
173  pn[pointI] = octantNormal[visI];
174  }
175  else
176  {
177  pn[pointI] = vector::zero;
178 
179  WarningIn
180  (
181  "Foam::meshTools::calcBoxPointNormals"
182  "(const primitivePatch& pp)"
183  ) << "No visible octant for point:" << pp.meshPoints()[pointI]
184  << " cooord:" << pp.points()[pp.meshPoints()[pointI]] << nl
185  << "Normal set to " << pn[pointI] << endl;
186  }
187  }
188  }
189 
190  return pn;
191 }
192 
193 
195 (
196  const primitiveMesh& mesh,
197  const label edgeI
198 )
199 {
200  vector eVec = mesh.edges()[edgeI].vec(mesh.points());
201 
202  eVec /= mag(eVec);
203 
204  return eVec;
205 }
206 
207 
209 (
210  Ostream& os,
211  const point& pt
212 )
213 {
214  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
215 }
216 
217 
219 (
220  Ostream& os,
221  const faceList& faces,
222  const pointField& points,
223  const labelList& faceLabels
224 )
225 {
226  Map<label> foamToObj(4*faceLabels.size());
227 
228  label vertI = 0;
229 
230  forAll(faceLabels, i)
231  {
232  const face& f = faces[faceLabels[i]];
233 
234  forAll(f, fp)
235  {
236  if (foamToObj.insert(f[fp], vertI))
237  {
238  writeOBJ(os, points[f[fp]]);
239  vertI++;
240  }
241  }
242 
243  os << 'l';
244  forAll(f, fp)
245  {
246  os << ' ' << foamToObj[f[fp]]+1;
247  }
248  os << ' ' << foamToObj[f[0]]+1 << endl;
249  }
250 }
251 
252 
254 (
255  Ostream& os,
256  const faceList& faces,
257  const pointField& points
258 )
259 {
260  labelList allFaces(faces.size());
261  forAll(allFaces, i)
262  {
263  allFaces[i] = i;
264  }
265  writeOBJ(os, faces, points, allFaces);
266 }
267 
268 
270 (
271  Ostream& os,
272  const cellList& cells,
273  const faceList& faces,
274  const pointField& points,
275  const labelList& cellLabels
276 )
277 {
278  labelHashSet usedFaces(4*cellLabels.size());
279 
280  forAll(cellLabels, i)
281  {
282  const cell& cFaces = cells[cellLabels[i]];
283 
284  forAll(cFaces, j)
285  {
286  usedFaces.insert(cFaces[j]);
287  }
288  }
289 
290  writeOBJ(os, faces, points, usedFaces.toc());
291 }
292 
293 
295 (
296  const primitiveMesh& mesh,
297  const label cellI,
298  const label edgeI
299 )
300 {
301  return findIndex(mesh.edgeCells(edgeI), cellI) != -1;
302 }
303 
304 
306 (
307  const primitiveMesh& mesh,
308  const label faceI,
309  const label edgeI
310 )
311 {
312  return findIndex(mesh.faceEdges(faceI), edgeI) != -1;
313 }
314 
315 
316 // Return true if faceI part of cellI
318 (
319  const primitiveMesh& mesh,
320  const label cellI,
321  const label faceI
322 )
323 {
324  if (mesh.isInternalFace(faceI))
325  {
326  if
327  (
328  (mesh.faceOwner()[faceI] == cellI)
329  || (mesh.faceNeighbour()[faceI] == cellI)
330  )
331  {
332  return true;
333  }
334  }
335  else
336  {
337  if (mesh.faceOwner()[faceI] == cellI)
338  {
339  return true;
340  }
341  }
342  return false;
343 }
344 
345 
346 Foam::label Foam::meshTools::findEdge
347 (
348  const edgeList& edges,
349  const labelList& candidates,
350  const label v0,
351  const label v1
352 )
353 {
354  forAll(candidates, i)
355  {
356  label edgeI = candidates[i];
357 
358  const edge& e = edges[edgeI];
359 
360  if ((e[0] == v0 && e[1] == v1) || (e[0] == v1 && e[1] == v0))
361  {
362  return edgeI;
363  }
364  }
365  return -1;
366 }
367 
368 
369 Foam::label Foam::meshTools::findEdge
370 (
371  const primitiveMesh& mesh,
372  const label v0,
373  const label v1
374 )
375 {
376  const edgeList& edges = mesh.edges();
377 
378  const labelList& v0Edges = mesh.pointEdges()[v0];
379 
380  forAll(v0Edges, i)
381  {
382  label edgeI = v0Edges[i];
383 
384  const edge& e = edges[edgeI];
385 
386  if ((e.start() == v1) || (e.end() == v1))
387  {
388  return edgeI;
389  }
390  }
391  return -1;
392 }
393 
394 
396 (
397  const primitiveMesh& mesh,
398  const label f0,
399  const label f1
400 )
401 {
402  const labelList& f0Edges = mesh.faceEdges()[f0];
403  const labelList& f1Edges = mesh.faceEdges()[f1];
404 
405  forAll(f0Edges, f0EdgeI)
406  {
407  label edge0 = f0Edges[f0EdgeI];
408 
409  forAll(f1Edges, f1EdgeI)
410  {
411  label edge1 = f1Edges[f1EdgeI];
412 
413  if (edge0 == edge1)
414  {
415  return edge0;
416  }
417  }
418  }
420  (
421  "meshTools::getSharedEdge(const primitiveMesh&, const label"
422  ", const label)"
423  ) << "Faces " << f0 << " and " << f1 << " do not share an edge"
424  << abort(FatalError);
425 
426  return -1;
427 
428 }
429 
430 
432 (
433  const primitiveMesh& mesh,
434  const label cell0I,
435  const label cell1I
436 )
437 {
438  const cell& cFaces = mesh.cells()[cell0I];
439 
440  forAll(cFaces, cFaceI)
441  {
442  label faceI = cFaces[cFaceI];
443 
444  if
445  (
446  mesh.isInternalFace(faceI)
447  && (
448  mesh.faceOwner()[faceI] == cell1I
449  || mesh.faceNeighbour()[faceI] == cell1I
450  )
451  )
452  {
453  return faceI;
454  }
455  }
456 
457 
459  (
460  "meshTools::getSharedFace(const primitiveMesh&, const label"
461  ", const label)"
462  ) << "No common face for"
463  << " cell0I:" << cell0I << " faces:" << cFaces
464  << " cell1I:" << cell1I << " faces:"
465  << mesh.cells()[cell1I]
466  << abort(FatalError);
467 
468  return -1;
469 }
470 
471 
472 // Get the two faces on cellI using edgeI.
474 (
475  const primitiveMesh& mesh,
476  const label cellI,
477  const label edgeI,
478  label& face0,
479  label& face1
480 )
481 {
482  const labelList& eFaces = mesh.edgeFaces(edgeI);
483 
484  face0 = -1;
485  face1 = -1;
486 
487  forAll(eFaces, eFaceI)
488  {
489  label faceI = eFaces[eFaceI];
490 
491  if (faceOnCell(mesh, cellI, faceI))
492  {
493  if (face0 == -1)
494  {
495  face0 = faceI;
496  }
497  else
498  {
499  face1 = faceI;
500 
501  return;
502  }
503  }
504  }
505 
506  if ((face0 == -1) || (face1 == -1))
507  {
509  (
510  "meshTools::getEdgeFaces(const primitiveMesh&, const label"
511  ", const label, label&, label&"
512  ) << "Can not find faces using edge " << mesh.edges()[edgeI]
513  << " on cell " << cellI << abort(FatalError);
514  }
515 }
516 
517 
518 // Return label of other edge connected to vertex
519 Foam::label Foam::meshTools::otherEdge
520 (
521  const primitiveMesh& mesh,
522  const labelList& edgeLabels,
523  const label thisEdgeI,
524  const label thisVertI
525 )
526 {
527  forAll(edgeLabels, edgeLabelI)
528  {
529  label edgeI = edgeLabels[edgeLabelI];
530 
531  if (edgeI != thisEdgeI)
532  {
533  const edge& e = mesh.edges()[edgeI];
534 
535  if ((e.start() == thisVertI) || (e.end() == thisVertI))
536  {
537  return edgeI;
538  }
539  }
540  }
541 
543  (
544  "meshTools::otherEdge(const primitiveMesh&, const labelList&"
545  ", const label, const label)"
546  ) << "Can not find edge in "
547  << UIndirectList<edge>(mesh.edges(), edgeLabels)()
548  << " connected to edge "
549  << thisEdgeI << " with vertices " << mesh.edges()[thisEdgeI]
550  << " on side " << thisVertI << abort(FatalError);
551 
552  return -1;
553 }
554 
555 
556 // Return face on other side of edgeI
557 Foam::label Foam::meshTools::otherFace
558 (
559  const primitiveMesh& mesh,
560  const label cellI,
561  const label faceI,
562  const label edgeI
563 )
564 {
565  label face0;
566  label face1;
567 
568  getEdgeFaces(mesh, cellI, edgeI, face0, face1);
569 
570  if (face0 == faceI)
571  {
572  return face1;
573  }
574  else
575  {
576  return face0;
577  }
578 }
579 
580 
581 // Return face on other side of edgeI
582 Foam::label Foam::meshTools::otherCell
583 (
584  const primitiveMesh& mesh,
585  const label otherCellI,
586  const label faceI
587 )
588 {
589  if (!mesh.isInternalFace(faceI))
590  {
592  (
593  "meshTools::otherCell(const primitiveMesh&, const label"
594  ", const label)"
595  ) << "Face " << faceI << " is not internal"
596  << abort(FatalError);
597  }
598 
599  label newCellI = mesh.faceOwner()[faceI];
600 
601  if (newCellI == otherCellI)
602  {
603  newCellI = mesh.faceNeighbour()[faceI];
604  }
605  return newCellI;
606 }
607 
608 
609 // Returns label of edge nEdges away from startEdge (in the direction of
610 // startVertI)
611 Foam::label Foam::meshTools::walkFace
612 (
613  const primitiveMesh& mesh,
614  const label faceI,
615  const label startEdgeI,
616  const label startVertI,
617  const label nEdges
618 )
619 {
620  const labelList& fEdges = mesh.faceEdges(faceI);
621 
622  label edgeI = startEdgeI;
623 
624  label vertI = startVertI;
625 
626  for (label iter = 0; iter < nEdges; iter++)
627  {
628  edgeI = otherEdge(mesh, fEdges, edgeI, vertI);
629 
630  vertI = mesh.edges()[edgeI].otherVertex(vertI);
631  }
632 
633  return edgeI;
634 }
635 
636 
638 (
639  const polyMesh& mesh,
640  point& pt
641 )
642 {
643  const Vector<label>& dirs = mesh.geometricD();
644 
645  const point& min = mesh.bounds().min();
646  const point& max = mesh.bounds().max();
647 
648  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
649  {
650  if (dirs[cmpt] == -1)
651  {
652  pt[cmpt] = 0.5*(min[cmpt] + max[cmpt]);
653  }
654  }
655 }
656 
657 
659 (
660  const polyMesh& mesh,
661  pointField& pts
662 )
663 {
664  const Vector<label>& dirs = mesh.geometricD();
665 
666  const point& min = mesh.bounds().min();
667  const point& max = mesh.bounds().max();
668 
669  bool isConstrained = false;
670  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
671  {
672  if (dirs[cmpt] == -1)
673  {
674  isConstrained = true;
675  break;
676  }
677  }
678 
679  if (isConstrained)
680  {
681  forAll(pts, i)
682  {
683  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
684  {
685  if (dirs[cmpt] == -1)
686  {
687  pts[i][cmpt] = 0.5*(min[cmpt] + max[cmpt]);
688  }
689  }
690  }
691  }
692 }
693 
694 
695 //- Set the constrained components of directions/velocity to zero
697 (
698  const polyMesh& mesh,
699  const Vector<label>& dirs,
700  vector& d
701 )
702 {
703  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
704  {
705  if (dirs[cmpt] == -1)
706  {
707  d[cmpt] = 0.0;
708  }
709  }
710 }
711 
712 
714 (
715  const polyMesh& mesh,
716  const Vector<label>& dirs,
717  vectorField& d
718 )
719 {
720  bool isConstrained = false;
721  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
722  {
723  if (dirs[cmpt] == -1)
724  {
725  isConstrained = true;
726  break;
727  }
728  }
729 
730  if (isConstrained)
731  {
732  forAll(d, i)
733  {
734  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
735  {
736  if (dirs[cmpt] == -1)
737  {
738  d[i][cmpt] = 0.0;
739  }
740  }
741  }
742  }
743 }
744 
745 
747 (
748  const primitiveMesh& mesh,
749  const label cellI,
750  const label e0,
751  label& e1,
752  label& e2,
753  label& e3
754 )
755 {
756  // Go to any face using e0
757  label faceI = meshTools::otherFace(mesh, cellI, -1, e0);
758 
759  // Opposite edge on face
760  e1 = meshTools::walkFace(mesh, faceI, e0, mesh.edges()[e0].end(), 2);
761 
762  faceI = meshTools::otherFace(mesh, cellI, faceI, e1);
763 
764  e2 = meshTools::walkFace(mesh, faceI, e1, mesh.edges()[e1].end(), 2);
765 
766  faceI = meshTools::otherFace(mesh, cellI, faceI, e2);
767 
768  e3 = meshTools::walkFace(mesh, faceI, e2, mesh.edges()[e2].end(), 2);
769 }
770 
771 
773 (
774  const primitiveMesh& mesh,
775  const label cellI,
776  const label startEdgeI
777 )
778 {
779  if (!hexMatcher().isA(mesh, cellI))
780  {
782  (
783  "Foam::meshTools::getCutDir(const label, const label)"
784  ) << "Not a hex : cell:" << cellI << abort(FatalError);
785  }
786 
787 
788  vector avgVec(normEdgeVec(mesh, startEdgeI));
789 
790  label edgeI = startEdgeI;
791 
792  label faceI = -1;
793 
794  for (label i = 0; i < 3; i++)
795  {
796  // Step to next face, next edge
797  faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
798 
799  vector eVec(normEdgeVec(mesh, edgeI));
800 
801  if ((eVec & avgVec) > 0)
802  {
803  avgVec += eVec;
804  }
805  else
806  {
807  avgVec -= eVec;
808  }
809 
810  label vertI = mesh.edges()[edgeI].end();
811 
812  edgeI = meshTools::walkFace(mesh, faceI, edgeI, vertI, 2);
813  }
814 
815  avgVec /= mag(avgVec) + VSMALL;
816 
817  return avgVec;
818 }
819 
820 
821 // Find edges most aligned with cutDir
823 (
824  const primitiveMesh& mesh,
825  const label cellI,
826  const vector& cutDir
827 )
828 {
829  if (!hexMatcher().isA(mesh, cellI))
830  {
832  (
833  "Foam::meshTools::getCutDir(const label, const vector&)"
834  ) << "Not a hex : cell:" << cellI << abort(FatalError);
835  }
836 
837  const labelList& cEdges = mesh.cellEdges()[cellI];
838 
839  labelHashSet doneEdges(2*cEdges.size());
840 
841  scalar maxCos = -GREAT;
842  label maxEdgeI = -1;
843 
844  for (label i = 0; i < 4; i++)
845  {
846  forAll(cEdges, cEdgeI)
847  {
848  label e0 = cEdges[cEdgeI];
849 
850  if (!doneEdges.found(e0))
851  {
852  vector avgDir(edgeToCutDir(mesh, cellI, e0));
853 
854  scalar cosAngle = mag(avgDir & cutDir);
855 
856  if (cosAngle > maxCos)
857  {
858  maxCos = cosAngle;
859  maxEdgeI = e0;
860  }
861 
862  // Mark off edges in cEdges.
863  label e1, e2, e3;
864  getParallelEdges(mesh, cellI, e0, e1, e2, e3);
865 
866  doneEdges.insert(e0);
867  doneEdges.insert(e1);
868  doneEdges.insert(e2);
869  doneEdges.insert(e3);
870  }
871  }
872  }
873 
874  forAll(cEdges, cEdgeI)
875  {
876  if (!doneEdges.found(cEdges[cEdgeI]))
877  {
879  (
880  "meshTools::cutDirToEdge(const label, const vector&)"
881  ) << "Cell:" << cellI << " edges:" << cEdges << endl
882  << "Edge:" << cEdges[cEdgeI] << " not yet handled"
883  << abort(FatalError);
884  }
885  }
886 
887  if (maxEdgeI == -1)
888  {
890  (
891  "meshTools::cutDirToEdge(const label, const vector&)"
892  ) << "Problem : did not find edge aligned with " << cutDir
893  << " on cell " << cellI << abort(FatalError);
894  }
895 
896  return maxEdgeI;
897 }
898 
899 
900 // ************************ vim: set sw=4 sts=4 et: ************************ //