repatchMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "repatchMesh.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(repatchMesh, 0);
43 
44 // Normal along which to divide faces into categories (used in getNearest)
45 const vector repatchMesh::splitNormal_(3, 2, 1);
46 
47 // Distance to face tolerance for getNearest
48 const scalar repatchMesh::distanceTol_ = 1e-2;
49 }
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 // Returns number of feature edges connected to pointi
55 Foam::label Foam::repatchMesh::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::repatchMesh::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::repatchMesh::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  {
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 // Gets labels of changed faces and propagates them to the edges. Returns
193 // labels of edges changed.
194 Foam::labelList Foam::repatchMesh::faceToEdge
195 (
196  const boolList& regionEdge,
197  const label region,
198  const labelList& changedFaces,
199  labelList& edgeRegion
200 ) const
201 {
202  labelList changedEdges(mesh().nEdges(), -1);
203  label changedI = 0;
204 
205  forAll(changedFaces, i)
206  {
207  label facei = changedFaces[i];
208 
209  const labelList& fEdges = mesh().faceEdges()[facei];
210 
211  forAll(fEdges, fEdgeI)
212  {
213  label edgeI = fEdges[fEdgeI];
214 
215  if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
216  {
217  edgeRegion[edgeI] = region;
218 
219  changedEdges[changedI++] = edgeI;
220  }
221  }
222  }
223 
224  changedEdges.setSize(changedI);
225 
226  return changedEdges;
227 }
228 
229 
230 // Reverse of faceToEdge: gets edges and returns faces
231 Foam::labelList Foam::repatchMesh::edgeToFace
232 (
233  const label region,
234  const labelList& changedEdges,
235  labelList& faceRegion
236 ) const
237 {
238  labelList changedFaces(mesh().size(), -1);
239  label changedI = 0;
240 
241  forAll(changedEdges, i)
242  {
243  label edgeI = changedEdges[i];
244 
245  const labelList& eFaces = mesh().edgeFaces()[edgeI];
246 
247  forAll(eFaces, eFacei)
248  {
249  label facei = eFaces[eFacei];
250 
251  if (faceRegion[facei] == -1)
252  {
253  faceRegion[facei] = region;
254 
255  changedFaces[changedI++] = facei;
256  }
257  }
258  }
259 
260  changedFaces.setSize(changedI);
261 
262  return changedFaces;
263 }
264 
265 
266 // Finds area, starting at facei, delimited by borderEdge
267 void Foam::repatchMesh::markZone
268 (
269  const boolList& borderEdge,
270  label facei,
271  label currentZone,
272  labelList& faceZone
273 ) const
274 {
275  faceZone[facei] = currentZone;
276 
277  // List of faces whose faceZone has been set.
278  labelList changedFaces(1, facei);
279  // List of edges whose faceZone has been set.
280  labelList changedEdges;
281 
282  // Zones on all edges.
283  labelList edgeZone(mesh().nEdges(), -1);
284 
285  while (true)
286  {
287  changedEdges = faceToEdge
288  (
289  borderEdge,
290  currentZone,
291  changedFaces,
292  edgeZone
293  );
294 
295  if (debug)
296  {
297  Pout<< "From changedFaces:" << changedFaces.size()
298  << " to changedEdges:" << changedEdges.size()
299  << endl;
300  }
301 
302  if (changedEdges.empty())
303  {
304  break;
305  }
306 
307  changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
308 
309  if (debug)
310  {
311  Pout<< "From changedEdges:" << changedEdges.size()
312  << " to changedFaces:" << changedFaces.size()
313  << endl;
314  }
315 
316  if (changedFaces.empty())
317  {
318  break;
319  }
320  }
321 }
322 
323 
324 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
325 
326 // Null constructor
328 :
329  meshPtr_(nullptr),
330  patches_(),
331  meshFace_(),
332  featurePoints_(),
333  featureEdges_(),
334  featureToEdge_(),
335  edgeToFeature_(),
336  featureSegments_()
337 {}
338 
339 
340 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
341 
343 {}
344 
345 
346 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
347 
349 {
350  patches_.clear();
351 
352  patches_.setSize(mesh.boundaryMesh().size());
353 
354  // Number of boundary faces
355  label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
356 
357  faceList bFaces(nBFaces);
358 
359  meshFace_.setSize(nBFaces);
360 
361  label bFacei = 0;
362 
363  // Collect all boundary faces.
364  forAll(mesh.boundaryMesh(), patchi)
365  {
366  const polyPatch& pp = mesh.boundaryMesh()[patchi];
367 
368  patches_.set
369  (
370  patchi,
371  new repatchPatch
372  (
373  pp.name(),
374  patchi,
375  pp.size(),
376  bFacei,
377  pp.type()
378  )
379  );
380 
381  // Collect all faces in global numbering.
382  forAll(pp, patchFacei)
383  {
384  meshFace_[bFacei] = pp.start() + patchFacei;
385 
386  bFaces[bFacei] = pp[patchFacei];
387 
388  bFacei++;
389  }
390  }
391 
392 
393  if (debug)
394  {
395  Pout<< "read : patches now:" << endl;
396 
397  forAll(patches_, patchi)
398  {
399  const repatchPatch& bp = patches_[patchi];
400 
401  Pout<< " name : " << bp.name() << endl
402  << " size : " << bp.size() << endl
403  << " start : " << bp.start() << endl
404  << " type : " << bp.physicalType() << endl
405  << endl;
406  }
407  }
408 
409  //
410  // Construct single patch for all of boundary
411  //
412 
413  // Temporary primitivePatch to calculate compact points & faces.
415  (
416  bFaces,
417  mesh.points()
418  );
419 
420  // Store in local(compact) addressing
421  meshPtr_ = new rMesh(globalPatch.localFaces(), globalPatch.localPoints());
422 
423 
424  if (debug & 2)
425  {
426  const rMesh& msh = *meshPtr_;
427 
428  Pout<< "** Start of Faces **" << endl;
429 
430  forAll(msh, facei)
431  {
432  const face& f = msh[facei];
433 
434  point ctr(Zero);
435 
436  forAll(f, fp)
437  {
438  ctr += msh.points()[f[fp]];
439  }
440  ctr /= f.size();
441 
442  Pout<< " " << facei
443  << " ctr:" << ctr
444  << " verts:" << f
445  << endl;
446  }
447 
448  Pout<< "** End of Faces **" << endl;
449 
450  Pout<< "** Start of Points **" << endl;
451 
452  forAll(msh.points(), pointi)
453  {
454  Pout<< " " << pointi
455  << " coord:" << msh.points()[pointi]
456  << endl;
457  }
458 
459  Pout<< "** End of Points **" << endl;
460  }
461 
462  // Clear edge storage
463  featurePoints_.setSize(0);
464  featureEdges_.setSize(0);
465 
466  featureToEdge_.setSize(0);
467  edgeToFeature_.setSize(meshPtr_->nEdges());
468  edgeToFeature_ = -1;
469 
470  featureSegments_.setSize(0);
471 }
472 
473 
475 {
476  triSurface surf(fName);
477 
478  if (surf.empty())
479  {
480  return;
481  }
482 
483  // Sort according to region
484  SortableList<label> regions(surf.size());
485 
486  forAll(surf, triI)
487  {
488  regions[triI] = surf[triI].region();
489  }
490  regions.sort();
491 
492  // Determine region mapping.
493  Map<label> regionToBoundaryPatch;
494 
495  label oldRegion = -1111;
496  label boundPatch = 0;
497 
498  forAll(regions, i)
499  {
500  if (regions[i] != oldRegion)
501  {
502  regionToBoundaryPatch.insert(regions[i], boundPatch);
503 
504  oldRegion = regions[i];
505  boundPatch++;
506  }
507  }
508 
509  const geometricSurfacePatchList& surfPatches = surf.patches();
510 
511  patches_.clear();
512 
513  if (surfPatches.size() == regionToBoundaryPatch.size())
514  {
515  // There are as many surface patches as region numbers in triangles
516  // so use the surface patches
517 
518  patches_.setSize(surfPatches.size());
519 
520  // Take over patches, setting size to 0 for now.
521  forAll(surfPatches, patchi)
522  {
523  const geometricSurfacePatch& surfPatch = surfPatches[patchi];
524 
525  patches_.set
526  (
527  patchi,
528  new repatchPatch
529  (
530  surfPatch.name(),
531  patchi,
532  0,
533  0,
534  surfPatch.geometricType()
535  )
536  );
537  }
538  }
539  else
540  {
541  // There are not enough surface patches. Make up my own.
542 
543  patches_.setSize(regionToBoundaryPatch.size());
544 
545  forAll(patches_, patchi)
546  {
547  patches_.set
548  (
549  patchi,
550  new repatchPatch
551  (
552  "patch" + name(patchi),
553  patchi,
554  0,
555  0,
556  "empty"
557  )
558  );
559  }
560  }
561 
562  //
563  // Copy according into bFaces according to regions
564  //
565 
566  const labelList& indices = regions.indices();
567 
568  faceList bFaces(surf.size());
569 
570  meshFace_.setSize(surf.size());
571 
572  label bFacei = 0;
573 
574  // Current region number
575  label surfRegion = regions[0];
576  label foamRegion = regionToBoundaryPatch[surfRegion];
577 
578  Pout<< "Surface region " << surfRegion << " becomes boundary patch "
579  << foamRegion << " with name " << patches_[foamRegion].name() << endl;
580 
581 
582  // Index in bFaces of start of current patch
583  label startFacei = 0;
584 
585  forAll(indices, indexI)
586  {
587  label triI = indices[indexI];
588 
589  const labelledTri& tri = surf.localFaces()[triI];
590 
591  if (tri.region() != surfRegion)
592  {
593  // Change of region. We now know the size of the previous one.
594  repatchPatch& bp = patches_[foamRegion];
595 
596  bp.size() = bFacei - startFacei;
597  bp.start() = startFacei;
598 
599  surfRegion = tri.region();
600  foamRegion = regionToBoundaryPatch[surfRegion];
601 
602  Pout<< "Surface region " << surfRegion << " becomes boundary patch "
603  << foamRegion << " with name " << patches_[foamRegion].name()
604  << endl;
605 
606  startFacei = bFacei;
607  }
608 
609  meshFace_[bFacei] = triI;
610 
611  bFaces[bFacei++] = face(tri);
612  }
613 
614  // Final region
615  repatchPatch& bp = patches_[foamRegion];
616 
617  bp.size() = bFacei - startFacei;
618  bp.start() = startFacei;
619 
620  //
621  // Construct single primitivePatch for all of boundary
622  //
623 
624  // Store compact.
625  meshPtr_ = new rMesh(bFaces, surf.localPoints());
626 
627  // Clear edge storage
628  featurePoints_.setSize(0);
629  featureEdges_.setSize(0);
630 
631  featureToEdge_.setSize(0);
632  edgeToFeature_.setSize(meshPtr_->nEdges());
633  edgeToFeature_ = -1;
634 
635  featureSegments_.setSize(0);
636 }
637 
638 
639 // Get index in this (repatchMesh) of face nearest to each boundary face in
640 // pMesh.
641 // Originally all triangles/faces of repatchMesh would be bunged into
642 // one big octree. Problem was that faces on top of each other, differing
643 // only in sign of normal, could not be found separately. It would always
644 // find only one. We could detect that it was probably finding the wrong one
645 // (based on normal) but could not 'tell' the octree to retrieve the other
646 // one (since they occupy exactly the same space)
647 // So now faces get put into different octrees depending on normal.
648 // !It still will not be possible to differentiate between two faces on top
649 // of each other having the same normal
651 (
652  const primitiveMesh& pMesh,
653  const vector& searchSpan
654 ) const
655 {
656 
657  // Divide faces into two bins acc. to normal
658  // - left of splitNormal
659  // - right ,,
660  DynamicList<label> leftFaces(mesh().size()/2);
661  DynamicList<label> rightFaces(mesh().size()/2);
662 
663  forAll(mesh(), bFacei)
664  {
665  scalar sign = mesh().faceNormals()[bFacei] & splitNormal_;
666 
667  if (sign > -1e-5)
668  {
669  rightFaces.append(bFacei);
670  }
671  if (sign < 1e-5)
672  {
673  leftFaces.append(bFacei);
674  }
675  }
676 
677  leftFaces.shrink();
678  rightFaces.shrink();
679 
680  if (debug)
681  {
682  Pout<< "getNearest :"
683  << " rightBin:" << rightFaces.size()
684  << " leftBin:" << leftFaces.size()
685  << endl;
686  }
687 
688  uindirectPrimitivePatch leftPatch
689  (
690  UIndirectList<face>(mesh(), leftFaces),
691  mesh().points()
692  );
693  uindirectPrimitivePatch rightPatch
694  (
695  UIndirectList<face>(mesh(), rightFaces),
696  mesh().points()
697  );
698 
699 
700  // Overall bb
701  treeBoundBox overallBb(mesh().localPoints());
702 
703  // Extend domain slightly (also makes it 3D if was 2D)
704  // Note asymmetry to avoid having faces align with octree cubes.
705  scalar tol = 1e-6 * overallBb.avgDim();
706 
707  point& bbMin = overallBb.min();
708  bbMin.x() -= tol;
709  bbMin.y() -= tol;
710  bbMin.z() -= tol;
711 
712  point& bbMax = overallBb.max();
713  bbMax.x() += 2*tol;
714  bbMax.y() += 2*tol;
715  bbMax.z() += 2*tol;
716 
717  const scalar planarTol =
719  perturbTol();
720 
721 
722  // Create the octrees
724  <
726  > leftTree
727  (
729  (
730  false, // cacheBb
731  leftPatch,
732  planarTol
733  ),
734  overallBb,
735  10, // maxLevel
736  10, // leafSize
737  3.0 // duplicity
738  );
740  <
742  > rightTree
743  (
745  (
746  false, // cacheBb
747  rightPatch,
748  planarTol
749  ),
750  overallBb,
751  10, // maxLevel
752  10, // leafSize
753  3.0 // duplicity
754  );
755 
756  if (debug)
757  {
758  Pout<< "getNearest : built trees" << endl;
759  }
760 
761 
762  const vectorField& ns = mesh().faceNormals();
763 
764 
765  //
766  // Search nearest triangle centre for every polyMesh boundary face
767  //
768 
769  labelList nearestBFacei(pMesh.nFaces() - pMesh.nInternalFaces());
770 
771  treeBoundBox tightest;
772 
773  const scalar searchDimSqr = magSqr(searchSpan);
774 
775  forAll(nearestBFacei, patchFacei)
776  {
777  label meshFacei = pMesh.nInternalFaces() + patchFacei;
778 
779  const point& ctr = pMesh.faceCentres()[meshFacei];
780 
781  if (debug && (patchFacei % 1000) == 0)
782  {
783  Pout<< "getNearest : patchFace:" << patchFacei
784  << " meshFacei:" << meshFacei << " ctr:" << ctr << endl;
785  }
786 
787 
788  // Get normal from area vector
789  vector n = pMesh.faceAreas()[meshFacei];
790  scalar area = mag(n);
791  n /= area;
792 
793  scalar typDim = -great;
794  const face& f = pMesh.faces()[meshFacei];
795 
796  forAll(f, fp)
797  {
798  typDim = max(typDim, mag(pMesh.points()[f[fp]] - ctr));
799  }
800 
801  // Search right tree
802  pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
803 
804  // Search left tree. Note: could start from rightDist bounding box
805  // instead of starting from top.
806  pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
807 
808  if (rightInfo.hit())
809  {
810  if (leftInfo.hit())
811  {
812  // Found in both trees. Compare normals.
813  label rightFacei = rightFaces[rightInfo.index()];
814  label leftFacei = leftFaces[leftInfo.index()];
815 
816  label rightDist = mag(rightInfo.hitPoint()-ctr);
817  label leftDist = mag(leftInfo.hitPoint()-ctr);
818 
819  scalar rightSign = n & ns[rightFacei];
820  scalar leftSign = n & ns[leftFacei];
821 
822  if
823  (
824  (rightSign > 0 && leftSign > 0)
825  || (rightSign < 0 && leftSign < 0)
826  )
827  {
828  // Both same sign. Choose nearest.
829  if (rightDist < leftDist)
830  {
831  nearestBFacei[patchFacei] = rightFacei;
832  }
833  else
834  {
835  nearestBFacei[patchFacei] = leftFacei;
836  }
837  }
838  else
839  {
840  // Differing sign.
841  // - if both near enough choose one with correct sign
842  // - otherwise choose nearest.
843 
844  // Get local dimension as max of distance between ctr and
845  // any face vertex.
846 
847  typDim *= distanceTol_;
848 
849  if (rightDist < typDim && leftDist < typDim)
850  {
851  // Different sign and nearby. Choosing matching normal
852  if (rightSign > 0)
853  {
854  nearestBFacei[patchFacei] = rightFacei;
855  }
856  else
857  {
858  nearestBFacei[patchFacei] = leftFacei;
859  }
860  }
861  else
862  {
863  // Different sign but faraway. Choosing nearest.
864  if (rightDist < leftDist)
865  {
866  nearestBFacei[patchFacei] = rightFacei;
867  }
868  else
869  {
870  nearestBFacei[patchFacei] = leftFacei;
871  }
872  }
873  }
874  }
875  else
876  {
877  // Found in right but not in left. Choose right regardless if
878  // correct sign. Note: do we want this?
879  label rightFacei = rightFaces[rightInfo.index()];
880  nearestBFacei[patchFacei] = rightFacei;
881  }
882  }
883  else
884  {
885  // No face found in right tree.
886 
887  if (leftInfo.hit())
888  {
889  // Found in left but not in right. Choose left regardless if
890  // correct sign. Note: do we want this?
891  nearestBFacei[patchFacei] = leftFaces[leftInfo.index()];
892  }
893  else
894  {
895  // No face found in left tree.
896  nearestBFacei[patchFacei] = -1;
897  }
898  }
899  }
900 
901  return nearestBFacei;
902 }
903 
904 
905 void Foam::repatchMesh::setFeatureEdges(const scalar minCos)
906 {
907  edgeToFeature_.setSize(mesh().nEdges());
908 
909  edgeToFeature_ = -1;
910 
911  // 1. Mark feature edges
912 
913  // Storage for edge labels that are features. Trim later.
914  featureToEdge_.setSize(mesh().nEdges());
915 
916  label featureI = 0;
917 
918  if (minCos >= 0.9999)
919  {
920  // Select everything
921  forAll(mesh().edges(), edgeI)
922  {
923  edgeToFeature_[edgeI] = featureI;
924  featureToEdge_[featureI++] = edgeI;
925  }
926  }
927  else
928  {
929  forAll(mesh().edges(), edgeI)
930  {
931  const labelList& eFaces = mesh().edgeFaces()[edgeI];
932 
933  if (eFaces.size() == 2)
934  {
935  label face0I = eFaces[0];
936 
937  label face1I = eFaces[1];
938 
941  // if (whichPatch(face0I) != whichPatch(face1I))
942  //{
943  // edgeToFeature_[edgeI] = featureI;
944  // featureToEdge_[featureI++] = edgeI;
945  //}
946  // else
947  {
948  const vector& n0 = mesh().faceNormals()[face0I];
949 
950  const vector& n1 = mesh().faceNormals()[face1I];
951 
952  float cosAng = n0 & n1;
953 
954  if (cosAng < minCos)
955  {
956  edgeToFeature_[edgeI] = featureI;
957  featureToEdge_[featureI++] = edgeI;
958  }
959  }
960  }
961  else
962  {
963  // Should not occur: 0 or more than two faces
964  edgeToFeature_[edgeI] = featureI;
965  featureToEdge_[featureI++] = edgeI;
966  }
967  }
968  }
969 
970  // Trim featureToEdge_ to actual number of edges.
971  featureToEdge_.setSize(featureI);
972 
973  //
974  // Compact edges i.e. relabel vertices.
975  //
976 
977  featureEdges_.setSize(featureI);
978  featurePoints_.setSize(mesh().nPoints());
979 
980  labelList featToMeshPoint(mesh().nPoints(), -1);
981 
982  label featPtI = 0;
983 
984  forAll(featureToEdge_, fEdgeI)
985  {
986  label edgeI = featureToEdge_[fEdgeI];
987 
988  const edge& e = mesh().edges()[edgeI];
989 
990  label start = featToMeshPoint[e.start()];
991 
992  if (start == -1)
993  {
994  featToMeshPoint[e.start()] = featPtI;
995 
996  featurePoints_[featPtI] = mesh().points()[e.start()];
997 
998  start = featPtI;
999 
1000  featPtI++;
1001  }
1002 
1003  label end = featToMeshPoint[e.end()];
1004 
1005  if (end == -1)
1006  {
1007  featToMeshPoint[e.end()] = featPtI;
1008 
1009  featurePoints_[featPtI] = mesh().points()[e.end()];
1010 
1011  end = featPtI;
1012 
1013  featPtI++;
1014  }
1015 
1016  // Store with renumbered vertices.
1017  featureEdges_[fEdgeI] = edge(start, end);
1018  }
1019 
1020  // Compact points
1021  featurePoints_.setSize(featPtI);
1022 
1023 
1024  //
1025  // 2. Mark endpoints of feature segments. These are points with
1026  // != 2 feature edges connected.
1027  // Note: can add geometric constraint here as well that if 2 feature
1028  // edges the angle between them should be less than xxx.
1029  //
1030 
1031  boolList isFeaturePoint(mesh().nPoints(), false);
1032 
1033  forAll(featureToEdge_, featI)
1034  {
1035  label edgeI = featureToEdge_[featI];
1036 
1037  const edge& e = mesh().edges()[edgeI];
1038 
1039  if (nFeatureEdges(e.start()) != 2)
1040  {
1041  isFeaturePoint[e.start()] = true;
1042  }
1043 
1044  if (nFeatureEdges(e.end()) != 2)
1045  {
1046  isFeaturePoint[e.end()] = true;
1047  }
1048  }
1049 
1050 
1051  //
1052  // 3: Split feature edges into segments:
1053  // find point with not 2 feature edges -> start of feature segment
1054  //
1055 
1056  DynamicList<labelList> segments;
1057 
1058 
1059  boolList featVisited(featureToEdge_.size(), false);
1060 
1061  do
1062  {
1063  label startFeatI = -1;
1064 
1065  forAll(featVisited, featI)
1066  {
1067  if (!featVisited[featI])
1068  {
1069  startFeatI = featI;
1070 
1071  break;
1072  }
1073  }
1074 
1075  if (startFeatI == -1)
1076  {
1077  // No feature lines left.
1078  break;
1079  }
1080 
1081  segments.append
1082  (
1083  collectSegment
1084  (
1085  isFeaturePoint,
1086  featureToEdge_[startFeatI],
1087  featVisited
1088  )
1089  );
1090  }
1091  while (true);
1092 
1093 
1094  //
1095  // Store in *this
1096  //
1097  featureSegments_.setSize(segments.size());
1098 
1099  forAll(featureSegments_, segmentI)
1100  {
1101  featureSegments_[segmentI] = segments[segmentI];
1102  }
1103 }
1104 
1105 
1107 {
1108  forAll(patches_, patchi)
1109  {
1110  const repatchPatch& bp = patches_[patchi];
1111 
1112  if ((facei >= bp.start()) && (facei < (bp.start() + bp.size())))
1113  {
1114  return patchi;
1115  }
1116  }
1117 
1119  << "Cannot find face " << facei << " in list of repatchPatches "
1120  << patches_
1121  << abort(FatalError);
1122 
1123  return -1;
1124 }
1125 
1126 
1128 {
1129  forAll(patches_, patchi)
1130  {
1131  if (patches_[patchi].name() == patchName)
1132  {
1133  return patchi;
1134  }
1135  }
1136 
1137  return -1;
1138 }
1139 
1140 
1141 void Foam::repatchMesh::addPatch(const word& patchName)
1142 {
1143  patches_.setSize(patches_.size() + 1);
1144 
1145  // Add empty patch at end of patch list.
1146 
1147  label patchi = patches_.size()-1;
1148 
1149  repatchPatch* bpPtr = new repatchPatch
1150  (
1151  patchName,
1152  patchi,
1153  0,
1154  mesh().size(),
1155  "empty"
1156  );
1157 
1158  patches_.set(patchi, bpPtr);
1159 
1160  if (debug)
1161  {
1162  Pout<< "addPatch : patches now:" << endl;
1163 
1164  forAll(patches_, patchi)
1165  {
1166  const repatchPatch& bp = patches_[patchi];
1167 
1168  Pout<< " name : " << bp.name() << endl
1169  << " size : " << bp.size() << endl
1170  << " start : " << bp.start() << endl
1171  << " type : " << bp.physicalType() << endl
1172  << endl;
1173  }
1174  }
1175 }
1176 
1177 
1178 void Foam::repatchMesh::deletePatch(const word& patchName)
1179 {
1180  const label delPatchi = findPatchID(patchName);
1181 
1182  if (delPatchi == -1)
1183  {
1185  << "Can't find patch named " << patchName
1186  << abort(FatalError);
1187  }
1188 
1189  if (patches_[delPatchi].size())
1190  {
1192  << "Trying to delete non-empty patch " << patchName
1193  << endl << "Current size:" << patches_[delPatchi].size()
1194  << abort(FatalError);
1195  }
1196 
1197  PtrList<repatchPatch> newPatches(patches_.size() - 1);
1198 
1199  for (label patchi = 0; patchi < delPatchi; patchi++)
1200  {
1201  newPatches.set(patchi, patches_[patchi].clone());
1202  }
1203 
1204  // Move patches down, starting from delPatchi.
1205 
1206  for (label patchi = delPatchi + 1; patchi < patches_.size(); patchi++)
1207  {
1208  newPatches.set(patchi - 1, patches_[patchi].clone());
1209  }
1210 
1211  patches_.clear();
1212 
1213  patches_ = newPatches;
1214 
1215  if (debug)
1216  {
1217  Pout<< "deletePatch : patches now:" << endl;
1218 
1219  forAll(patches_, patchi)
1220  {
1221  const repatchPatch& bp = patches_[patchi];
1222 
1223  Pout<< " name : " << bp.name() << endl
1224  << " size : " << bp.size() << endl
1225  << " start : " << bp.start() << endl
1226  << " type : " << bp.physicalType() << endl
1227  << endl;
1228  }
1229  }
1230 }
1231 
1232 
1235  const word& patchName,
1236  const word& patchType
1237 )
1238 {
1239  const label changeI = findPatchID(patchName);
1240 
1241  if (changeI == -1)
1242  {
1244  << "Can't find patch named " << patchName
1245  << abort(FatalError);
1246  }
1247 
1248 
1249  // Cause we can't reassign to individual PtrList elems ;-(
1250  // work on copy
1251 
1252 
1253  PtrList<repatchPatch> newPatches(patches_.size());
1254 
1255  forAll(patches_, patchi)
1256  {
1257  if (patchi == changeI)
1258  {
1259  // Create copy but for type
1260  const repatchPatch& oldBp = patches_[patchi];
1261 
1262  repatchPatch* bpPtr = new repatchPatch
1263  (
1264  oldBp.name(),
1265  oldBp.index(),
1266  oldBp.size(),
1267  oldBp.start(),
1268  patchType
1269  );
1270 
1271  newPatches.set(patchi, bpPtr);
1272  }
1273  else
1274  {
1275  // Create copy
1276  newPatches.set(patchi, patches_[patchi].clone());
1277  }
1278  }
1279 
1280  patches_ = newPatches;
1281 }
1282 
1283 
1286  const labelList& protectedEdges,
1287  const label seedFacei,
1288  boolList& visited
1289 ) const
1290 {
1291  boolList protectedEdge(mesh().nEdges(), false);
1292 
1293  forAll(protectedEdges, i)
1294  {
1295  protectedEdge[protectedEdges[i]] = true;
1296  }
1297 
1298 
1299  // Initialise zone for all faces to -1
1300  labelList currentZone(mesh().size(), -1);
1301 
1302  // Mark with 0 all faces reachable from seedFacei
1303  markZone(protectedEdge, seedFacei, 0, currentZone);
1304 
1305  // Set in visited all reached ones.
1306  visited.setSize(mesh().size());
1307 
1308  forAll(currentZone, facei)
1309  {
1310  if (currentZone[facei] == 0)
1311  {
1312  visited[facei] = true;
1313  }
1314  else
1315  {
1316  visited[facei] = false;
1317  }
1318  }
1319 }
1320 
1321 
1322 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
dimensionedScalar sign(const dimensionedScalar &ds)
const rMesh & mesh() const
Access the boundary mesh.
Definition: repatchMesh.H:180
void deletePatch(const word &patchName)
Delete from patch list.
Definition: repatchMesh.C:1178
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
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
const word & name() const
Return name.
A class for handling file names.
Definition: fileName.H:79
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
const word & name() const
Return name.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const labelListList & faceEdges() const
void setFeatureEdges(const scalar minCos)
Set featureEdges, edgeToFeature, featureSegments according.
Definition: repatchMesh.C:905
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
label nInternalFaces() const
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:80
const labelListList & pointEdges() const
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
void read(const polyMesh &)
Read from repatchMesh of polyMesh.
Definition: repatchMesh.C:348
label nFaces() const
void changePatchType(const word &patchName, const word &type)
Change patch.
Definition: repatchMesh.C:1234
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const word & geometricType() const
Return the type of the patch.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label size() const
Definition: repatchPatch.H:101
virtual const pointField & points() const =0
Return mesh points.
const Cmpt & z() const
Definition: VectorI.H:87
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const Field< PointType > & localPoints() const
Return pointField of points in patch.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const Cmpt & y() const
Definition: VectorI.H:81
scalar avgDim() const
Average length/height/width dimension.
Definition: boundBoxI.H:114
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
A list of faces which address into the list of points.
const Point & hitPoint() const
Return hit point.
label region() const
Return region label.
Definition: labelledTriI.H:68
T clone(const T &t)
Definition: List.H:54
dynamicFvMesh & mesh
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A class for handling words, derived from string.
Definition: word.H:59
label nPoints
void readTriSurface(const fileName &)
Read from triSurface.
Definition: repatchMesh.C:474
Encapsulation of data needed to search on PrimitivePatches.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
Triangle with additional region number.
Definition: labelledTri.H:57
const Field< PointType > & faceNormals() const
Return face normals for patch.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
dimensioned< scalar > magSqr(const dimensioned< Type > &)
label findPatchID(const word &patchName) const
Get index of patch by name.
Definition: repatchMesh.C:1127
label nEdges() const
Return number of edges in patch.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
repatchMesh()
Construct null.
Definition: repatchMesh.C:327
void markFaces(const labelList &protectedEdges, const label facei, boolList &visited) const
Definition: repatchMesh.C:1285
const word & physicalType() const
Return the optional physical type of the patch.
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label whichPatch(const label facei) const
Get index of patch face is in.
Definition: repatchMesh.C:1106
PrimitivePatch< faceList, const pointField > rMesh
Type of the mesh.
Definition: repatchMesh.H:64
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
~repatchMesh()
Destructor.
Definition: repatchMesh.C:342
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
label patchi
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
void addPatch(const word &patchName)
Add to back of patch list.
Definition: repatchMesh.C:1141
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const vectorField & faceAreas() const
A List with indirect addressing.
Definition: fvMatrix.H:106
label index() const
Return the index of this patch in the boundaryMesh.
virtual const faceList & faces() const =0
Return faces.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
The geometricSurfacePatch is like patchIdentifier but for surfaces. Holds type, name and index...
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
labelList getNearest(const primitiveMesh &pMesh, const vector &searchSpan) const
Get rMesh index of nearest face for every boundary face in.
Definition: repatchMesh.C:651
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label start() const
Definition: repatchPatch.H:111
Triangulated surface description with patch information.
Definition: triSurface.H:66
label index() const
Return index.
const labelListList & edgeFaces() const
label start() const
Return start vertex label.
Definition: edgeI.H:81
Namespace for OpenFOAM.
Like polyPatch but without reference to mesh. patchIdentifier::index is not used. Used in repatchMesh...
Definition: repatchPatch.H:57