addPatchCellLayer.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-2018 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 "addPatchCellLayer.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "meshTools.H"
30 #include "mapPolyMesh.H"
31 #include "syncTools.H"
32 #include "polyAddPoint.H"
33 #include "polyAddFace.H"
34 #include "polyModifyFace.H"
35 #include "polyAddCell.H"
36 #include "globalIndex.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(addPatchCellLayer, 0);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 Foam::label Foam::addPatchCellLayer::nbrFace
49 (
50  const labelListList& edgeFaces,
51  const label edgeI,
52  const label facei
53 )
54 {
55  const labelList& eFaces = edgeFaces[edgeI];
56 
57  if (eFaces.size() == 2)
58  {
59  return (eFaces[0] != facei ? eFaces[0] : eFaces[1]);
60  }
61  else
62  {
63  return -1;
64  }
65 }
66 
67 
68 void Foam::addPatchCellLayer::addVertex
69 (
70  const label pointi,
71  face& f,
72  label& fp
73 )
74 {
75  if (fp == 0)
76  {
77  f[fp++] = pointi;
78  }
79  else
80  {
81  if (f[fp-1] != pointi && f[0] != pointi)
82  {
83  f[fp++] = pointi;
84  }
85  }
86 }
87 
88 
89 // Is edge to the same neighbour? (and needs extrusion and has not been
90 // dealt with already)
91 bool Foam::addPatchCellLayer::sameEdgeNeighbour
92 (
93  const indirectPrimitivePatch& pp,
94  const labelListList& globalEdgeFaces,
95  const boolList& doneEdge,
96  const label thisGlobalFacei,
97  const label nbrGlobalFacei,
98  const label edgeI
99 ) const
100 {
101  const edge& e = pp.edges()[edgeI];
102 
103  return
104  !doneEdge[edgeI] // not yet handled
105  && (
106  addedPoints_[e[0]].size() // is extruded
107  || addedPoints_[e[1]].size()
108  )
109  && (
110  nbrFace(globalEdgeFaces, edgeI, thisGlobalFacei)
111  == nbrGlobalFacei // is to same neighbour
112  );
113 }
114 
115 
116 // Collect consecutive string of edges that connects the same two
117 // (possibly coupled) faces. Returns -1 if no unvisited edge can be found.
118 // Otherwise returns start and end index in face.
119 Foam::labelPair Foam::addPatchCellLayer::getEdgeString
120 (
121  const indirectPrimitivePatch& pp,
122  const labelListList& globalEdgeFaces,
123  const boolList& doneEdge,
124  const label patchFacei,
125  const label globalFacei
126 ) const
127 {
128  const labelList& fEdges = pp.faceEdges()[patchFacei];
129 
130  label startFp = -1;
131  label endFp = -1;
132 
133  // Get edge that hasn't been done yet but needs extrusion
134  forAll(fEdges, fp)
135  {
136  label edgeI = fEdges[fp];
137  const edge& e = pp.edges()[edgeI];
138 
139  if
140  (
141  !doneEdge[edgeI]
142  && ( addedPoints_[e[0]].size() || addedPoints_[e[1]].size() )
143  )
144  {
145  startFp = fp;
146  break;
147  }
148  }
149 
150  if (startFp != -1)
151  {
152  // We found an edge that needs extruding but hasn't been done yet.
153  // Now find the face on the other side
154  label nbrGlobalFacei = nbrFace
155  (
156  globalEdgeFaces,
157  fEdges[startFp],
158  globalFacei
159  );
160 
161  if (nbrGlobalFacei == -1)
162  {
163  // Proper boundary edge. Only extrude single edge.
164  endFp = startFp;
165  }
166  else
167  {
168  // Search back for edge
169  // - which hasn't been handled yet
170  // - with same neighbour
171  // - that needs extrusion
172  while (true)
173  {
174  label prevFp = fEdges.rcIndex(startFp);
175 
176  if
177  (
178  !sameEdgeNeighbour
179  (
180  pp,
181  globalEdgeFaces,
182  doneEdge,
183  globalFacei,
184  nbrGlobalFacei,
185  fEdges[prevFp]
186  )
187  )
188  {
189  break;
190  }
191  startFp = prevFp;
192  }
193 
194  // Search forward for end of string
195  endFp = startFp;
196  while (true)
197  {
198  label nextFp = fEdges.fcIndex(endFp);
199 
200  if
201  (
202  !sameEdgeNeighbour
203  (
204  pp,
205  globalEdgeFaces,
206  doneEdge,
207  globalFacei,
208  nbrGlobalFacei,
209  fEdges[nextFp]
210  )
211  )
212  {
213  break;
214  }
215  endFp = nextFp;
216  }
217  }
218  }
219 
220  return labelPair(startFp, endFp);
221 }
222 
223 
224 // Adds a side face i.e. extrudes a patch edge.
225 Foam::label Foam::addPatchCellLayer::addSideFace
226 (
227  const indirectPrimitivePatch& pp,
228  const labelListList& addedCells, // per pp face the new extruded cell
229  const face& newFace,
230  const label newPatchID,
231 
232  const label ownFacei, // pp face that provides owner
233  const label nbrFacei,
234  const label meshEdgeI, // corresponding mesh edge
235  const label layerI, // layer
236  const label numEdgeFaces, // number of layers for edge
237  const labelList& meshFaces, // precalculated edgeFaces
238  polyTopoChange& meshMod
239 ) const
240 {
241  // Face or edge to 'inflate' from
242  label inflateEdgeI = -1;
243  label inflateFacei = -1;
244 
245  // Check mesh faces using edge
246  if (addToMesh_)
247  {
248  forAll(meshFaces, i)
249  {
250  if (mesh_.isInternalFace(meshFaces[i]))
251  {
252  // meshEdge uses internal faces so ok to inflate from it
253  inflateEdgeI = meshEdgeI;
254  break;
255  }
256  }
257  }
258 
259  // Zone info comes from any side patch face. Otherwise -1 since we
260  // don't know what to put it in - inherit from the extruded faces?
261  label zoneI = -1; // mesh_.faceZones().whichZone(meshFacei);
262  bool flip = false;
263 
264  label addedFacei = -1;
265 
266  // Is patch edge external edge of indirectPrimitivePatch?
267  if (nbrFacei == -1)
268  {
269  // External edge so external face.
270 
271  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
272 
273  // Loop over all faces connected to edge to inflate and
274  // see if we can find a face that is otherPatchID
275 
276  // Get my mesh face and its zone.
277  label meshFacei = pp.addressing()[ownFacei];
278 
279  forAll(meshFaces, k)
280  {
281  label facei = meshFaces[k];
282 
283  if
284  (
285  (facei != meshFacei)
286  && (patches.whichPatch(facei) == newPatchID)
287  )
288  {
289  // Found the patch face. Use it to inflate from
290  inflateEdgeI = -1;
291  inflateFacei = facei;
292 
293  zoneI = mesh_.faceZones().whichZone(facei);
294  if (zoneI != -1)
295  {
296  label index = mesh_.faceZones()[zoneI].whichFace(facei);
297  flip = mesh_.faceZones()[zoneI].flipMap()[index];
298  }
299  break;
300  }
301  }
302 
303  // Determine if different number of layer on owner and neighbour side
304  // (relevant only for coupled faces). See section for internal edge
305  // below.
306 
307  label layerOwn;
308 
309  if (addedCells[ownFacei].size() < numEdgeFaces)
310  {
311  label offset = numEdgeFaces - addedCells[ownFacei].size();
312  if (layerI <= offset)
313  {
314  layerOwn = 0;
315  }
316  else
317  {
318  layerOwn = layerI - offset;
319  }
320  }
321  else
322  {
323  layerOwn = layerI;
324  }
325 
326 
327  // Pout<< "Added boundary face:" << newFace
328  // << " own:" << addedCells[ownFacei][layerOwn]
329  // << " patch:" << newPatchID
330  // << endl;
331 
332  addedFacei = meshMod.setAction
333  (
334  polyAddFace
335  (
336  newFace, // face
337  addedCells[ownFacei][layerOwn], // owner
338  -1, // neighbour
339  -1, // master point
340  inflateEdgeI, // master edge
341  inflateFacei, // master face
342  false, // flux flip
343  newPatchID, // patch for face
344  zoneI, // zone for face
345  flip // face zone flip
346  )
347  );
348  }
349  else
350  {
351  // When adding side faces we need to modify neighbour and owners
352  // in region where layer mesh is stopped. Determine which side
353  // has max number of faces and make sure layers match closest to
354  // original pp if there are different number of layers.
355 
356  label layerNbr;
357  label layerOwn;
358 
359  if (addedCells[ownFacei].size() > addedCells[nbrFacei].size())
360  {
361  label offset =
362  addedCells[ownFacei].size() - addedCells[nbrFacei].size();
363 
364  layerOwn = layerI;
365 
366  if (layerI <= offset)
367  {
368  layerNbr = 0;
369  }
370  else
371  {
372  layerNbr = layerI - offset;
373  }
374  }
375  else if (addedCells[nbrFacei].size() > addedCells[ownFacei].size())
376  {
377  label offset =
378  addedCells[nbrFacei].size() - addedCells[ownFacei].size();
379 
380  layerNbr = layerI;
381 
382  if (layerI <= offset)
383  {
384  layerOwn = 0;
385  }
386  else
387  {
388  layerOwn = layerI - offset;
389  }
390  }
391  else
392  {
393  // Same number of layers on both sides.
394  layerNbr = layerI;
395  layerOwn = layerI;
396  }
397 
398  addedFacei = meshMod.setAction
399  (
400  polyAddFace
401  (
402  newFace, // face
403  addedCells[ownFacei][layerOwn], // owner
404  addedCells[nbrFacei][layerNbr], // neighbour
405  -1, // master point
406  inflateEdgeI, // master edge
407  -1, // master face
408  false, // flux flip
409  -1, // patch for face
410  zoneI, // zone for face
411  flip // face zone flip
412  )
413  );
414 
415  // Pout<< "Added internal face:" << newFace
416  // << " own:" << addedCells[ownFacei][layerOwn]
417  // << " nei:" << addedCells[nbrFacei][layerNbr]
418  // << endl;
419  }
420 
421  return addedFacei;
422 }
423 
424 
425 Foam::label Foam::addPatchCellLayer::findProcPatch
426 (
427  const polyMesh& mesh,
428  const label nbrProcID
429 )
430 {
431  const polyBoundaryMesh& patches = mesh.boundaryMesh();
432 
433  forAll(mesh.globalData().processorPatches(), i)
434  {
435  label patchi = mesh.globalData().processorPatches()[i];
436 
437  if
438  (
439  refCast<const processorPolyPatch>(patches[patchi]).neighbProcNo()
440  == nbrProcID
441  )
442  {
443  return patchi;
444  }
445  }
446  return -1;
447 }
448 
449 
450 void Foam::addPatchCellLayer::setFaceProps
451 (
452  const polyMesh& mesh,
453  const label facei,
454 
455  label& patchi,
456  label& zoneI,
457  bool& zoneFlip
458 )
459 {
460  patchi = mesh.boundaryMesh().whichPatch(facei);
461  zoneI = mesh.faceZones().whichZone(facei);
462  if (zoneI != -1)
463  {
464  label index = mesh.faceZones()[zoneI].whichFace(facei);
465  zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
466  }
467 }
468 
469 
470 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
471 
472 // Construct from mesh
474 (
475  const polyMesh& mesh,
476  const bool addToMesh
477 )
478 :
479  mesh_(mesh),
480  addToMesh_(addToMesh),
481  addedPoints_(0),
482  layerFaces_(0)
483 {}
484 
485 
486 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
487 
489 (
490  const polyMesh& mesh,
491  const labelListList& layerFaces
492 )
493 {
494  labelListList layerCells(layerFaces.size());
495 
496  forAll(layerFaces, patchFacei)
497  {
498  const labelList& faceLabels = layerFaces[patchFacei];
499 
500  if (faceLabels.size())
501  {
502  labelList& added = layerCells[patchFacei];
503  added.setSize(faceLabels.size()-1);
504 
505  for (label i = 0; i < faceLabels.size()-1; i++)
506  {
507  added[i] = mesh.faceNeighbour()[faceLabels[i]];
508  }
509  }
510  }
511  return layerCells;
512 }
513 
514 
516 {
517  return addedCells(mesh_, layerFaces_);
518 }
519 
520 
521 // Calculate global faces per pp edge.
523 (
524  const polyMesh& mesh,
525  const globalIndex& globalFaces,
526  const indirectPrimitivePatch& pp
527 )
528 {
529  // Precalculate mesh edges for pp.edges.
530  const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
531 
532  // From mesh edge to global face labels. Non-empty sublists only for
533  // pp edges.
534  labelListList globalEdgeFaces(mesh.nEdges());
535 
536  const labelListList& edgeFaces = pp.edgeFaces();
537 
538  forAll(edgeFaces, edgeI)
539  {
540  label meshEdgeI = meshEdges[edgeI];
541 
542  const labelList& eFaces = edgeFaces[edgeI];
543 
544  // Store face and processor as unique tag.
545  labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
546  globalEFaces.setSize(eFaces.size());
547  forAll(eFaces, i)
548  {
549  globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
550  }
551  }
552 
553  // Synchronise across coupled edges.
555  (
556  mesh,
557  globalEdgeFaces,
558  uniqueEqOp(),
559  labelList() // null value
560  );
561 
562  // Extract pp part
563  return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
564 }
565 
566 
568 (
569  const polyMesh& mesh,
570  const globalIndex& globalFaces,
571  const labelListList& globalEdgeFaces,
572  const indirectPrimitivePatch& pp,
573 
574  labelList& sidePatchID,
575  label& nPatches,
576  Map<label>& nbrProcToPatch,
577  Map<label>& patchToNbrProc
578 )
579 {
580  const polyBoundaryMesh& patches = mesh.boundaryMesh();
581 
582  // Precalculate mesh edges for pp.edges.
583  const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
584 
585  sidePatchID.setSize(pp.nEdges());
586  sidePatchID = -1;
587 
588  // These also get determined but not (yet) exported:
589  // - whether face is created from other face or edge
590  // - what zone&orientation face should have
591 
592  labelList inflateEdgeI(pp.nEdges(), -1);
593  labelList inflateFacei(pp.nEdges(), -1);
594  labelList sideZoneID(pp.nEdges(), -1);
595  boolList sideFlip(pp.nEdges(), false);
596 
597  nPatches = patches.size();
598 
599  forAll(globalEdgeFaces, edgeI)
600  {
601  const labelList& eGlobalFaces = globalEdgeFaces[edgeI];
602  if
603  (
604  eGlobalFaces.size() == 2
605  && pp.edgeFaces()[edgeI].size() == 1
606  )
607  {
608  // Locally but not globally a boundary edge. Hence a coupled
609  // edge. Find the patch to use if on different
610  // processors.
611 
612  label f0 = eGlobalFaces[0];
613  label f1 = eGlobalFaces[1];
614 
615  label otherProci = -1;
616  if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1))
617  {
618  otherProci = globalFaces.whichProcID(f1);
619  }
620  else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1))
621  {
622  otherProci = globalFaces.whichProcID(f0);
623  }
624 
625 
626  if (otherProci != -1)
627  {
628  sidePatchID[edgeI] = findProcPatch(mesh, otherProci);
629  if (sidePatchID[edgeI] == -1)
630  {
631  // Cannot find a patch to processor. See if already
632  // marked for addition
633  if (nbrProcToPatch.found(otherProci))
634  {
635  sidePatchID[edgeI] = nbrProcToPatch[otherProci];
636  }
637  else
638  {
639  sidePatchID[edgeI] = nPatches;
640  nbrProcToPatch.insert(otherProci, nPatches);
641  patchToNbrProc.insert(nPatches, otherProci);
642  nPatches++;
643  }
644  }
645  }
646  }
647  }
648 
649 
650 
651  // Determine face properties for all other boundary edges
652  // ------------------------------------------------------
653 
654  const labelListList& edgeFaces = pp.edgeFaces();
655 
656  DynamicList<label> dynMeshEdgeFaces;
657 
658  forAll(edgeFaces, edgeI)
659  {
660  if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
661  {
662  // Proper, uncoupled patch edge.
663 
664  label myFacei = pp.addressing()[edgeFaces[edgeI][0]];
665 
666  // Pick up any boundary face on this edge and use its properties
667  label meshEdgeI = meshEdges[edgeI];
668  const labelList& meshFaces = mesh.edgeFaces
669  (
670  meshEdgeI,
671  dynMeshEdgeFaces
672  );
673 
674  forAll(meshFaces, k)
675  {
676  label facei = meshFaces[k];
677 
678  if (facei != myFacei && !mesh.isInternalFace(facei))
679  {
680  setFaceProps
681  (
682  mesh,
683  facei,
684 
685  sidePatchID[edgeI],
686  sideZoneID[edgeI],
687  sideFlip[edgeI]
688  );
689  inflateFacei[edgeI] = facei;
690  inflateEdgeI[edgeI] = -1;
691 
692  break;
693  }
694  }
695  }
696  }
697 
698 
699 
700  // Now hopefully every boundary edge has a side patch. Check
701  if (debug)
702  {
703  forAll(edgeFaces, edgeI)
704  {
705  if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
706  {
707  const edge& e = pp.edges()[edgeI];
709  << "Have no sidePatchID for edge " << edgeI << " points "
710  << pp.points()[pp.meshPoints()[e[0]]]
711  << pp.points()[pp.meshPoints()[e[1]]]
712  << endl;
713  }
714  }
715  }
716 
717 
718 
719  // Now we have sidepatch see if we have patchface or edge to inflate
720  // from.
721  forAll(edgeFaces, edgeI)
722  {
723  if
724  (
725  edgeFaces[edgeI].size() == 1
726  && sidePatchID[edgeI] != -1
727  && inflateFacei[edgeI] == -1
728  )
729  {
730  // 1. Do we have a boundary face to inflate from
731 
732  label myFacei = pp.addressing()[edgeFaces[edgeI][0]];
733 
734  // Pick up any boundary face on this edge and use its properties
735  label meshEdgeI = meshEdges[edgeI];
736  const labelList& meshFaces = mesh.edgeFaces
737  (
738  meshEdgeI,
739  dynMeshEdgeFaces
740  );
741 
742  forAll(meshFaces, k)
743  {
744  label facei = meshFaces[k];
745 
746  if (facei != myFacei)
747  {
748  if (mesh.isInternalFace(facei))
749  {
750  inflateEdgeI[edgeI] = meshEdgeI;
751  }
752  else
753  {
754  if (patches.whichPatch(facei) == sidePatchID[edgeI])
755  {
756  setFaceProps
757  (
758  mesh,
759  facei,
760 
761  sidePatchID[edgeI],
762  sideZoneID[edgeI],
763  sideFlip[edgeI]
764  );
765  inflateFacei[edgeI] = facei;
766  inflateEdgeI[edgeI] = -1;
767 
768  break;
769  }
770  }
771  }
772  }
773  }
774  }
775 }
776 
777 
779 (
780  const globalIndex& globalFaces,
781  const labelListList& globalEdgeFaces,
782  const scalarField& expansionRatio,
783  const indirectPrimitivePatch& pp,
784  const labelList& sidePatchID,
785  const labelList& exposedPatchID,
786  const labelList& nFaceLayers,
787  const labelList& nPointLayers,
788  const vectorField& firstLayerDisp,
789  polyTopoChange& meshMod
790 )
791 {
792  if (debug)
793  {
794  Pout<< "addPatchCellLayer::setRefinement : Adding up to "
795  << gMax(nPointLayers)
796  << " layers of cells to indirectPrimitivePatch with "
797  << pp.nPoints() << " points" << endl;
798  }
799 
800  if
801  (
802  pp.nPoints() != firstLayerDisp.size()
803  || pp.nPoints() != nPointLayers.size()
804  || pp.size() != nFaceLayers.size()
805  )
806  {
808  << "Size of new points is not same as number of points used by"
809  << " the face subset" << endl
810  << " patch.nPoints:" << pp.nPoints()
811  << " displacement:" << firstLayerDisp.size()
812  << " nPointLayers:" << nPointLayers.size() << nl
813  << " patch.nFaces:" << pp.size()
814  << " nFaceLayers:" << nFaceLayers.size()
815  << abort(FatalError);
816  }
817 
818  forAll(nPointLayers, i)
819  {
820  if (nPointLayers[i] < 0)
821  {
823  << "Illegal number of layers " << nPointLayers[i]
824  << " at patch point " << i << abort(FatalError);
825  }
826  }
827  forAll(nFaceLayers, i)
828  {
829  if (nFaceLayers[i] < 0)
830  {
832  << "Illegal number of layers " << nFaceLayers[i]
833  << " at patch face " << i << abort(FatalError);
834  }
835  }
836 
837  forAll(globalEdgeFaces, edgeI)
838  {
839  if (globalEdgeFaces[edgeI].size() > 2)
840  {
841  const edge& e = pp.edges()[edgeI];
842 
843  if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
844  {
846  << "Trying to extrude edge "
847  << e.line(pp.localPoints())
848  << " which is non-manifold (has "
849  << globalEdgeFaces[edgeI].size()
850  << " faces using it)"
851  << abort(FatalError);
852  }
853  }
854  }
855 
856 
857  const labelList& meshPoints = pp.meshPoints();
858 
859  // Some storage for edge-face-addressing.
861 
862  // Precalculate mesh edges for pp.edges.
863  const labelList meshEdges(pp.meshEdges(mesh_.edges(), mesh_.pointEdges()));
864 
865  if (debug)
866  {
867  // Check synchronisation
868  // ~~~~~~~~~~~~~~~~~~~~~
869 
870  {
871  labelList n(mesh_.nPoints(), 0);
872  UIndirectList<label>(n, meshPoints) = nPointLayers;
874 
875  // Non-synced
876  forAll(meshPoints, i)
877  {
878  label meshPointi = meshPoints[i];
879 
880  if (n[meshPointi] != nPointLayers[i])
881  {
883  << "At mesh point:" << meshPointi
884  << " coordinate:" << mesh_.points()[meshPointi]
885  << " specified nLayers:" << nPointLayers[i] << endl
886  << "On coupled point a different nLayers:"
887  << n[meshPointi] << " was specified."
888  << abort(FatalError);
889  }
890  }
891 
892 
893  // Check that nPointLayers equals the max layers of connected faces
894  // (or 0). Anything else makes no sense.
895  labelList nFromFace(mesh_.nPoints(), 0);
896  forAll(nFaceLayers, i)
897  {
898  const face& f = pp[i];
899 
900  forAll(f, fp)
901  {
902  label pointi = f[fp];
903 
904  nFromFace[pointi] = max(nFromFace[pointi], nFaceLayers[i]);
905  }
906  }
908  (
909  mesh_,
910  nFromFace,
911  maxEqOp<label>(),
912  label(0)
913  );
914 
915  forAll(nPointLayers, i)
916  {
917  label meshPointi = meshPoints[i];
918 
919  if
920  (
921  nPointLayers[i] > 0
922  && nPointLayers[i] != nFromFace[meshPointi]
923  )
924  {
926  << "At mesh point:" << meshPointi
927  << " coordinate:" << mesh_.points()[meshPointi]
928  << " specified nLayers:" << nPointLayers[i] << endl
929  << "but the max nLayers of surrounding faces is:"
930  << nFromFace[meshPointi]
931  << abort(FatalError);
932  }
933  }
934  }
935 
936  {
937  pointField d(mesh_.nPoints(), vector::max);
938  UIndirectList<point>(d, meshPoints) = firstLayerDisp;
940  (
941  mesh_,
942  d,
943  minEqOp<vector>(),
945  );
946 
947  forAll(meshPoints, i)
948  {
949  label meshPointi = meshPoints[i];
950 
951  if (mag(d[meshPointi] - firstLayerDisp[i]) > small)
952  {
954  << "At mesh point:" << meshPointi
955  << " coordinate:" << mesh_.points()[meshPointi]
956  << " specified displacement:" << firstLayerDisp[i]
957  << endl
958  << "On coupled point a different displacement:"
959  << d[meshPointi] << " was specified."
960  << abort(FatalError);
961  }
962  }
963  }
964 
965  // Check that edges of pp (so ones that become boundary faces)
966  // connect to only one boundary face. Guarantees uniqueness of
967  // patch that they go into so if this is a coupled patch both
968  // sides decide the same.
969  // ~~~~~~~~~~~~~~~~~~~~~~
970 
971  for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
972  {
973  const edge& e = pp.edges()[edgeI];
974 
975  if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
976  {
977  // Edge is to become a face
978 
979  const labelList& eFaces = pp.edgeFaces()[edgeI];
980 
981  // First check: pp should be single connected.
982  if (eFaces.size() != 1)
983  {
985  << "boundary-edge-to-be-extruded:"
986  << pp.points()[meshPoints[e[0]]]
987  << pp.points()[meshPoints[e[1]]]
988  << " has more than two faces using it:" << eFaces
989  << abort(FatalError);
990  }
991 
992  label myFacei = pp.addressing()[eFaces[0]];
993 
994  label meshEdgeI = meshEdges[edgeI];
995 
996  // Mesh faces using edge
997  const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
998 
999  // Check that there is only one patchface using edge.
1000  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1001 
1002  label bFacei = -1;
1003 
1004  forAll(meshFaces, i)
1005  {
1006  label facei = meshFaces[i];
1007 
1008  if (facei != myFacei)
1009  {
1010  if (!mesh_.isInternalFace(facei))
1011  {
1012  if (bFacei == -1)
1013  {
1014  bFacei = facei;
1015  }
1016  else
1017  {
1019  << "boundary-edge-to-be-extruded:"
1020  << pp.points()[meshPoints[e[0]]]
1021  << pp.points()[meshPoints[e[1]]]
1022  << " has more than two boundary faces"
1023  << " using it:"
1024  << bFacei << " fc:"
1025  << mesh_.faceCentres()[bFacei]
1026  << " patch:" << patches.whichPatch(bFacei)
1027  << " and " << facei << " fc:"
1028  << mesh_.faceCentres()[facei]
1029  << " patch:" << patches.whichPatch(facei)
1030  << abort(FatalError);
1031  }
1032  }
1033  }
1034  }
1035  }
1036  }
1037  }
1038 
1039 
1040  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1041 
1042  // Precalculated patchID for each patch face
1043  labelList patchID(pp.size());
1044 
1045  forAll(pp, patchFacei)
1046  {
1047  label meshFacei = pp.addressing()[patchFacei];
1048 
1049  patchID[patchFacei] = patches.whichPatch(meshFacei);
1050  }
1051 
1052 
1053  // From master point (in patch point label) to added points (in mesh point
1054  // label)
1055  addedPoints_.setSize(pp.nPoints());
1056 
1057  // Mark points that do not get extruded by setting size of addedPoints_ to 0
1058  label nTruncated = 0;
1059 
1060  forAll(nPointLayers, patchPointi)
1061  {
1062  if (nPointLayers[patchPointi] > 0)
1063  {
1064  addedPoints_[patchPointi].setSize(nPointLayers[patchPointi]);
1065  }
1066  else
1067  {
1068  nTruncated++;
1069  }
1070  }
1071 
1072  if (debug)
1073  {
1074  Pout<< "Not adding points at " << nTruncated << " out of "
1075  << pp.nPoints() << " points" << endl;
1076  }
1077 
1078 
1079  //
1080  // Create new points
1081  //
1082 
1083  // If creating new mesh: copy existing patch points
1084  labelList copiedPatchPoints;
1085  if (!addToMesh_)
1086  {
1087  copiedPatchPoints.setSize(firstLayerDisp.size());
1088  forAll(firstLayerDisp, patchPointi)
1089  {
1090  if (addedPoints_[patchPointi].size())
1091  {
1092  label meshPointi = meshPoints[patchPointi];
1093  label zoneI = mesh_.pointZones().whichZone(meshPointi);
1094  copiedPatchPoints[patchPointi] = meshMod.setAction
1095  (
1096  polyAddPoint
1097  (
1098  mesh_.points()[meshPointi], // point
1099  -1, // master point
1100  zoneI, // zone for point
1101  true // supports a cell
1102  )
1103  );
1104  }
1105  }
1106  }
1107 
1108 
1109  // Create points for additional layers
1110  forAll(firstLayerDisp, patchPointi)
1111  {
1112  if (addedPoints_[patchPointi].size())
1113  {
1114  label meshPointi = meshPoints[patchPointi];
1115 
1116  label zoneI = mesh_.pointZones().whichZone(meshPointi);
1117 
1118  point pt = mesh_.points()[meshPointi];
1119 
1120  vector disp = firstLayerDisp[patchPointi];
1121 
1122  forAll(addedPoints_[patchPointi], i)
1123  {
1124  pt += disp;
1125 
1126  label addedVertI = meshMod.setAction
1127  (
1128  polyAddPoint
1129  (
1130  pt, // point
1131  (addToMesh_ ? meshPointi : -1), // master point
1132  zoneI, // zone for point
1133  true // supports a cell
1134  )
1135  );
1136 
1137  addedPoints_[patchPointi][i] = addedVertI;
1138 
1139  disp *= expansionRatio[patchPointi];
1140  }
1141  }
1142  }
1143 
1144 
1145  //
1146  // Add cells to all boundaryFaces
1147  //
1148 
1149  labelListList addedCells(pp.size());
1150 
1151  forAll(pp, patchFacei)
1152  {
1153  if (nFaceLayers[patchFacei] > 0)
1154  {
1155  addedCells[patchFacei].setSize(nFaceLayers[patchFacei]);
1156 
1157  label meshFacei = pp.addressing()[patchFacei];
1158 
1159  label ownZoneI = mesh_.cellZones().whichZone
1160  (
1161  mesh_.faceOwner()[meshFacei]
1162  );
1163 
1164  for (label i = 0; i < nFaceLayers[patchFacei]; i++)
1165  {
1166  // Note: add from cell (owner of patch face) or from face?
1167  // for now add from cell so we can map easily.
1168  addedCells[patchFacei][i] = meshMod.setAction
1169  (
1170  polyAddCell
1171  (
1172  -1, // master point
1173  -1, // master edge
1174  -1, // master face
1175  (addToMesh_ ? mesh_.faceOwner()[meshFacei] : -1),
1176  // master
1177  ownZoneI // zone for cell
1178  )
1179  );
1180  }
1181  }
1182  }
1183 
1184 
1185 
1186  // Create faces on top of the original patch faces.
1187  // These faces are created from original patch faces outwards so in order
1188  // of increasing cell number. So orientation should be same as original
1189  // patch face for them to have owner<neighbour.
1190 
1191  layerFaces_.setSize(pp.size());
1192 
1193  forAll(pp.localFaces(), patchFacei)
1194  {
1195  label meshFacei = pp.addressing()[patchFacei];
1196 
1197  if (addedCells[patchFacei].size())
1198  {
1199  layerFaces_[patchFacei].setSize(addedCells[patchFacei].size() + 1);
1200 
1201  // Get duplicated vertices on the patch face.
1202  const face& f = pp.localFaces()[patchFacei];
1203 
1204  face newFace(f.size());
1205 
1206  forAll(addedCells[patchFacei], i)
1207  {
1208  forAll(f, fp)
1209  {
1210  if (addedPoints_[f[fp]].empty())
1211  {
1212  // Keep original point
1213  newFace[fp] =
1214  (
1215  addToMesh_
1216  ? meshPoints[f[fp]]
1217  : copiedPatchPoints[f[fp]]
1218  );
1219  }
1220  else
1221  {
1222  // Get new outside point
1223  label offset =
1224  addedPoints_[f[fp]].size()
1225  - addedCells[patchFacei].size();
1226  newFace[fp] = addedPoints_[f[fp]][i+offset];
1227  }
1228  }
1229 
1230 
1231  // Get new neighbour
1232  label nei;
1233  label patchi;
1234  label zoneI = -1;
1235  bool flip = false;
1236 
1237 
1238  if (i == addedCells[patchFacei].size()-1)
1239  {
1240  // Top layer so is patch face.
1241  nei = -1;
1242  patchi = patchID[patchFacei];
1243  zoneI = mesh_.faceZones().whichZone(meshFacei);
1244  if (zoneI != -1)
1245  {
1246  const faceZone& fz = mesh_.faceZones()[zoneI];
1247  flip = fz.flipMap()[fz.whichFace(meshFacei)];
1248  }
1249  }
1250  else
1251  {
1252  // Internal face between layer i and i+1
1253  nei = addedCells[patchFacei][i+1];
1254  patchi = -1;
1255  }
1256 
1257 
1258  layerFaces_[patchFacei][i+1] = meshMod.setAction
1259  (
1260  polyAddFace
1261  (
1262  newFace, // face
1263  addedCells[patchFacei][i], // owner
1264  nei, // neighbour
1265  -1, // master point
1266  -1, // master edge
1267  (addToMesh_ ? meshFacei : -1), // master face
1268  false, // flux flip
1269  patchi, // patch for face
1270  zoneI, // zone for face
1271  flip // face zone flip
1272  )
1273  );
1274  }
1275  }
1276  }
1277 
1278  //
1279  // Modify old patch faces to be on the inside
1280  //
1281 
1282  if (addToMesh_)
1283  {
1284  forAll(pp, patchFacei)
1285  {
1286  if (addedCells[patchFacei].size())
1287  {
1288  label meshFacei = pp.addressing()[patchFacei];
1289 
1290  layerFaces_[patchFacei][0] = meshFacei;
1291 
1292  meshMod.setAction
1293  (
1295  (
1296  pp[patchFacei], // modified face
1297  meshFacei, // label of face
1298  mesh_.faceOwner()[meshFacei], // owner
1299  addedCells[patchFacei][0], // neighbour
1300  false, // face flip
1301  -1, // patch for face
1302  true, // false, // remove from zone
1303  -1, // zoneI, // zone for face
1304  false // face flip in zone
1305  )
1306  );
1307  }
1308  }
1309  }
1310  else
1311  {
1312  // If creating new mesh: reverse original faces and put them
1313  // in the exposed patch ID.
1314  forAll(pp, patchFacei)
1315  {
1316  if (nFaceLayers[patchFacei] > 0)
1317  {
1318  label meshFacei = pp.addressing()[patchFacei];
1319  label zoneI = mesh_.faceZones().whichZone(meshFacei);
1320  bool zoneFlip = false;
1321  if (zoneI != -1)
1322  {
1323  const faceZone& fz = mesh_.faceZones()[zoneI];
1324  zoneFlip = !fz.flipMap()[fz.whichFace(meshFacei)];
1325  }
1326 
1327  // Reverse and renumber old patch face.
1328  face f(pp.localFaces()[patchFacei].reverseFace());
1329  forAll(f, fp)
1330  {
1331  f[fp] = copiedPatchPoints[f[fp]];
1332  }
1333 
1334  layerFaces_[patchFacei][0] = meshMod.setAction
1335  (
1336  polyAddFace
1337  (
1338  f, // modified face
1339  addedCells[patchFacei][0], // owner
1340  -1, // neighbour
1341  -1, // masterPoint
1342  -1, // masterEdge
1343  -1, // masterFace
1344  true, // face flip
1345  exposedPatchID[patchFacei], // patch for face
1346  zoneI, // zone for face
1347  zoneFlip // face flip in zone
1348  )
1349  );
1350  }
1351  }
1352  }
1353 
1354 
1355 
1356  //
1357  // Create 'side' faces, one per edge that is being extended.
1358  //
1359 
1360  const labelListList& faceEdges = pp.faceEdges();
1361  const faceList& localFaces = pp.localFaces();
1362  const edgeList& edges = pp.edges();
1363 
1364  // Get number of layers per edge. This is 0 if edge is not extruded;
1365  // max of connected faces otherwise.
1366  labelList edgeLayers(pp.nEdges());
1367 
1368  {
1369  // Use list over mesh.nEdges() since syncTools does not yet support
1370  // partial list synchronisation.
1371  labelList meshEdgeLayers(mesh_.nEdges(), -1);
1372 
1373  forAll(meshEdges, edgeI)
1374  {
1375  const edge& e = edges[edgeI];
1376 
1377  label meshEdgeI = meshEdges[edgeI];
1378 
1379  if ((nPointLayers[e[0]] == 0) && (nPointLayers[e[1]] == 0))
1380  {
1381  meshEdgeLayers[meshEdgeI] = 0;
1382  }
1383  else
1384  {
1385  const labelList& eFaces = pp.edgeFaces()[edgeI];
1386 
1387  forAll(eFaces, i)
1388  {
1389  meshEdgeLayers[meshEdgeI] = max
1390  (
1391  nFaceLayers[eFaces[i]],
1392  meshEdgeLayers[meshEdgeI]
1393  );
1394  }
1395  }
1396  }
1397 
1399  (
1400  mesh_,
1401  meshEdgeLayers,
1402  maxEqOp<label>(),
1403  label(0) // initial value
1404  );
1405 
1406  forAll(meshEdges, edgeI)
1407  {
1408  edgeLayers[edgeI] = meshEdgeLayers[meshEdges[edgeI]];
1409  }
1410  }
1411 
1412 
1413  // Mark off which edges have been extruded
1414  boolList doneEdge(pp.nEdges(), false);
1415 
1416 
1417  // Create faces. Per face walk connected edges and find string of edges
1418  // between the same two faces and extrude string into a single face.
1419  forAll(pp, patchFacei)
1420  {
1421  const labelList& fEdges = faceEdges[patchFacei];
1422 
1423  forAll(fEdges, fp)
1424  {
1425  // Get string of edges that needs to be extruded as a single face.
1426  // Returned as indices in fEdges.
1427  labelPair indexPair
1428  (
1429  getEdgeString
1430  (
1431  pp,
1432  globalEdgeFaces,
1433  doneEdge,
1434  patchFacei,
1435  globalFaces.toGlobal(pp.addressing()[patchFacei])
1436  )
1437  );
1438 
1439  // Pout<< "Found unextruded edges in edges:" << fEdges
1440  // << " start:" << indexPair[0]
1441  // << " end:" << indexPair[1]
1442  // << endl;
1443 
1444  const label startFp = indexPair[0];
1445  const label endFp = indexPair[1];
1446 
1447  if (startFp != -1)
1448  {
1449  // Extrude edges from indexPair[0] up to indexPair[1]
1450  // (note indexPair = indices of edges. There is one more vertex
1451  // than edges)
1452  const face& f = localFaces[patchFacei];
1453 
1454  labelList stringedVerts;
1455  if (endFp >= startFp)
1456  {
1457  stringedVerts.setSize(endFp-startFp+2);
1458  }
1459  else
1460  {
1461  stringedVerts.setSize(endFp+f.size()-startFp+2);
1462  }
1463 
1464  label fp = startFp;
1465 
1466  for (label i = 0; i < stringedVerts.size()-1; i++)
1467  {
1468  stringedVerts[i] = f[fp];
1469  doneEdge[fEdges[fp]] = true;
1470  fp = f.fcIndex(fp);
1471  }
1472  stringedVerts.last() = f[fp];
1473 
1474 
1475  // Now stringedVerts contains the vertices in order of face f.
1476  // This is consistent with the order if f becomes the owner cell
1477  // and nbrFacei the neighbour cell. Note that the cells get
1478  // added in order of pp so we can just use face ordering and
1479  // because we loop in incrementing order as well we will
1480  // always have nbrFacei > patchFacei.
1481 
1482  label startEdgeI = fEdges[startFp];
1483 
1484  label meshEdgeI = meshEdges[startEdgeI];
1485 
1486  label numEdgeSideFaces = edgeLayers[startEdgeI];
1487 
1488  for (label i = 0; i < numEdgeSideFaces; i++)
1489  {
1490  label vEnd = stringedVerts.last();
1491  label vStart = stringedVerts[0];
1492 
1493  // calculate number of points making up a face
1494  label newFp = 2*stringedVerts.size();
1495 
1496  if (i == 0)
1497  {
1498  // layer 0 gets all the truncation of neighbouring
1499  // faces with more layers.
1500  if (addedPoints_[vEnd].size())
1501  {
1502  newFp +=
1503  addedPoints_[vEnd].size() - numEdgeSideFaces;
1504  }
1505  if (addedPoints_[vStart].size())
1506  {
1507  newFp +=
1508  addedPoints_[vStart].size() - numEdgeSideFaces;
1509  }
1510  }
1511 
1512  face newFace(newFp);
1513 
1514  newFp = 0;
1515 
1516  // For layer 0 get pp points, for all other layers get
1517  // points of layer-1.
1518  if (i == 0)
1519  {
1520  forAll(stringedVerts, stringedI)
1521  {
1522  label v = stringedVerts[stringedI];
1523  addVertex
1524  (
1525  (
1526  addToMesh_
1527  ? meshPoints[v]
1528  : copiedPatchPoints[v]
1529  ),
1530  newFace,
1531  newFp
1532  );
1533  }
1534  }
1535  else
1536  {
1537  forAll(stringedVerts, stringedI)
1538  {
1539  label v = stringedVerts[stringedI];
1540  if (addedPoints_[v].size())
1541  {
1542  label offset =
1543  addedPoints_[v].size() - numEdgeSideFaces;
1544  addVertex
1545  (
1546  addedPoints_[v][i+offset-1],
1547  newFace,
1548  newFp
1549  );
1550  }
1551  else
1552  {
1553  addVertex
1554  (
1555  (
1556  addToMesh_
1557  ? meshPoints[v]
1558  : copiedPatchPoints[v]
1559  ),
1560  newFace,
1561  newFp
1562  );
1563  }
1564  }
1565  }
1566 
1567  // add points between stringed vertices (end)
1568  if (numEdgeSideFaces < addedPoints_[vEnd].size())
1569  {
1570  if (i == 0 && addedPoints_[vEnd].size())
1571  {
1572  label offset =
1573  addedPoints_[vEnd].size() - numEdgeSideFaces;
1574  for (label ioff = 0; ioff < offset; ioff++)
1575  {
1576  addVertex
1577  (
1578  addedPoints_[vEnd][ioff],
1579  newFace,
1580  newFp
1581  );
1582  }
1583  }
1584  }
1585 
1586  forAllReverse(stringedVerts, stringedI)
1587  {
1588  label v = stringedVerts[stringedI];
1589  if (addedPoints_[v].size())
1590  {
1591  label offset =
1592  addedPoints_[v].size() - numEdgeSideFaces;
1593  addVertex
1594  (
1595  addedPoints_[v][i+offset],
1596  newFace,
1597  newFp
1598  );
1599  }
1600  else
1601  {
1602  addVertex
1603  (
1604  (
1605  addToMesh_
1606  ? meshPoints[v]
1607  : copiedPatchPoints[v]
1608  ),
1609  newFace,
1610  newFp
1611  );
1612  }
1613  }
1614 
1615 
1616  // add points between stringed vertices (start)
1617  if (numEdgeSideFaces < addedPoints_[vStart].size())
1618  {
1619  if (i == 0 && addedPoints_[vStart].size())
1620  {
1621  label offset =
1622  addedPoints_[vStart].size() - numEdgeSideFaces;
1623  for (label ioff = offset-1; ioff >= 0; ioff--)
1624  {
1625  addVertex
1626  (
1627  addedPoints_[vStart][ioff],
1628  newFace,
1629  newFp
1630  );
1631  }
1632  }
1633  }
1634 
1635  if (newFp >= 3)
1636  {
1637  // Add face in between faces patchFacei and nbrFacei
1638  // (possibly -1 for external edges)
1639 
1640  newFace.setSize(newFp);
1641 
1642  if (debug)
1643  {
1644  labelHashSet verts(2*newFace.size());
1645  forAll(newFace, fp)
1646  {
1647  if (!verts.insert(newFace[fp]))
1648  {
1650  << "Duplicate vertex in face"
1651  << " to be added." << nl
1652  << "newFace:" << newFace << nl
1653  << "points:"
1655  (
1656  meshMod.points(),
1657  newFace
1658  ) << nl
1659  << "Layer:" << i
1660  << " out of:" << numEdgeSideFaces << nl
1661  << "ExtrudeEdge:" << meshEdgeI
1662  << " at:"
1663  << mesh_.edges()[meshEdgeI].line
1664  (
1665  mesh_.points()
1666  ) << nl
1667  << "string:" << stringedVerts
1668  << "stringpoints:"
1670  (
1671  pp.localPoints(),
1672  stringedVerts
1673  ) << nl
1674  << "stringNLayers:"
1676  (
1677  nPointLayers,
1678  stringedVerts
1679  ) << nl
1680  << abort(FatalError);
1681  }
1682  }
1683  }
1684 
1685  label nbrFacei = nbrFace
1686  (
1687  pp.edgeFaces(),
1688  startEdgeI,
1689  patchFacei
1690  );
1691 
1692  const labelList& meshFaces = mesh_.edgeFaces
1693  (
1694  meshEdgeI,
1695  ef
1696  );
1697 
1698  addSideFace
1699  (
1700  pp,
1701  addedCells,
1702 
1703  newFace, // vertices of new face
1704  sidePatchID[startEdgeI],// -1 or patch for face
1705 
1706  patchFacei,
1707  nbrFacei,
1708  meshEdgeI, // (mesh) edge to inflate
1709  i, // layer
1710  numEdgeSideFaces, // num layers
1711  meshFaces, // edgeFaces
1712  meshMod
1713  );
1714  }
1715  }
1716  }
1717  }
1718  }
1719 }
1720 
1721 
1724  const mapPolyMesh& morphMap,
1725  const labelList& faceMap, // new to old patch faces
1726  const labelList& pointMap // new to old patch points
1727 )
1728 {
1729  {
1730  labelListList newAddedPoints(pointMap.size());
1731 
1732  forAll(newAddedPoints, newPointi)
1733  {
1734  label oldPointi = pointMap[newPointi];
1735 
1736  const labelList& added = addedPoints_[oldPointi];
1737 
1738  labelList& newAdded = newAddedPoints[newPointi];
1739  newAdded.setSize(added.size());
1740  label newI = 0;
1741 
1742  forAll(added, i)
1743  {
1744  label newPointi = morphMap.reversePointMap()[added[i]];
1745 
1746  if (newPointi >= 0)
1747  {
1748  newAdded[newI++] = newPointi;
1749  }
1750  }
1751  newAdded.setSize(newI);
1752  }
1753  addedPoints_.transfer(newAddedPoints);
1754  }
1755 
1756  {
1757  labelListList newLayerFaces(faceMap.size());
1758 
1759  forAll(newLayerFaces, newFacei)
1760  {
1761  label oldFacei = faceMap[newFacei];
1762 
1763  const labelList& added = layerFaces_[oldFacei];
1764 
1765  labelList& newAdded = newLayerFaces[newFacei];
1766  newAdded.setSize(added.size());
1767  label newI = 0;
1768 
1769  forAll(added, i)
1770  {
1771  label newFacei = morphMap.reverseFaceMap()[added[i]];
1772 
1773  if (newFacei >= 0)
1774  {
1775  newAdded[newI++] = newFacei;
1776  }
1777  }
1778  newAdded.setSize(newI);
1779  }
1780  layerFaces_.transfer(newLayerFaces);
1781  }
1782 }
1783 
1784 
1785 // ************************************************************************* //
label nPatches
Definition: readKivaGrid.H:402
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
bool isLocal(const label i) const
Is on local processor.
Definition: globalIndexI.H:95
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
label nPoints() const
Return number of points supporting patch faces.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:458
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
Class describing modification of a face.
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 > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
const labelListList & pointEdges() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:248
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
label nInternalEdges() const
Number of internal edges.
Class containing data for point addition.
Definition: polyAddPoint.H:47
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:306
label k
Boltzmann constant.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
scalar f1
Definition: createFields.H:28
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
A face addition data class. A face can be inflated either from a point or from another face and can e...
Definition: polyAddFace.H:51
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
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
Class containing data for cell addition.
Definition: polyAddCell.H:46
A list of faces which address into the list of points.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
linePointRef line(const pointField &) const
Return edge line.
Definition: edgeI.H:187
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
const Field< PointType > & points() const
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
label whichProcID(const label i) const
Which processor does global come from? Binary search.
Definition: globalIndexI.H:123
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
label nEdges() const
Return number of edges in patch.
static const char nl
Definition: Ostream.H:260
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
label nEdges() const
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void updateMesh(const mapPolyMesh &, const labelList &faceMap, const labelList &pointMap)
Update any locally stored mesh information. Gets additional.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void setSize(const label)
Reset size of List.
Definition: List.C:281
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
label newPointi
Definition: readKivaGrid.H:501
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
A List with indirect addressing.
Definition: fvMatrix.H:106
Direct mesh changes based on v1.3 polyTopoChange syntax.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
addPatchCellLayer(const polyMesh &, const bool addToMesh=true)
Construct from mesh.
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
static labelListList globalEdgeFaces(const polyMesh &, const globalIndex &globalFaces, const indirectPrimitivePatch &pp)
Per patch edge the pp faces (in global indices) using it. Uses.
label n
labelListList addedCells() const
Added cells given current mesh & layerfaces.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
T & last()
Return the last element of the list.
Definition: UListI.H:128
void setRefinement(const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const scalarField &expansionRatio, const indirectPrimitivePatch &pp, const labelList &sidePatchID, const labelList &exposedPatchID, const labelList &nFaceLayers, const labelList &nPointLayers, const vectorField &firstLayerDisp, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layers on top.
static void calcSidePatch(const polyMesh &, const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const indirectPrimitivePatch &pp, labelList &sidePatchID, label &nPatches, Map< label > &nbrProcToPatch, Map< label > &patchToNbrProc)
Boundary edges get extruded into boundary faces. Determine patch.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
const labelListList & edgeFaces() const
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:490
Namespace for OpenFOAM.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.