addCellLayer.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "layerAdditionRemoval.H"
27 #include "polyMesh.H"
28 #include "primitiveMesh.H"
29 #include "polyTopoChange.H"
30 #include "polyTopoChanger.H"
31 #include "polyAddPoint.H"
32 #include "polyAddCell.H"
33 #include "polyAddFace.H"
34 #include "polyModifyFace.H"
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
39 {
40  const polyMesh& mesh = topoChanger().mesh();
41  const primitiveFacePatch& masterFaceLayer =
42  mesh.faceZones()[faceZoneID_.index()]();
43 
44  const pointField& points = mesh.points();
45  const labelList& mp = masterFaceLayer.meshPoints();
46 
47  tmp<vectorField> textrusionDir(new vectorField(mp.size()));
48  vectorField& extrusionDir = textrusionDir.ref();
49 
50  if (setLayerPairing())
51  {
52  if (debug)
53  {
54  Pout<< "void layerAdditionRemoval::extrusionDir() const "
55  << " for object " << name() << " : "
56  << "Using edges for point insertion" << endl;
57  }
58 
59  // Detected a valid layer. Grab the point and face collapse mapping
60  const labelList& ptc = pointsPairing();
61 
62  forAll(extrusionDir, mpI)
63  {
64  extrusionDir[mpI] = points[ptc[mpI]] - points[mp[mpI]];
65  }
66  }
67  else
68  {
69  if (debug)
70  {
71  Pout<< "void layerAdditionRemoval::extrusionDir() const "
72  << " for object " << name() << " : "
73  << "A valid layer could not be found in front of "
74  << "the addition face layer. Using face-based "
75  << "point normals for point addition"
76  << endl;
77  }
78 
79  extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
80  }
81 
82  return textrusionDir;
83 }
84 
85 
86 void Foam::layerAdditionRemoval::addCellLayer
87 (
88  polyTopoChange& ref
89 ) const
90 {
91  // Insert the layer addition instructions into the topological change
92 
93  // Algorithm:
94  // 1. For every point in the master zone add a new point
95  // 2. For every face in the master zone add a cell
96  // 3. Go through all the edges of the master zone. For all internal edges,
97  // add a face with the correct orientation and make it internal.
98  // For all boundary edges, find the patch it belongs to and add the face
99  // to this patch
100  // 4. Create a set of new faces on the patch image and assign them to be
101  // between the old master cells and the newly created cells
102  // 5. Modify all the faces in the patch such that they are located
103  // between the old slave cells and newly created cells
104 
105  if (debug)
106  {
107  Pout<< "void layerAdditionRemoval::addCellLayer("
108  << "polyTopoChange& ref) const for object " << name() << " : "
109  << "Adding cell layer" << endl;
110  }
111 
112  // Create the points
113 
114  const polyMesh& mesh = topoChanger().mesh();
115  const primitiveFacePatch& masterFaceLayer =
116  mesh.faceZones()[faceZoneID_.index()]();
117 
118  const pointField& points = mesh.points();
119  const labelList& mp = masterFaceLayer.meshPoints();
120 
121  // Get the extrusion direction for the added points
122 
123  const vectorField pointOffsets(extrusionDir());
124 
125  // Add the new points
126  labelList addedPoints(mp.size());
127 
128  forAll(mp, pointi)
129  {
130  // Add the point nominal thickness away from the master zone point
131  // and grab the label
132  addedPoints[pointi] =
133  ref.setAction
134  (
135  polyAddPoint
136  (
137  points[mp[pointi]] // point
138  + addDelta_*pointOffsets[pointi],
139  mp[pointi], // master point
140  -1, // zone for point
141  true // supports a cell
142  )
143  );
144  }
145 
146  if (debug > 1)
147  {
148  Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
149  }
150 
151  // Create the cells
152 
153  const labelList& mc =
154  mesh.faceZones()[faceZoneID_.index()].masterCells();
155  const labelList& sc =
156  mesh.faceZones()[faceZoneID_.index()].slaveCells();
157 
158  const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
159  const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
160 
161  const labelList& own = mesh.faceOwner();
162  const labelList& nei = mesh.faceNeighbour();
163 
164  labelList addedCells(mf.size());
165 
166  forAll(mf, facei)
167  {
168  label celli = mc[facei];
169  label zoneI = mesh.cellZones().whichZone(celli);
170 
171  addedCells[facei] =
172  ref.setAction
173  (
174  polyAddCell
175  (
176  -1, // master point
177  -1, // master edge
178  mf[facei], // master face
179  -1, // master cell
180  zoneI // zone for cell
181  )
182  );
183  }
184 
185  // Create the new faces
186 
187  // Grab the local faces from the master zone
188  const faceList& zoneFaces = masterFaceLayer.localFaces();
189 
190  // Create the new faces by copying the master zone. All the added
191  // faces need to point into the newly created cells, which means
192  // all the faces from the master layer are flipped. The flux flip
193  // is determined from the original master layer face and the face
194  // owner: if the master cell is equal to the face owner the flux
195  // remains the same; otherwise it is flipped
196 
197  forAll(zoneFaces, facei)
198  {
199  const face oldFace = zoneFaces[facei].reverseFace();
200 
201  face newFace(oldFace.size());
202 
203  forAll(oldFace, pointi)
204  {
205  newFace[pointi] = addedPoints[oldFace[pointi]];
206  }
207 
208  bool flipFaceFlux = false;
209 
210  // Flip the face as necessary
211  if
212  (
213  !mesh.isInternalFace(mf[facei])
214  || mc[facei] == nei[mf[facei]]
215  )
216  {
217  flipFaceFlux = true;
218  }
219 
220  // Add the face
221  ref.setAction
222  (
223  polyAddFace
224  (
225  newFace, // face
226  mc[facei], // owner
227  addedCells[facei], // neighbour
228  -1, // master point
229  -1, // master edge
230  mf[facei], // master face for addition
231  flipFaceFlux, // flux flip
232  -1, // patch for face
233  -1, // zone for face
234  false // face zone flip
235  )
236  );
237 
238  if (debug > 1)
239  {
240  Pout<< "adding face: " << newFace
241  << " own: " << mc[facei]
242  << " nei: " << addedCells[facei]
243  << endl;
244  }
245  }
246 
247  // Modify the faces from the master zone for the new neighbour
248 
249  const faceList& faces = mesh.faces();
250 
251  // Pout<< "mfFlip: " << mfFlip << endl;
252 
253  forAll(mf, facei)
254  {
255  const label curfaceID = mf[facei];
256 
257  // If the face is internal, modify its owner to be the newly
258  // created cell. No flip is necessary
259  if (!mesh.isInternalFace(curfaceID))
260  {
261  ref.setAction
262  (
263  polyModifyFace
264  (
265  faces[curfaceID], // modified face
266  curfaceID, // label of face being modified
267  addedCells[facei], // owner
268  -1, // neighbour
269  false, // face flip
270  mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
271  false, // remove from zone
272  faceZoneID_.index(), // zone for face
273  mfFlip[facei] // face flip in zone
274  )
275  );
276 
277  if (debug > 1)
278  {
279  Pout<< "Modifying a boundary face. Face: " << curfaceID
280  << " flip: " << mfFlip[facei]
281  << endl;
282  }
283  }
284 
285  // If slave cell is owner, the face remains the same (but with
286  // a new neighbour - the newly created cell). Otherwise, the
287  // face is flipped.
288  else if (sc[facei] == own[curfaceID])
289  {
290  // Orientation is good, just change neighbour
291  ref.setAction
292  (
293  polyModifyFace
294  (
295  faces[curfaceID], // modified face
296  curfaceID, // label of face being modified
297  own[curfaceID], // owner
298  addedCells[facei], // neighbour
299  false, // face flip
300  mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
301  false, // remove from zone
302  faceZoneID_.index(), // zone for face
303  mfFlip[facei] // face flip in zone
304  )
305  );
306 
307  if (debug > 1)
308  {
309  Pout<< "modify face, no flip " << curfaceID
310  << " own: " << own[curfaceID]
311  << " nei: " << addedCells[facei]
312  << endl;
313  }
314  }
315  else
316  {
317  // Change in orientation; flip face
318  ref.setAction
319  (
320  polyModifyFace
321  (
322  faces[curfaceID].reverseFace(), // modified face
323  curfaceID, // label of face being modified
324  nei[curfaceID], // owner
325  addedCells[facei], // neighbour
326  true, // face flip
327  mesh.boundaryMesh().whichPatch(curfaceID), // patch for face
328  false, // remove from zone
329  faceZoneID_.index(), // zone for face
330  !mfFlip[facei] // face flip in zone
331  )
332  );
333 
334  if (debug > 1)
335  {
336  Pout<< "modify face, with flip " << curfaceID
337  << " own: " << own[curfaceID]
338  << " nei: " << addedCells[facei]
339  << endl;
340  }
341  }
342  }
343 
344  // Create new faces from edges
345 
346  const edgeList& zoneLocalEdges = masterFaceLayer.edges();
347 
348  const labelListList& edgeFaces = masterFaceLayer.edgeFaces();
349 
350  label nInternalEdges = masterFaceLayer.nInternalEdges();
351 
352  const labelList& meshEdges =
353  mesh.faceZones()[faceZoneID_.index()].meshEdges();
354 
355  // Do all internal edges
356  for (label curEdgeID = 0; curEdgeID < nInternalEdges; curEdgeID++)
357  {
358  face newFace(4);
359 
360  newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
361  newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
362  newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
363  newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
364 
365  ref.setAction
366  (
367  polyAddFace
368  (
369  newFace, // face
370  addedCells[edgeFaces[curEdgeID][0]], // owner
371  addedCells[edgeFaces[curEdgeID][1]], // neighbour
372  -1, // master point
373  meshEdges[curEdgeID], // master edge
374  -1, // master face
375  false, // flip flux
376  -1, // patch for face
377  -1, // zone for face
378  false // face zone flip
379  )
380  );
381 
382  if (debug > 1)
383  {
384  Pout<< "Add internal face off edge: " << newFace
385  << " own: " << addedCells[edgeFaces[curEdgeID][0]]
386  << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
387  << endl;
388  }
389  }
390 
391  // Prepare creation of faces from boundary edges.
392  // Note:
393  // A tricky part of the algorithm is finding the patch into which the
394  // newly created face will be added. For this purpose, take the edge
395  // and grab all the faces coming from it. From the set of faces
396  // eliminate the internal faces and find the boundary face from the rest.
397  // If there is more than one boundary face (excluding the ones in
398  // the master zone), the patch with the lower label is selected.
399 
400  const labelListList& meshEdgeFaces = mesh.edgeFaces();
401 
402  const meshFaceZones& meshZones = mesh.faceZones();
403 
404  // Do all boundary edges
405 
406  for
407  (
408  label curEdgeID = nInternalEdges;
409  curEdgeID < zoneLocalEdges.size();
410  curEdgeID++
411  )
412  {
413  face newFace(4);
414  newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
415  newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
416  newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
417  newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
418 
419  // Determine the patch for insertion
420  const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
421 
422  label patchID = -1;
423  label zoneID = -1;
424 
425  forAll(curFaces, facei)
426  {
427  const label cf = curFaces[facei];
428 
429  if (!mesh.isInternalFace(cf))
430  {
431  // Face not internal. Check if it is in the zone
432  if (meshZones.whichZone(cf) != faceZoneID_.index())
433  {
434  // Found the face in a boundary patch which is not in zone
435  patchID = mesh.boundaryMesh().whichPatch(cf);
436  zoneID = mesh.faceZones().whichZone(cf);
437 
438  break;
439  }
440  }
441  }
442 
443  if (patchID < 0)
444  {
446  << "Cannot find patch for edge " << meshEdges[curEdgeID]
447  << ". Edge: " << mesh.edges()[meshEdges[curEdgeID]]
448  << abort(FatalError);
449  }
450 
451  ref.setAction
452  (
453  polyAddFace
454  (
455  newFace, // face
456  addedCells[edgeFaces[curEdgeID][0]], // owner
457  -1, // neighbour
458  -1, // master point
459  meshEdges[curEdgeID], // master edge
460  -1, // master face
461  false, // flip flux
462  patchID, // patch for face
463  zoneID, // zone
464  false // zone face flip
465  )
466  );
467 
468  if (debug > 1)
469  {
470  Pout<< "add boundary face: " << newFace
471  << " into patch " << patchID
472  << " own: " << addedCells[edgeFaces[curEdgeID][0]]
473  << endl;
474  }
475  }
476 
477  // Modify the remaining faces of the master cells to reconnect to the new
478  // layer of faces.
479  // Algorithm: Go through all the cells of the master zone and make
480  // a map of faces to avoid duplicates. Do not insert the faces in
481  // the master patch (as they have already been dealt with). Make
482  // a master layer point renumbering map, which for every point in
483  // the master layer gives its new label. Loop through all faces in
484  // the map and attempt to renumber them using the master layer
485  // point renumbering map. Once the face is renumbered, compare it
486  // with the original face; if they are the same, the face has not
487  // changed; if not, modify the face but keep all of its old
488  // attributes (apart from the vertex numbers).
489 
490  // Create the map of faces in the master cell layer
491  labelHashSet masterCellFaceMap(primitiveMesh::facesPerCell_*mc.size());
492 
493  const cellList& cells = mesh.cells();
494 
495  forAll(mc, celli)
496  {
497  const labelList& curFaces = cells[mc[celli]];
498 
499  forAll(curFaces, facei)
500  {
501  // Check if the face belongs to the master zone; if not add it
502  if (meshZones.whichZone(curFaces[facei]) != faceZoneID_.index())
503  {
504  masterCellFaceMap.insert(curFaces[facei]);
505  }
506  }
507  }
508 
509  // Create the master layer point map
510  Map<label> masterLayerPointMap(2*mp.size());
511 
512  forAll(mp, pointi)
513  {
514  masterLayerPointMap.insert
515  (
516  mp[pointi],
517  addedPoints[pointi]
518  );
519  }
520 
521  // Grab the list of faces of the master layer
522  const labelList masterCellFaces = masterCellFaceMap.toc();
523 
524  forAll(masterCellFaces, facei)
525  {
526  // Attempt to renumber the face using the masterLayerPointMap.
527  // Missing point remain the same
528 
529  const label curFaceID = masterCellFaces[facei];
530 
531  const face& oldFace = faces[curFaceID];
532 
533  face newFace(oldFace.size());
534 
535  bool changed = false;
536 
537  forAll(oldFace, pointi)
538  {
539  if (masterLayerPointMap.found(oldFace[pointi]))
540  {
541  changed = true;
542 
543  newFace[pointi] = masterLayerPointMap.find(oldFace[pointi])();
544  }
545  else
546  {
547  newFace[pointi] = oldFace[pointi];
548  }
549  }
550 
551  // If the face has changed, create a modification entry
552  if (changed)
553  {
554  // Get face zone and its flip
555  label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
556  bool modifiedFaceZoneFlip = false;
557 
558  if (modifiedFaceZone >= 0)
559  {
560  modifiedFaceZoneFlip =
561  mesh.faceZones()[modifiedFaceZone].flipMap()
562  [
563  mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
564  ];
565  }
566 
567  if (mesh.isInternalFace(curFaceID))
568  {
569  ref.setAction
570  (
571  polyModifyFace
572  (
573  newFace, // modified face
574  curFaceID, // label of face being modified
575  own[curFaceID], // owner
576  nei[curFaceID], // neighbour
577  false, // face flip
578  -1, // patch for face
579  false, // remove from zone
580  modifiedFaceZone, // zone for face
581  modifiedFaceZoneFlip // face flip in zone
582  )
583  );
584 
585  if (debug > 1)
586  {
587  Pout<< "modifying stick-out face. Internal Old face: "
588  << oldFace
589  << " new face: " << newFace
590  << " own: " << own[curFaceID]
591  << " nei: " << nei[curFaceID]
592  << endl;
593  }
594  }
595  else
596  {
597  ref.setAction
598  (
599  polyModifyFace
600  (
601  newFace, // modified face
602  curFaceID, // label of face being modified
603  own[curFaceID], // owner
604  -1, // neighbour
605  false, // face flip
606  mesh.boundaryMesh().whichPatch(curFaceID),
607  // patch for face
608  false, // remove from zone
609  modifiedFaceZone, // zone for face
610  modifiedFaceZoneFlip // face flip in zone
611  )
612  );
613 
614  if (debug > 1)
615  {
616  Pout<< "modifying stick-out face. Boundary Old face: "
617  << oldFace
618  << " new face: " << newFace
619  << " own: " << own[curFaceID]
620  << " patch: "
621  << mesh.boundaryMesh().whichPatch(curFaceID)
622  << endl;
623  }
624  }
625  }
626  }
627 
628  if (debug)
629  {
630  Pout<< "void layerAdditionRemoval::addCellLayer(polyTopoChange&) const "
631  << " for object " << name() << ": "
632  << "Finished adding cell layer" << endl;
633  }
634 }
635 
636 
637 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
static const unsigned facesPerCell_
Estimated number of faces per cell.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< edge > edgeList
Definition: edgeList.H:38
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
MeshZones< faceZone, polyMesh > meshFaceZones
A MeshZones with the type faceZone.
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const word & name() const
Return name of this modifier.
Field< vector > vectorField
Specialisation of Field<T> for vector.
A class for managing temporary objects.
Definition: PtrList.H:53
const polyMesh & mesh() const
Return the mesh reference.
List< cell > cellList
list of cells
Definition: cellList.H:42