FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/polyMesh.H>
28 #include "polyTopoChange.H"
29 #include <meshTools/meshTools.H>
30 #include <OpenFOAM/mapPolyMesh.H>
31 #include <OpenFOAM/syncTools.H>
36 #include <meshTools/wallPoint.H>
37 #include <OpenFOAM/globalIndex.H>
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 // Calculate global faces per pp edge.
47 Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
48 (
49  const polyMesh& mesh,
50  const globalIndex& globalFaces,
51  const indirectPrimitivePatch& pp,
52  const labelList& meshEdges
53 )
54 {
57  //
58  //PackedBoolList isCoupledEdge(mesh.nEdges());
59  //
60  //const polyBoundaryMesh& patches = mesh.boundaryMesh();
61  //
62  //forAll(patches, patchI)
63  //{
64  // const polyPatch& pp = patches[patchI];
65  //
66  // if (pp.coupled())
67  // {
68  // const labelList& meshEdges = pp.meshEdges();
69  //
70  // forAll(meshEdges, i)
71  // {
72  // isCoupledEdge.set(meshEdges[i], 1);
73  // }
74  // }
75  //}
76 
77  // From mesh edge to global face labels. Sized only for pp edges.
78  labelListList globalEdgeFaces(mesh.nEdges());
79 
80  const labelListList& edgeFaces = pp.edgeFaces();
81 
82  forAll(edgeFaces, edgeI)
83  {
84  label meshEdgeI = meshEdges[edgeI];
85 
86  //if (isCoupledEdge.get(meshEdgeI) == 1)
87  {
88  const labelList& eFaces = edgeFaces[edgeI];
89 
90  // Store face and processor as unique tag.
91  labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
92  globalEFaces.setSize(eFaces.size());
93  forAll(eFaces, i)
94  {
95  globalEFaces[i] =
96  globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
97  }
98  }
99  }
100 
101  // Synchronise across coupled edges.
103  (
104  mesh,
105  globalEdgeFaces,
106  uniqueEqOp(),
107  labelList(), // null value
108  false // no separation
109  );
110 
111  // Extract pp part
112  return UIndirectList<labelList>(globalEdgeFaces, meshEdges);
113 }
114 
115 
116 Foam::label Foam::addPatchCellLayer::nbrFace
117 (
118  const labelListList& edgeFaces,
119  const label edgeI,
120  const label faceI
121 )
122 {
123  const labelList& eFaces = edgeFaces[edgeI];
124 
125  if (eFaces.size() == 2)
126  {
127  return (eFaces[0] != faceI ? eFaces[0] : eFaces[1]);
128  }
129  else
130  {
131  return -1;
132  }
133 }
134 
135 
136 void Foam::addPatchCellLayer::addVertex
137 (
138  const label pointI,
139  face& f,
140  label& fp
141 )
142 {
143  if (fp == 0)
144  {
145  f[fp++] = pointI;
146  }
147  else
148  {
149  if (f[fp-1] != pointI && f[0] != pointI)
150  {
151  f[fp++] = pointI;
152  }
153  }
154 
155  // Check for duplicates.
156  if (debug)
157  {
158  label n = 0;
159  for (label i = 0; i < fp; i++)
160  {
161  if (f[i] == pointI)
162  {
163  n++;
164 
165  if (n == 2)
166  {
167  f.setSize(fp);
169  (
170  "addPatchCellLayer::addVertex(const label, face&"
171  ", label&)"
172  ) << "Point " << pointI << " already present in face "
173  << f << abort(FatalError);
174  }
175  }
176  }
177  }
178 }
179 
180 
181 // Is edge to the same neighbour? (and needs extrusion and has not been
182 // dealt with already)
183 bool Foam::addPatchCellLayer::sameEdgeNeighbour
184 (
185  const indirectPrimitivePatch& pp,
186  const labelListList& globalEdgeFaces,
187  const boolList& doneEdge,
188  const label thisGlobalFaceI,
189  const label nbrGlobalFaceI,
190  const label edgeI
191 ) const
192 {
193  const edge& e = pp.edges()[edgeI];
194 
195  return
196  !doneEdge[edgeI] // not yet handled
197  && (
198  addedPoints_[e[0]].size() // is extruded
199  || addedPoints_[e[1]].size()
200  )
201  && (
202  nbrFace(globalEdgeFaces, edgeI, thisGlobalFaceI)
203  == nbrGlobalFaceI // is to same neighbour
204  );
205 }
206 
207 
208 // Collect consecutive string of edges that connects the same two
209 // (possibly coupled) faces. Returns -1 if no unvisited edge can be found.
210 // Otherwise returns start and end index in face.
211 Foam::labelPair Foam::addPatchCellLayer::getEdgeString
212 (
213  const indirectPrimitivePatch& pp,
214  const labelListList& globalEdgeFaces,
215  const boolList& doneEdge,
216  const label patchFaceI,
217  const label globalFaceI
218 ) const
219 {
220  const labelList& fEdges = pp.faceEdges()[patchFaceI];
221 
222  label startFp = -1;
223  label endFp = -1;
224 
225  // Get edge that hasn't been done yet but needs extrusion
226  forAll(fEdges, fp)
227  {
228  label edgeI = fEdges[fp];
229  const edge& e = pp.edges()[edgeI];
230 
231  if
232  (
233  !doneEdge[edgeI]
234  && ( addedPoints_[e[0]].size() || addedPoints_[e[1]].size() )
235  )
236  {
237  startFp = fp;
238  break;
239  }
240  }
241 
242  if (startFp != -1)
243  {
244  // We found an edge that needs extruding but hasn't been done yet.
245  // Now find the face on the other side
246  label nbrGlobalFaceI = nbrFace
247  (
248  globalEdgeFaces,
249  fEdges[startFp],
250  globalFaceI
251  );
252 
253  if (nbrGlobalFaceI == -1)
254  {
255  // Proper boundary edge. Only extrude single edge.
256  endFp = startFp;
257  }
258  else
259  {
260  // Search back for edge
261  // - which hasn't been handled yet
262  // - with same neighbour
263  // - that needs extrusion
264  while(true)
265  {
266  label prevFp = fEdges.rcIndex(startFp);
267 
268  if
269  (
270  !sameEdgeNeighbour
271  (
272  pp,
273  globalEdgeFaces,
274  doneEdge,
275  globalFaceI,
276  nbrGlobalFaceI,
277  fEdges[prevFp]
278  )
279  )
280  {
281  break;
282  }
283  startFp = prevFp;
284  }
285 
286  // Search forward for end of string
287  endFp = startFp;
288  while(true)
289  {
290  label nextFp = fEdges.fcIndex(endFp);
291 
292  if
293  (
294  !sameEdgeNeighbour
295  (
296  pp,
297  globalEdgeFaces,
298  doneEdge,
299  globalFaceI,
300  nbrGlobalFaceI,
301  fEdges[nextFp]
302  )
303  )
304  {
305  break;
306  }
307  endFp = nextFp;
308  }
309  }
310  }
311 
312  return labelPair(startFp, endFp);
313 }
314 
315 
316 // Adds a side face i.e. extrudes a patch edge.
317 Foam::label Foam::addPatchCellLayer::addSideFace
318 (
319  const indirectPrimitivePatch& pp,
320  const labelList& patchID, // prestored patch per pp face
321  const labelListList& addedCells, // per pp face the new extruded cell
322  const face& newFace,
323  const label ownFaceI, // pp face that provides owner
324  const label nbrFaceI,
325  const label patchEdgeI, // edge to add to
326  const label meshEdgeI, // corresponding mesh edge
327  const label layerI, // layer
328  const label numEdgeFaces, // number of layers for edge
329  const labelList& meshFaces, // precalculated edgeFaces
330  polyTopoChange& meshMod
331 ) const
332 {
333  // Edge to 'inflate' from
334  label inflateEdgeI = -1;
335 
336  // Check mesh faces using edge
337  forAll(meshFaces, i)
338  {
339  if (mesh_.isInternalFace(meshFaces[i]))
340  {
341  // meshEdge uses internal faces so ok to inflate from it
342  inflateEdgeI = meshEdgeI;
343  break;
344  }
345  }
346 
347 
348  // Get my mesh face and its zone.
349  label meshFaceI = pp.addressing()[ownFaceI];
350  // Zone info comes from any side patch face. Otherwise -1 since we
351  // don't know what to put it in - inherit from the extruded faces?
352  label zoneI = -1; //mesh_.faceZones().whichZone(meshFaceI);
353  bool flip = false;
354 
355  label addedFaceI = -1;
356 
357  // Is patch edge external edge of indirectPrimitivePatch?
358  if (nbrFaceI == -1)
359  {
360  // External edge so external face. Patch id is obtained from
361  // any other patch connected to edge.
362 
363  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
364 
365  // Loop over all faces connected to edge to inflate and
366  // see if any boundary face (but not meshFaceI)
367  label otherPatchID = patchID[ownFaceI];
368 
369  forAll(meshFaces, k)
370  {
371  label faceI = meshFaces[k];
372 
373  if
374  (
375  faceI != meshFaceI
376  && !mesh_.isInternalFace(faceI)
377  )
378  {
379  otherPatchID = patches.whichPatch(faceI);
380  zoneI = mesh_.faceZones().whichZone(faceI);
381  if (zoneI != -1)
382  {
383  label index = mesh_.faceZones()[zoneI].whichFace(faceI);
384  flip = mesh_.faceZones()[zoneI].flipMap()[index];
385  }
386  break;
387  }
388  }
389 
390  // Determine if different number of layer on owner and neighbour side
391  // (relevant only for coupled faces). See section for internal edge
392  // below.
393 
394  label layerOwn;
395 
396  if (addedCells[ownFaceI].size() < numEdgeFaces)
397  {
398  label offset = numEdgeFaces - addedCells[ownFaceI].size();
399  if (layerI <= offset)
400  {
401  layerOwn = 0;
402  }
403  else
404  {
405  layerOwn = layerI - offset;
406  }
407  }
408  else
409  {
410  layerOwn = layerI;
411  }
412 
413 
414  //Pout<< "Added boundary face:" << newFace
415  // << " own:" << addedCells[ownFaceI][layerOwn]
416  // << " patch:" << otherPatchID
417  // << endl;
418 
419  addedFaceI = meshMod.setAction
420  (
421  polyAddFace
422  (
423  newFace, // face
424  addedCells[ownFaceI][layerOwn], // owner
425  -1, // neighbour
426  -1, // master point
427  inflateEdgeI, // master edge
428  -1, // master face
429  false, // flux flip
430  otherPatchID, // patch for face
431  zoneI, // zone for face
432  flip // face zone flip
433  )
434  );
435  }
436  else
437  {
438  // When adding side faces we need to modify neighbour and owners
439  // in region where layer mesh is stopped. Determine which side
440  // has max number of faces and make sure layers match closest to
441  // original pp if there are different number of layers.
442 
443  label layerNbr;
444  label layerOwn;
445 
446  if (addedCells[ownFaceI].size() > addedCells[nbrFaceI].size())
447  {
448  label offset =
449  addedCells[ownFaceI].size() - addedCells[nbrFaceI].size();
450 
451  layerOwn = layerI;
452 
453  if (layerI <= offset)
454  {
455  layerNbr = 0;
456  }
457  else
458  {
459  layerNbr = layerI - offset;
460  }
461  }
462  else if (addedCells[nbrFaceI].size() > addedCells[ownFaceI].size())
463  {
464  label offset =
465  addedCells[nbrFaceI].size() - addedCells[ownFaceI].size();
466 
467  layerNbr = layerI;
468 
469  if (layerI <= offset)
470  {
471  layerOwn = 0;
472  }
473  else
474  {
475  layerOwn = layerI - offset;
476  }
477  }
478  else
479  {
480  // Same number of layers on both sides.
481  layerNbr = layerI;
482  layerOwn = layerI;
483  }
484 
485  addedFaceI = meshMod.setAction
486  (
487  polyAddFace
488  (
489  newFace, // face
490  addedCells[ownFaceI][layerOwn], // owner
491  addedCells[nbrFaceI][layerNbr], // neighbour
492  -1, // master point
493  inflateEdgeI, // master edge
494  -1, // master face
495  false, // flux flip
496  -1, // patch for face
497  zoneI, // zone for face
498  flip // face zone flip
499  )
500  );
501 
502  //Pout<< "Added internal face:" << newFace
503  // << " own:" << addedCells[ownFaceI][layerOwn]
504  // << " nei:" << addedCells[nbrFaceI][layerNbr]
505  // << endl;
506  }
507 
508  return addedFaceI;
509 }
510 
511 
512 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
513 
514 // Construct from mesh
515 Foam::addPatchCellLayer::addPatchCellLayer(const polyMesh& mesh)
516 :
517  mesh_(mesh),
518  addedPoints_(0),
519  layerFaces_(0)
520 {}
521 
522 
523 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
524 
526 (
527  const polyMesh& mesh,
528  const labelListList& layerFaces
529 )
530 {
531  labelListList layerCells(layerFaces.size());
532 
533  forAll(layerFaces, patchFaceI)
534  {
535  const labelList& faceLabels = layerFaces[patchFaceI];
536 
537  if (faceLabels.size())
538  {
539  labelList& added = layerCells[patchFaceI];
540  added.setSize(faceLabels.size()-1);
541 
542  for (label i = 0; i < faceLabels.size()-1; i++)
543  {
544  added[i] = mesh.faceNeighbour()[faceLabels[i]];
545  }
546  }
547  }
548  return layerCells;
549 }
550 
551 
553 {
554  return addedCells(mesh_, layerFaces_);
555 }
556 
557 
559 (
560  const scalarField& expansionRatio,
561  const indirectPrimitivePatch& pp,
562  const labelList& nFaceLayers,
563  const labelList& nPointLayers,
564  const vectorField& firstLayerDisp,
565  polyTopoChange& meshMod
566 )
567 {
568  if (debug)
569  {
570  Pout<< "addPatchCellLayer::setRefinement : Adding up to "
571  << max(nPointLayers)
572  << " layers of cells to indirectPrimitivePatch with "
573  << pp.nPoints() << " points" << endl;
574  }
575 
576  if
577  (
578  pp.nPoints() != firstLayerDisp.size()
579  || pp.nPoints() != nPointLayers.size()
580  || pp.size() != nFaceLayers.size()
581  )
582  {
584  (
585  "addPatchCellLayer::setRefinement"
586  "(const scalar, const indirectPrimitivePatch&"
587  ", const labelList&, const vectorField&, polyTopoChange&)"
588  ) << "Size of new points is not same as number of points used by"
589  << " the face subset" << endl
590  << " patch.nPoints:" << pp.nPoints()
591  << " displacement:" << firstLayerDisp.size()
592  << " nPointLayers:" << nPointLayers.size() << nl
593  << " patch.nFaces:" << pp.size()
594  << " nFaceLayers:" << nFaceLayers.size()
595  << abort(FatalError);
596  }
597 
598  forAll(nPointLayers, i)
599  {
600  if (nPointLayers[i] < 0)
601  {
603  (
604  "addPatchCellLayer::setRefinement"
605  "(const scalar, const indirectPrimitivePatch&"
606  ", const labelList&, const vectorField&, polyTopoChange&)"
607  ) << "Illegal number of layers " << nPointLayers[i]
608  << " at patch point " << i << abort(FatalError);
609  }
610  }
611  forAll(nFaceLayers, i)
612  {
613  if (nFaceLayers[i] < 0)
614  {
616  (
617  "addPatchCellLayer::setRefinement"
618  "(const scalar, const indirectPrimitivePatch&"
619  ", const labelList&, const vectorField&, polyTopoChange&)"
620  ) << "Illegal number of layers " << nFaceLayers[i]
621  << " at patch face " << i << abort(FatalError);
622  }
623  }
624 
625  const labelList& meshPoints = pp.meshPoints();
626 
627  // Some storage for edge-face-addressing.
629 
630  // Precalculate mesh edges for pp.edges.
631  labelList meshEdges(calcMeshEdges(mesh_, pp));
632 
633  if (debug)
634  {
635  // Check synchronisation
636  // ~~~~~~~~~~~~~~~~~~~~~
637 
638  {
639  labelList n(mesh_.nPoints(), 0);
640  UIndirectList<label>(n, meshPoints) = nPointLayers;
641  syncTools::syncPointList(mesh_, n, maxEqOp<label>(), 0, false);
642 
643  // Non-synced
644  forAll(meshPoints, i)
645  {
646  label meshPointI = meshPoints[i];
647 
648  if (n[meshPointI] != nPointLayers[i])
649  {
651  (
652  "addPatchCellLayer::setRefinement"
653  "(const scalar, const indirectPrimitivePatch&"
654  ", const labelList&, const vectorField&"
655  ", polyTopoChange&)"
656  ) << "At mesh point:" << meshPointI
657  << " coordinate:" << mesh_.points()[meshPointI]
658  << " specified nLayers:" << nPointLayers[i] << endl
659  << "On coupled point a different nLayers:"
660  << n[meshPointI] << " was specified."
661  << abort(FatalError);
662  }
663  }
664 
665 
666  // Check that nPointLayers equals the max layers of connected faces
667  // (or 0). Anything else makes no sense.
668  labelList nFromFace(mesh_.nPoints(), 0);
669  forAll(nFaceLayers, i)
670  {
671  const face& f = pp[i];
672 
673  forAll(f, fp)
674  {
675  label pointI = f[fp];
676 
677  nFromFace[pointI] = max(nFromFace[pointI], nFaceLayers[i]);
678  }
679  }
681  (
682  mesh_,
683  nFromFace,
684  maxEqOp<label>(),
685  0,
686  false
687  );
688 
689  forAll(nPointLayers, i)
690  {
691  label meshPointI = meshPoints[i];
692 
693  if
694  (
695  nPointLayers[i] > 0
696  && nPointLayers[i] != nFromFace[meshPointI]
697  )
698  {
700  (
701  "addPatchCellLayer::setRefinement"
702  "(const scalar, const indirectPrimitivePatch&"
703  ", const labelList&, const vectorField&"
704  ", polyTopoChange&)"
705  ) << "At mesh point:" << meshPointI
706  << " coordinate:" << mesh_.points()[meshPointI]
707  << " specified nLayers:" << nPointLayers[i] << endl
708  << "but the max nLayers of surrounding faces is:"
709  << nFromFace[meshPointI]
710  << abort(FatalError);
711  }
712  }
713  }
714 
715  {
716  pointField d(mesh_.nPoints(), wallPoint::greatPoint);
717  UIndirectList<point>(d, meshPoints) = firstLayerDisp;
719  (
720  mesh_,
721  d,
722  minEqOp<vector>(),
724  false
725  );
726 
727  forAll(meshPoints, i)
728  {
729  label meshPointI = meshPoints[i];
730 
731  if (mag(d[meshPointI] - firstLayerDisp[i]) > SMALL)
732  {
734  (
735  "addPatchCellLayer::setRefinement"
736  "(const scalar, const indirectPrimitivePatch&"
737  ", const labelList&, const vectorField&"
738  ", polyTopoChange&)"
739  ) << "At mesh point:" << meshPointI
740  << " coordinate:" << mesh_.points()[meshPointI]
741  << " specified displacement:" << firstLayerDisp[i]
742  << endl
743  << "On coupled point a different displacement:"
744  << d[meshPointI] << " was specified."
745  << abort(FatalError);
746  }
747  }
748  }
749 
750  // Check that edges of pp (so ones that become boundary faces)
751  // connect to only one boundary face. Guarantees uniqueness of
752  // patch that they go into so if this is a coupled patch both
753  // sides decide the same.
754  // ~~~~~~~~~~~~~~~~~~~~~~
755 
756  for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
757  {
758  const edge& e = pp.edges()[edgeI];
759 
760  if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
761  {
762  // Edge is to become a face
763 
764  const labelList& eFaces = pp.edgeFaces()[edgeI];
765 
766  // First check: pp should be single connected.
767  if (eFaces.size() != 1)
768  {
770  (
771  "addPatchCellLayer::setRefinement"
772  "(const scalar, const indirectPrimitivePatch&"
773  ", const labelList&, const vectorField&"
774  ", polyTopoChange&)"
775  ) << "boundary-edge-to-be-extruded:"
776  << pp.points()[meshPoints[e[0]]]
777  << pp.points()[meshPoints[e[1]]]
778  << " has more than two faces using it:" << eFaces
779  << abort(FatalError);
780  }
781 
782  label myFaceI = pp.addressing()[eFaces[0]];
783 
784  label meshEdgeI = meshEdges[edgeI];
785 
786  // Mesh faces using edge
787 
788  // Mesh faces using edge
789  const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
790 
791  // Check that there is only one patchface using edge.
792  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
793 
794  label bFaceI = -1;
795 
796  forAll(meshFaces, i)
797  {
798  label faceI = meshFaces[i];
799 
800  if (faceI != myFaceI)
801  {
802  if (!mesh_.isInternalFace(faceI))
803  {
804  if (bFaceI == -1)
805  {
806  bFaceI = faceI;
807  }
808  else
809  {
811  (
812  "addPatchCellLayer::setRefinement"
813  "(const scalar"
814  ", const indirectPrimitivePatch&"
815  ", const labelList&, const vectorField&"
816  ", polyTopoChange&)"
817  ) << "boundary-edge-to-be-extruded:"
818  << pp.points()[meshPoints[e[0]]]
819  << pp.points()[meshPoints[e[1]]]
820  << " has more than two boundary faces"
821  << " using it:"
822  << bFaceI << " fc:"
823  << mesh_.faceCentres()[bFaceI]
824  << " patch:" << patches.whichPatch(bFaceI)
825  << " and " << faceI << " fc:"
826  << mesh_.faceCentres()[faceI]
827  << " patch:" << patches.whichPatch(faceI)
828  << abort(FatalError);
829  }
830  }
831  }
832  }
833  }
834  }
835  }
836 
837 
838  // From master point (in patch point label) to added points (in mesh point
839  // label)
840  addedPoints_.setSize(pp.nPoints());
841 
842  // Mark points that do not get extruded by setting size of addedPoints_ to 0
843  label nTruncated = 0;
844 
845  forAll(nPointLayers, patchPointI)
846  {
847  if (nPointLayers[patchPointI] > 0)
848  {
849  addedPoints_[patchPointI].setSize(nPointLayers[patchPointI]);
850  }
851  else
852  {
853  nTruncated++;
854  }
855  }
856 
857  if (debug)
858  {
859  Pout<< "Not adding points at " << nTruncated << " out of "
860  << pp.nPoints() << " points" << endl;
861  }
862 
863 
864  //
865  // Create new points
866  //
867 
868  forAll(firstLayerDisp, patchPointI)
869  {
870  if (addedPoints_[patchPointI].size())
871  {
872  label meshPointI = meshPoints[patchPointI];
873 
874  label zoneI = mesh_.pointZones().whichZone(meshPointI);
875 
876  point pt = mesh_.points()[meshPointI];
877 
878  vector disp = firstLayerDisp[patchPointI];
879 
880  for (label i = 0; i < addedPoints_[patchPointI].size(); i++)
881  {
882  pt += disp;
883 
884  label addedVertI = meshMod.setAction
885  (
887  (
888  pt, // point
889  meshPointI, // master point
890  zoneI, // zone for point
891  true // supports a cell
892  )
893  );
894 
895  addedPoints_[patchPointI][i] = addedVertI;
896 
897  disp *= expansionRatio[patchPointI];
898  }
899  }
900  }
901 
902 
903  //
904  // Add cells to all boundaryFaces
905  //
906 
907  labelListList addedCells(pp.size());
908 
909  forAll(pp, patchFaceI)
910  {
911  if (nFaceLayers[patchFaceI] > 0)
912  {
913  addedCells[patchFaceI].setSize(nFaceLayers[patchFaceI]);
914 
915  label meshFaceI = pp.addressing()[patchFaceI];
916 
917  label ownZoneI = mesh_.cellZones().whichZone
918  (
919  mesh_.faceOwner()[meshFaceI]
920  );
921 
922  for (label i = 0; i < nFaceLayers[patchFaceI]; i++)
923  {
924  // Note: add from cell (owner of patch face) or from face?
925  // for now add from cell so we can map easily.
926  addedCells[patchFaceI][i] = meshMod.setAction
927  (
929  (
930  -1, // master point
931  -1, // master edge
932  -1, // master face
933  mesh_.faceOwner()[meshFaceI], // master cell id
934  ownZoneI // zone for cell
935  )
936  );
937 
938  //Pout<< "For patchFace:" << patchFaceI
939  // << " meshFace:" << pp.addressing()[patchFaceI]
940  // << " layer:" << i << " added cell:"
941  // << addedCells[patchFaceI][i]
942  // << endl;
943  }
944  }
945  }
946 
947 
948  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
949 
950  // Precalculated patchID for each patch face
951  labelList patchID(pp.size());
952 
953  forAll(pp, patchFaceI)
954  {
955  label meshFaceI = pp.addressing()[patchFaceI];
956 
957  patchID[patchFaceI] = patches.whichPatch(meshFaceI);
958  }
959 
960 
961 
962  // Create faces on top of the original patch faces.
963  // These faces are created from original patch faces outwards so in order
964  // of increasing cell number. So orientation should be same as original
965  // patch face for them to have owner<neighbour.
966 
967  layerFaces_.setSize(pp.size());
968 
969  forAll(pp.localFaces(), patchFaceI)
970  {
971  label meshFaceI = pp.addressing()[patchFaceI];
972 
973  if (addedCells[patchFaceI].size())
974  {
975  layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1);
976  layerFaces_[patchFaceI][0] = meshFaceI;
977 
978  // Get duplicated vertices on the patch face.
979  const face& f = pp.localFaces()[patchFaceI];
980 
981  face newFace(f.size());
982 
983  for (label i = 0; i < addedCells[patchFaceI].size(); i++)
984  {
985  forAll(f, fp)
986  {
987  if (addedPoints_[f[fp]].empty())
988  {
989  // Keep original point
990  newFace[fp] = meshPoints[f[fp]];
991  }
992  else
993  {
994  // Get new outside point
995  label offset =
996  addedPoints_[f[fp]].size()
997  - addedCells[patchFaceI].size();
998  newFace[fp] = addedPoints_[f[fp]][i+offset];
999  }
1000  }
1001 
1002 
1003  // Get new neighbour
1004  label nei;
1005  label patchI;
1006  label zoneI = -1;
1007  bool flip = false;
1008 
1009 
1010  if (i == addedCells[patchFaceI].size()-1)
1011  {
1012  // Top layer so is patch face.
1013  nei = -1;
1014  patchI = patchID[patchFaceI];
1015  zoneI = mesh_.faceZones().whichZone(meshFaceI);
1016  if (zoneI != -1)
1017  {
1018  const faceZone& fz = mesh_.faceZones()[zoneI];
1019  flip = fz.flipMap()[fz.whichFace(meshFaceI)];
1020  }
1021  }
1022  else
1023  {
1024  // Internal face between layer i and i+1
1025  nei = addedCells[patchFaceI][i+1];
1026  patchI = -1;
1027  }
1028 
1029 
1030  layerFaces_[patchFaceI][i+1] = meshMod.setAction
1031  (
1032  polyAddFace
1033  (
1034  newFace, // face
1035  addedCells[patchFaceI][i], // owner
1036  nei, // neighbour
1037  -1, // master point
1038  -1, // master edge
1039  meshFaceI, // master face for addition
1040  false, // flux flip
1041  patchI, // patch for face
1042  zoneI, // zone for face
1043  flip // face zone flip
1044  )
1045  );
1046  //Pout<< "Added inbetween face " << newFace
1047  // << " own:" << addedCells[patchFaceI][i]
1048  // << " nei:" << nei
1049  // << " patch:" << patchI
1050  // << endl;
1051  }
1052  }
1053  }
1054 
1055  //
1056  // Modify old patch faces to be on the inside
1057  //
1058  forAll(pp, patchFaceI)
1059  {
1060  if (addedCells[patchFaceI].size())
1061  {
1062  label meshFaceI = pp.addressing()[patchFaceI];
1063 
1064  meshMod.setAction
1065  (
1067  (
1068  pp[patchFaceI], // modified face
1069  meshFaceI, // label of face
1070  mesh_.faceOwner()[meshFaceI], // owner
1071  addedCells[patchFaceI][0], // neighbour
1072  false, // face flip
1073  -1, // patch for face
1074  true, //false, // remove from zone
1075  -1, //zoneI, // zone for face
1076  false // face flip in zone
1077  )
1078  );
1079  //Pout<< "Modified old patch face " << meshFaceI
1080  // << " own:" << mesh_.faceOwner()[meshFaceI]
1081  // << " nei:" << addedCells[patchFaceI][0]
1082  // << endl;
1083  }
1084  }
1085 
1086 
1087  //
1088  // Create 'side' faces, one per edge that is being extended.
1089  //
1090 
1091  const labelListList& faceEdges = pp.faceEdges();
1092  const faceList& localFaces = pp.localFaces();
1093  const edgeList& edges = pp.edges();
1094 
1095  // Get number of layers per edge. This is 0 if edge is not extruded;
1096  // max of connected faces otherwise.
1097  labelList edgeLayers(pp.nEdges());
1098 
1099  {
1100  // Use list over mesh.nEdges() since syncTools does not yet support
1101  // partial list synchronisation.
1102  labelList meshEdgeLayers(mesh_.nEdges(), -1);
1103 
1104  forAll(meshEdges, edgeI)
1105  {
1106  const edge& e = edges[edgeI];
1107 
1108  label meshEdgeI = meshEdges[edgeI];
1109 
1110  if ((nPointLayers[e[0]] == 0) && (nPointLayers[e[1]] == 0))
1111  {
1112  meshEdgeLayers[meshEdgeI] = 0;
1113  }
1114  else
1115  {
1116  const labelList& eFaces = pp.edgeFaces()[edgeI];
1117 
1118  forAll(eFaces, i)
1119  {
1120  meshEdgeLayers[meshEdgeI] = max
1121  (
1122  nFaceLayers[eFaces[i]],
1123  meshEdgeLayers[meshEdgeI]
1124  );
1125  }
1126  }
1127  }
1128 
1130  (
1131  mesh_,
1132  meshEdgeLayers,
1133  maxEqOp<label>(),
1134  0, // initial value
1135  false // no separation
1136  );
1137 
1138  forAll(meshEdges, edgeI)
1139  {
1140  edgeLayers[edgeI] = meshEdgeLayers[meshEdges[edgeI]];
1141  }
1142  }
1143 
1144 
1145  // Global indices engine
1146  const globalIndex globalFaces(mesh_.nFaces());
1147 
1148  // Get for all pp edgeFaces a unique faceID
1149  labelListList globalEdgeFaces
1150  (
1151  calcGlobalEdgeFaces
1152  (
1153  mesh_,
1154  globalFaces,
1155  pp,
1156  meshEdges
1157  )
1158  );
1159 
1160 
1161  // Mark off which edges have been extruded
1162  boolList doneEdge(pp.nEdges(), false);
1163 
1164 
1165  // Create faces. Per face walk connected edges and find string of edges
1166  // between the same two faces and extrude string into a single face.
1167  forAll(pp, patchFaceI)
1168  {
1169  const labelList& fEdges = faceEdges[patchFaceI];
1170 
1171  forAll(fEdges, fp)
1172  {
1173  // Get string of edges that needs to be extruded as a single face.
1174  // Returned as indices in fEdges.
1175  labelPair indexPair
1176  (
1177  getEdgeString
1178  (
1179  pp,
1180  globalEdgeFaces,
1181  doneEdge,
1182  patchFaceI,
1183  globalFaces.toGlobal(pp.addressing()[patchFaceI])
1184  )
1185  );
1186 
1187  //Pout<< "Found unextruded edges in edges:" << fEdges
1188  // << " start:" << indexPair[0]
1189  // << " end:" << indexPair[1]
1190  // << endl;
1191 
1192  const label startFp = indexPair[0];
1193  const label endFp = indexPair[1];
1194 
1195  if (startFp != -1)
1196  {
1197  // Extrude edges from indexPair[0] up to indexPair[1]
1198  // (note indexPair = indices of edges. There is one more vertex
1199  // than edges)
1200  const face& f = localFaces[patchFaceI];
1201 
1202  labelList stringedVerts;
1203  if (endFp >= startFp)
1204  {
1205  stringedVerts.setSize(endFp-startFp+2);
1206  }
1207  else
1208  {
1209  stringedVerts.setSize(endFp+f.size()-startFp+2);
1210  }
1211 
1212  label fp = startFp;
1213 
1214  for (label i = 0; i < stringedVerts.size()-1; i++)
1215  {
1216  stringedVerts[i] = f[fp];
1217  doneEdge[fEdges[fp]] = true;
1218  fp = f.fcIndex(fp);
1219  }
1220  stringedVerts[stringedVerts.size()-1] = f[fp];
1221 
1222 
1223  // Now stringedVerts contains the vertices in order of face f.
1224  // This is consistent with the order if f becomes the owner cell
1225  // and nbrFaceI the neighbour cell. Note that the cells get
1226  // added in order of pp so we can just use face ordering and
1227  // because we loop in incrementing order as well we will
1228  // always have nbrFaceI > patchFaceI.
1229 
1230  label startEdgeI = fEdges[startFp];
1231 
1232  label meshEdgeI = meshEdges[startEdgeI];
1233 
1234  label numEdgeSideFaces = edgeLayers[startEdgeI];
1235 
1236  for (label i = 0; i < numEdgeSideFaces; i++)
1237  {
1238  label vEnd = stringedVerts[stringedVerts.size()-1];
1239  label vStart = stringedVerts[0];
1240 
1241  // calculate number of points making up a face
1242  label newFp = 2*stringedVerts.size();
1243 
1244  if (i == 0)
1245  {
1246  // layer 0 gets all the truncation of neighbouring
1247  // faces with more layers.
1248  if (addedPoints_[vEnd].size())
1249  {
1250  newFp +=
1251  addedPoints_[vEnd].size() - numEdgeSideFaces;
1252  }
1253  if (addedPoints_[vStart].size())
1254  {
1255  newFp +=
1256  addedPoints_[vStart].size() - numEdgeSideFaces;
1257  }
1258  }
1259 
1260  face newFace(newFp);
1261 
1262  newFp = 0;
1263 
1264  // For layer 0 get pp points, for all other layers get
1265  // points of layer-1.
1266  if (i == 0)
1267  {
1268  forAll(stringedVerts, stringedI)
1269  {
1270  label v = stringedVerts[stringedI];
1271  addVertex(meshPoints[v], newFace, newFp);
1272  }
1273  }
1274  else
1275  {
1276  forAll(stringedVerts, stringedI)
1277  {
1278  label v = stringedVerts[stringedI];
1279  if (addedPoints_[v].size())
1280  {
1281  label offset =
1282  addedPoints_[v].size() - numEdgeSideFaces;
1283  addVertex
1284  (
1285  addedPoints_[v][i+offset-1],
1286  newFace,
1287  newFp
1288  );
1289  }
1290  else
1291  {
1292  addVertex(meshPoints[v], newFace, newFp);
1293  }
1294  }
1295  }
1296 
1297  // add points between stringed vertices (end)
1298  if (numEdgeSideFaces < addedPoints_[vEnd].size())
1299  {
1300  if (i == 0 && addedPoints_[vEnd].size())
1301  {
1302  label offset =
1303  addedPoints_[vEnd].size() - numEdgeSideFaces;
1304  for (label ioff = 0; ioff < offset; ioff++)
1305  {
1306  addVertex
1307  (
1308  addedPoints_[vEnd][ioff],
1309  newFace,
1310  newFp
1311  );
1312  }
1313  }
1314  }
1315 
1316  forAllReverse(stringedVerts, stringedI)
1317  {
1318  label v = stringedVerts[stringedI];
1319  if (addedPoints_[v].size())
1320  {
1321  label offset =
1322  addedPoints_[v].size() - numEdgeSideFaces;
1323  addVertex
1324  (
1325  addedPoints_[v][i+offset],
1326  newFace,
1327  newFp
1328  );
1329  }
1330  else
1331  {
1332  addVertex(meshPoints[v], newFace, newFp);
1333  }
1334  }
1335 
1336 
1337  // add points between stringed vertices (start)
1338  if (numEdgeSideFaces < addedPoints_[vStart].size())
1339  {
1340  if (i == 0 && addedPoints_[vStart].size())
1341  {
1342  label offset =
1343  addedPoints_[vStart].size() - numEdgeSideFaces;
1344  for (label ioff = offset; ioff > 0; ioff--)
1345  {
1346  addVertex
1347  (
1348  addedPoints_[vStart][ioff-1],
1349  newFace,
1350  newFp
1351  );
1352  }
1353  }
1354  }
1355 
1356  if (newFp >= 3)
1357  {
1358  // Add face inbetween faces patchFaceI and nbrFaceI
1359  // (possibly -1 for external edges)
1360 
1361  newFace.setSize(newFp);
1362 
1363  label nbrFaceI = nbrFace
1364  (
1365  pp.edgeFaces(),
1366  startEdgeI,
1367  patchFaceI
1368  );
1369 
1370  const labelList& meshFaces = mesh_.edgeFaces
1371  (
1372  meshEdgeI,
1373  ef
1374  );
1375 
1376  addSideFace
1377  (
1378  pp,
1379  patchID,
1380  addedCells,
1381  newFace,
1382  patchFaceI,
1383  nbrFaceI,
1384  startEdgeI, // edge to inflate from
1385  meshEdgeI, // corresponding mesh edge
1386  i,
1387  numEdgeSideFaces,
1388  meshFaces,
1389  meshMod
1390  );
1391  }
1392  }
1393  }
1394  }
1395  }
1396 }
1397 
1398 
1401  const mapPolyMesh& morphMap,
1402  const labelList& faceMap, // new to old patch faces
1403  const labelList& pointMap // new to old patch points
1404 )
1405 {
1406  {
1407  labelListList newAddedPoints(pointMap.size());
1408 
1409  forAll(newAddedPoints, newPointI)
1410  {
1411  label oldPointI = pointMap[newPointI];
1412 
1413  const labelList& added = addedPoints_[oldPointI];
1414 
1415  labelList& newAdded = newAddedPoints[newPointI];
1416  newAdded.setSize(added.size());
1417  label newI = 0;
1418 
1419  forAll(added, i)
1420  {
1421  label newPointI = morphMap.reversePointMap()[added[i]];
1422 
1423  if (newPointI >= 0)
1424  {
1425  newAdded[newI++] = newPointI;
1426  }
1427  }
1428  newAdded.setSize(newI);
1429  }
1430  addedPoints_.transfer(newAddedPoints);
1431  }
1432 
1433  {
1434  labelListList newLayerFaces(faceMap.size());
1435 
1436  forAll(newLayerFaces, newFaceI)
1437  {
1438  label oldFaceI = faceMap[newFaceI];
1439 
1440  const labelList& added = layerFaces_[oldFaceI];
1441 
1442  labelList& newAdded = newLayerFaces[newFaceI];
1443  newAdded.setSize(added.size());
1444  label newI = 0;
1445 
1446  forAll(added, i)
1447  {
1448  label newFaceI = morphMap.reverseFaceMap()[added[i]];
1449 
1450  if (newFaceI >= 0)
1451  {
1452  newAdded[newI++] = newFaceI;
1453  }
1454  }
1455  newAdded.setSize(newI);
1456  }
1457  layerFaces_.transfer(newLayerFaces);
1458  }
1459 }
1460 
1461 
1464  const primitiveMesh& mesh,
1465  const indirectPrimitivePatch& pp
1466 )
1467 {
1468  labelList meshEdges(pp.nEdges());
1469 
1470  forAll(meshEdges, patchEdgeI)
1471  {
1472  const edge& e = pp.edges()[patchEdgeI];
1473 
1474  label v0 = pp.meshPoints()[e[0]];
1475  label v1 = pp.meshPoints()[e[1]];
1476  meshEdges[patchEdgeI] = meshTools::findEdge
1477  (
1478  mesh.edges(),
1479  mesh.pointEdges()[v0],
1480  v0,
1481  v1
1482  );
1483  }
1484  return meshEdges;
1485 }
1486 
1487 
1488 // ************************ vim: set sw=4 sts=4 et: ************************ //