FreeFOAM The Cross-Platform CFD Toolkit
boundaryMesh.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 "boundaryMesh.H"
27 #include <OpenFOAM/Time.H>
28 #include <OpenFOAM/polyMesh.H>
30 #include <OpenFOAM/faceList.H>
31 #include <meshTools/octree.H>
32 #include "octreeDataFaceList.H"
33 #include <triSurface/triSurface.H>
34 #include <OpenFOAM/SortableList.H>
35 #include <OpenFOAM/OFstream.H>
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
40 
41 // Normal along which to divide faces into categories (used in getNearest)
42 const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
43 
44 // Distance to face tolerance for getNearest
45 const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1E-2;
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 // Returns number of feature edges connected to pointI
51 Foam::label Foam::boundaryMesh::nFeatureEdges(label pointI) const
52 {
53  label nFeats = 0;
54 
55  const labelList& pEdges = mesh().pointEdges()[pointI];
56 
57  forAll(pEdges, pEdgeI)
58  {
59  label edgeI = pEdges[pEdgeI];
60 
61  if (edgeToFeature_[edgeI] != -1)
62  {
63  nFeats++;
64  }
65  }
66  return nFeats;
67 }
68 
69 
70 // Returns next feature edge connected to pointI
71 Foam::label Foam::boundaryMesh::nextFeatureEdge
72 (
73  const label edgeI,
74  const label vertI
75 ) const
76 {
77  const labelList& pEdges = mesh().pointEdges()[vertI];
78 
79  forAll(pEdges, pEdgeI)
80  {
81  label nbrEdgeI = pEdges[pEdgeI];
82 
83  if (nbrEdgeI != edgeI)
84  {
85  label featI = edgeToFeature_[nbrEdgeI];
86 
87  if (featI != -1)
88  {
89  return nbrEdgeI;
90  }
91  }
92  }
93 
94  return -1;
95 }
96 
97 
98 // Finds connected feature edges, starting from startPointI and returns
99 // feature labels (not edge labels). Marks feature edges handled in
100 // featVisited.
101 Foam::labelList Foam::boundaryMesh::collectSegment
102 (
103  const boolList& isFeaturePoint,
104  const label startEdgeI,
105  boolList& featVisited
106 ) const
107 {
108  // Find starting feature point on edge.
109 
110  label edgeI = startEdgeI;
111 
112  const edge& e = mesh().edges()[edgeI];
113 
114  label vertI = e.start();
115 
116  while (!isFeaturePoint[vertI])
117  {
118  // Step to next feature edge
119 
120  edgeI = nextFeatureEdge(edgeI, vertI);
121 
122  if ((edgeI == -1) || (edgeI == startEdgeI))
123  {
124  break;
125  }
126 
127  // Step to next vertex on edge
128 
129  const edge& e = mesh().edges()[edgeI];
130 
131  vertI = e.otherVertex(vertI);
132  }
133 
134  //
135  // Now we have:
136  // edgeI : first edge on this segment
137  // vertI : one of the endpoints of this segment
138  //
139  // Start walking other way and storing edges as we go along.
140  //
141 
142  // Untrimmed storage for current segment
143  labelList featLabels(featureEdges_.size());
144 
145  label featLabelI = 0;
146 
147  label initEdgeI = edgeI;
148 
149  do
150  {
151  // Mark edge as visited
152  label featI = edgeToFeature_[edgeI];
153 
154  if (featI == -1)
155  {
156  FatalErrorIn("boundaryMesh::collectSegment")
157  << "Problem" << abort(FatalError);
158  }
159  featLabels[featLabelI++] = featI;
160 
161  featVisited[featI] = true;
162 
163  // Step to next vertex on edge
164 
165  const edge& e = mesh().edges()[edgeI];
166 
167  vertI = e.otherVertex(vertI);
168 
169  // Step to next feature edge
170 
171  edgeI = nextFeatureEdge(edgeI, vertI);
172 
173  if ((edgeI == -1) || (edgeI == initEdgeI))
174  {
175  break;
176  }
177  }
178  while (!isFeaturePoint[vertI]);
179 
180 
181  // Trim to size
182  featLabels.setSize(featLabelI);
183 
184  return featLabels;
185 }
186 
187 
188 void Foam::boundaryMesh::markEdges
189 (
190  const label maxDistance,
191  const label edgeI,
192  const label distance,
193  labelList& minDistance,
194  DynamicList<label>& visited
195 ) const
196 {
197  if (distance < maxDistance)
198  {
199  // Don't do anything if reached beyond maxDistance.
200 
201  if (minDistance[edgeI] == -1)
202  {
203  // First visit of edge. Store edge label.
204  visited.append(edgeI);
205  }
206  else if (minDistance[edgeI] <= distance)
207  {
208  // Already done this edge
209  return;
210  }
211 
212  minDistance[edgeI] = distance;
213 
214  const edge& e = mesh().edges()[edgeI];
215 
216  // Do edges connected to e.start
217  const labelList& startEdges = mesh().pointEdges()[e.start()];
218 
219  forAll(startEdges, pEdgeI)
220  {
221  markEdges
222  (
223  maxDistance,
224  startEdges[pEdgeI],
225  distance+1,
226  minDistance,
227  visited
228  );
229  }
230 
231  // Do edges connected to e.end
232  const labelList& endEdges = mesh().pointEdges()[e.end()];
233 
234  forAll(endEdges, pEdgeI)
235  {
236  markEdges
237  (
238  maxDistance,
239  endEdges[pEdgeI],
240  distance+1,
241  minDistance,
242  visited
243  );
244  }
245  }
246 }
247 
248 
249 Foam::label Foam::boundaryMesh::findPatchID
250 (
251  const polyPatchList& patches,
252  const word& patchName
253 ) const
254 {
255  forAll(patches, patchI)
256  {
257  if (patches[patchI].name() == patchName)
258  {
259  return patchI;
260  }
261  }
262 
263  return -1;
264 }
265 
266 
268 {
269  wordList names(patches_.size());
270 
271  forAll(patches_, patchI)
272  {
273  names[patchI] = patches_[patchI].name();
274  }
275  return names;
276 }
277 
278 
279 Foam::label Foam::boundaryMesh::whichPatch
280 (
281  const polyPatchList& patches,
282  const label faceI
283 ) const
284 {
285  forAll(patches, patchI)
286  {
287  const polyPatch& pp = patches[patchI];
288 
289  if ((faceI >= pp.start()) && (faceI < (pp.start() + pp.size())))
290  {
291  return patchI;
292  }
293  }
294  return -1;
295 }
296 
297 
298 // Gets labels of changed faces and propagates them to the edges. Returns
299 // labels of edges changed.
300 Foam::labelList Foam::boundaryMesh::faceToEdge
301 (
302  const boolList& regionEdge,
303  const label region,
304  const labelList& changedFaces,
305  labelList& edgeRegion
306 ) const
307 {
308  labelList changedEdges(mesh().nEdges(), -1);
309  label changedI = 0;
310 
311  forAll(changedFaces, i)
312  {
313  label faceI = changedFaces[i];
314 
315  const labelList& fEdges = mesh().faceEdges()[faceI];
316 
317  forAll(fEdges, fEdgeI)
318  {
319  label edgeI = fEdges[fEdgeI];
320 
321  if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
322  {
323  edgeRegion[edgeI] = region;
324 
325  changedEdges[changedI++] = edgeI;
326  }
327  }
328  }
329 
330  changedEdges.setSize(changedI);
331 
332  return changedEdges;
333 }
334 
335 
336 // Reverse of faceToEdge: gets edges and returns faces
337 Foam::labelList Foam::boundaryMesh::edgeToFace
338 (
339  const label region,
340  const labelList& changedEdges,
341  labelList& faceRegion
342 ) const
343 {
344  labelList changedFaces(mesh().size(), -1);
345  label changedI = 0;
346 
347  forAll(changedEdges, i)
348  {
349  label edgeI = changedEdges[i];
350 
351  const labelList& eFaces = mesh().edgeFaces()[edgeI];
352 
353  forAll(eFaces, eFaceI)
354  {
355  label faceI = eFaces[eFaceI];
356 
357  if (faceRegion[faceI] == -1)
358  {
359  faceRegion[faceI] = region;
360 
361  changedFaces[changedI++] = faceI;
362  }
363  }
364  }
365 
366  changedFaces.setSize(changedI);
367 
368  return changedFaces;
369 }
370 
371 
372 // Finds area, starting at faceI, delimited by borderEdge
373 void Foam::boundaryMesh::markZone
374 (
375  const boolList& borderEdge,
376  label faceI,
377  label currentZone,
378  labelList& faceZone
379 ) const
380 {
381  faceZone[faceI] = currentZone;
382 
383  // List of faces whose faceZone has been set.
384  labelList changedFaces(1, faceI);
385  // List of edges whose faceZone has been set.
386  labelList changedEdges;
387 
388  // Zones on all edges.
389  labelList edgeZone(mesh().nEdges(), -1);
390 
391  while(1)
392  {
393  changedEdges = faceToEdge
394  (
395  borderEdge,
396  currentZone,
397  changedFaces,
398  edgeZone
399  );
400 
401  if (debug)
402  {
403  Pout<< "From changedFaces:" << changedFaces.size()
404  << " to changedEdges:" << changedEdges.size()
405  << endl;
406  }
407 
408  if (changedEdges.empty())
409  {
410  break;
411  }
412 
413  changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
414 
415  if (debug)
416  {
417  Pout<< "From changedEdges:" << changedEdges.size()
418  << " to changedFaces:" << changedFaces.size()
419  << endl;
420  }
421 
422  if (changedFaces.empty())
423  {
424  break;
425  }
426  }
427 }
428 
429 
430 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
431 
432 // Null constructor
434 :
435  meshPtr_(NULL),
436  patches_(),
437  meshFace_(),
438  featurePoints_(),
439  featureEdges_(),
440  featureToEdge_(),
441  edgeToFeature_(),
442  featureSegments_(),
443  extraEdges_()
444 {}
445 
446 
447 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
448 
450 {
451  clearOut();
452 }
453 
454 
456 {
457  if (meshPtr_)
458  {
459  delete meshPtr_;
460 
461  meshPtr_ = NULL;
462  }
463 }
464 
465 
466 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
467 
469 {
470  patches_.clear();
471 
472  patches_.setSize(mesh.boundaryMesh().size());
473 
474  // Number of boundary faces
475  label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
476 
477  faceList bFaces(nBFaces);
478 
479  meshFace_.setSize(nBFaces);
480 
481  label bFaceI = 0;
482 
483  // Collect all boundary faces.
484  forAll(mesh.boundaryMesh(), patchI)
485  {
486  const polyPatch& pp = mesh.boundaryMesh()[patchI];
487 
488  patches_.set
489  (
490  patchI,
491  new boundaryPatch
492  (
493  pp.name(),
494  patchI,
495  pp.size(),
496  bFaceI,
497  pp.type()
498  )
499  );
500 
501  // Collect all faces in global numbering.
502  forAll(pp, patchFaceI)
503  {
504  meshFace_[bFaceI] = pp.start() + patchFaceI;
505 
506  bFaces[bFaceI] = pp[patchFaceI];
507 
508  bFaceI++;
509  }
510  }
511 
512 
513  if (debug)
514  {
515  Pout<< "read : patches now:" << endl;
516 
517  forAll(patches_, patchI)
518  {
519  const boundaryPatch& bp = patches_[patchI];
520 
521  Pout<< " name : " << bp.name() << endl
522  << " size : " << bp.size() << endl
523  << " start : " << bp.start() << endl
524  << " type : " << bp.physicalType() << endl
525  << endl;
526  }
527  }
528 
529  //
530  // Construct single patch for all of boundary
531  //
532 
533  // Temporary primitivePatch to calculate compact points & faces.
535  (
536  bFaces,
537  mesh.points()
538  );
539 
540  // Store in local(compact) addressing
541  clearOut();
542 
543  meshPtr_ = new bMesh(globalPatch.localFaces(), globalPatch.localPoints());
544 
545 
546  if (debug & 2)
547  {
548  const bMesh& msh = *meshPtr_;
549 
550  Pout<< "** Start of Faces **" << endl;
551 
552  forAll(msh, faceI)
553  {
554  const face& f = msh[faceI];
555 
556  point ctr(vector::zero);
557 
558  forAll(f, fp)
559  {
560  ctr += msh.points()[f[fp]];
561  }
562  ctr /= f.size();
563 
564  Pout<< " " << faceI
565  << " ctr:" << ctr
566  << " verts:" << f
567  << endl;
568  }
569 
570  Pout<< "** End of Faces **" << endl;
571 
572  Pout<< "** Start of Points **" << endl;
573 
574  forAll(msh.points(), pointI)
575  {
576  Pout<< " " << pointI
577  << " coord:" << msh.points()[pointI]
578  << endl;
579  }
580 
581  Pout<< "** End of Points **" << endl;
582  }
583 
584  // Clear edge storage
585  featurePoints_.setSize(0);
586  featureEdges_.setSize(0);
587 
588  featureToEdge_.setSize(0);
589  edgeToFeature_.setSize(meshPtr_->nEdges());
590  edgeToFeature_ = -1;
591 
592  featureSegments_.setSize(0);
593 
594  extraEdges_.setSize(0);
595 }
596 
597 
599 {
600  triSurface surf(fName);
601 
602  if (surf.empty())
603  {
604  return;
605  }
606 
607  // Sort according to region
608  SortableList<label> regions(surf.size());
609 
610  forAll(surf, triI)
611  {
612  regions[triI] = surf[triI].region();
613  }
614  regions.sort();
615 
616  // Determine region mapping.
617  Map<label> regionToBoundaryPatch;
618 
619  label oldRegion = -1111;
620  label boundPatch = 0;
621 
622  forAll(regions, i)
623  {
624  if (regions[i] != oldRegion)
625  {
626  regionToBoundaryPatch.insert(regions[i], boundPatch);
627 
628  oldRegion = regions[i];
629  boundPatch++;
630  }
631  }
632 
633  const geometricSurfacePatchList& surfPatches = surf.patches();
634 
635  patches_.clear();
636 
637  if (surfPatches.size() == regionToBoundaryPatch.size())
638  {
639  // There are as many surface patches as region numbers in triangles
640  // so use the surface patches
641 
642  patches_.setSize(surfPatches.size());
643 
644  // Take over patches, setting size to 0 for now.
645  forAll(surfPatches, patchI)
646  {
647  const geometricSurfacePatch& surfPatch = surfPatches[patchI];
648 
649  patches_.set
650  (
651  patchI,
652  new boundaryPatch
653  (
654  surfPatch.name(),
655  patchI,
656  0,
657  0,
658  surfPatch.geometricType()
659  )
660  );
661  }
662  }
663  else
664  {
665  // There are not enough surface patches. Make up my own.
666 
667  patches_.setSize(regionToBoundaryPatch.size());
668 
669  forAll(patches_, patchI)
670  {
671  patches_.set
672  (
673  patchI,
674  new boundaryPatch
675  (
676  "patch" + name(patchI),
677  patchI,
678  0,
679  0,
680  "empty"
681  )
682  );
683  }
684  }
685 
686  //
687  // Copy according into bFaces according to regions
688  //
689 
690  const labelList& indices = regions.indices();
691 
692  faceList bFaces(surf.size());
693 
694  meshFace_.setSize(surf.size());
695 
696  label bFaceI = 0;
697 
698  // Current region number
699  label surfRegion = regions[0];
700  label foamRegion = regionToBoundaryPatch[surfRegion];
701 
702  Pout<< "Surface region " << surfRegion << " becomes boundary patch "
703  << foamRegion << " with name " << patches_[foamRegion].name() << endl;
704 
705 
706  // Index in bFaces of start of current patch
707  label startFaceI = 0;
708 
709  forAll(indices, indexI)
710  {
711  label triI = indices[indexI];
712 
713  const labelledTri& tri = surf.localFaces()[triI];
714 
715  if (tri.region() != surfRegion)
716  {
717  // Change of region. We now know the size of the previous one.
718  boundaryPatch& bp = patches_[foamRegion];
719 
720  bp.size() = bFaceI - startFaceI;
721  bp.start() = startFaceI;
722 
723  surfRegion = tri.region();
724  foamRegion = regionToBoundaryPatch[surfRegion];
725 
726  Pout<< "Surface region " << surfRegion << " becomes boundary patch "
727  << foamRegion << " with name " << patches_[foamRegion].name()
728  << endl;
729 
730  startFaceI = bFaceI;
731  }
732 
733  meshFace_[bFaceI] = triI;
734 
735  bFaces[bFaceI++] = face(tri);
736  }
737 
738  // Final region
739  boundaryPatch& bp = patches_[foamRegion];
740 
741  bp.size() = bFaceI - startFaceI;
742  bp.start() = startFaceI;
743 
744  //
745  // Construct single primitivePatch for all of boundary
746  //
747 
748  clearOut();
749 
750  // Store compact.
751  meshPtr_ = new bMesh(bFaces, surf.localPoints());
752 
753  // Clear edge storage
754  featurePoints_.setSize(0);
755  featureEdges_.setSize(0);
756 
757  featureToEdge_.setSize(0);
758  edgeToFeature_.setSize(meshPtr_->nEdges());
759  edgeToFeature_ = -1;
760 
761  featureSegments_.setSize(0);
762 
763  extraEdges_.setSize(0);
764 }
765 
766 
768 {
769  geometricSurfacePatchList surfPatches(patches_.size());
770 
771  forAll(patches_, patchI)
772  {
773  const boundaryPatch& bp = patches_[patchI];
774 
775  surfPatches[patchI] =
777  (
778  bp.physicalType(),
779  bp.name(),
780  patchI
781  );
782  }
783 
784  //
785  // Simple triangulation.
786  //
787 
788  // Get number of triangles per face
789  labelList nTris(mesh().size());
790 
791  label totalNTris = getNTris(0, mesh().size(), nTris);
792 
793  // Determine per face the starting triangle.
794  labelList startTri(mesh().size());
795 
796  label triI = 0;
797 
798  forAll(mesh(), faceI)
799  {
800  startTri[faceI] = triI;
801 
802  triI += nTris[faceI];
803  }
804 
805  // Triangulate
806  labelList triVerts(3*totalNTris);
807 
808  triangulate(0, mesh().size(), totalNTris, triVerts);
809 
810  // Convert to labelledTri
811 
812  List<labelledTri> tris(totalNTris);
813 
814  triI = 0;
815 
816  forAll(patches_, patchI)
817  {
818  const boundaryPatch& bp = patches_[patchI];
819 
820  forAll(bp, patchFaceI)
821  {
822  label faceI = bp.start() + patchFaceI;
823 
824  label triVertI = 3*startTri[faceI];
825 
826  for (label faceTriI = 0; faceTriI < nTris[faceI]; faceTriI++)
827  {
828  label v0 = triVerts[triVertI++];
829  label v1 = triVerts[triVertI++];
830  label v2 = triVerts[triVertI++];
831 
832  tris[triI++] = labelledTri(v0, v1, v2, patchI);
833  }
834  }
835  }
836 
837  triSurface surf(tris, surfPatches, mesh().points());
838 
839  OFstream surfStream(fName);
840 
841  surf.write(surfStream);
842 }
843 
844 
845 // Get index in this (boundaryMesh) of face nearest to each boundary face in
846 // pMesh.
847 // Origininally all triangles/faces of boundaryMesh would be bunged into
848 // one big octree. Problem was that faces on top of each other, differing
849 // only in sign of normal, could not be found separately. It would always
850 // find only one. We could detect that it was probably finding the wrong one
851 // (based on normal) but could not 'tell' the octree to retrieve the other
852 // one (since they occupy exactly the same space)
853 // So now faces get put into different octrees depending on normal.
854 // !It still will not be possible to differentiate between two faces on top
855 // of each other having the same normal
857 (
858  const primitiveMesh& pMesh,
859  const vector& searchSpan
860 ) const
861 {
862 
863  // Divide faces into two bins acc. to normal
864  // - left of splitNormal
865  // - right ,,
866  DynamicList<label> leftFaces(mesh().size()/2);
867  DynamicList<label> rightFaces(mesh().size()/2);
868 
869  forAll(mesh(), bFaceI)
870  {
871  scalar sign = mesh().faceNormals()[bFaceI] & splitNormal_;
872 
873  if (sign > -1E-5)
874  {
875  rightFaces.append(bFaceI);
876  }
877  if (sign < 1E-5)
878  {
879  leftFaces.append(bFaceI);
880  }
881  }
882 
883  leftFaces.shrink();
884  rightFaces.shrink();
885 
886  if (debug)
887  {
888  Pout<< "getNearest :"
889  << " rightBin:" << rightFaces.size()
890  << " leftBin:" << leftFaces.size()
891  << endl;
892  }
893 
894 
895  // Overall bb
896  treeBoundBox overallBb(mesh().localPoints());
897 
898  // Extend domain slightly (also makes it 3D if was 2D)
899  // Note asymmetry to avoid having faces align with octree cubes.
900  scalar tol = 1E-6 * overallBb.avgDim();
901 
902  point& bbMin = overallBb.min();
903  bbMin.x() -= tol;
904  bbMin.y() -= tol;
905  bbMin.z() -= tol;
906 
907  point& bbMax = overallBb.max();
908  bbMax.x() += 2*tol;
909  bbMax.y() += 2*tol;
910  bbMax.z() += 2*tol;
911 
912  // Create the octrees
914  (
915  overallBb,
917  (
918  mesh(),
919  leftFaces
920  ),
921  1,
922  20,
923  10
924  );
926  (
927  overallBb,
929  (
930  mesh(),
931  rightFaces
932  ),
933  1,
934  20,
935  10
936  );
937 
938  if (debug)
939  {
940  Pout<< "getNearest : built trees" << endl;
941  }
942 
943 
944  const vectorField& ns = mesh().faceNormals();
945 
946 
947  //
948  // Search nearest triangle centre for every polyMesh boundary face
949  //
950 
951  labelList nearestBFaceI(pMesh.nFaces() - pMesh.nInternalFaces());
952 
953  treeBoundBox tightest;
954 
955  const scalar searchDim = mag(searchSpan);
956 
957  forAll(nearestBFaceI, patchFaceI)
958  {
959  label meshFaceI = pMesh.nInternalFaces() + patchFaceI;
960 
961  const point& ctr = pMesh.faceCentres()[meshFaceI];
962 
963  if (debug && (patchFaceI % 1000) == 0)
964  {
965  Pout<< "getNearest : patchFace:" << patchFaceI
966  << " meshFaceI:" << meshFaceI << " ctr:" << ctr << endl;
967  }
968 
969 
970  // Get normal from area vector
971  vector n = pMesh.faceAreas()[meshFaceI];
972  scalar area = mag(n);
973  n /= area;
974 
975  scalar typDim = -GREAT;
976  const face& f = pMesh.faces()[meshFaceI];
977 
978  forAll(f, fp)
979  {
980  typDim = max(typDim, mag(pMesh.points()[f[fp]] - ctr));
981  }
982 
983  // Search right tree
984  tightest.min() = ctr - searchSpan;
985  tightest.max() = ctr + searchSpan;
986  scalar rightDist = searchDim;
987  label rightI = rightTree.findNearest(ctr, tightest, rightDist);
988 
989 
990  // Search left tree. Note: could start from rightDist bounding box
991  // instead of starting from top.
992  tightest.min() = ctr - searchSpan;
993  tightest.max() = ctr + searchSpan;
994  scalar leftDist = searchDim;
995  label leftI = leftTree.findNearest(ctr, tightest, leftDist);
996 
997 
998  if (rightI == -1)
999  {
1000  // No face found in right tree.
1001 
1002  if (leftI == -1)
1003  {
1004  // No face found in left tree.
1005  nearestBFaceI[patchFaceI] = -1;
1006  }
1007  else
1008  {
1009  // Found in left but not in right. Choose left regardless if
1010  // correct sign. Note: do we want this?
1011  nearestBFaceI[patchFaceI] = leftFaces[leftI];
1012  }
1013  }
1014  else
1015  {
1016  if (leftI == -1)
1017  {
1018  // Found in right but not in left. Choose right regardless if
1019  // correct sign. Note: do we want this?
1020  nearestBFaceI[patchFaceI] = rightFaces[rightI];
1021  }
1022  else
1023  {
1024  // Found in both trees. Compare normals.
1025 
1026  scalar rightSign = n & ns[rightFaces[rightI]];
1027  scalar leftSign = n & ns[leftFaces[leftI]];
1028 
1029  if
1030  (
1031  (rightSign > 0 && leftSign > 0)
1032  || (rightSign < 0 && leftSign < 0)
1033  )
1034  {
1035  // Both same sign. Choose nearest.
1036  if (rightDist < leftDist)
1037  {
1038  nearestBFaceI[patchFaceI] = rightFaces[rightI];
1039  }
1040  else
1041  {
1042  nearestBFaceI[patchFaceI] = leftFaces[leftI];
1043  }
1044  }
1045  else
1046  {
1047  // Differing sign.
1048  // - if both near enough choose one with correct sign
1049  // - otherwise choose nearest.
1050 
1051  // Get local dimension as max of distance between ctr and
1052  // any face vertex.
1053 
1054  typDim *= distanceTol_;
1055 
1056  if (rightDist < typDim && leftDist < typDim)
1057  {
1058  // Different sign and nearby. Choosing matching normal
1059  if (rightSign > 0)
1060  {
1061  nearestBFaceI[patchFaceI] = rightFaces[rightI];
1062  }
1063  else
1064  {
1065  nearestBFaceI[patchFaceI] = leftFaces[leftI];
1066  }
1067  }
1068  else
1069  {
1070  // Different sign but faraway. Choosing nearest.
1071  if (rightDist < leftDist)
1072  {
1073  nearestBFaceI[patchFaceI] = rightFaces[rightI];
1074  }
1075  else
1076  {
1077  nearestBFaceI[patchFaceI] = leftFaces[leftI];
1078  }
1079  }
1080  }
1081  }
1082  }
1083  }
1084 
1085  return nearestBFaceI;
1086 }
1087 
1088 
1091  const labelList& nearest,
1092  const polyBoundaryMesh& oldPatches,
1093  polyMesh& newMesh
1094 ) const
1095 {
1096 
1097  // 2 cases to be handled:
1098  // A- patches in boundaryMesh patches_
1099  // B- patches not in boundaryMesh patches_ but in polyMesh
1100 
1101  // Create maps from patch name to new patch index.
1102  HashTable<label> nameToIndex(2*patches_.size());
1103 
1104  Map<word> indexToName(2*patches_.size());
1105 
1106 
1107  label nNewPatches = patches_.size();
1108 
1109  forAll(oldPatches, oldPatchI)
1110  {
1111  const polyPatch& patch = oldPatches[oldPatchI];
1112 
1113  label newPatchI = findPatchID(patch.name());
1114 
1115  if (newPatchI != -1)
1116  {
1117  nameToIndex.insert(patch.name(), newPatchI);
1118 
1119  indexToName.insert(newPatchI, patch.name());
1120  }
1121  }
1122 
1123  // Include all boundaryPatches not yet in nameToIndex (i.e. not in old
1124  // patches)
1125  forAll(patches_, bPatchI)
1126  {
1127  const boundaryPatch& bp = patches_[bPatchI];
1128 
1129  if (!nameToIndex.found(bp.name()))
1130  {
1131  nameToIndex.insert(bp.name(), bPatchI);
1132 
1133  indexToName.insert(bPatchI, bp.name());
1134  }
1135  }
1136 
1137  // Pass1:
1138  // Copy names&type of patches (with zero size) from old mesh as far as
1139  // possible. First patch created gets all boundary faces; all others get
1140  // zero faces (repatched later on). Exception is coupled patches which
1141  // keep their size.
1142 
1143  List<polyPatch*> newPatchPtrList(nNewPatches);
1144 
1145  label meshFaceI = newMesh.nInternalFaces();
1146 
1147  // First patch gets all non-coupled faces
1148  label facesToBeDone = newMesh.nFaces() - newMesh.nInternalFaces();
1149 
1150  forAll(patches_, bPatchI)
1151  {
1152  const boundaryPatch& bp = patches_[bPatchI];
1153 
1154  label newPatchI = nameToIndex[bp.name()];
1155 
1156  // Find corresponding patch in polyMesh
1157  label oldPatchI = findPatchID(oldPatches, bp.name());
1158 
1159  if (oldPatchI == -1)
1160  {
1161  // Newly created patch. Gets all or zero faces.
1162  if (debug)
1163  {
1164  Pout<< "patchify : Creating new polyPatch:" << bp.name()
1165  << " type:" << bp.physicalType() << endl;
1166  }
1167 
1168  newPatchPtrList[newPatchI] = polyPatch::New
1169  (
1170  bp.physicalType(),
1171  bp.name(),
1172  facesToBeDone,
1173  meshFaceI,
1174  newPatchI,
1175  newMesh.boundaryMesh()
1176  ).ptr();
1177 
1178  meshFaceI += facesToBeDone;
1179 
1180  // first patch gets all boundary faces; all others get 0.
1181  facesToBeDone = 0;
1182  }
1183  else
1184  {
1185  // Existing patch. Gets all or zero faces.
1186  const polyPatch& oldPatch = oldPatches[oldPatchI];
1187 
1188  if (debug)
1189  {
1190  Pout<< "patchify : Cloning existing polyPatch:"
1191  << oldPatch.name() << endl;
1192  }
1193 
1194  newPatchPtrList[newPatchI] = oldPatch.clone
1195  (
1196  newMesh.boundaryMesh(),
1197  newPatchI,
1198  facesToBeDone,
1199  meshFaceI
1200  ).ptr();
1201 
1202  meshFaceI += facesToBeDone;
1203 
1204  // first patch gets all boundary faces; all others get 0.
1205  facesToBeDone = 0;
1206  }
1207  }
1208 
1209 
1210  if (debug)
1211  {
1212  Pout<< "Patchify : new polyPatch list:" << endl;
1213 
1214  forAll(newPatchPtrList, patchI)
1215  {
1216  const polyPatch& newPatch = *newPatchPtrList[patchI];
1217 
1218  if (debug)
1219  {
1220  Pout<< "polyPatch:" << newPatch.name() << endl
1221  << " type :" << newPatch.typeName << endl
1222  << " size :" << newPatch.size() << endl
1223  << " start:" << newPatch.start() << endl
1224  << " index:" << patchI << endl;
1225  }
1226  }
1227  }
1228 
1229  // Actually add new list of patches
1230  repatchPolyTopoChanger polyMeshRepatcher(newMesh);
1231  polyMeshRepatcher.changePatches(newPatchPtrList);
1232 
1233 
1234  // Pass2:
1235  // Change patch type for face
1236 
1237  if (newPatchPtrList.size())
1238  {
1239  List<DynamicList<label> > patchFaces(nNewPatches);
1240 
1241  // Give reasonable estimate for size of patches
1242  label nAvgFaces =
1243  (newMesh.nFaces() - newMesh.nInternalFaces())
1244  / nNewPatches;
1245 
1246  forAll(patchFaces, newPatchI)
1247  {
1248  patchFaces[newPatchI].setCapacity(nAvgFaces);
1249  }
1250 
1251  //
1252  // Sort faces acc. to new patch index. Can loop over all old patches
1253  // since will contain all faces.
1254  //
1255 
1256  forAll(oldPatches, oldPatchI)
1257  {
1258  const polyPatch& patch = oldPatches[oldPatchI];
1259 
1260  forAll(patch, patchFaceI)
1261  {
1262  // Put face into region given by nearest boundary face
1263 
1264  label meshFaceI = patch.start() + patchFaceI;
1265 
1266  label bFaceI = meshFaceI - newMesh.nInternalFaces();
1267 
1268  patchFaces[whichPatch(nearest[bFaceI])].append(meshFaceI);
1269  }
1270  }
1271 
1272  forAll(patchFaces, newPatchI)
1273  {
1274  patchFaces[newPatchI].shrink();
1275  }
1276 
1277 
1278  // Change patch > 0. (since above we put all faces into the zeroth
1279  // patch)
1280 
1281  for (label newPatchI = 1; newPatchI < patchFaces.size(); newPatchI++)
1282  {
1283  const labelList& pFaces = patchFaces[newPatchI];
1284 
1285  forAll(pFaces, pFaceI)
1286  {
1287  polyMeshRepatcher.changePatchID(pFaces[pFaceI], newPatchI);
1288  }
1289  }
1290 
1291  polyMeshRepatcher.repatch();
1292  }
1293 }
1294 
1295 
1296 void Foam::boundaryMesh::setFeatureEdges(const scalar minCos)
1297 {
1298  edgeToFeature_.setSize(mesh().nEdges());
1299 
1300  edgeToFeature_ = -1;
1301 
1302  // 1. Mark feature edges
1303 
1304  // Storage for edge labels that are features. Trim later.
1305  featureToEdge_.setSize(mesh().nEdges());
1306 
1307  label featureI = 0;
1308 
1309  if (minCos >= 0.9999)
1310  {
1311  // Select everything
1312  forAll(mesh().edges(), edgeI)
1313  {
1314  edgeToFeature_[edgeI] = featureI;
1315  featureToEdge_[featureI++] = edgeI;
1316  }
1317  }
1318  else
1319  {
1320  forAll(mesh().edges(), edgeI)
1321  {
1322  const labelList& eFaces = mesh().edgeFaces()[edgeI];
1323 
1324  if (eFaces.size() == 2)
1325  {
1326  label face0I = eFaces[0];
1327 
1328  label face1I = eFaces[1];
1329 
1332  //if (whichPatch(face0I) != whichPatch(face1I))
1333  //{
1334  // edgeToFeature_[edgeI] = featureI;
1335  // featureToEdge_[featureI++] = edgeI;
1336  //}
1337  //else
1338  {
1339  const vector& n0 = mesh().faceNormals()[face0I];
1340 
1341  const vector& n1 = mesh().faceNormals()[face1I];
1342 
1343  float cosAng = n0 & n1;
1344 
1345  if (cosAng < minCos)
1346  {
1347  edgeToFeature_[edgeI] = featureI;
1348  featureToEdge_[featureI++] = edgeI;
1349  }
1350  }
1351  }
1352  else
1353  {
1354  //Should not occur: 0 or more than two faces
1355  edgeToFeature_[edgeI] = featureI;
1356  featureToEdge_[featureI++] = edgeI;
1357  }
1358  }
1359  }
1360 
1361  // Trim featureToEdge_ to actual number of edges.
1362  featureToEdge_.setSize(featureI);
1363 
1364  //
1365  // Compact edges i.e. relabel vertices.
1366  //
1367 
1368  featureEdges_.setSize(featureI);
1369  featurePoints_.setSize(mesh().nPoints());
1370 
1371  labelList featToMeshPoint(mesh().nPoints(), -1);
1372 
1373  label featPtI = 0;
1374 
1375  forAll(featureToEdge_, fEdgeI)
1376  {
1377  label edgeI = featureToEdge_[fEdgeI];
1378 
1379  const edge& e = mesh().edges()[edgeI];
1380 
1381  label start = featToMeshPoint[e.start()];
1382 
1383  if (start == -1)
1384  {
1385  featToMeshPoint[e.start()] = featPtI;
1386 
1387  featurePoints_[featPtI] = mesh().points()[e.start()];
1388 
1389  start = featPtI;
1390 
1391  featPtI++;
1392  }
1393 
1394  label end = featToMeshPoint[e.end()];
1395 
1396  if (end == -1)
1397  {
1398  featToMeshPoint[e.end()] = featPtI;
1399 
1400  featurePoints_[featPtI] = mesh().points()[e.end()];
1401 
1402  end = featPtI;
1403 
1404  featPtI++;
1405  }
1406 
1407  // Store with renumbered vertices.
1408  featureEdges_[fEdgeI] = edge(start, end);
1409  }
1410 
1411  // Compact points
1412  featurePoints_.setSize(featPtI);
1413 
1414 
1415  //
1416  // 2. Mark endpoints of feature segments. These are points with
1417  // != 2 feature edges connected.
1418  // Note: can add geometric constraint here as well that if 2 feature
1419  // edges the angle between them should be less than xxx.
1420  //
1421 
1422  boolList isFeaturePoint(mesh().nPoints(), false);
1423 
1424  forAll(featureToEdge_, featI)
1425  {
1426  label edgeI = featureToEdge_[featI];
1427 
1428  const edge& e = mesh().edges()[edgeI];
1429 
1430  if (nFeatureEdges(e.start()) != 2)
1431  {
1432  isFeaturePoint[e.start()] = true;
1433  }
1434 
1435  if (nFeatureEdges(e.end()) != 2)
1436  {
1437  isFeaturePoint[e.end()] = true;
1438  }
1439  }
1440 
1441 
1442  //
1443  // 3: Split feature edges into segments:
1444  // find point with not 2 feature edges -> start of feature segment
1445  //
1446 
1447  DynamicList<labelList> segments;
1448 
1449 
1450  boolList featVisited(featureToEdge_.size(), false);
1451 
1452  do
1453  {
1454  label startFeatI = -1;
1455 
1456  forAll(featVisited, featI)
1457  {
1458  if (!featVisited[featI])
1459  {
1460  startFeatI = featI;
1461 
1462  break;
1463  }
1464  }
1465 
1466  if (startFeatI == -1)
1467  {
1468  // No feature lines left.
1469  break;
1470  }
1471 
1472  segments.append
1473  (
1474  collectSegment
1475  (
1476  isFeaturePoint,
1477  featureToEdge_[startFeatI],
1478  featVisited
1479  )
1480  );
1481  }
1482  while (true);
1483 
1484 
1485  //
1486  // Store in *this
1487  //
1488  featureSegments_.setSize(segments.size());
1489 
1490  forAll(featureSegments_, segmentI)
1491  {
1492  featureSegments_[segmentI] = segments[segmentI];
1493  }
1494 }
1495 
1496 
1497 void Foam::boundaryMesh::setExtraEdges(const label edgeI)
1498 {
1499  labelList minDistance(mesh().nEdges(), -1);
1500 
1501  // All edge labels encountered
1502  DynamicList<label> visitedEdges;
1503 
1504  // Floodfill from edgeI starting from distance 0. Stop at distance.
1505  markEdges(8, edgeI, 0, minDistance, visitedEdges);
1506 
1507  // Set edge labels to display
1508  extraEdges_.transfer(visitedEdges);
1509 }
1510 
1511 
1512 Foam::label Foam::boundaryMesh::whichPatch(const label faceI) const
1513 {
1514  forAll(patches_, patchI)
1515  {
1516  const boundaryPatch& bp = patches_[patchI];
1517 
1518  if ((faceI >= bp.start()) && (faceI < (bp.start() + bp.size())))
1519  {
1520  return patchI;
1521  }
1522  }
1523 
1524  FatalErrorIn("boundaryMesh::whichPatch(const label)")
1525  << "Cannot find face " << faceI << " in list of boundaryPatches "
1526  << patches_
1527  << abort(FatalError);
1528 
1529  return -1;
1530 }
1531 
1532 
1533 Foam::label Foam::boundaryMesh::findPatchID(const word& patchName) const
1534 {
1535  forAll(patches_, patchI)
1536  {
1537  if (patches_[patchI].name() == patchName)
1538  {
1539  return patchI;
1540  }
1541  }
1542 
1543  return -1;
1544 }
1545 
1546 
1547 void Foam::boundaryMesh::addPatch(const word& patchName)
1548 {
1549  patches_.setSize(patches_.size() + 1);
1550 
1551  // Add empty patch at end of patch list.
1552 
1553  label patchI = patches_.size()-1;
1554 
1555  boundaryPatch* bpPtr = new boundaryPatch
1556  (
1557  patchName,
1558  patchI,
1559  0,
1560  mesh().size(),
1561  "empty"
1562  );
1563 
1564  patches_.set(patchI, bpPtr);
1565 
1566  if (debug)
1567  {
1568  Pout<< "addPatch : patches now:" << endl;
1569 
1570  forAll(patches_, patchI)
1571  {
1572  const boundaryPatch& bp = patches_[patchI];
1573 
1574  Pout<< " name : " << bp.name() << endl
1575  << " size : " << bp.size() << endl
1576  << " start : " << bp.start() << endl
1577  << " type : " << bp.physicalType() << endl
1578  << endl;
1579  }
1580  }
1581 }
1582 
1583 
1585 {
1586  label delPatchI = findPatchID(patchName);
1587 
1588  if (delPatchI == -1)
1589  {
1590  FatalErrorIn("boundaryMesh::deletePatch(const word&")
1591  << "Can't find patch named " << patchName
1592  << abort(FatalError);
1593  }
1594 
1595  if (patches_[delPatchI].size())
1596  {
1597  FatalErrorIn("boundaryMesh::deletePatch(const word&")
1598  << "Trying to delete non-empty patch " << patchName
1599  << endl << "Current size:" << patches_[delPatchI].size()
1600  << abort(FatalError);
1601  }
1602 
1603  PtrList<boundaryPatch> newPatches(patches_.size() - 1);
1604 
1605  for (label patchI = 0; patchI < delPatchI; patchI++)
1606  {
1607  newPatches.set(patchI, patches_[patchI].clone());
1608  }
1609 
1610  // Move patches down, starting from delPatchI.
1611 
1612  for (label patchI = delPatchI + 1; patchI < patches_.size(); patchI++)
1613  {
1614  newPatches.set(patchI - 1, patches_[patchI].clone());
1615  }
1616 
1617  patches_.clear();
1618 
1619  patches_ = newPatches;
1620 
1621  if (debug)
1622  {
1623  Pout<< "deletePatch : patches now:" << endl;
1624 
1625  forAll(patches_, patchI)
1626  {
1627  const boundaryPatch& bp = patches_[patchI];
1628 
1629  Pout<< " name : " << bp.name() << endl
1630  << " size : " << bp.size() << endl
1631  << " start : " << bp.start() << endl
1632  << " type : " << bp.physicalType() << endl
1633  << endl;
1634  }
1635  }
1636 }
1637 
1638 
1641  const word& patchName,
1642  const word& patchType
1643 )
1644 {
1645  label changeI = findPatchID(patchName);
1646 
1647  if (changeI == -1)
1648  {
1649  FatalErrorIn("boundaryMesh::changePatchType(const word&, const word&)")
1650  << "Can't find patch named " << patchName
1651  << abort(FatalError);
1652  }
1653 
1654 
1655  // Cause we can't reassign to individual PtrList elems ;-(
1656  // work on copy
1657 
1658 
1659  PtrList<boundaryPatch> newPatches(patches_.size());
1660 
1661  forAll(patches_, patchI)
1662  {
1663  if (patchI == changeI)
1664  {
1665  // Create copy but for type
1666  const boundaryPatch& oldBp = patches_[patchI];
1667 
1668  boundaryPatch* bpPtr = new boundaryPatch
1669  (
1670  oldBp.name(),
1671  oldBp.index(),
1672  oldBp.size(),
1673  oldBp.start(),
1674  patchType
1675  );
1676 
1677  newPatches.set(patchI, bpPtr);
1678  }
1679  else
1680  {
1681  // Create copy
1682  newPatches.set(patchI, patches_[patchI].clone());
1683  }
1684  }
1685 
1686  patches_ = newPatches;
1687 }
1688 
1689 
1692  const labelList& patchIDs,
1693  labelList& oldToNew
1694 )
1695 {
1696  if (patchIDs.size() != mesh().size())
1697  {
1698  FatalErrorIn("boundaryMesh::changeFaces(const labelList& patchIDs)")
1699  << "List of patchIDs not equal to number of faces." << endl
1700  << "PatchIDs size:" << patchIDs.size()
1701  << " nFaces::" << mesh().size()
1702  << abort(FatalError);
1703  }
1704 
1705  // Count number of faces for each patch
1706 
1707  labelList nFaces(patches_.size(), 0);
1708 
1709  forAll(patchIDs, faceI)
1710  {
1711  label patchID = patchIDs[faceI];
1712 
1713  if (patchID < 0 || patchID >= patches_.size())
1714  {
1715  FatalErrorIn("boundaryMesh::changeFaces(const labelList&)")
1716  << "PatchID " << patchID << " out of range"
1717  << abort(FatalError);
1718  }
1719  nFaces[patchID]++;
1720  }
1721 
1722 
1723  // Determine position in faces_ for each patch
1724 
1725  labelList startFace(patches_.size());
1726 
1727  startFace[0] = 0;
1728 
1729  for (label patchI = 1; patchI < patches_.size(); patchI++)
1730  {
1731  startFace[patchI] = startFace[patchI-1] + nFaces[patchI-1];
1732  }
1733 
1734  // Update patch info
1735  PtrList<boundaryPatch> newPatches(patches_.size());
1736 
1737  forAll(patches_, patchI)
1738  {
1739  const boundaryPatch& bp = patches_[patchI];
1740 
1741  newPatches.set
1742  (
1743  patchI,
1744  new boundaryPatch
1745  (
1746  bp.name(),
1747  patchI,
1748  nFaces[patchI],
1749  startFace[patchI],
1750  bp.physicalType()
1751  )
1752  );
1753  }
1754  patches_ = newPatches;
1755 
1756  if (debug)
1757  {
1758  Pout<< "changeFaces : patches now:" << endl;
1759 
1760  forAll(patches_, patchI)
1761  {
1762  const boundaryPatch& bp = patches_[patchI];
1763 
1764  Pout<< " name : " << bp.name() << endl
1765  << " size : " << bp.size() << endl
1766  << " start : " << bp.start() << endl
1767  << " type : " << bp.physicalType() << endl
1768  << endl;
1769  }
1770  }
1771 
1772 
1773  // Construct face mapping array
1774  oldToNew.setSize(patchIDs.size());
1775 
1776  forAll(patchIDs, faceI)
1777  {
1778  int patchID = patchIDs[faceI];
1779 
1780  oldToNew[faceI] = startFace[patchID]++;
1781  }
1782 
1783  // Copy faces into correct position and maintain label of original face
1784  faceList newFaces(mesh().size());
1785 
1786  labelList newMeshFace(mesh().size());
1787 
1788  forAll(oldToNew, faceI)
1789  {
1790  newFaces[oldToNew[faceI]] = mesh()[faceI];
1791  newMeshFace[oldToNew[faceI]] = meshFace_[faceI];
1792  }
1793 
1794  // Reconstruct 'mesh' from new faces and (copy of) existing points.
1795  bMesh* newMeshPtr_ = new bMesh(newFaces, mesh().points());
1796 
1797  // Reset meshFace_ to new ordering.
1798  meshFace_.transfer(newMeshFace);
1799 
1800 
1801  // Remove old PrimitivePatch on meshPtr_.
1802  clearOut();
1803 
1804  // And insert new 'mesh'.
1805  meshPtr_ = newMeshPtr_;
1806 }
1807 
1808 
1809 Foam::label Foam::boundaryMesh::getNTris(const label faceI) const
1810 {
1811  const face& f = mesh()[faceI];
1812 
1813  return f.nTriangles(mesh().points());
1814 }
1815 
1816 
1817 Foam::label Foam::boundaryMesh::getNTris
1819  const label startFaceI,
1820  const label nFaces,
1821  labelList& nTris
1822 ) const
1823 {
1824  label totalNTris = 0;
1825 
1826  nTris.setSize(nFaces);
1827 
1828  for (label i = 0; i < nFaces; i++)
1829  {
1830  label faceNTris = getNTris(startFaceI + i);
1831 
1832  nTris[i] = faceNTris;
1833 
1834  totalNTris += faceNTris;
1835  }
1836  return totalNTris;
1837 }
1838 
1839 
1840 // Simple triangulation of face subset. Stores vertices in tris[] as three
1841 // consecutive vertices per triangle.
1844  const label startFaceI,
1845  const label nFaces,
1846  const label totalNTris,
1847  labelList& triVerts
1848 ) const
1849 {
1850  // Triangulate faces.
1851  triVerts.setSize(3*totalNTris);
1852 
1853  label vertI = 0;
1854 
1855  for (label i = 0; i < nFaces; i++)
1856  {
1857  label faceI = startFaceI + i;
1858 
1859  const face& f = mesh()[faceI];
1860 
1861  // Have face triangulate itself (results in faceList)
1862  faceList triFaces(f.nTriangles(mesh().points()));
1863 
1864  label nTri = 0;
1865 
1866  f.triangles(mesh().points(), nTri, triFaces);
1867 
1868  // Copy into triVerts
1869 
1870  forAll(triFaces, triFaceI)
1871  {
1872  const face& triF = triFaces[triFaceI];
1873 
1874  triVerts[vertI++] = triF[0];
1875  triVerts[vertI++] = triF[1];
1876  triVerts[vertI++] = triF[2];
1877  }
1878  }
1879 }
1880 
1881 
1882 // Number of local points in subset.
1885  const label startFaceI,
1886  const label nFaces
1887 ) const
1888 {
1889  SubList<face> patchFaces(mesh(), nFaces, startFaceI);
1890 
1891  primitivePatch patch(patchFaces, mesh().points());
1892 
1893  return patch.nPoints();
1894 }
1895 
1896 
1897 // Triangulation of face subset in local coords.
1900  const label startFaceI,
1901  const label nFaces,
1902  const label totalNTris,
1903  labelList& triVerts,
1904  labelList& localToGlobal
1905 ) const
1906 {
1907  SubList<face> patchFaces(mesh(), nFaces, startFaceI);
1908 
1909  primitivePatch patch(patchFaces, mesh().points());
1910 
1911  // Map from local to mesh().points() addressing
1912  localToGlobal = patch.meshPoints();
1913 
1914  // Triangulate patch faces.
1915  triVerts.setSize(3*totalNTris);
1916 
1917  label vertI = 0;
1918 
1919  for (label i = 0; i < nFaces; i++)
1920  {
1921  // Local face
1922  const face& f = patch.localFaces()[i];
1923 
1924  // Have face triangulate itself (results in faceList)
1925  faceList triFaces(f.nTriangles(patch.localPoints()));
1926 
1927  label nTri = 0;
1928 
1929  f.triangles(patch.localPoints(), nTri, triFaces);
1930 
1931  // Copy into triVerts
1932 
1933  forAll(triFaces, triFaceI)
1934  {
1935  const face& triF = triFaces[triFaceI];
1936 
1937  triVerts[vertI++] = triF[0];
1938  triVerts[vertI++] = triF[1];
1939  triVerts[vertI++] = triF[2];
1940  }
1941  }
1942 }
1943 
1944 
1947  const labelList& protectedEdges,
1948  const label seedFaceI,
1949  boolList& visited
1950 ) const
1951 {
1952  boolList protectedEdge(mesh().nEdges(), false);
1953 
1954  forAll(protectedEdges, i)
1955  {
1956  protectedEdge[protectedEdges[i]] = true;
1957  }
1958 
1959 
1960  // Initialize zone for all faces to -1
1961  labelList currentZone(mesh().size(), -1);
1962 
1963  // Mark with 0 all faces reachable from seedFaceI
1964  markZone(protectedEdge, seedFaceI, 0, currentZone);
1965 
1966  // Set in visited all reached ones.
1967  visited.setSize(mesh().size());
1968 
1969  forAll(currentZone, faceI)
1970  {
1971  if (currentZone[faceI] == 0)
1972  {
1973  visited[faceI] = true;
1974  }
1975  else
1976  {
1977  visited[faceI] = false;
1978  }
1979  }
1980 }
1981 
1982 
1983 // ************************ vim: set sw=4 sts=4 et: ************************ //