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