meshTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "polyMesh.H"
28 #include "hexMatcher.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
33 (
34  const vector& n,
35  const vectorField& faceNormals,
36  const labelList& faceLabels
37 )
38 {
39  forAll(faceLabels, i)
40  {
41  if ((faceNormals[faceLabels[i]] & n) < small)
42  {
43  // Found normal in different direction from n.
44  return false;
45  }
46  }
47 
48  return true;
49 }
50 
51 
53 {
54  vectorField octantNormal(8);
55  octantNormal[mXmYmZ] = vector(-1, -1, -1);
56  octantNormal[pXmYmZ] = vector(1, -1, -1);
57  octantNormal[mXpYmZ] = vector(-1, 1, -1);
58  octantNormal[pXpYmZ] = vector(1, 1, -1);
59 
60  octantNormal[mXmYpZ] = vector(-1, -1, 1);
61  octantNormal[pXmYpZ] = vector(1, -1, 1);
62  octantNormal[mXpYpZ] = vector(-1, 1, 1);
63  octantNormal[pXpYpZ] = vector(1, 1, 1);
64 
65  octantNormal /= mag(octantNormal);
66 
67 
68  vectorField pn(pp.nPoints());
69 
70  const vectorField& faceNormals = pp.faceNormals();
71  const vectorField& pointNormals = pp.pointNormals();
72  const labelListList& pointFaces = pp.pointFaces();
73 
74  forAll(pointFaces, pointi)
75  {
76  const labelList& pFaces = pointFaces[pointi];
77 
78  if (visNormal(pointNormals[pointi], faceNormals, pFaces))
79  {
80  pn[pointi] = pointNormals[pointi];
81  }
82  else
83  {
85  << "Average point normal not visible for point:"
86  << pp.meshPoints()[pointi] << endl;
87 
88  label visOctant =
90  | pXmYmZMask
91  | mXpYmZMask
92  | pXpYmZMask
93  | mXmYpZMask
94  | pXmYpZMask
95  | mXpYpZMask
96  | pXpYpZMask;
97 
98  forAll(pFaces, i)
99  {
100  const vector& n = faceNormals[pFaces[i]];
101 
102  if (n.x() > small)
103  {
104  // All -x octants become invisible
105  visOctant &= ~mXmYmZMask;
106  visOctant &= ~mXmYpZMask;
107  visOctant &= ~mXpYmZMask;
108  visOctant &= ~mXpYpZMask;
109  }
110  else if (n.x() < -small)
111  {
112  // All +x octants become invisible
113  visOctant &= ~pXmYmZMask;
114  visOctant &= ~pXmYpZMask;
115  visOctant &= ~pXpYmZMask;
116  visOctant &= ~pXpYpZMask;
117  }
118 
119  if (n.y() > small)
120  {
121  visOctant &= ~mXmYmZMask;
122  visOctant &= ~mXmYpZMask;
123  visOctant &= ~pXmYmZMask;
124  visOctant &= ~pXmYpZMask;
125  }
126  else if (n.y() < -small)
127  {
128  visOctant &= ~mXpYmZMask;
129  visOctant &= ~mXpYpZMask;
130  visOctant &= ~pXpYmZMask;
131  visOctant &= ~pXpYpZMask;
132  }
133 
134  if (n.z() > small)
135  {
136  visOctant &= ~mXmYmZMask;
137  visOctant &= ~mXpYmZMask;
138  visOctant &= ~pXmYmZMask;
139  visOctant &= ~pXpYmZMask;
140  }
141  else if (n.z() < -small)
142  {
143  visOctant &= ~mXmYpZMask;
144  visOctant &= ~mXpYpZMask;
145  visOctant &= ~pXmYpZMask;
146  visOctant &= ~pXpYpZMask;
147  }
148  }
149 
150  label visI = -1;
151 
152  label mask = 1;
153 
154  for (label octant = 0; octant < 8; octant++)
155  {
156  if (visOctant & mask)
157  {
158  // Take first visible octant found. Note: could use
159  // combination of visible octants here.
160  visI = octant;
161 
162  break;
163  }
164  mask <<= 1;
165  }
166 
167  if (visI != -1)
168  {
169  // Take a vector in this octant.
170  pn[pointi] = octantNormal[visI];
171  }
172  else
173  {
174  pn[pointi] = Zero;
175 
177  << "No visible octant for point:" << pp.meshPoints()[pointi]
178  << " cooord:" << pp.points()[pp.meshPoints()[pointi]] << nl
179  << "Normal set to " << pn[pointi] << endl;
180  }
181  }
182  }
183 
184  return pn;
185 }
186 
187 
189 (
190  const primitiveMesh& mesh,
191  const label edgeI
192 )
193 {
194  vector eVec = mesh.edges()[edgeI].vec(mesh.points());
195 
196  eVec /= mag(eVec);
197 
198  return eVec;
199 }
200 
201 
203 (
204  Ostream& os,
205  const point& pt
206 )
207 {
208  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
209 }
210 
211 
213 (
214  Ostream& os,
215  const triad& t,
216  const point& pt
217 )
218 {
219  forAll(t, dirI)
220  {
221  writeOBJ(os, pt, pt + t[dirI]);
222  }
223 }
224 
225 
227 (
228  Ostream& os,
229  const point& p1,
230  const point& p2,
231  label& count
232 )
233 {
234  os << "v" << ' ' << p1.x() << ' ' << p1.y() << ' ' << p1.z() << endl;
235  os << "v" << ' ' << p2.x() << ' ' << p2.y() << ' ' << p2.z() << endl;
236 
237  os << "l" << " " << (count + 1) << " " << (count + 2) << endl;
238 
239  count += 2;
240 }
241 
242 
244 (
245  Ostream& os,
246  const point& p1,
247  const point& p2
248 )
249 {
250  os << "v" << ' ' << p1.x() << ' ' << p1.y() << ' ' << p1.z() << endl;
251 
252  os << "vn"
253  << ' ' << p2.x() - p1.x()
254  << ' ' << p2.y() - p1.y()
255  << ' ' << p2.z() - p1.z() << endl;
256 }
257 
258 
260 (
261  Ostream& os,
262  const cellList& cells,
263  const faceList& faces,
264  const pointField& points,
265  const labelList& cellLabels
266 )
267 {
268  labelHashSet usedFaces(4*cellLabels.size());
269 
270  forAll(cellLabels, i)
271  {
272  const cell& cFaces = cells[cellLabels[i]];
273 
274  forAll(cFaces, j)
275  {
276  usedFaces.insert(cFaces[j]);
277  }
278  }
279 
280  writeOBJ(os, faces, points, usedFaces.toc());
281 }
282 
283 
285 (
286  const primitiveMesh& mesh,
287  const label celli,
288  const label edgeI
289 )
290 {
291  return findIndex(mesh.edgeCells(edgeI), celli) != -1;
292 }
293 
294 
296 (
297  const primitiveMesh& mesh,
298  const label facei,
299  const label edgeI
300 )
301 {
302  return findIndex(mesh.faceEdges(facei), edgeI) != -1;
303 }
304 
305 
306 // Return true if facei part of celli
308 (
309  const primitiveMesh& mesh,
310  const label celli,
311  const label facei
312 )
313 {
314  if (mesh.isInternalFace(facei))
315  {
316  if
317  (
318  (mesh.faceOwner()[facei] == celli)
319  || (mesh.faceNeighbour()[facei] == celli)
320  )
321  {
322  return true;
323  }
324  }
325  else
326  {
327  if (mesh.faceOwner()[facei] == celli)
328  {
329  return true;
330  }
331  }
332  return false;
333 }
334 
335 
337 (
338  const edgeList& edges,
339  const labelList& candidates,
340  const label v0,
341  const label v1
342 )
343 {
344  forAll(candidates, i)
345  {
346  label edgeI = candidates[i];
347 
348  const edge& e = edges[edgeI];
349 
350  if ((e[0] == v0 && e[1] == v1) || (e[0] == v1 && e[1] == v0))
351  {
352  return edgeI;
353  }
354  }
355  return -1;
356 }
357 
358 
360 (
361  const primitiveMesh& mesh,
362  const label v0,
363  const label v1
364 )
365 {
366  const edgeList& edges = mesh.edges();
367 
368  const labelList& v0Edges = mesh.pointEdges()[v0];
369 
370  forAll(v0Edges, i)
371  {
372  label edgeI = v0Edges[i];
373 
374  const edge& e = edges[edgeI];
375 
376  if ((e.start() == v1) || (e.end() == v1))
377  {
378  return edgeI;
379  }
380  }
381  return -1;
382 }
383 
384 
386 (
387  const primitiveMesh& mesh,
388  const label f0,
389  const label f1
390 )
391 {
392  const labelList& f0Edges = mesh.faceEdges()[f0];
393  const labelList& f1Edges = mesh.faceEdges()[f1];
394 
395  forAll(f0Edges, f0EdgeI)
396  {
397  label edge0 = f0Edges[f0EdgeI];
398 
399  forAll(f1Edges, f1EdgeI)
400  {
401  label edge1 = f1Edges[f1EdgeI];
402 
403  if (edge0 == edge1)
404  {
405  return edge0;
406  }
407  }
408  }
410  << "Faces " << f0 << " and " << f1 << " do not share an edge"
411  << abort(FatalError);
412 
413  return -1;
414 
415 }
416 
417 
419 (
420  const primitiveMesh& mesh,
421  const label cell0I,
422  const label cell1I
423 )
424 {
425  const cell& cFaces = mesh.cells()[cell0I];
426 
427  forAll(cFaces, cFacei)
428  {
429  label facei = cFaces[cFacei];
430 
431  if
432  (
433  mesh.isInternalFace(facei)
434  && (
435  mesh.faceOwner()[facei] == cell1I
436  || mesh.faceNeighbour()[facei] == cell1I
437  )
438  )
439  {
440  return facei;
441  }
442  }
443 
444 
446  << "No common face for"
447  << " cell0I:" << cell0I << " faces:" << cFaces
448  << " cell1I:" << cell1I << " faces:"
449  << mesh.cells()[cell1I]
450  << abort(FatalError);
451 
452  return -1;
453 }
454 
455 
456 // Get the two faces on celli using edgeI.
458 (
459  const primitiveMesh& mesh,
460  const label celli,
461  const label edgeI,
462  label& face0,
463  label& face1
464 )
465 {
466  const labelList& eFaces = mesh.edgeFaces(edgeI);
467 
468  face0 = -1;
469  face1 = -1;
470 
471  forAll(eFaces, eFacei)
472  {
473  label facei = eFaces[eFacei];
474 
475  if (faceOnCell(mesh, celli, facei))
476  {
477  if (face0 == -1)
478  {
479  face0 = facei;
480  }
481  else
482  {
483  face1 = facei;
484 
485  return;
486  }
487  }
488  }
489 
490  if ((face0 == -1) || (face1 == -1))
491  {
493  << "Can not find faces using edge " << mesh.edges()[edgeI]
494  << " on cell " << celli << abort(FatalError);
495  }
496 }
497 
498 
499 // Return label of other edge connected to vertex
501 (
502  const primitiveMesh& mesh,
503  const labelList& edgeLabels,
504  const label thisEdgeI,
505  const label thisVertI
506 )
507 {
508  forAll(edgeLabels, edgeLabelI)
509  {
510  label edgeI = edgeLabels[edgeLabelI];
511 
512  if (edgeI != thisEdgeI)
513  {
514  const edge& e = mesh.edges()[edgeI];
515 
516  if ((e.start() == thisVertI) || (e.end() == thisVertI))
517  {
518  return edgeI;
519  }
520  }
521  }
522 
524  << "Can not find edge in "
525  << UIndirectList<edge>(mesh.edges(), edgeLabels)()
526  << " connected to edge "
527  << thisEdgeI << " with vertices " << mesh.edges()[thisEdgeI]
528  << " on side " << thisVertI << abort(FatalError);
529 
530  return -1;
531 }
532 
533 
534 // Return face on other side of edgeI
536 (
537  const primitiveMesh& mesh,
538  const label celli,
539  const label facei,
540  const label edgeI
541 )
542 {
543  label face0;
544  label face1;
545 
546  getEdgeFaces(mesh, celli, edgeI, face0, face1);
547 
548  if (face0 == facei)
549  {
550  return face1;
551  }
552  else
553  {
554  return face0;
555  }
556 }
557 
558 
559 // Return face on other side of edgeI
561 (
562  const primitiveMesh& mesh,
563  const label otherCelli,
564  const label facei
565 )
566 {
567  if (!mesh.isInternalFace(facei))
568  {
570  << "Face " << facei << " is not internal"
571  << abort(FatalError);
572  }
573 
574  label newCelli = mesh.faceOwner()[facei];
575 
576  if (newCelli == otherCelli)
577  {
578  newCelli = mesh.faceNeighbour()[facei];
579  }
580  return newCelli;
581 }
582 
583 
584 // Returns label of edge nEdges away from startEdge (in the direction of
585 // startVertI)
587 (
588  const primitiveMesh& mesh,
589  const label facei,
590  const label startEdgeI,
591  const label startVertI,
592  const label nEdges
593 )
594 {
595  const labelList& fEdges = mesh.faceEdges(facei);
596 
597  label edgeI = startEdgeI;
598 
599  label vertI = startVertI;
600 
601  for (label iter = 0; iter < nEdges; iter++)
602  {
603  edgeI = otherEdge(mesh, fEdges, edgeI, vertI);
604 
605  vertI = mesh.edges()[edgeI].otherVertex(vertI);
606  }
607 
608  return edgeI;
609 }
610 
611 
613 (
614  const polyMesh& mesh,
615  point& pt
616 )
617 {
618  const Vector<label>& dirs = mesh.geometricD();
619 
620  const point& min = mesh.bounds().min();
621  const point& max = mesh.bounds().max();
622 
623  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
624  {
625  if (dirs[cmpt] == -1)
626  {
627  pt[cmpt] = 0.5*(min[cmpt] + max[cmpt]);
628  }
629  }
630 }
631 
632 
634 (
635  const polyMesh& mesh,
636  pointField& pts
637 )
638 {
639  const Vector<label>& dirs = mesh.geometricD();
640 
641  const point& min = mesh.bounds().min();
642  const point& max = mesh.bounds().max();
643 
644  bool isConstrained = false;
645  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
646  {
647  if (dirs[cmpt] == -1)
648  {
649  isConstrained = true;
650  break;
651  }
652  }
653 
654  if (isConstrained)
655  {
656  forAll(pts, i)
657  {
658  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
659  {
660  if (dirs[cmpt] == -1)
661  {
662  pts[i][cmpt] = 0.5*(min[cmpt] + max[cmpt]);
663  }
664  }
665  }
666  }
667 }
668 
669 
671 (
672  const polyMesh& mesh,
673  const Vector<label>& dirs,
674  vector& d
675 )
676 {
677  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
678  {
679  if (dirs[cmpt] == -1)
680  {
681  d[cmpt] = 0.0;
682  }
683  }
684 }
685 
686 
688 (
689  const polyMesh& mesh,
690  const Vector<label>& dirs,
691  vectorField& d
692 )
693 {
694  bool isConstrained = false;
695  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
696  {
697  if (dirs[cmpt] == -1)
698  {
699  isConstrained = true;
700  break;
701  }
702  }
703 
704  if (isConstrained)
705  {
706  forAll(d, i)
707  {
708  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
709  {
710  if (dirs[cmpt] == -1)
711  {
712  d[i][cmpt] = 0.0;
713  }
714  }
715  }
716  }
717 }
718 
719 
721 (
722  const primitiveMesh& mesh,
723  const label celli,
724  const label e0,
725  label& e1,
726  label& e2,
727  label& e3
728 )
729 {
730  // Go to any face using e0
731  label facei = meshTools::otherFace(mesh, celli, -1, e0);
732 
733  // Opposite edge on face
734  e1 = meshTools::walkFace(mesh, facei, e0, mesh.edges()[e0].end(), 2);
735 
736  facei = meshTools::otherFace(mesh, celli, facei, e1);
737 
738  e2 = meshTools::walkFace(mesh, facei, e1, mesh.edges()[e1].end(), 2);
739 
740  facei = meshTools::otherFace(mesh, celli, facei, e2);
741 
742  e3 = meshTools::walkFace(mesh, facei, e2, mesh.edges()[e2].end(), 2);
743 }
744 
745 
747 (
748  const primitiveMesh& mesh,
749  const label celli,
750  const label startEdgeI
751 )
752 {
753  if (!hexMatcher().isA(mesh, celli))
754  {
756  << "Not a hex : cell:" << celli << abort(FatalError);
757  }
758 
759 
760  vector avgVec(normEdgeVec(mesh, startEdgeI));
761 
762  label edgeI = startEdgeI;
763 
764  label facei = -1;
765 
766  for (label i = 0; i < 3; i++)
767  {
768  // Step to next face, next edge
769  facei = meshTools::otherFace(mesh, celli, facei, edgeI);
770 
771  vector eVec(normEdgeVec(mesh, edgeI));
772 
773  if ((eVec & avgVec) > 0)
774  {
775  avgVec += eVec;
776  }
777  else
778  {
779  avgVec -= eVec;
780  }
781 
782  label vertI = mesh.edges()[edgeI].end();
783 
784  edgeI = meshTools::walkFace(mesh, facei, edgeI, vertI, 2);
785  }
786 
787  avgVec /= mag(avgVec) + vSmall;
788 
789  return avgVec;
790 }
791 
792 
793 // Find edges most aligned with cutDir
795 (
796  const primitiveMesh& mesh,
797  const label celli,
798  const vector& cutDir
799 )
800 {
801  if (!hexMatcher().isA(mesh, celli))
802  {
804  << "Not a hex : cell:" << celli << abort(FatalError);
805  }
806 
807  const labelList& cEdges = mesh.cellEdges()[celli];
808 
809  labelHashSet doneEdges(2*cEdges.size());
810 
811  scalar maxCos = -great;
812  label maxEdgeI = -1;
813 
814  for (label i = 0; i < 4; i++)
815  {
816  forAll(cEdges, cEdgeI)
817  {
818  label e0 = cEdges[cEdgeI];
819 
820  if (!doneEdges.found(e0))
821  {
822  vector avgDir(edgeToCutDir(mesh, celli, e0));
823 
824  scalar cosAngle = mag(avgDir & cutDir);
825 
826  if (cosAngle > maxCos)
827  {
828  maxCos = cosAngle;
829  maxEdgeI = e0;
830  }
831 
832  // Mark off edges in cEdges.
833  label e1, e2, e3;
834  getParallelEdges(mesh, celli, e0, e1, e2, e3);
835 
836  doneEdges.insert(e0);
837  doneEdges.insert(e1);
838  doneEdges.insert(e2);
839  doneEdges.insert(e3);
840  }
841  }
842  }
843 
844  forAll(cEdges, cEdgeI)
845  {
846  if (!doneEdges.found(cEdges[cEdgeI]))
847  {
849  << "Cell:" << celli << " edges:" << cEdges << endl
850  << "Edge:" << cEdges[cEdgeI] << " not yet handled"
851  << abort(FatalError);
852  }
853  }
854 
855  if (maxEdgeI == -1)
856  {
858  << "Problem : did not find edge aligned with " << cutDir
859  << " on cell " << celli << abort(FatalError);
860  }
861 
862  return maxEdgeI;
863 }
864 
865 
866 // ************************************************************************* //
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
Definition: typeInfo.H:134
vector edgeToCutDir(const primitiveMesh &, const label celli, const label edgeI)
Given edge on hex find all &#39;parallel&#39; (i.e. non-connected)
Definition: meshTools.C:747
label nPoints() const
Return number of points supporting patch faces.
const labelListList & cellEdges() const
bool edgeOnFace(const primitiveMesh &, const label facei, const label edgeI)
Is edge used by face.
Definition: meshTools.C:296
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
vector normEdgeVec(const primitiveMesh &, const label edgeI)
Normalised edge vector.
Definition: meshTools.C:189
label otherEdge(const primitiveMesh &, const labelList &edgeLabels, const label edgeI, const label vertI)
Return label of other edge (out of candidates edgeLabels)
Definition: meshTools.C:501
error FatalError
const labelListList & faceEdges() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
bool edgeOnCell(const primitiveMesh &, const label celli, const label edgeI)
Is edge used by cell.
Definition: meshTools.C:285
const labelListList & pointEdges() const
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label otherCell(const primitiveMesh &, const label celli, const label facei)
Return cell on other side of face. Throws error.
Definition: meshTools.C:561
virtual const pointField & points() const =0
Return mesh points.
label findEdge(const edgeList &edges, const labelList &candidates, const label v0, const label v1)
Return edge among candidates that uses the two vertices.
Definition: meshTools.C:337
const Cmpt & z() const
Definition: VectorI.H:87
label getSharedEdge(const primitiveMesh &, const label f0, const label f1)
Return edge shared by two faces. Throws error if no edge found.
Definition: meshTools.C:386
static const label mXmYmZMask
Definition: meshTools.H:73
A cellMatcher for hex cells.
Definition: hexMatcher.H:51
const cellList & cells() const
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
static const label pXpYmZ
Definition: meshTools.H:66
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:224
static const label mXpYmZ
Definition: meshTools.H:65
static const label pXmYpZ
Definition: meshTools.H:69
bool faceOnCell(const primitiveMesh &, const label celli, const label facei)
Is face used by cell.
Definition: meshTools.C:308
static const label mXpYpZ
Definition: meshTools.H:70
scalar f1
Definition: createFields.H:15
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
const Cmpt & y() const
Definition: VectorI.H:81
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:894
label getSharedFace(const primitiveMesh &, const label cell0, const label cell1)
Return face shared by two cells. Throws error if none found.
Definition: meshTools.C:419
const labelListList & edgeCells() const
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
static const label pXpYmZMask
Definition: meshTools.H:76
A list of faces which address into the list of points.
static const label mXmYpZ
Definition: meshTools.H:68
label cutDirToEdge(const primitiveMesh &, const label celli, const vector &cutDir)
Reverse of edgeToCutDir: given direction find edge bundle and.
Definition: meshTools.C:795
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static const label pXmYmZMask
Definition: meshTools.H:74
const Field< PointType > & points() const
Return reference to global points.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const zero Zero
Definition: zero.H:97
const Field< PointType > & faceNormals() const
Return face normals for patch.
bool visNormal(const vector &n, const vectorField &faceNormals, const labelList &faceLabels)
Check if n is in same direction as normals of all faceLabels.
Definition: meshTools.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
static const label pXpYpZ
Definition: meshTools.H:71
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const labelListList & pointFaces() const
Return point-face addressing.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
vectorField calcBoxPointNormals(const primitivePatch &pp)
Calculate point normals on a &#39;box&#39; mesh (all edges aligned with.
Definition: meshTools.C:52
static const char nl
Definition: Ostream.H:260
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Representation of a 3D Cartesian coordinate system as a Vector of vectors.
Definition: triad.H:65
const Field< PointType > & pointNormals() const
Return point normals for patch.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
static const label mXpYpZMask
Definition: meshTools.H:80
static const label mXmYpZMask
Definition: meshTools.H:78
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:459
static const label pXmYpZMask
Definition: meshTools.H:79
#define WarningInFunction
Report a warning using Foam::Warning.
label end() const
Return end vertex label.
Definition: edgeI.H:92
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A List with indirect addressing.
Definition: fvMatrix.H:106
static const label pXmYmZ
Definition: meshTools.H:64
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:671
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
label walkFace(const primitiveMesh &, const label facei, const label startEdgeI, const label startVertI, const label nEdges)
Returns label of edge nEdges away from startEdge (in the direction.
Definition: meshTools.C:587
static const label pXpYpZMask
Definition: meshTools.H:81
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:613
void getParallelEdges(const primitiveMesh &, const label celli, const label e0, label &, label &, label &)
Given edge on hex find other &#39;parallel&#39;, non-connected edges.
Definition: meshTools.C:721
static const label mXpYmZMask
Definition: meshTools.H:75
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
void getEdgeFaces(const primitiveMesh &, const label celli, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
Definition: meshTools.C:458
static const label mXmYmZ
Definition: meshTools.H:63
const labelListList & edgeFaces() const
label start() const
Return start vertex label.
Definition: edgeI.H:81
label otherFace(const primitiveMesh &, const label celli, const label facei, const label edgeI)
Return face on cell using edgeI but not facei. Throws error.
Definition: meshTools.C:536