createShellMesh.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019-2022 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "createShellMesh.H"
30 #include "polyTopoChange.H"
31 #include "meshTools.H"
32 #include "mapPolyMesh.H"
33 #include "polyAddPoint.H"
34 #include "polyAddFace.H"
35 #include "polyModifyFace.H"
36 #include "polyAddCell.H"
37 #include "labelPair.H"
38 #include "indirectPrimitivePatch.H"
39 #include "mapDistribute.H"
40 #include "globalMeshData.H"
41 #include "PatchTools.H"
42 #include "globalIndex.H"
43 
44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
45 
46 namespace Foam
47 {
49 
50 template<>
51 class minEqOp<labelPair>
52 {
53 public:
54  void operator()(labelPair& x, const labelPair& y) const
55  {
56  x[0] = min(x[0], y[0]);
57  x[1] = min(x[1], y[1]);
58  }
59 };
60 }
61 
62 
63 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
64 
65 // Synchronise edges
66 void Foam::createShellMesh::syncEdges
67 (
68  const globalMeshData& globalData,
69 
70  const labelList& patchEdges,
71  const labelList& coupledEdges,
72  const bitSet& sameEdgeOrientation,
73  const bool syncNonCollocated,
74 
75  bitSet& isChangedEdge,
76  DynamicList<label>& changedEdges,
77  labelPairList& allEdgeData
78 )
79 {
80  const mapDistribute& map = globalData.globalEdgeSlavesMap();
81  const bitSet& cppOrientation = globalData.globalEdgeOrientation();
82 
83  // Convert patch-edge data into cpp-edge data
84  labelPairList cppEdgeData
85  (
86  map.constructSize(),
88  );
89 
90  forAll(patchEdges, i)
91  {
92  label patchEdgeI = patchEdges[i];
93  label coupledEdgeI = coupledEdges[i];
94 
95  if (isChangedEdge[patchEdgeI])
96  {
97  const labelPair& data = allEdgeData[patchEdgeI];
98 
99  // Patch-edge data needs to be converted into coupled-edge data
100  // (optionally flipped) and consistent in orientation with
101  // other coupled edge (optionally flipped)
102  if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
103  {
104  cppEdgeData[coupledEdgeI] = data;
105  }
106  else
107  {
108  cppEdgeData[coupledEdgeI] = labelPair(data[1], data[0]);
109  }
110  }
111  }
112 
113  // Synchronise
114  globalData.syncData
115  (
116  cppEdgeData,
117  globalData.globalEdgeSlaves(),
118  (
119  syncNonCollocated
120  ? globalData.globalEdgeTransformedSlaves() // transformed elems
121  : labelListList(globalData.globalEdgeSlaves().size()) //no transformed
122  ),
123  map,
124  minEqOp<labelPair>()
125  );
126 
127  // Back from cpp-edge to patch-edge data
128  forAll(patchEdges, i)
129  {
130  label patchEdgeI = patchEdges[i];
131  label coupledEdgeI = coupledEdges[i];
132 
133  if (cppEdgeData[coupledEdgeI] != labelPair(labelMax, labelMax))
134  {
135  const labelPair& data = cppEdgeData[coupledEdgeI];
136 
137  if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
138  {
139  allEdgeData[patchEdgeI] = data;
140  }
141  else
142  {
143  allEdgeData[patchEdgeI] = labelPair(data[1], data[0]);
144  }
145 
146  if (!isChangedEdge[patchEdgeI])
147  {
148  changedEdges.append(patchEdgeI);
149  isChangedEdge.set(patchEdgeI);
150  }
151  }
152  }
153 }
154 
155 
157 (
158  const globalMeshData& globalData,
159  const primitiveFacePatch& patch,
160  const bitSet& nonManifoldEdge,
161  const bool syncNonCollocated,
162 
163  faceList& pointGlobalRegions,
164  faceList& pointLocalRegions,
165  labelList& localToGlobalRegion
166 )
167 {
168  const indirectPrimitivePatch& cpp = globalData.coupledPatch();
169 
170  // Calculate correspondence between patch and globalData.coupledPatch.
171  labelList patchEdges;
172  labelList coupledEdges;
173  bitSet sameEdgeOrientation;
175  (
176  cpp,
177  patch,
178 
179  coupledEdges,
180  patchEdges,
181  sameEdgeOrientation
182  );
183 
184 
185  // Initial unique regions
186  // ~~~~~~~~~~~~~~~~~~~~~~
187  // These get merged later on across connected edges.
188 
189  // 1. Count
190  label nMaxRegions = 0;
191  forAll(patch.localFaces(), facei)
192  {
193  const face& f = patch.localFaces()[facei];
194  nMaxRegions += f.size();
195  }
196 
197  const globalIndex globalRegions(nMaxRegions);
198 
199  // 2. Assign unique regions
200  label nRegions = 0;
201 
202  pointGlobalRegions.setSize(patch.size());
203  forAll(pointGlobalRegions, facei)
204  {
205  const face& f = patch.localFaces()[facei];
206  labelList& pRegions = pointGlobalRegions[facei];
207  pRegions.setSize(f.size());
208  forAll(pRegions, fp)
209  {
210  pRegions[fp] = globalRegions.toGlobal(nRegions++);
211  }
212  }
213 
214 
215  DynamicList<label> changedEdges(patch.nEdges());
216  labelPairList allEdgeData(patch.nEdges(), labelPair(labelMax, labelMax));
217  bitSet isChangedEdge(patch.nEdges());
218 
219 
220  // Fill initial seed
221  // ~~~~~~~~~~~~~~~~~
222 
223  forAll(patch.edgeFaces(), edgeI)
224  {
225  if (!nonManifoldEdge[edgeI])
226  {
227  // Take over value from one face only.
228  const edge& e = patch.edges()[edgeI];
229  label facei = patch.edgeFaces()[edgeI][0];
230  const face& f = patch.localFaces()[facei];
231 
232  label fp0 = f.find(e[0]);
233  label fp1 = f.find(e[1]);
234  allEdgeData[edgeI] = labelPair
235  (
236  pointGlobalRegions[facei][fp0],
237  pointGlobalRegions[facei][fp1]
238  );
239  if (!isChangedEdge[edgeI])
240  {
241  changedEdges.append(edgeI);
242  isChangedEdge.set(edgeI);
243  }
244  }
245  }
246 
247 
248  syncEdges
249  (
250  globalData,
251 
252  patchEdges,
253  coupledEdges,
254  sameEdgeOrientation,
255  syncNonCollocated,
256 
257  isChangedEdge,
258  changedEdges,
259  allEdgeData
260  );
261 
262 
263  // Edge-Face-Edge walk across patch
264  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
265  // Across edge minimum regions win
266 
267  while (true)
268  {
269  // From edge to face
270  // ~~~~~~~~~~~~~~~~~
271 
272  DynamicList<label> changedFaces(patch.size());
273  bitSet isChangedFace(patch.size());
274 
275  forAll(changedEdges, changedI)
276  {
277  label edgeI = changedEdges[changedI];
278  const labelPair& edgeData = allEdgeData[edgeI];
279 
280  const edge& e = patch.edges()[edgeI];
281  const labelList& eFaces = patch.edgeFaces()[edgeI];
282 
283  forAll(eFaces, i)
284  {
285  label facei = eFaces[i];
286  const face& f = patch.localFaces()[facei];
287 
288  // Combine edgeData with face data
289  label fp0 = f.find(e[0]);
290  if (pointGlobalRegions[facei][fp0] > edgeData[0])
291  {
292  pointGlobalRegions[facei][fp0] = edgeData[0];
293  if (!isChangedFace[facei])
294  {
295  isChangedFace.set(facei);
296  changedFaces.append(facei);
297  }
298  }
299 
300  label fp1 = f.find(e[1]);
301  if (pointGlobalRegions[facei][fp1] > edgeData[1])
302  {
303  pointGlobalRegions[facei][fp1] = edgeData[1];
304  if (!isChangedFace[facei])
305  {
306  isChangedFace.set(facei);
307  changedFaces.append(facei);
308  }
309  }
310  }
311  }
312 
313 
314  if (returnReduceAnd(changedFaces.empty()))
315  {
316  break;
317  }
318 
319 
320  // From face to edge
321  // ~~~~~~~~~~~~~~~~~
322 
323  isChangedEdge = false;
324  changedEdges.clear();
325 
326  forAll(changedFaces, i)
327  {
328  label facei = changedFaces[i];
329  const face& f = patch.localFaces()[facei];
330  const labelList& fEdges = patch.faceEdges()[facei];
331 
332  forAll(fEdges, fp)
333  {
334  label edgeI = fEdges[fp];
335 
336  if (!nonManifoldEdge[edgeI])
337  {
338  const edge& e = patch.edges()[edgeI];
339  label fp0 = f.find(e[0]);
340  label region0 = pointGlobalRegions[facei][fp0];
341  label fp1 = f.find(e[1]);
342  label region1 = pointGlobalRegions[facei][fp1];
343 
344  if
345  (
346  (allEdgeData[edgeI][0] > region0)
347  || (allEdgeData[edgeI][1] > region1)
348  )
349  {
350  allEdgeData[edgeI] = labelPair(region0, region1);
351  if (!isChangedEdge[edgeI])
352  {
353  changedEdges.append(edgeI);
354  isChangedEdge.set(edgeI);
355  }
356  }
357  }
358  }
359  }
360 
361  syncEdges
362  (
363  globalData,
364 
365  patchEdges,
366  coupledEdges,
367  sameEdgeOrientation,
368  syncNonCollocated,
369 
370  isChangedEdge,
371  changedEdges,
372  allEdgeData
373  );
374 
375 
376  if (returnReduceAnd(changedEdges.empty()))
377  {
378  break;
379  }
380  }
381 
382 
383 
384  // Assign local regions
385  // ~~~~~~~~~~~~~~~~~~~~
386 
387  // Calculate addressing from global region back to local region
388  pointLocalRegions.setSize(patch.size());
389  Map<label> globalToLocalRegion(globalRegions.localSize()/4);
390  DynamicList<label> dynLocalToGlobalRegion(globalToLocalRegion.size());
391  forAll(patch.localFaces(), facei)
392  {
393  const face& f = patch.localFaces()[facei];
394  face& pRegions = pointLocalRegions[facei];
395  pRegions.setSize(f.size());
396  forAll(f, fp)
397  {
398  label globalRegionI = pointGlobalRegions[facei][fp];
399 
400  const auto fnd = globalToLocalRegion.cfind(globalRegionI);
401 
402  if (fnd.good())
403  {
404  // Already encountered this global region. Assign same local one
405  pRegions[fp] = fnd();
406  }
407  else
408  {
409  // Region not yet seen. Create new one
410  label localRegionI = globalToLocalRegion.size();
411  pRegions[fp] = localRegionI;
412  globalToLocalRegion.insert(globalRegionI, localRegionI);
413  dynLocalToGlobalRegion.append(globalRegionI);
414  }
415  }
416  }
417  localToGlobalRegion.transfer(dynLocalToGlobalRegion);
418 }
419 
420 
421 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
422 
423 Foam::createShellMesh::createShellMesh
424 (
425  const primitiveFacePatch& patch,
426  const faceList& pointRegions,
427  const labelList& regionPoints
428 )
429 :
430  patch_(patch),
431  pointRegions_(pointRegions),
432  regionPoints_(regionPoints)
433 {
434  if (pointRegions_.size() != patch_.size())
435  {
437  << "nFaces:" << patch_.size()
438  << " pointRegions:" << pointRegions.size()
439  << exit(FatalError);
440  }
441 }
442 
443 
444 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
445 
447 (
448  const pointField& firstLayerDisp,
449  const scalar expansionRatio,
450  const label nLayers,
451  const labelList& topPatchID,
452  const labelList& bottomPatchID,
453  const labelListList& extrudeEdgePatches,
454  polyTopoChange& meshMod
455 )
456 {
457  if (firstLayerDisp.size() != regionPoints_.size())
458  {
460  << "nRegions:" << regionPoints_.size()
461  << " firstLayerDisp:" << firstLayerDisp.size()
462  << exit(FatalError);
463  }
464 
465  if
466  (
467  topPatchID.size() != patch_.size()
468  && bottomPatchID.size() != patch_.size()
469  )
470  {
472  << "nFaces:" << patch_.size()
473  << " topPatchID:" << topPatchID.size()
474  << " bottomPatchID:" << bottomPatchID.size()
475  << exit(FatalError);
476  }
477 
478  if (extrudeEdgePatches.size() != patch_.nEdges())
479  {
481  << "nEdges:" << patch_.nEdges()
482  << " extrudeEdgePatches:" << extrudeEdgePatches.size()
483  << exit(FatalError);
484  }
485 
486 
487 
488  // From cell to patch (trivial)
489  DynamicList<label> cellToFaceMap(nLayers*patch_.size());
490  // From face to patch+turning index
491  DynamicList<label> faceToFaceMap
492  (
493  (nLayers+1)*(patch_.size()+patch_.nEdges())
494  );
495  // From face to patch edge index
496  DynamicList<label> faceToEdgeMap(nLayers*(patch_.nEdges()+patch_.nEdges()));
497  // From point to patch point index
498  DynamicList<label> pointToPointMap((nLayers+1)*patch_.nPoints());
499 
500 
501  // Introduce new cell for every face
502  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
503 
504  labelList addedCells(nLayers*patch_.size());
505  forAll(patch_, facei)
506  {
507  for (label layerI = 0; layerI < nLayers; layerI++)
508  {
509  addedCells[nLayers*facei+layerI] = meshMod.addCell
510  (
511  -1, // masterPointID
512  -1, // masterEdgeID
513  -1, // masterFaceID
514  cellToFaceMap.size(), // masterCellID
515  -1 // zoneID
516  );
517  cellToFaceMap.append(facei);
518  }
519  }
520 
521 
522  // Introduce original points
523  // ~~~~~~~~~~~~~~~~~~~~~~~~~
524 
525  // Original point numbers in local point ordering so no need to store.
526  forAll(patch_.localPoints(), pointi)
527  {
528  //label addedPointi =
529  meshMod.addPoint
530  (
531  patch_.localPoints()[pointi], // point
532  pointToPointMap.size(), // masterPointID
533  -1, // zoneID
534  true // inCell
535  );
536  pointToPointMap.append(pointi);
537 
538  //Pout<< "Added bottom point " << addedPointi
539  // << " at " << patch_.localPoints()[pointi]
540  // << " from point " << pointi
541  // << endl;
542  }
543 
544 
545  // Introduce new points (one for every region)
546  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547 
548  labelList addedPoints(nLayers*regionPoints_.size());
549  forAll(regionPoints_, regionI)
550  {
551  label pointi = regionPoints_[regionI];
552 
553  point pt = patch_.localPoints()[pointi];
554  point disp = firstLayerDisp[regionI];
555  for (label layerI = 0; layerI < nLayers; layerI++)
556  {
557  pt += disp;
558 
559  addedPoints[nLayers*regionI+layerI] = meshMod.addPoint
560  (
561  pt, // point
562  pointToPointMap.size(), // masterPointID - used only addressing
563  -1, // zoneID
564  true // inCell
565  );
566  pointToPointMap.append(pointi);
567 
568  disp *= expansionRatio;
569  }
570  }
571 
572 
573  // Add face on bottom side
574  forAll(patch_.localFaces(), facei)
575  {
576  meshMod.addFace
577  (
578  patch_.localFaces()[facei].reverseFace(),// vertices
579  addedCells[nLayers*facei], // own
580  -1, // nei
581  -1, // masterPointID
582  -1, // masterEdgeID
583  faceToFaceMap.size(), // masterFaceID : current facei
584  true, // flipFaceFlux
585  bottomPatchID[facei], // patchID
586  -1, // zoneID
587  false // zoneFlip
588  );
589  faceToFaceMap.append(-facei-1); // points to flipped original face
590  faceToEdgeMap.append(-1);
591 
592  //const face newF(patch_.localFaces()[facei].reverseFace());
593  //Pout<< "Added bottom face "
594  // << newF
595  // << " coords:" << UIndirectList<point>(meshMod.points(), newF)
596  // << " own " << addedCells[facei]
597  // << " patch:" << bottomPatchID[facei]
598  // << " at " << patch_.faceCentres()[facei]
599  // << endl;
600  }
601 
602  // Add inbetween faces and face on top
603  forAll(patch_.localFaces(), facei)
604  {
605  // Get face in original ordering
606  const face& f = patch_.localFaces()[facei];
607 
608  face newF(f.size());
609 
610  for (label layerI = 0; layerI < nLayers; layerI++)
611  {
612  // Pick up point based on region and layer
613  forAll(f, fp)
614  {
615  label region = pointRegions_[facei][fp];
616  newF[fp] = addedPoints[region*nLayers+layerI];
617  }
618 
619  label own = addedCells[facei*nLayers+layerI];
620  label nei;
621  label patchi;
622  if (layerI == nLayers-1)
623  {
624  nei = -1;
625  patchi = topPatchID[facei];
626  }
627  else
628  {
629  nei = addedCells[facei*nLayers+layerI+1];
630  patchi = -1;
631  }
632 
633  meshMod.addFace
634  (
635  newF, // vertices
636  own, // own
637  nei, // nei
638  -1, // masterPointID
639  -1, // masterEdgeID
640  faceToFaceMap.size(), // masterFaceID : current facei
641  false, // flipFaceFlux
642  patchi, // patchID
643  -1, // zoneID
644  false // zoneFlip
645  );
646  faceToFaceMap.append(facei+1); // unflipped
647  faceToEdgeMap.append(-1);
648 
649  //Pout<< "Added inbetween face " << newF
650  // << " coords:" << UIndirectList<point>(meshMod.points(), newF)
651  // << " at layer " << layerI
652  // << " own " << own
653  // << " nei " << nei
654  // << " at " << patch_.faceCentres()[facei]
655  // << endl;
656  }
657  }
658 
659 
660  // Add side faces
661  // ~~~~~~~~~~~~~~
662  // Note that we loop over edges multiple times so for edges with
663  // two cyclic faces they get added in two passes (for correct ordering)
664 
665  // Pass1. Internal edges and first face of other edges
666  forAll(extrudeEdgePatches, edgeI)
667  {
668  const labelList& eFaces = patch_.edgeFaces()[edgeI];
669  const labelList& ePatches = extrudeEdgePatches[edgeI];
670 
671  if (ePatches.size() == 0)
672  {
673  // Internal face
674  if (eFaces.size() != 2)
675  {
677  << "edge:" << edgeI
678  << " not internal but does not have side-patches defined."
679  << exit(FatalError);
680  }
681  }
682  else
683  {
684  if (eFaces.size() != ePatches.size())
685  {
687  << "external/feature edge:" << edgeI
688  << " has " << eFaces.size() << " connected extruded faces "
689  << " but only " << ePatches.size()
690  << " boundary faces defined." << exit(FatalError);
691  }
692  }
693 
694 
695 
696  // Make face pointing in to eFaces[0] so out of new master face
697  const face& f = patch_.localFaces()[eFaces[0]];
698  const edge& e = patch_.edges()[edgeI];
699 
700  label fp0 = f.find(e[0]);
701  label fp1 = f.fcIndex(fp0);
702 
703  if (f[fp1] != e[1])
704  {
705  fp1 = fp0;
706  fp0 = f.rcIndex(fp1);
707  }
708 
709  face newF(4);
710 
711  for (label layerI = 0; layerI < nLayers; layerI++)
712  {
713  label region0 = pointRegions_[eFaces[0]][fp0];
714  label region1 = pointRegions_[eFaces[0]][fp1];
715 
716  // Pick up points with correct normal
717  if (layerI == 0)
718  {
719  newF[0] = f[fp0];
720  newF[1] = f[fp1];
721  newF[2] = addedPoints[nLayers*region1+layerI];
722  newF[3] = addedPoints[nLayers*region0+layerI];
723  }
724  else
725  {
726  newF[0] = addedPoints[nLayers*region0+layerI-1];
727  newF[1] = addedPoints[nLayers*region1+layerI-1];
728  newF[2] = addedPoints[nLayers*region1+layerI];
729  newF[3] = addedPoints[nLayers*region0+layerI];
730  }
731 
732  // Optionally rotate so e[0] is always 0th vertex. Note that
733  // this normally is automatically done by coupled face ordering
734  // but with NOORDERING we have to do it ourselves.
735  if (f[fp0] != e[0])
736  {
737  // rotate one back to get newF[1] (originating from e[0])
738  // into newF[0]
739  label v0 = newF[0];
740  for (label i = 0; i < newF.size()-1; i++)
741  {
742  newF[i] = newF[newF.fcIndex(i)];
743  }
744  newF.last() = v0;
745  }
746 
747 
748  label minCelli = addedCells[nLayers*eFaces[0]+layerI];
749  label maxCelli;
750  label patchi;
751  if (ePatches.size() == 0)
752  {
753  maxCelli = addedCells[nLayers*eFaces[1]+layerI];
754  if (minCelli > maxCelli)
755  {
756  // Swap
757  std::swap(minCelli, maxCelli);
758  newF.flip();
759  }
760  patchi = -1;
761  }
762  else
763  {
764  maxCelli = -1;
765  patchi = ePatches[0];
766  }
767 
768  //{
769  // Pout<< "Adding from face:" << patch_.faceCentres()[eFaces[0]]
770  // << " from edge:"
771  // << patch_.localPoints()[f[fp0]]
772  // << patch_.localPoints()[f[fp1]]
773  // << " at layer:" << layerI
774  // << " with new points:" << newF
775  // << " locations:"
776  // << UIndirectList<point>(meshMod.points(), newF)
777  // << " own:" << minCelli
778  // << " nei:" << maxCelli
779  // << endl;
780  //}
781 
782 
783  // newF already outwards pointing.
784  meshMod.addFace
785  (
786  newF, // vertices
787  minCelli, // own
788  maxCelli, // nei
789  -1, // masterPointID
790  -1, // masterEdgeID
791  faceToFaceMap.size(), // masterFaceID
792  false, // flipFaceFlux
793  patchi, // patchID
794  -1, // zoneID
795  false // zoneFlip
796  );
797  faceToFaceMap.append(0);
798  faceToEdgeMap.append(edgeI);
799  }
800  }
801 
802  // Pass2. Other faces of boundary edges
803  forAll(extrudeEdgePatches, edgeI)
804  {
805  const labelList& eFaces = patch_.edgeFaces()[edgeI];
806  const labelList& ePatches = extrudeEdgePatches[edgeI];
807 
808  if (ePatches.size() >= 2)
809  {
810  for (label i = 1; i < ePatches.size(); i++)
811  {
812  // Extrude eFaces[i]
813  label minFacei = eFaces[i];
814 
815  // Make face pointing in to eFaces[0] so out of new master face
816  const face& f = patch_.localFaces()[minFacei];
817 
818  const edge& e = patch_.edges()[edgeI];
819  label fp0 = f.find(e[0]);
820  label fp1 = f.fcIndex(fp0);
821 
822  if (f[fp1] != e[1])
823  {
824  fp1 = fp0;
825  fp0 = f.rcIndex(fp1);
826  }
827 
828  face newF(4);
829  for (label layerI = 0; layerI < nLayers; layerI++)
830  {
831  label region0 = pointRegions_[minFacei][fp0];
832  label region1 = pointRegions_[minFacei][fp1];
833 
834  if (layerI == 0)
835  {
836  newF[0] = f[fp0];
837  newF[1] = f[fp1];
838  newF[2] = addedPoints[nLayers*region1+layerI];
839  newF[3] = addedPoints[nLayers*region0+layerI];
840  }
841  else
842  {
843  newF[0] = addedPoints[nLayers*region0+layerI-1];
844  newF[1] = addedPoints[nLayers*region1+layerI-1];
845  newF[2] = addedPoints[nLayers*region1+layerI];
846  newF[3] = addedPoints[nLayers*region0+layerI];
847  }
848 
849 
850  // Optionally rotate so e[0] is always 0th vertex. Note that
851  // this normally is automatically done by coupled face
852  // ordering but with NOORDERING we have to do it ourselves.
853  if (f[fp0] != e[0])
854  {
855  // rotate one back to get newF[1] (originating
856  // from e[0]) into newF[0].
857  label v0 = newF[0];
858  for (label i = 0; i < newF.size()-1; i++)
859  {
860  newF[i] = newF[newF.fcIndex(i)];
861  }
862  newF.last() = v0;
863  }
865  //{
866  // Pout<< "Adding from MULTI face:"
867  // << patch_.faceCentres()[minFacei]
868  // << " from edge:"
869  // << patch_.localPoints()[f[fp0]]
870  // << patch_.localPoints()[f[fp1]]
871  // << " at layer:" << layerI
872  // << " with new points:" << newF
873  // << " locations:"
874  // << UIndirectList<point>(meshMod.points(), newF)
875  // << endl;
876  //}
877 
878  // newF already outwards pointing.
879  meshMod.addFace
880  (
881  newF, // vertices
882  addedCells[nLayers*minFacei+layerI], // own
883  -1, // nei
884  -1, // masterPointID
885  -1, // masterEdgeID
886  faceToFaceMap.size(), // masterFaceID
887  false, // flipFaceFlux
888  ePatches[i], // patchID
889  -1, // zoneID
890  false // zoneFlip
891  );
892  faceToFaceMap.append(0);
893  faceToEdgeMap.append(edgeI);
894  }
895  }
896  }
897  }
898 
899 
900  cellToFaceMap_.transfer(cellToFaceMap);
901  faceToFaceMap_.transfer(faceToFaceMap);
902  faceToEdgeMap_.transfer(faceToEdgeMap);
903  pointToPointMap_.transfer(pointToPointMap);
904 }
905 
906 
908 {
909  inplaceReorder(map.reverseCellMap(), cellToFaceMap_);
910  inplaceReorder(map.reverseFaceMap(), faceToFaceMap_);
911  inplaceReorder(map.reverseFaceMap(), faceToEdgeMap_);
912  inplaceReorder(map.reversePointMap(), pointToPointMap_);
913 }
914 
915 
916 // ************************************************************************* //
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:579
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:504
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598
const mapDistribute & globalEdgeSlavesMap() const
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
static void matchEdges(const PrimitivePatch< FaceList1, PointField1 > &p1, const PrimitivePatch< FaceList2, PointField2 > &p2, labelList &p1EdgeLabels, labelList &p2EdgeLabels, bitSet &sameOrientation)
Find corresponding edges on patches sharing the same points.
bool empty() const noexcept
True if List is empty (ie, size() is zero)
Definition: UList.H:666
const bitSet & globalEdgeOrientation() const
Is my edge same orientation as master edge.
label constructSize() const noexcept
Constructed data size.
List< labelPair > labelPairList
List of labelPair.
Definition: labelPair.H:33
const labelListList & globalEdgeSlaves() const
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
List< labelList > labelListList
List of labelList.
Definition: labelList.H:38
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:157
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:421
static void syncData(List< Type > &elems, const labelListList &slaves, const labelListList &transformedSlaves, const mapDistribute &slavesMap, const globalIndexAndTransform &, const CombineOp &cop, const TransformOp &top)
Helper: synchronise data with transforms.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
Definition: UListI.H:97
List< face > faceList
List of faces.
Definition: faceListFwd.H:39
scalar y
A list of faces which address into the list of points.
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
void setSize(const label n)
Alias for resize()
Definition: List.H:316
Creates mesh by extruding a patch.
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicList.H:584
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:26
label find(const T &val) const
Find index of the first occurrence of the value.
Definition: UList.C:173
static void calcPointRegions(const globalMeshData &globalData, const primitiveFacePatch &patch, const bitSet &nonManifoldEdge, const bool syncNonCollocated, faceList &pointGlobalRegions, faceList &pointLocalRegions, labelList &localToGlobalRegion)
Helper: calculate point regions. The point region is the.
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:653
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:51
void updateMesh(const mapPolyMesh &)
Update any locally stored mesh information.
defineTypeNameAndDebug(combustionModel, 0)
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:405
labelList f(nPoints)
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
T & last()
Access last element of the list, position [size()-1].
Definition: UList.H:867
void operator()(T &x, const T &y) const
Definition: ops.H:76
Class containing processor-to-processor mapping information.
vector point
Point is a vector.
Definition: point.H:37
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:59
const std::string patch
OpenFOAM patch number as a std::string.
constexpr label labelMax
Definition: label.H:55
void setRefinement(const pointField &firstLayerThickness, const scalar expansionRatio, const label nLayers, const labelList &topPatchID, const labelList &bottomPatchID, const labelListList &extrudeEdgePatches, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layer mesh.
List< label > labelList
A List of labels.
Definition: List.H:62
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
Definition: UListI.H:104
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:617
Namespace for OpenFOAM.
const labelListList & globalEdgeTransformedSlaves() const