polyDualMesh.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) 2011-2016 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 InClass
25  polyDualMesh
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "polyDualMesh.H"
30 #include "meshTools.H"
31 #include "OFstream.H"
32 #include "Time.H"
33 #include "SortableList.H"
34 #include "pointSet.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(polyDualMesh, 0);
41 }
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 // Determine order for faces:
46 // - upper-triangular order for internal faces
47 // - external faces after internal faces (were already so)
48 Foam::labelList Foam::polyDualMesh::getFaceOrder
49 (
50  const labelList& faceOwner,
51  const labelList& faceNeighbour,
52  const cellList& cells,
53  label& nInternalFaces
54 )
55 {
56  labelList oldToNew(faceOwner.size(), -1);
57 
58  // First unassigned face
59  label newFacei = 0;
60 
61  forAll(cells, celli)
62  {
63  const labelList& cFaces = cells[celli];
64 
65  SortableList<label> nbr(cFaces.size());
66 
67  forAll(cFaces, i)
68  {
69  label facei = cFaces[i];
70 
71  label nbrCelli = faceNeighbour[facei];
72 
73  if (nbrCelli != -1)
74  {
75  // Internal face. Get cell on other side.
76  if (nbrCelli == celli)
77  {
78  nbrCelli = faceOwner[facei];
79  }
80 
81  if (celli < nbrCelli)
82  {
83  // Celli is master
84  nbr[i] = nbrCelli;
85  }
86  else
87  {
88  // nbrCell is master. Let it handle this face.
89  nbr[i] = -1;
90  }
91  }
92  else
93  {
94  // External face. Do later.
95  nbr[i] = -1;
96  }
97  }
98 
99  nbr.sort();
100 
101  forAll(nbr, i)
102  {
103  if (nbr[i] != -1)
104  {
105  oldToNew[cFaces[nbr.indices()[i]]] = newFacei++;
106  }
107  }
108  }
109 
110  nInternalFaces = newFacei;
111 
112  Pout<< "nInternalFaces:" << nInternalFaces << endl;
113  Pout<< "nFaces:" << faceOwner.size() << endl;
114  Pout<< "nCells:" << cells.size() << endl;
115 
116  // Leave patch faces intact.
117  for (label facei = newFacei; facei < faceOwner.size(); facei++)
118  {
119  oldToNew[facei] = facei;
120  }
121 
122 
123  // Check done all faces.
124  forAll(oldToNew, facei)
125  {
126  if (oldToNew[facei] == -1)
127  {
129  << "Did not determine new position"
130  << " for face " << facei
131  << abort(FatalError);
132  }
133  }
134 
135  return oldToNew;
136 }
137 
138 
139 // Get the two edges on facei using pointi. Returns them such that the order
140 // - otherVertex of e0
141 // - pointi
142 // - otherVertex(pointi) of e1
143 // is in face order
144 void Foam::polyDualMesh::getPointEdges
145 (
146  const primitivePatch& patch,
147  const label facei,
148  const label pointi,
149  label& e0,
150  label& e1
151 )
152 {
153  const labelList& fEdges = patch.faceEdges()[facei];
154  const face& f = patch.localFaces()[facei];
155 
156  e0 = -1;
157  e1 = -1;
158 
159  forAll(fEdges, i)
160  {
161  label edgeI = fEdges[i];
162 
163  const edge& e = patch.edges()[edgeI];
164 
165  if (e[0] == pointi)
166  {
167  // One of the edges using pointi. Check which one.
168  label index = findIndex(f, pointi);
169 
170  if (f.nextLabel(index) == e[1])
171  {
172  e1 = edgeI;
173  }
174  else
175  {
176  e0 = edgeI;
177  }
178 
179  if (e0 != -1 && e1 != -1)
180  {
181  return;
182  }
183  }
184  else if (e[1] == pointi)
185  {
186  // One of the edges using pointi. Check which one.
187  label index = findIndex(f, pointi);
188 
189  if (f.nextLabel(index) == e[0])
190  {
191  e1 = edgeI;
192  }
193  else
194  {
195  e0 = edgeI;
196  }
197 
198  if (e0 != -1 && e1 != -1)
199  {
200  return;
201  }
202  }
203  }
204 
206  << " vertices:" << patch.localFaces()[facei]
207  << " that uses point:" << pointi
208  << abort(FatalError);
209 }
210 
211 
212 // Collect the face on an external point of the patch.
213 Foam::labelList Foam::polyDualMesh::collectPatchSideFace
214 (
215  const polyPatch& patch,
216  const label patchToDualOffset,
217  const labelList& edgeToDualPoint,
218  const labelList& pointToDualPoint,
219  const label pointi,
220 
221  label& edgeI
222 )
223 {
224  // Construct face by walking around e[eI] starting from
225  // patchEdgeI
226 
227  label meshPointi = patch.meshPoints()[pointi];
228  const labelList& pFaces = patch.pointFaces()[pointi];
229 
230  DynamicList<label> dualFace;
231 
232  if (pointToDualPoint[meshPointi] >= 0)
233  {
234  // Number of pFaces + 2 boundary edge + feature point
235  dualFace.setCapacity(pFaces.size()+2+1);
236  // Store dualVertex for feature edge
237  dualFace.append(pointToDualPoint[meshPointi]);
238  }
239  else
240  {
241  dualFace.setCapacity(pFaces.size()+2);
242  }
243 
244  // Store dual vertex for starting edge.
245  if (edgeToDualPoint[patch.meshEdges()[edgeI]] < 0)
246  {
248  << abort(FatalError);
249  }
250 
251  dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
252 
253  label facei = patch.edgeFaces()[edgeI][0];
254 
255  // Check order of vertices of edgeI in face to see if we need to reverse.
256  bool reverseFace;
257 
258  label e0, e1;
259  getPointEdges(patch, facei, pointi, e0, e1);
260 
261  if (e0 == edgeI)
262  {
263  reverseFace = true;
264  }
265  else
266  {
267  reverseFace = false;
268  }
269 
270  while (true)
271  {
272  // Store dual vertex for facei.
273  dualFace.append(facei + patchToDualOffset);
274 
275  // Cross face to other edge on pointi
276  label e0, e1;
277  getPointEdges(patch, facei, pointi, e0, e1);
278 
279  if (e0 == edgeI)
280  {
281  edgeI = e1;
282  }
283  else
284  {
285  edgeI = e0;
286  }
287 
288  if (edgeToDualPoint[patch.meshEdges()[edgeI]] >= 0)
289  {
290  // Feature edge. Insert dual vertex for edge.
291  dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
292  }
293 
294  const labelList& eFaces = patch.edgeFaces()[edgeI];
295 
296  if (eFaces.size() == 1)
297  {
298  // Reached other edge of patch
299  break;
300  }
301 
302  // Cross edge to other face.
303  if (eFaces[0] == facei)
304  {
305  facei = eFaces[1];
306  }
307  else
308  {
309  facei = eFaces[0];
310  }
311  }
312 
313  dualFace.shrink();
314 
315  if (reverseFace)
316  {
317  reverse(dualFace);
318  }
319 
320  return dualFace;
321 }
322 
323 
324 // Collect face around pointi which is not on the outside of the patch.
325 // Returns the vertices of the face and the indices in these vertices of
326 // any points which are on feature edges.
327 void Foam::polyDualMesh::collectPatchInternalFace
328 (
329  const polyPatch& patch,
330  const label patchToDualOffset,
331  const labelList& edgeToDualPoint,
332 
333  const label pointi,
334  const label startEdgeI,
335 
336  labelList& dualFace2,
337  labelList& featEdgeIndices2
338 )
339 {
340  // Construct face by walking around pointi starting from startEdgeI
341  const labelList& meshEdges = patch.meshEdges();
342  const labelList& pFaces = patch.pointFaces()[pointi];
343 
344  // Vertices of dualFace
345  DynamicList<label> dualFace(pFaces.size());
346  // Indices in dualFace of vertices added because of feature edge
347  DynamicList<label> featEdgeIndices(dualFace.size());
348 
349 
350  label edgeI = startEdgeI;
351  label facei = patch.edgeFaces()[edgeI][0];
352 
353  // Check order of vertices of edgeI in face to see if we need to reverse.
354  bool reverseFace;
355 
356  label e0, e1;
357  getPointEdges(patch, facei, pointi, e0, e1);
358 
359  if (e0 == edgeI)
360  {
361  reverseFace = true;
362  }
363  else
364  {
365  reverseFace = false;
366  }
367 
368  while (true)
369  {
370  // Insert dual vertex for face
371  dualFace.append(facei + patchToDualOffset);
372 
373  // Cross face to other edge on pointi
374  label e0, e1;
375  getPointEdges(patch, facei, pointi, e0, e1);
376 
377  if (e0 == edgeI)
378  {
379  edgeI = e1;
380  }
381  else
382  {
383  edgeI = e0;
384  }
385 
386  if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
387  {
388  // Feature edge. Insert dual vertex for edge.
389  dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
390  featEdgeIndices.append(dualFace.size()-1);
391  }
392 
393  if (edgeI == startEdgeI)
394  {
395  break;
396  }
397 
398  // Cross edge to other face.
399  const labelList& eFaces = patch.edgeFaces()[edgeI];
400 
401  if (eFaces[0] == facei)
402  {
403  facei = eFaces[1];
404  }
405  else
406  {
407  facei = eFaces[0];
408  }
409  }
410 
411  dualFace2.transfer(dualFace);
412 
413  featEdgeIndices2.transfer(featEdgeIndices);
414 
415  if (reverseFace)
416  {
417  reverse(dualFace2);
418 
419  // Correct featEdgeIndices for change in dualFace2
420  forAll(featEdgeIndices2, i)
421  {
422  featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
423  }
424  // Reverse indices (might not be necessary but do anyway)
425  reverse(featEdgeIndices2);
426  }
427 }
428 
429 
430 void Foam::polyDualMesh::splitFace
431 (
432  const polyPatch& patch,
433  const labelList& pointToDualPoint,
434 
435  const label pointi,
436  const labelList& dualFace,
437  const labelList& featEdgeIndices,
438 
439  DynamicList<face>& dualFaces,
440  DynamicList<label>& dualOwner,
441  DynamicList<label>& dualNeighbour,
442  DynamicList<label>& dualRegion
443 )
444 {
445 
446  // Split face because of feature edges/points
447  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
448  label meshPointi = patch.meshPoints()[pointi];
449 
450  if (pointToDualPoint[meshPointi] >= 0)
451  {
452  // Feature point. Do face-centre decomposition.
453 
454  if (featEdgeIndices.size() < 2)
455  {
456  // Feature point but no feature edges. Not handled at the moment
457  dualFaces.append(face(dualFace));
458  dualOwner.append(meshPointi);
459  dualNeighbour.append(-1);
460  dualRegion.append(patch.index());
461  }
462  else
463  {
464  // Do 'face-centre' decomposition. Start from first feature
465  // edge create face up until next feature edge.
466 
467  forAll(featEdgeIndices, i)
468  {
469  label startFp = featEdgeIndices[i];
470 
471  label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
472 
473  // Collect face from startFp to endFp
474  label sz = 0;
475 
476  if (endFp > startFp)
477  {
478  sz = endFp - startFp + 2;
479  }
480  else
481  {
482  sz = endFp + dualFace.size() - startFp + 2;
483  }
484  face subFace(sz);
485 
486  // feature point becomes face centre.
487  subFace[0] = pointToDualPoint[patch.meshPoints()[pointi]];
488 
489  // Copy from startFp up to endFp.
490  for (label subFp = 1; subFp < subFace.size(); subFp++)
491  {
492  subFace[subFp] = dualFace[startFp];
493 
494  startFp = (startFp+1) % dualFace.size();
495  }
496 
497  dualFaces.append(face(subFace));
498  dualOwner.append(meshPointi);
499  dualNeighbour.append(-1);
500  dualRegion.append(patch.index());
501  }
502  }
503  }
504  else
505  {
506  // No feature point. Check if feature edges.
507  if (featEdgeIndices.size() < 2)
508  {
509  // Not enough feature edges. No split.
510  dualFaces.append(face(dualFace));
511  dualOwner.append(meshPointi);
512  dualNeighbour.append(-1);
513  dualRegion.append(patch.index());
514  }
515  else
516  {
517  // Multiple feature edges but no feature point.
518  // Split face along feature edges. Gives weird result if
519  // number of feature edges > 2.
520 
521  // Storage for new face
522  DynamicList<label> subFace(dualFace.size());
523 
524  forAll(featEdgeIndices, featI)
525  {
526  label startFp = featEdgeIndices[featI];
527  label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
528 
529  label fp = startFp;
530 
531  while (true)
532  {
533  label vertI = dualFace[fp];
534 
535  subFace.append(vertI);
536 
537  if (fp == endFp)
538  {
539  break;
540  }
541 
542  fp = dualFace.fcIndex(fp);
543  }
544 
545  if (subFace.size() > 2)
546  {
547  // Enough vertices to create a face from.
548  subFace.shrink();
549 
550  dualFaces.append(face(subFace));
551  dualOwner.append(meshPointi);
552  dualNeighbour.append(-1);
553  dualRegion.append(patch.index());
554 
555  subFace.clear();
556  }
557  }
558  // Check final face.
559  if (subFace.size() > 2)
560  {
561  // Enough vertices to create a face from.
562  subFace.shrink();
563 
564  dualFaces.append(face(subFace));
565  dualOwner.append(meshPointi);
566  dualNeighbour.append(-1);
567  dualRegion.append(patch.index());
568 
569  subFace.clear();
570  }
571  }
572  }
573 }
574 
575 
576 // Create boundary face for every point in patch
577 void Foam::polyDualMesh::dualPatch
578 (
579  const polyPatch& patch,
580  const label patchToDualOffset,
581  const labelList& edgeToDualPoint,
582  const labelList& pointToDualPoint,
583 
584  const pointField& dualPoints,
585 
586  DynamicList<face>& dualFaces,
587  DynamicList<label>& dualOwner,
588  DynamicList<label>& dualNeighbour,
589  DynamicList<label>& dualRegion
590 )
591 {
592  const labelList& meshEdges = patch.meshEdges();
593 
594  // Whether edge has been done.
595  // 0 : not
596  // 1 : done e.start()
597  // 2 : done e.end()
598  // 3 : done both
599  labelList doneEdgeSide(meshEdges.size(), 0);
600 
601  boolList donePoint(patch.nPoints(), false);
602 
603 
604  // Do points on edge of patch
605  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
606 
607  forAll(doneEdgeSide, patchEdgeI)
608  {
609  const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
610 
611  if (eFaces.size() == 1)
612  {
613  const edge& e = patch.edges()[patchEdgeI];
614 
615  forAll(e, eI)
616  {
617  label bitMask = 1<<eI;
618 
619  if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
620  {
621  // Construct face by walking around e[eI] starting from
622  // patchEdgeI
623 
624  label pointi = e[eI];
625 
626  label edgeI = patchEdgeI;
627  labelList dualFace
628  (
629  collectPatchSideFace
630  (
631  patch,
632  patchToDualOffset,
633  edgeToDualPoint,
634  pointToDualPoint,
635 
636  pointi,
637  edgeI
638  )
639  );
640 
641  // Now edgeI is end edge. Mark as visited
642  if (patch.edges()[edgeI][0] == pointi)
643  {
644  doneEdgeSide[edgeI] |= 1;
645  }
646  else
647  {
648  doneEdgeSide[edgeI] |= 2;
649  }
650 
651  dualFaces.append(face(dualFace));
652  dualOwner.append(patch.meshPoints()[pointi]);
653  dualNeighbour.append(-1);
654  dualRegion.append(patch.index());
655 
656  doneEdgeSide[patchEdgeI] |= bitMask;
657  donePoint[pointi] = true;
658  }
659  }
660  }
661  }
662 
663 
664 
665  // Do patch-internal points
666  // ~~~~~~~~~~~~~~~~~~~~~~~~
667 
668  forAll(donePoint, pointi)
669  {
670  if (!donePoint[pointi])
671  {
672  labelList dualFace, featEdgeIndices;
673 
674  collectPatchInternalFace
675  (
676  patch,
677  patchToDualOffset,
678  edgeToDualPoint,
679  pointi,
680  patch.pointEdges()[pointi][0], // Arbitrary starting edge
681 
682  dualFace,
683  featEdgeIndices
684  );
685 
686  //- Either keep in one piece or split face according to feature.
687 
689  //dualFaces.append(face(dualFace));
690  //dualOwner.append(patch.meshPoints()[pointi]);
691  //dualNeighbour.append(-1);
692  //dualRegion.append(patch.index());
693 
694  splitFace
695  (
696  patch,
697  pointToDualPoint,
698  pointi,
699  dualFace,
700  featEdgeIndices,
701 
702  dualFaces,
703  dualOwner,
704  dualNeighbour,
705  dualRegion
706  );
707 
708  donePoint[pointi] = true;
709  }
710  }
711 }
712 
713 
714 void Foam::polyDualMesh::calcDual
715 (
716  const polyMesh& mesh,
717  const labelList& featureEdges,
718  const labelList& featurePoints
719 )
720 {
721  const polyBoundaryMesh& patches = mesh.boundaryMesh();
722  const label nIntFaces = mesh.nInternalFaces();
723 
724 
725  // Get patch for all of outside
726  primitivePatch allBoundary
727  (
728  SubList<face>
729  (
730  mesh.faces(),
731  mesh.nFaces() - nIntFaces,
732  nIntFaces
733  ),
734  mesh.points()
735  );
736 
737  // Correspondence from patch edge to mesh edge.
738  labelList meshEdges
739  (
740  allBoundary.meshEdges
741  (
742  mesh.edges(),
743  mesh.pointEdges()
744  )
745  );
746 
747  {
748  pointSet nonManifoldPoints
749  (
750  mesh,
751  "nonManifoldPoints",
752  meshEdges.size() / 100
753  );
754 
755  allBoundary.checkPointManifold(true, &nonManifoldPoints);
756 
757  if (nonManifoldPoints.size())
758  {
759  nonManifoldPoints.write();
760 
762  << "There are " << nonManifoldPoints.size() << " points where"
763  << " the outside of the mesh is non-manifold." << nl
764  << "Such a mesh cannot be converted to a dual." << nl
765  << "Writing points at non-manifold sites to pointSet "
766  << nonManifoldPoints.name()
767  << exit(FatalError);
768  }
769  }
770 
771 
772  // Assign points
773  // ~~~~~~~~~~~~~
774 
775  // mesh label dualMesh vertex
776  // ---------- ---------------
777  // celli celli
778  // facei nCells+facei-nIntFaces
779  // featEdgeI nCells+nFaces-nIntFaces+featEdgeI
780  // featPointi nCells+nFaces-nIntFaces+nFeatEdges+featPointi
781 
782  pointField dualPoints
783  (
784  mesh.nCells() // cell centres
785  + mesh.nFaces() - nIntFaces // boundary face centres
786  + featureEdges.size() // additional boundary edges
787  + featurePoints.size() // additional boundary points
788  );
789 
790  label dualPointi = 0;
791 
792 
793  // Cell centres.
794  const pointField& cellCentres = mesh.cellCentres();
795 
796  cellPoint_.setSize(cellCentres.size());
797 
798  forAll(cellCentres, celli)
799  {
800  cellPoint_[celli] = dualPointi;
801  dualPoints[dualPointi++] = cellCentres[celli];
802  }
803 
804  // Boundary faces centres
805  const pointField& faceCentres = mesh.faceCentres();
806 
807  boundaryFacePoint_.setSize(mesh.nFaces() - nIntFaces);
808 
809  for (label facei = nIntFaces; facei < mesh.nFaces(); facei++)
810  {
811  boundaryFacePoint_[facei - nIntFaces] = dualPointi;
812  dualPoints[dualPointi++] = faceCentres[facei];
813  }
814 
815  // Edge status:
816  // >0 : dual point label (edge is feature edge)
817  // -1 : is boundary edge.
818  // -2 : is internal edge.
819  labelList edgeToDualPoint(mesh.nEdges(), -2);
820 
821  forAll(meshEdges, patchEdgeI)
822  {
823  label edgeI = meshEdges[patchEdgeI];
824 
825  edgeToDualPoint[edgeI] = -1;
826  }
827 
828  forAll(featureEdges, i)
829  {
830  label edgeI = featureEdges[i];
831 
832  const edge& e = mesh.edges()[edgeI];
833 
834  edgeToDualPoint[edgeI] = dualPointi;
835  dualPoints[dualPointi++] = e.centre(mesh.points());
836  }
837 
838 
839 
840  // Point status:
841  // >0 : dual point label (point is feature point)
842  // -1 : is point on edge of patch
843  // -2 : is point on patch (but not on edge)
844  // -3 : is internal point.
845  labelList pointToDualPoint(mesh.nPoints(), -3);
846 
847  forAll(patches, patchi)
848  {
849  const labelList& meshPoints = patches[patchi].meshPoints();
850 
851  forAll(meshPoints, i)
852  {
853  pointToDualPoint[meshPoints[i]] = -2;
854  }
855 
856  const labelListList& loops = patches[patchi].edgeLoops();
857 
858  forAll(loops, i)
859  {
860  const labelList& loop = loops[i];
861 
862  forAll(loop, j)
863  {
864  pointToDualPoint[meshPoints[loop[j]]] = -1;
865  }
866  }
867  }
868 
869  forAll(featurePoints, i)
870  {
871  label pointi = featurePoints[i];
872 
873  pointToDualPoint[pointi] = dualPointi;
874  dualPoints[dualPointi++] = mesh.points()[pointi];
875  }
876 
877 
878  // Storage for new faces.
879  // Dynamic sized since we don't know size.
880 
881  DynamicList<face> dynDualFaces(mesh.nEdges());
882  DynamicList<label> dynDualOwner(mesh.nEdges());
883  DynamicList<label> dynDualNeighbour(mesh.nEdges());
884  DynamicList<label> dynDualRegion(mesh.nEdges());
885 
886 
887  // Generate faces from edges on the boundary
888  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
889 
890  forAll(meshEdges, patchEdgeI)
891  {
892  label edgeI = meshEdges[patchEdgeI];
893 
894  const edge& e = mesh.edges()[edgeI];
895 
896  label owner = -1;
897  label neighbour = -1;
898 
899  if (e[0] < e[1])
900  {
901  owner = e[0];
902  neighbour = e[1];
903  }
904  else
905  {
906  owner = e[1];
907  neighbour = e[0];
908  }
909 
910  // Find the boundary faces using the edge.
911  const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
912 
913  if (patchFaces.size() != 2)
914  {
916  << "Cannot handle edges with " << patchFaces.size()
917  << " connected boundary faces."
918  << abort(FatalError);
919  }
920 
921  label face0 = patchFaces[0] + nIntFaces;
922  const face& f0 = mesh.faces()[face0];
923 
924  label face1 = patchFaces[1] + nIntFaces;
925 
926  // We want to start walking from patchFaces[0] or patchFaces[1],
927  // depending on which one uses owner,neighbour in the right order.
928 
929  label startFacei = -1;
930  label endFacei = -1;
931 
932  label index = findIndex(f0, neighbour);
933 
934  if (f0.nextLabel(index) == owner)
935  {
936  startFacei = face0;
937  endFacei = face1;
938  }
939  else
940  {
941  startFacei = face1;
942  endFacei = face0;
943  }
944 
945  // Now walk from startFacei to cell to face to cell etc. until endFacei
946 
947  DynamicList<label> dualFace;
948 
949  if (edgeToDualPoint[edgeI] >= 0)
950  {
951  // Number of cells + 2 boundary faces + feature edge point
952  dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
953  // Store dualVertex for feature edge
954  dualFace.append(edgeToDualPoint[edgeI]);
955  }
956  else
957  {
958  dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
959  }
960 
961  // Store dual vertex for starting face.
962  dualFace.append(mesh.nCells() + startFacei - nIntFaces);
963 
964  label celli = mesh.faceOwner()[startFacei];
965  label facei = startFacei;
966 
967  while (true)
968  {
969  // Store dual vertex from celli.
970  dualFace.append(celli);
971 
972  // Cross cell to other face on edge.
973  label f0, f1;
974  meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
975 
976  if (f0 == facei)
977  {
978  facei = f1;
979  }
980  else
981  {
982  facei = f0;
983  }
984 
985  // Cross face to other cell.
986  if (facei == endFacei)
987  {
988  break;
989  }
990 
991  if (mesh.faceOwner()[facei] == celli)
992  {
993  celli = mesh.faceNeighbour()[facei];
994  }
995  else
996  {
997  celli = mesh.faceOwner()[facei];
998  }
999  }
1000 
1001  // Store dual vertex for endFace.
1002  dualFace.append(mesh.nCells() + endFacei - nIntFaces);
1003 
1004  dynDualFaces.append(face(dualFace.shrink()));
1005  dynDualOwner.append(owner);
1006  dynDualNeighbour.append(neighbour);
1007  dynDualRegion.append(-1);
1008 
1009  {
1010  // Check orientation.
1011  const face& f = dynDualFaces.last();
1012  vector n = f.normal(dualPoints);
1013  if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1014  {
1016  << " on boundary edge:" << edgeI
1017  << mesh.points()[mesh.edges()[edgeI][0]]
1018  << mesh.points()[mesh.edges()[edgeI][1]]
1019  << endl;
1020  }
1021  }
1022  }
1023 
1024 
1025  // Generate faces from internal edges
1026  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1027 
1028  forAll(edgeToDualPoint, edgeI)
1029  {
1030  if (edgeToDualPoint[edgeI] == -2)
1031  {
1032  // Internal edge.
1033 
1034  // Find dual owner, neighbour
1035 
1036  const edge& e = mesh.edges()[edgeI];
1037 
1038  label owner = -1;
1039  label neighbour = -1;
1040 
1041  if (e[0] < e[1])
1042  {
1043  owner = e[0];
1044  neighbour = e[1];
1045  }
1046  else
1047  {
1048  owner = e[1];
1049  neighbour = e[0];
1050  }
1051 
1052  // Get a starting cell
1053  const labelList& eCells = mesh.edgeCells()[edgeI];
1054 
1055  label celli = eCells[0];
1056 
1057  // Get the two faces on the cell and edge.
1058  label face0, face1;
1059  meshTools::getEdgeFaces(mesh, celli, edgeI, face0, face1);
1060 
1061  // Find the starting face by looking at the order in which
1062  // the face uses the owner, neighbour
1063  const face& f0 = mesh.faces()[face0];
1064 
1065  label index = findIndex(f0, neighbour);
1066 
1067  bool f0OrderOk = (f0.nextLabel(index) == owner);
1068 
1069  label startFacei = -1;
1070 
1071  if (f0OrderOk == (mesh.faceOwner()[face0] == celli))
1072  {
1073  startFacei = face0;
1074  }
1075  else
1076  {
1077  startFacei = face1;
1078  }
1079 
1080  // Walk face-cell-face until starting face reached.
1081  DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
1082 
1083  label facei = startFacei;
1084 
1085  while (true)
1086  {
1087  // Store dual vertex from celli.
1088  dualFace.append(celli);
1089 
1090  // Cross cell to other face on edge.
1091  label f0, f1;
1092  meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
1093 
1094  if (f0 == facei)
1095  {
1096  facei = f1;
1097  }
1098  else
1099  {
1100  facei = f0;
1101  }
1102 
1103  // Cross face to other cell.
1104  if (facei == startFacei)
1105  {
1106  break;
1107  }
1108 
1109  if (mesh.faceOwner()[facei] == celli)
1110  {
1111  celli = mesh.faceNeighbour()[facei];
1112  }
1113  else
1114  {
1115  celli = mesh.faceOwner()[facei];
1116  }
1117  }
1118 
1119  dynDualFaces.append(face(dualFace.shrink()));
1120  dynDualOwner.append(owner);
1121  dynDualNeighbour.append(neighbour);
1122  dynDualRegion.append(-1);
1123 
1124  {
1125  // Check orientation.
1126  const face& f = dynDualFaces.last();
1127  vector n = f.normal(dualPoints);
1128  if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1129  {
1131  << " on internal edge:" << edgeI
1132  << mesh.points()[mesh.edges()[edgeI][0]]
1133  << mesh.points()[mesh.edges()[edgeI][1]]
1134  << endl;
1135  }
1136  }
1137  }
1138  }
1139 
1140  // Dump faces.
1141  if (debug)
1142  {
1143  dynDualFaces.shrink();
1144  dynDualOwner.shrink();
1145  dynDualNeighbour.shrink();
1146  dynDualRegion.shrink();
1147 
1148  OFstream str("dualInternalFaces.obj");
1149  Pout<< "polyDualMesh::calcDual : dumping internal faces to "
1150  << str.name() << endl;
1151 
1152  forAll(dualPoints, dualPointi)
1153  {
1154  meshTools::writeOBJ(str, dualPoints[dualPointi]);
1155  }
1156  forAll(dynDualFaces, dualFacei)
1157  {
1158  const face& f = dynDualFaces[dualFacei];
1159 
1160  str<< 'f';
1161  forAll(f, fp)
1162  {
1163  str<< ' ' << f[fp]+1;
1164  }
1165  str<< nl;
1166  }
1167  }
1168 
1169  const label nInternalFaces = dynDualFaces.size();
1170 
1171  // Outside faces
1172  // ~~~~~~~~~~~~~
1173 
1174  forAll(patches, patchi)
1175  {
1176  const polyPatch& pp = patches[patchi];
1177 
1178  dualPatch
1179  (
1180  pp,
1181  mesh.nCells() + pp.start() - nIntFaces,
1182  edgeToDualPoint,
1183  pointToDualPoint,
1184 
1185  dualPoints,
1186 
1187  dynDualFaces,
1188  dynDualOwner,
1189  dynDualNeighbour,
1190  dynDualRegion
1191  );
1192  }
1193 
1194 
1195  // Transfer face info to straight lists
1196  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1197  faceList dualFaces(dynDualFaces.shrink(), true);
1198  dynDualFaces.clear();
1199 
1200  labelList dualOwner(dynDualOwner.shrink(), true);
1201  dynDualOwner.clear();
1202 
1203  labelList dualNeighbour(dynDualNeighbour.shrink(), true);
1204  dynDualNeighbour.clear();
1205 
1206  labelList dualRegion(dynDualRegion.shrink(), true);
1207  dynDualRegion.clear();
1208 
1209 
1210 
1211  // Dump faces.
1212  if (debug)
1213  {
1214  OFstream str("dualFaces.obj");
1215  Pout<< "polyDualMesh::calcDual : dumping all faces to "
1216  << str.name() << endl;
1217 
1218  forAll(dualPoints, dualPointi)
1219  {
1220  meshTools::writeOBJ(str, dualPoints[dualPointi]);
1221  }
1222  forAll(dualFaces, dualFacei)
1223  {
1224  const face& f = dualFaces[dualFacei];
1225 
1226  str<< 'f';
1227  forAll(f, fp)
1228  {
1229  str<< ' ' << f[fp]+1;
1230  }
1231  str<< nl;
1232  }
1233  }
1234 
1235  // Create cells.
1236  cellList dualCells(mesh.nPoints());
1237 
1238  forAll(dualCells, celli)
1239  {
1240  dualCells[celli].setSize(0);
1241  }
1242 
1243  forAll(dualOwner, facei)
1244  {
1245  label celli = dualOwner[facei];
1246 
1247  labelList& cFaces = dualCells[celli];
1248 
1249  label sz = cFaces.size();
1250  cFaces.setSize(sz+1);
1251  cFaces[sz] = facei;
1252  }
1253  forAll(dualNeighbour, facei)
1254  {
1255  label celli = dualNeighbour[facei];
1256 
1257  if (celli != -1)
1258  {
1259  labelList& cFaces = dualCells[celli];
1260 
1261  label sz = cFaces.size();
1262  cFaces.setSize(sz+1);
1263  cFaces[sz] = facei;
1264  }
1265  }
1266 
1267 
1268  // Do upper-triangular face ordering. Determines face reordering map and
1269  // number of internal faces.
1270  label dummy;
1271 
1272  labelList oldToNew
1273  (
1274  getFaceOrder
1275  (
1276  dualOwner,
1277  dualNeighbour,
1278  dualCells,
1279  dummy
1280  )
1281  );
1282 
1283  // Reorder faces.
1284  inplaceReorder(oldToNew, dualFaces);
1285  inplaceReorder(oldToNew, dualOwner);
1286  inplaceReorder(oldToNew, dualNeighbour);
1287  inplaceReorder(oldToNew, dualRegion);
1288  forAll(dualCells, celli)
1289  {
1290  inplaceRenumber(oldToNew, dualCells[celli]);
1291  }
1292 
1293 
1294  // Create patches
1295  labelList patchSizes(patches.size(), 0);
1296 
1297  forAll(dualRegion, facei)
1298  {
1299  if (dualRegion[facei] >= 0)
1300  {
1301  patchSizes[dualRegion[facei]]++;
1302  }
1303  }
1304 
1305  labelList patchStarts(patches.size(), 0);
1306 
1307  label facei = nInternalFaces;
1308 
1309  forAll(patches, patchi)
1310  {
1311  patchStarts[patchi] = facei;
1312 
1313  facei += patchSizes[patchi];
1314  }
1315 
1316 
1317  Pout<< "nFaces:" << dualFaces.size()
1318  << " patchSizes:" << patchSizes
1319  << " patchStarts:" << patchStarts
1320  << endl;
1321 
1322 
1323  // Add patches. First add zero sized (since mesh still 0 size)
1324  List<polyPatch*> dualPatches(patches.size());
1325 
1326  forAll(patches, patchi)
1327  {
1328  const polyPatch& pp = patches[patchi];
1329 
1330  dualPatches[patchi] = pp.clone
1331  (
1332  boundaryMesh(),
1333  patchi,
1334  0, //patchSizes[patchi],
1335  0 //patchStarts[patchi]
1336  ).ptr();
1337  }
1338  addPatches(dualPatches);
1339 
1340  // Assign to mesh.
1341  resetPrimitives
1342  (
1343  xferMove(dualPoints),
1344  xferMove(dualFaces),
1345  xferMove(dualOwner),
1346  xferMove(dualNeighbour),
1347  patchSizes,
1348  patchStarts
1349  );
1350 }
1351 
1352 
1353 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1354 
1355 // Construct from components
1356 Foam::polyDualMesh::polyDualMesh(const IOobject& io)
1357 :
1358  polyMesh(io),
1359  cellPoint_
1360  (
1361  IOobject
1362  (
1363  "cellPoint",
1364  time().findInstance(meshDir(), "cellPoint"),
1365  meshSubDir,
1366  *this,
1367  IOobject::MUST_READ,
1368  IOobject::NO_WRITE
1369  )
1370  ),
1371  boundaryFacePoint_
1372  (
1373  IOobject
1374  (
1375  "boundaryFacePoint",
1376  time().findInstance(meshDir(), "boundaryFacePoint"),
1377  meshSubDir,
1378  *this,
1379  IOobject::MUST_READ,
1380  IOobject::NO_WRITE
1381  )
1382  )
1383 {}
1384 
1385 
1386 // Construct from polyMesh
1387 Foam::polyDualMesh::polyDualMesh
1389  const polyMesh& mesh,
1390  const labelList& featureEdges,
1391  const labelList& featurePoints
1392 )
1393 :
1394  polyMesh
1395  (
1396  mesh,
1397  xferCopy(pointField()),// to prevent any warnings "points not allocated"
1398  xferCopy(faceList()), // to prevent any warnings "faces not allocated"
1399  xferCopy(cellList())
1400  ),
1401  cellPoint_
1402  (
1403  IOobject
1404  (
1405  "cellPoint",
1406  time().findInstance(meshDir(), "faces"),
1407  meshSubDir,
1408  *this,
1411  ),
1412  labelList(mesh.nCells(), -1)
1413  ),
1414  boundaryFacePoint_
1415  (
1416  IOobject
1417  (
1418  "boundaryFacePoint",
1419  time().findInstance(meshDir(), "faces"),
1420  meshSubDir,
1421  *this,
1424  ),
1425  labelList(mesh.nFaces() - mesh.nInternalFaces())
1426  )
1427 {
1428  calcDual(mesh, featureEdges, featurePoints);
1429 }
1430 
1431 
1432 // Construct from polyMesh and feature angle
1433 Foam::polyDualMesh::polyDualMesh
1435  const polyMesh& mesh,
1436  const scalar featureCos
1437 )
1438 :
1439  polyMesh
1440  (
1441  mesh,
1442  xferCopy(pointField()),// to prevent any warnings "points not allocated"
1443  xferCopy(faceList()), // to prevent any warnings "faces not allocated"
1444  xferCopy(cellList())
1445  ),
1446  cellPoint_
1447  (
1448  IOobject
1449  (
1450  "cellPoint",
1451  time().findInstance(meshDir(), "faces"),
1452  meshSubDir,
1453  *this,
1456  ),
1457  labelList(mesh.nCells(), -1)
1458  ),
1459  boundaryFacePoint_
1460  (
1461  IOobject
1462  (
1463  "boundaryFacePoint",
1464  time().findInstance(meshDir(), "faces"),
1465  meshSubDir,
1466  *this,
1469  ),
1470  labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
1471  )
1472 {
1473  labelList featureEdges, featurePoints;
1474 
1475  calcFeatures(mesh, featureCos, featureEdges, featurePoints);
1476  calcDual(mesh, featureEdges, featurePoints);
1477 }
1478 
1479 
1482  const polyMesh& mesh,
1483  const scalar featureCos,
1484  labelList& featureEdges,
1485  labelList& featurePoints
1486 )
1487 {
1488  // Create big primitivePatch for all outside.
1489  primitivePatch allBoundary
1490  (
1492  (
1493  mesh.faces(),
1494  mesh.nFaces() - mesh.nInternalFaces(),
1495  mesh.nInternalFaces()
1496  ),
1497  mesh.points()
1498  );
1499 
1500  // For ease of use store patch number per face in allBoundary.
1501  labelList allRegion(allBoundary.size());
1502 
1503  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1504 
1505  forAll(patches, patchi)
1506  {
1507  const polyPatch& pp = patches[patchi];
1508 
1509  forAll(pp, i)
1510  {
1511  allRegion[i + pp.start() - mesh.nInternalFaces()] = patchi;
1512  }
1513  }
1514 
1515 
1516  // Calculate patch/feature edges
1517  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1518 
1519  const labelListList& edgeFaces = allBoundary.edgeFaces();
1520  const vectorField& faceNormals = allBoundary.faceNormals();
1521  const labelList& meshPoints = allBoundary.meshPoints();
1522 
1523  boolList isFeatureEdge(edgeFaces.size(), false);
1524 
1525  forAll(edgeFaces, edgeI)
1526  {
1527  const labelList& eFaces = edgeFaces[edgeI];
1528 
1529  if (eFaces.size() != 2)
1530  {
1531  // Non-manifold. Problem?
1532  const edge& e = allBoundary.edges()[edgeI];
1533 
1535  << meshPoints[e[0]] << ' ' << meshPoints[e[1]]
1536  << " coords:" << mesh.points()[meshPoints[e[0]]]
1537  << mesh.points()[meshPoints[e[1]]]
1538  << " has more than 2 faces connected to it:"
1539  << eFaces.size() << endl;
1540 
1541  isFeatureEdge[edgeI] = true;
1542  }
1543  else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1544  {
1545  isFeatureEdge[edgeI] = true;
1546  }
1547  else if
1548  (
1549  (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1550  < featureCos
1551  )
1552  {
1553  isFeatureEdge[edgeI] = true;
1554  }
1555  }
1556 
1557 
1558  // Calculate feature points
1559  // ~~~~~~~~~~~~~~~~~~~~~~~~
1560 
1561  const labelListList& pointEdges = allBoundary.pointEdges();
1562 
1563  DynamicList<label> allFeaturePoints(pointEdges.size());
1564 
1565  forAll(pointEdges, pointi)
1566  {
1567  const labelList& pEdges = pointEdges[pointi];
1568 
1569  label nFeatEdges = 0;
1570 
1571  forAll(pEdges, i)
1572  {
1573  if (isFeatureEdge[pEdges[i]])
1574  {
1575  nFeatEdges++;
1576  }
1577  }
1578  if (nFeatEdges > 2)
1579  {
1580  // Store in mesh vertex label
1581  allFeaturePoints.append(allBoundary.meshPoints()[pointi]);
1582  }
1583  }
1584  featurePoints.transfer(allFeaturePoints);
1585 
1586  if (debug)
1587  {
1588  OFstream str("featurePoints.obj");
1589 
1590  Pout<< "polyDualMesh::calcFeatures : dumping feature points to "
1591  << str.name() << endl;
1592 
1593  forAll(featurePoints, i)
1594  {
1595  label pointi = featurePoints[i];
1596  meshTools::writeOBJ(str, mesh.points()[pointi]);
1597  }
1598  }
1599 
1600 
1601  // Get all feature edges.
1602  labelList meshEdges
1603  (
1604  allBoundary.meshEdges
1605  (
1606  mesh.edges(),
1607  mesh.cellEdges(),
1609  (
1610  mesh.faceOwner(),
1611  allBoundary.size(),
1612  mesh.nInternalFaces()
1613  )
1614  )
1615  );
1616 
1617  DynamicList<label> allFeatureEdges(isFeatureEdge.size());
1618  forAll(isFeatureEdge, edgeI)
1619  {
1620  if (isFeatureEdge[edgeI])
1621  {
1622  // Store in mesh edge label.
1623  allFeatureEdges.append(meshEdges[edgeI]);
1624  }
1625  }
1626  featureEdges.transfer(allFeatureEdges);
1627 }
1628 
1629 
1630 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1631 
1633 {}
1634 
1635 
1636 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
const labelListList & cellEdges() const
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nInternalFaces() const
const labelListList & pointEdges() const
label nFaces() const
Output to file stream.
Definition: OFstream.H:82
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:312
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
static void calcFeatures(const polyMesh &, const scalar featureCos, labelList &featureEdges, labelList &featurePoints)
Helper function to create feature edges and points based on.
List< face > faceList
Definition: faceListFwd.H:43
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Xfer< T > xferCopy(const T &)
Construct by copying the contents of the arg.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:120
scalar f1
Definition: createFields.H:28
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1011
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
A list of faces which address into the list of points.
A List obtained as a section of another List.
Definition: SubList.H:53
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1049
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1036
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
static const char nl
Definition: Ostream.H:262
const Time & time() const
Return time.
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
~polyDualMesh()
Destructor.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
fileName meshDir() const
Return the local mesh directory (dbDir()/meshSubDir)
Definition: polyMesh.C:789
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
List< cell > cellList
list of cells
Definition: cellList.H:42
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
const labelListList & edgeFaces() const
Namespace for OpenFOAM.