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