addPatchCellLayer.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "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
473 Foam::addPatchCellLayer::addPatchCellLayer
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];
708  //FatalErrorIn("addPatchCellLayer::calcSidePatch(..)")
709  WarningIn("addPatchCellLayer::calcSidePatch(..)")
710  << "Have no sidePatchID for edge " << edgeI << " points "
711  << pp.points()[pp.meshPoints()[e[0]]]
712  << pp.points()[pp.meshPoints()[e[1]]]
713  //<< abort(FatalError);
714  << endl;
715  }
716  }
717  }
718 
719 
720 
721  // Now we have sidepatch see if we have patchface or edge to inflate
722  // from.
723  forAll(edgeFaces, edgeI)
724  {
725  if
726  (
727  edgeFaces[edgeI].size() == 1
728  && sidePatchID[edgeI] != -1
729  && inflateFaceI[edgeI] == -1
730  )
731  {
732  // 1. Do we have a boundary face to inflate from
733 
734  label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
735 
736  // Pick up any boundary face on this edge and use its properties
737  label meshEdgeI = meshEdges[edgeI];
738  const labelList& meshFaces = mesh.edgeFaces
739  (
740  meshEdgeI,
741  dynMeshEdgeFaces
742  );
743 
744  forAll(meshFaces, k)
745  {
746  label faceI = meshFaces[k];
747 
748  if (faceI != myFaceI)
749  {
750  if (mesh.isInternalFace(faceI))
751  {
752  inflateEdgeI[edgeI] = meshEdgeI;
753  }
754  else
755  {
756  if (patches.whichPatch(faceI) == sidePatchID[edgeI])
757  {
758  setFaceProps
759  (
760  mesh,
761  faceI,
762 
763  sidePatchID[edgeI],
764  sideZoneID[edgeI],
765  sideFlip[edgeI]
766  );
767  inflateFaceI[edgeI] = faceI;
768  inflateEdgeI[edgeI] = -1;
769 
770  break;
771  }
772  }
773  }
774  }
775  }
776  }
777 }
778 
779 
781 (
782  const globalIndex& globalFaces,
783  const labelListList& globalEdgeFaces,
784  const scalarField& expansionRatio,
785  const indirectPrimitivePatch& pp,
786  const labelList& sidePatchID,
787  const labelList& exposedPatchID,
788  const labelList& nFaceLayers,
789  const labelList& nPointLayers,
790  const vectorField& firstLayerDisp,
791  polyTopoChange& meshMod
792 )
793 {
794  if (debug)
795  {
796  Pout<< "addPatchCellLayer::setRefinement : Adding up to "
797  << gMax(nPointLayers)
798  << " layers of cells to indirectPrimitivePatch with "
799  << pp.nPoints() << " points" << endl;
800  }
801 
802  if
803  (
804  pp.nPoints() != firstLayerDisp.size()
805  || pp.nPoints() != nPointLayers.size()
806  || pp.size() != nFaceLayers.size()
807  )
808  {
810  (
811  "addPatchCellLayer::setRefinement"
812  "(const scalar, const indirectPrimitivePatch&"
813  ", const labelList&, const vectorField&, polyTopoChange&)"
814  ) << "Size of new points is not same as number of points used by"
815  << " the face subset" << endl
816  << " patch.nPoints:" << pp.nPoints()
817  << " displacement:" << firstLayerDisp.size()
818  << " nPointLayers:" << nPointLayers.size() << nl
819  << " patch.nFaces:" << pp.size()
820  << " nFaceLayers:" << nFaceLayers.size()
821  << abort(FatalError);
822  }
823 
824  forAll(nPointLayers, i)
825  {
826  if (nPointLayers[i] < 0)
827  {
829  (
830  "addPatchCellLayer::setRefinement"
831  "(const scalar, const indirectPrimitivePatch&"
832  ", const labelList&, const vectorField&, polyTopoChange&)"
833  ) << "Illegal number of layers " << nPointLayers[i]
834  << " at patch point " << i << abort(FatalError);
835  }
836  }
837  forAll(nFaceLayers, i)
838  {
839  if (nFaceLayers[i] < 0)
840  {
842  (
843  "addPatchCellLayer::setRefinement"
844  "(const scalar, const indirectPrimitivePatch&"
845  ", const labelList&, const vectorField&, polyTopoChange&)"
846  ) << "Illegal number of layers " << nFaceLayers[i]
847  << " at patch face " << i << abort(FatalError);
848  }
849  }
850 
851  forAll(globalEdgeFaces, edgeI)
852  {
853  if (globalEdgeFaces[edgeI].size() > 2)
854  {
855  const edge& e = pp.edges()[edgeI];
856 
857  if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
858  {
860  (
861  "addPatchCellLayer::setRefinement"
862  "(const scalar, const indirectPrimitivePatch&"
863  ", const labelList&, const vectorField&, polyTopoChange&)"
864  ) << "Trying to extrude edge "
865  << e.line(pp.localPoints())
866  << " which is non-manifold (has "
867  << globalEdgeFaces[edgeI].size()
868  << " faces using it)"
869  << abort(FatalError);
870  }
871  }
872  }
873 
874 
875  const labelList& meshPoints = pp.meshPoints();
876 
877  // Some storage for edge-face-addressing.
879 
880  // Precalculate mesh edges for pp.edges.
881  const labelList meshEdges(pp.meshEdges(mesh_.edges(), mesh_.pointEdges()));
882 
883  if (debug)
884  {
885  // Check synchronisation
886  // ~~~~~~~~~~~~~~~~~~~~~
887 
888  {
889  labelList n(mesh_.nPoints(), 0);
890  UIndirectList<label>(n, meshPoints) = nPointLayers;
892 
893  // Non-synced
894  forAll(meshPoints, i)
895  {
896  label meshPointI = meshPoints[i];
897 
898  if (n[meshPointI] != nPointLayers[i])
899  {
901  (
902  "addPatchCellLayer::setRefinement"
903  "(const scalar, const indirectPrimitivePatch&"
904  ", const labelList&, const vectorField&"
905  ", polyTopoChange&)"
906  ) << "At mesh point:" << meshPointI
907  << " coordinate:" << mesh_.points()[meshPointI]
908  << " specified nLayers:" << nPointLayers[i] << endl
909  << "On coupled point a different nLayers:"
910  << n[meshPointI] << " was specified."
911  << abort(FatalError);
912  }
913  }
914 
915 
916  // Check that nPointLayers equals the max layers of connected faces
917  // (or 0). Anything else makes no sense.
918  labelList nFromFace(mesh_.nPoints(), 0);
919  forAll(nFaceLayers, i)
920  {
921  const face& f = pp[i];
922 
923  forAll(f, fp)
924  {
925  label pointI = f[fp];
926 
927  nFromFace[pointI] = max(nFromFace[pointI], nFaceLayers[i]);
928  }
929  }
931  (
932  mesh_,
933  nFromFace,
934  maxEqOp<label>(),
935  label(0)
936  );
937 
938  forAll(nPointLayers, i)
939  {
940  label meshPointI = meshPoints[i];
941 
942  if
943  (
944  nPointLayers[i] > 0
945  && nPointLayers[i] != nFromFace[meshPointI]
946  )
947  {
949  (
950  "addPatchCellLayer::setRefinement"
951  "(const scalar, const indirectPrimitivePatch&"
952  ", const labelList&, const vectorField&"
953  ", polyTopoChange&)"
954  ) << "At mesh point:" << meshPointI
955  << " coordinate:" << mesh_.points()[meshPointI]
956  << " specified nLayers:" << nPointLayers[i] << endl
957  << "but the max nLayers of surrounding faces is:"
958  << nFromFace[meshPointI]
959  << abort(FatalError);
960  }
961  }
962  }
963 
964  {
965  pointField d(mesh_.nPoints(), vector::max);
966  UIndirectList<point>(d, meshPoints) = firstLayerDisp;
968  (
969  mesh_,
970  d,
971  minEqOp<vector>(),
973  );
974 
975  forAll(meshPoints, i)
976  {
977  label meshPointI = meshPoints[i];
978 
979  if (mag(d[meshPointI] - firstLayerDisp[i]) > SMALL)
980  {
982  (
983  "addPatchCellLayer::setRefinement"
984  "(const scalar, const indirectPrimitivePatch&"
985  ", const labelList&, const vectorField&"
986  ", polyTopoChange&)"
987  ) << "At mesh point:" << meshPointI
988  << " coordinate:" << mesh_.points()[meshPointI]
989  << " specified displacement:" << firstLayerDisp[i]
990  << endl
991  << "On coupled point a different displacement:"
992  << d[meshPointI] << " was specified."
993  << abort(FatalError);
994  }
995  }
996  }
997 
998  // Check that edges of pp (so ones that become boundary faces)
999  // connect to only one boundary face. Guarantees uniqueness of
1000  // patch that they go into so if this is a coupled patch both
1001  // sides decide the same.
1002  // ~~~~~~~~~~~~~~~~~~~~~~
1003 
1004  for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
1005  {
1006  const edge& e = pp.edges()[edgeI];
1007 
1008  if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
1009  {
1010  // Edge is to become a face
1011 
1012  const labelList& eFaces = pp.edgeFaces()[edgeI];
1013 
1014  // First check: pp should be single connected.
1015  if (eFaces.size() != 1)
1016  {
1017  FatalErrorIn
1018  (
1019  "addPatchCellLayer::setRefinement"
1020  "(const scalar, const indirectPrimitivePatch&"
1021  ", const labelList&, const vectorField&"
1022  ", polyTopoChange&)"
1023  ) << "boundary-edge-to-be-extruded:"
1024  << pp.points()[meshPoints[e[0]]]
1025  << pp.points()[meshPoints[e[1]]]
1026  << " has more than two faces using it:" << eFaces
1027  << abort(FatalError);
1028  }
1029 
1030  label myFaceI = pp.addressing()[eFaces[0]];
1031 
1032  label meshEdgeI = meshEdges[edgeI];
1033 
1034  // Mesh faces using edge
1035  const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
1036 
1037  // Check that there is only one patchface using edge.
1038  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1039 
1040  label bFaceI = -1;
1041 
1042  forAll(meshFaces, i)
1043  {
1044  label faceI = meshFaces[i];
1045 
1046  if (faceI != myFaceI)
1047  {
1048  if (!mesh_.isInternalFace(faceI))
1049  {
1050  if (bFaceI == -1)
1051  {
1052  bFaceI = faceI;
1053  }
1054  else
1055  {
1056  FatalErrorIn
1057  (
1058  "addPatchCellLayer::setRefinement"
1059  "(const scalar"
1060  ", const indirectPrimitivePatch&"
1061  ", const labelList&, const vectorField&"
1062  ", polyTopoChange&)"
1063  ) << "boundary-edge-to-be-extruded:"
1064  << pp.points()[meshPoints[e[0]]]
1065  << pp.points()[meshPoints[e[1]]]
1066  << " has more than two boundary faces"
1067  << " using it:"
1068  << bFaceI << " fc:"
1069  << mesh_.faceCentres()[bFaceI]
1070  << " patch:" << patches.whichPatch(bFaceI)
1071  << " and " << faceI << " fc:"
1072  << mesh_.faceCentres()[faceI]
1073  << " patch:" << patches.whichPatch(faceI)
1074  << abort(FatalError);
1075  }
1076  }
1077  }
1078  }
1079  }
1080  }
1081  }
1082 
1083 
1084  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1085 
1086  // Precalculated patchID for each patch face
1087  labelList patchID(pp.size());
1088 
1089  forAll(pp, patchFaceI)
1090  {
1091  label meshFaceI = pp.addressing()[patchFaceI];
1092 
1093  patchID[patchFaceI] = patches.whichPatch(meshFaceI);
1094  }
1095 
1096 
1097  // From master point (in patch point label) to added points (in mesh point
1098  // label)
1099  addedPoints_.setSize(pp.nPoints());
1100 
1101  // Mark points that do not get extruded by setting size of addedPoints_ to 0
1102  label nTruncated = 0;
1103 
1104  forAll(nPointLayers, patchPointI)
1105  {
1106  if (nPointLayers[patchPointI] > 0)
1107  {
1108  addedPoints_[patchPointI].setSize(nPointLayers[patchPointI]);
1109  }
1110  else
1111  {
1112  nTruncated++;
1113  }
1114  }
1115 
1116  if (debug)
1117  {
1118  Pout<< "Not adding points at " << nTruncated << " out of "
1119  << pp.nPoints() << " points" << endl;
1120  }
1121 
1122 
1123  //
1124  // Create new points
1125  //
1126 
1127  // If creating new mesh: copy existing patch points
1128  labelList copiedPatchPoints;
1129  if (!addToMesh_)
1130  {
1131  copiedPatchPoints.setSize(firstLayerDisp.size());
1132  forAll(firstLayerDisp, patchPointI)
1133  {
1134  if (addedPoints_[patchPointI].size())
1135  {
1136  label meshPointI = meshPoints[patchPointI];
1137  label zoneI = mesh_.pointZones().whichZone(meshPointI);
1138  copiedPatchPoints[patchPointI] = meshMod.setAction
1139  (
1140  polyAddPoint
1141  (
1142  mesh_.points()[meshPointI], // point
1143  -1, // master point
1144  zoneI, // zone for point
1145  true // supports a cell
1146  )
1147  );
1148  }
1149  }
1150  }
1151 
1152 
1153  // Create points for additional layers
1154  forAll(firstLayerDisp, patchPointI)
1155  {
1156  if (addedPoints_[patchPointI].size())
1157  {
1158  label meshPointI = meshPoints[patchPointI];
1159 
1160  label zoneI = mesh_.pointZones().whichZone(meshPointI);
1161 
1162  point pt = mesh_.points()[meshPointI];
1163 
1164  vector disp = firstLayerDisp[patchPointI];
1165 
1166  forAll(addedPoints_[patchPointI], i)
1167  {
1168  pt += disp;
1169 
1170  label addedVertI = meshMod.setAction
1171  (
1172  polyAddPoint
1173  (
1174  pt, // point
1175  (addToMesh_ ? meshPointI : -1), // master point
1176  zoneI, // zone for point
1177  true // supports a cell
1178  )
1179  );
1180 
1181  addedPoints_[patchPointI][i] = addedVertI;
1182 
1183  disp *= expansionRatio[patchPointI];
1184  }
1185  }
1186  }
1187 
1188 
1189  //
1190  // Add cells to all boundaryFaces
1191  //
1192 
1193  labelListList addedCells(pp.size());
1194 
1195  forAll(pp, patchFaceI)
1196  {
1197  if (nFaceLayers[patchFaceI] > 0)
1198  {
1199  addedCells[patchFaceI].setSize(nFaceLayers[patchFaceI]);
1200 
1201  label meshFaceI = pp.addressing()[patchFaceI];
1202 
1203  label ownZoneI = mesh_.cellZones().whichZone
1204  (
1205  mesh_.faceOwner()[meshFaceI]
1206  );
1207 
1208  for (label i = 0; i < nFaceLayers[patchFaceI]; i++)
1209  {
1210  // Note: add from cell (owner of patch face) or from face?
1211  // for now add from cell so we can map easily.
1212  addedCells[patchFaceI][i] = meshMod.setAction
1213  (
1214  polyAddCell
1215  (
1216  -1, // master point
1217  -1, // master edge
1218  -1, // master face
1219  (addToMesh_ ? mesh_.faceOwner()[meshFaceI] : -1),
1220  //master
1221  ownZoneI // zone for cell
1222  )
1223  );
1224  }
1225  }
1226  }
1227 
1228 
1229 
1230  // Create faces on top of the original patch faces.
1231  // These faces are created from original patch faces outwards so in order
1232  // of increasing cell number. So orientation should be same as original
1233  // patch face for them to have owner<neighbour.
1234 
1235  layerFaces_.setSize(pp.size());
1236 
1237  forAll(pp.localFaces(), patchFaceI)
1238  {
1239  label meshFaceI = pp.addressing()[patchFaceI];
1240 
1241  if (addedCells[patchFaceI].size())
1242  {
1243  layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1);
1244 
1245  // Get duplicated vertices on the patch face.
1246  const face& f = pp.localFaces()[patchFaceI];
1247 
1248  face newFace(f.size());
1249 
1250  forAll(addedCells[patchFaceI], i)
1251  {
1252  forAll(f, fp)
1253  {
1254  if (addedPoints_[f[fp]].empty())
1255  {
1256  // Keep original point
1257  newFace[fp] =
1258  (
1259  addToMesh_
1260  ? meshPoints[f[fp]]
1261  : copiedPatchPoints[f[fp]]
1262  );
1263  }
1264  else
1265  {
1266  // Get new outside point
1267  label offset =
1268  addedPoints_[f[fp]].size()
1269  - addedCells[patchFaceI].size();
1270  newFace[fp] = addedPoints_[f[fp]][i+offset];
1271  }
1272  }
1273 
1274 
1275  // Get new neighbour
1276  label nei;
1277  label patchI;
1278  label zoneI = -1;
1279  bool flip = false;
1280 
1281 
1282  if (i == addedCells[patchFaceI].size()-1)
1283  {
1284  // Top layer so is patch face.
1285  nei = -1;
1286  patchI = patchID[patchFaceI];
1287  zoneI = mesh_.faceZones().whichZone(meshFaceI);
1288  if (zoneI != -1)
1289  {
1290  const faceZone& fz = mesh_.faceZones()[zoneI];
1291  flip = fz.flipMap()[fz.whichFace(meshFaceI)];
1292  }
1293  }
1294  else
1295  {
1296  // Internal face between layer i and i+1
1297  nei = addedCells[patchFaceI][i+1];
1298  patchI = -1;
1299  }
1300 
1301 
1302  layerFaces_[patchFaceI][i+1] = meshMod.setAction
1303  (
1304  polyAddFace
1305  (
1306  newFace, // face
1307  addedCells[patchFaceI][i], // owner
1308  nei, // neighbour
1309  -1, // master point
1310  -1, // master edge
1311  (addToMesh_ ? meshFaceI : -1), // master face
1312  false, // flux flip
1313  patchI, // patch for face
1314  zoneI, // zone for face
1315  flip // face zone flip
1316  )
1317  );
1318  }
1319  }
1320  }
1321 
1322  //
1323  // Modify old patch faces to be on the inside
1324  //
1325 
1326  if (addToMesh_)
1327  {
1328  forAll(pp, patchFaceI)
1329  {
1330  if (addedCells[patchFaceI].size())
1331  {
1332  label meshFaceI = pp.addressing()[patchFaceI];
1333 
1334  layerFaces_[patchFaceI][0] = meshFaceI;
1335 
1336  meshMod.setAction
1337  (
1339  (
1340  pp[patchFaceI], // modified face
1341  meshFaceI, // label of face
1342  mesh_.faceOwner()[meshFaceI], // owner
1343  addedCells[patchFaceI][0], // neighbour
1344  false, // face flip
1345  -1, // patch for face
1346  true, //false, // remove from zone
1347  -1, //zoneI, // zone for face
1348  false // face flip in zone
1349  )
1350  );
1351  }
1352  }
1353  }
1354  else
1355  {
1356  // If creating new mesh: reverse original faces and put them
1357  // in the exposed patch ID.
1358  forAll(pp, patchFaceI)
1359  {
1360  if (nFaceLayers[patchFaceI] > 0)
1361  {
1362  label meshFaceI = pp.addressing()[patchFaceI];
1363  label zoneI = mesh_.faceZones().whichZone(meshFaceI);
1364  bool zoneFlip = false;
1365  if (zoneI != -1)
1366  {
1367  const faceZone& fz = mesh_.faceZones()[zoneI];
1368  zoneFlip = !fz.flipMap()[fz.whichFace(meshFaceI)];
1369  }
1370 
1371  // Reverse and renumber old patch face.
1372  face f(pp.localFaces()[patchFaceI].reverseFace());
1373  forAll(f, fp)
1374  {
1375  f[fp] = copiedPatchPoints[f[fp]];
1376  }
1377 
1378  layerFaces_[patchFaceI][0] = meshMod.setAction
1379  (
1380  polyAddFace
1381  (
1382  f, // modified face
1383  addedCells[patchFaceI][0], // owner
1384  -1, // neighbour
1385  -1, // masterPoint
1386  -1, // masterEdge
1387  -1, // masterFace
1388  true, // face flip
1389  exposedPatchID[patchFaceI], // patch for face
1390  zoneI, // zone for face
1391  zoneFlip // face flip in zone
1392  )
1393  );
1394  }
1395  }
1396  }
1397 
1398 
1399 
1400  //
1401  // Create 'side' faces, one per edge that is being extended.
1402  //
1403 
1404  const labelListList& faceEdges = pp.faceEdges();
1405  const faceList& localFaces = pp.localFaces();
1406  const edgeList& edges = pp.edges();
1407 
1408  // Get number of layers per edge. This is 0 if edge is not extruded;
1409  // max of connected faces otherwise.
1410  labelList edgeLayers(pp.nEdges());
1411 
1412  {
1413  // Use list over mesh.nEdges() since syncTools does not yet support
1414  // partial list synchronisation.
1415  labelList meshEdgeLayers(mesh_.nEdges(), -1);
1416 
1417  forAll(meshEdges, edgeI)
1418  {
1419  const edge& e = edges[edgeI];
1420 
1421  label meshEdgeI = meshEdges[edgeI];
1422 
1423  if ((nPointLayers[e[0]] == 0) && (nPointLayers[e[1]] == 0))
1424  {
1425  meshEdgeLayers[meshEdgeI] = 0;
1426  }
1427  else
1428  {
1429  const labelList& eFaces = pp.edgeFaces()[edgeI];
1430 
1431  forAll(eFaces, i)
1432  {
1433  meshEdgeLayers[meshEdgeI] = max
1434  (
1435  nFaceLayers[eFaces[i]],
1436  meshEdgeLayers[meshEdgeI]
1437  );
1438  }
1439  }
1440  }
1441 
1443  (
1444  mesh_,
1445  meshEdgeLayers,
1446  maxEqOp<label>(),
1447  label(0) // initial value
1448  );
1449 
1450  forAll(meshEdges, edgeI)
1451  {
1452  edgeLayers[edgeI] = meshEdgeLayers[meshEdges[edgeI]];
1453  }
1454  }
1455 
1456 
1457  // Mark off which edges have been extruded
1458  boolList doneEdge(pp.nEdges(), false);
1459 
1460 
1461  // Create faces. Per face walk connected edges and find string of edges
1462  // between the same two faces and extrude string into a single face.
1463  forAll(pp, patchFaceI)
1464  {
1465  const labelList& fEdges = faceEdges[patchFaceI];
1466 
1467  forAll(fEdges, fp)
1468  {
1469  // Get string of edges that needs to be extruded as a single face.
1470  // Returned as indices in fEdges.
1471  labelPair indexPair
1472  (
1473  getEdgeString
1474  (
1475  pp,
1476  globalEdgeFaces,
1477  doneEdge,
1478  patchFaceI,
1479  globalFaces.toGlobal(pp.addressing()[patchFaceI])
1480  )
1481  );
1482 
1483  //Pout<< "Found unextruded edges in edges:" << fEdges
1484  // << " start:" << indexPair[0]
1485  // << " end:" << indexPair[1]
1486  // << endl;
1487 
1488  const label startFp = indexPair[0];
1489  const label endFp = indexPair[1];
1490 
1491  if (startFp != -1)
1492  {
1493  // Extrude edges from indexPair[0] up to indexPair[1]
1494  // (note indexPair = indices of edges. There is one more vertex
1495  // than edges)
1496  const face& f = localFaces[patchFaceI];
1497 
1498  labelList stringedVerts;
1499  if (endFp >= startFp)
1500  {
1501  stringedVerts.setSize(endFp-startFp+2);
1502  }
1503  else
1504  {
1505  stringedVerts.setSize(endFp+f.size()-startFp+2);
1506  }
1507 
1508  label fp = startFp;
1509 
1510  for (label i = 0; i < stringedVerts.size()-1; i++)
1511  {
1512  stringedVerts[i] = f[fp];
1513  doneEdge[fEdges[fp]] = true;
1514  fp = f.fcIndex(fp);
1515  }
1516  stringedVerts.last() = f[fp];
1517 
1518 
1519  // Now stringedVerts contains the vertices in order of face f.
1520  // This is consistent with the order if f becomes the owner cell
1521  // and nbrFaceI the neighbour cell. Note that the cells get
1522  // added in order of pp so we can just use face ordering and
1523  // because we loop in incrementing order as well we will
1524  // always have nbrFaceI > patchFaceI.
1525 
1526  label startEdgeI = fEdges[startFp];
1527 
1528  label meshEdgeI = meshEdges[startEdgeI];
1529 
1530  label numEdgeSideFaces = edgeLayers[startEdgeI];
1531 
1532  for (label i = 0; i < numEdgeSideFaces; i++)
1533  {
1534  label vEnd = stringedVerts.last();
1535  label vStart = stringedVerts[0];
1536 
1537  // calculate number of points making up a face
1538  label newFp = 2*stringedVerts.size();
1539 
1540  if (i == 0)
1541  {
1542  // layer 0 gets all the truncation of neighbouring
1543  // faces with more layers.
1544  if (addedPoints_[vEnd].size())
1545  {
1546  newFp +=
1547  addedPoints_[vEnd].size() - numEdgeSideFaces;
1548  }
1549  if (addedPoints_[vStart].size())
1550  {
1551  newFp +=
1552  addedPoints_[vStart].size() - numEdgeSideFaces;
1553  }
1554  }
1555 
1556  face newFace(newFp);
1557 
1558  newFp = 0;
1559 
1560  // For layer 0 get pp points, for all other layers get
1561  // points of layer-1.
1562  if (i == 0)
1563  {
1564  forAll(stringedVerts, stringedI)
1565  {
1566  label v = stringedVerts[stringedI];
1567  addVertex
1568  (
1569  (
1570  addToMesh_
1571  ? meshPoints[v]
1572  : copiedPatchPoints[v]
1573  ),
1574  newFace,
1575  newFp
1576  );
1577  }
1578  }
1579  else
1580  {
1581  forAll(stringedVerts, stringedI)
1582  {
1583  label v = stringedVerts[stringedI];
1584  if (addedPoints_[v].size())
1585  {
1586  label offset =
1587  addedPoints_[v].size() - numEdgeSideFaces;
1588  addVertex
1589  (
1590  addedPoints_[v][i+offset-1],
1591  newFace,
1592  newFp
1593  );
1594  }
1595  else
1596  {
1597  addVertex
1598  (
1599  (
1600  addToMesh_
1601  ? meshPoints[v]
1602  : copiedPatchPoints[v]
1603  ),
1604  newFace,
1605  newFp
1606  );
1607  }
1608  }
1609  }
1610 
1611  // add points between stringed vertices (end)
1612  if (numEdgeSideFaces < addedPoints_[vEnd].size())
1613  {
1614  if (i == 0 && addedPoints_[vEnd].size())
1615  {
1616  label offset =
1617  addedPoints_[vEnd].size() - numEdgeSideFaces;
1618  for (label ioff = 0; ioff < offset; ioff++)
1619  {
1620  addVertex
1621  (
1622  addedPoints_[vEnd][ioff],
1623  newFace,
1624  newFp
1625  );
1626  }
1627  }
1628  }
1629 
1630  forAllReverse(stringedVerts, stringedI)
1631  {
1632  label v = stringedVerts[stringedI];
1633  if (addedPoints_[v].size())
1634  {
1635  label offset =
1636  addedPoints_[v].size() - numEdgeSideFaces;
1637  addVertex
1638  (
1639  addedPoints_[v][i+offset],
1640  newFace,
1641  newFp
1642  );
1643  }
1644  else
1645  {
1646  addVertex
1647  (
1648  (
1649  addToMesh_
1650  ? meshPoints[v]
1651  : copiedPatchPoints[v]
1652  ),
1653  newFace,
1654  newFp
1655  );
1656  }
1657  }
1658 
1659 
1660  // add points between stringed vertices (start)
1661  if (numEdgeSideFaces < addedPoints_[vStart].size())
1662  {
1663  if (i == 0 && addedPoints_[vStart].size())
1664  {
1665  label offset =
1666  addedPoints_[vStart].size() - numEdgeSideFaces;
1667  for (label ioff = offset-1; ioff >= 0; ioff--)
1668  {
1669  addVertex
1670  (
1671  addedPoints_[vStart][ioff],
1672  newFace,
1673  newFp
1674  );
1675  }
1676  }
1677  }
1678 
1679  if (newFp >= 3)
1680  {
1681  // Add face inbetween faces patchFaceI and nbrFaceI
1682  // (possibly -1 for external edges)
1683 
1684  newFace.setSize(newFp);
1685 
1686  if (debug)
1687  {
1688  labelHashSet verts(2*newFace.size());
1689  forAll(newFace, fp)
1690  {
1691  if (!verts.insert(newFace[fp]))
1692  {
1693  FatalErrorIn
1694  (
1695  "addPatchCellLayer::setRefinement(..)"
1696  ) << "Duplicate vertex in face"
1697  << " to be added." << nl
1698  << "newFace:" << newFace << nl
1699  << "points:"
1701  (
1702  meshMod.points(),
1703  newFace
1704  ) << nl
1705  << "Layer:" << i
1706  << " out of:" << numEdgeSideFaces << nl
1707  << "ExtrudeEdge:" << meshEdgeI
1708  << " at:"
1709  << mesh_.edges()[meshEdgeI].line
1710  (
1711  mesh_.points()
1712  ) << nl
1713  << "string:" << stringedVerts
1714  << "stringpoints:"
1716  (
1717  pp.localPoints(),
1718  stringedVerts
1719  ) << nl
1720  << "stringNLayers:"
1722  (
1723  nPointLayers,
1724  stringedVerts
1725  ) << nl
1726  << abort(FatalError);
1727  }
1728  }
1729  }
1730 
1731  label nbrFaceI = nbrFace
1732  (
1733  pp.edgeFaces(),
1734  startEdgeI,
1735  patchFaceI
1736  );
1737 
1738  const labelList& meshFaces = mesh_.edgeFaces
1739  (
1740  meshEdgeI,
1741  ef
1742  );
1743 
1744  addSideFace
1745  (
1746  pp,
1747  addedCells,
1748 
1749  newFace, // vertices of new face
1750  sidePatchID[startEdgeI],// -1 or patch for face
1751 
1752  patchFaceI,
1753  nbrFaceI,
1754  meshEdgeI, // (mesh) edge to inflate
1755  i, // layer
1756  numEdgeSideFaces, // num layers
1757  meshFaces, // edgeFaces
1758  meshMod
1759  );
1760  }
1761  }
1762  }
1763  }
1764  }
1765 }
1766 
1767 
1770  const mapPolyMesh& morphMap,
1771  const labelList& faceMap, // new to old patch faces
1772  const labelList& pointMap // new to old patch points
1773 )
1774 {
1775  {
1776  labelListList newAddedPoints(pointMap.size());
1777 
1778  forAll(newAddedPoints, newPointI)
1779  {
1780  label oldPointI = pointMap[newPointI];
1781 
1782  const labelList& added = addedPoints_[oldPointI];
1783 
1784  labelList& newAdded = newAddedPoints[newPointI];
1785  newAdded.setSize(added.size());
1786  label newI = 0;
1787 
1788  forAll(added, i)
1789  {
1790  label newPointI = morphMap.reversePointMap()[added[i]];
1791 
1792  if (newPointI >= 0)
1793  {
1794  newAdded[newI++] = newPointI;
1795  }
1796  }
1797  newAdded.setSize(newI);
1798  }
1799  addedPoints_.transfer(newAddedPoints);
1800  }
1801 
1802  {
1803  labelListList newLayerFaces(faceMap.size());
1804 
1805  forAll(newLayerFaces, newFaceI)
1806  {
1807  label oldFaceI = faceMap[newFaceI];
1808 
1809  const labelList& added = layerFaces_[oldFaceI];
1810 
1811  labelList& newAdded = newLayerFaces[newFaceI];
1812  newAdded.setSize(added.size());
1813  label newI = 0;
1814 
1815  forAll(added, i)
1816  {
1817  label newFaceI = morphMap.reverseFaceMap()[added[i]];
1818 
1819  if (newFaceI >= 0)
1820  {
1821  newAdded[newI++] = newFaceI;
1822  }
1823  }
1824  newAdded.setSize(newI);
1825  }
1826  layerFaces_.transfer(newLayerFaces);
1827  }
1828 }
1829 
1830 
1831 // ************************************************************************* //
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
bool isLocal(const label i) const
Is on local processor.
Definition: globalIndexI.H:95
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
label size() const
Return the number of elements in the PtrList.
Definition: PtrListI.H:32
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:497
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
dimensioned< scalar > mag(const dimensioned< Type > &)
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & pointEdges() const
labelList f(nPoints)
labelListList addedCells() const
Added cells given current mesh & layerfaces.
const labelListList & edgeFaces() const
Return edge-face addressing.
static labelListList globalEdgeFaces(const polyMesh &, const globalIndex &globalFaces, const indirectPrimitivePatch &pp)
Per patch edge the pp faces (in global indices) using it. Uses.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:465
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
T & last()
Return the last element of the list.
Definition: UListI.H:131
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
scalar f1
Definition: createFields.H:28
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
static const Vector max
Definition: Vector.H:82
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
label nEdges() const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
label nPoints() const
Return number of points supporting patch faces.
Namespace for OpenFOAM.
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:317
Class containing data for cell addition.
Definition: polyAddCell.H:46
label n
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Field< PointType > & localPoints() const
Return pointField of points in patch.
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
A List with indirect addressing.
Definition: fvMatrix.H:106
Foam::polyBoundaryMesh.
linePointRef line(const pointField &) const
Return edge line.
Definition: edgeI.H:169
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
#define forAll(list, i)
Definition: UList.H:421
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
label k
Boltzmann constant.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Class describing modification of a face.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A list of faces which address into the list of points.
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.
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.
#define forAllReverse(list, i)
Definition: UList.H:424
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
label nInternalEdges() const
Number of internal edges.
error FatalError
label nPatches
Definition: readKivaGrid.H:402
const labelListList & edgeFaces() const
List< label > labelList
A List of labels.
Definition: labelList.H:56
const labelListList & faceEdges() const
Return face-edge addressing.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
Direct mesh changes based on v1.3 polyTopoChange syntax.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
Class containing data for point addition.
Definition: polyAddPoint.H:47
label nEdges() const
Return number of edges in patch.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
void updateMesh(const mapPolyMesh &, const labelList &faceMap, const labelList &pointMap)
Update any locally stored mesh information. Gets additional.
Type gMax(const FieldField< Field, Type > &f)
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
const Field< PointType > & points() const
Return reference to global points.
label whichProcID(const label i) const
Which processor does global come from? Binary search.
Definition: globalIndexI.H:123
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50