FreeFOAM The Cross-Platform CFD Toolkit
addCellLayer.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 "layerAdditionRemoval.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/primitiveMesh.H>
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 //- Calculate the offset to the next layer
39 Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
40 {
41  const polyMesh& mesh = topoChanger().mesh();
42  const primitiveFacePatch& masterFaceLayer =
43  mesh.faceZones()[faceZoneID_.index()]();
44 
45  const pointField& points = mesh.points();
46  const labelList& mp = masterFaceLayer.meshPoints();
47 
48  tmp<vectorField> textrusionDir(new vectorField(mp.size()));
49  vectorField& extrusionDir = textrusionDir();
50 
51  if (setLayerPairing())
52  {
53  if (debug)
54  {
55  Pout<< "void layerAdditionRemoval::extrusionDir() const "
56  << " for object " << name() << " : "
57  << "Using edges for point insertion" << endl;
58  }
59 
60  // Detected a valid layer. Grab the point and face collapse mapping
61  const labelList& ptc = pointsPairing();
62 
63  forAll (extrusionDir, mpI)
64  {
65  extrusionDir[mpI] = points[ptc[mpI]] - points[mp[mpI]];
66  }
67  }
68  else
69  {
70  if (debug)
71  {
72  Pout<< "void layerAdditionRemoval::extrusionDir() const "
73  << " for object " << name() << " : "
74  << "A valid layer could not be found in front of "
75  << "the addition face layer. Using face-based "
76  << "point normals for point addition"
77  << endl;
78  }
79 
80  extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
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  tmp<vectorField> tpointOffsets = 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_*tpointOffsets()[pointI],
139  mp[pointI], // master point
140  -1, // zone for point
141  true // supports a cell
142  )
143  );
144  }
145 
146 // Pout << "mp: " << mp << " addedPoints: " << addedPoints << endl;
147  // Create the cells
148 
149  const labelList& mc =
150  mesh.faceZones()[faceZoneID_.index()].masterCells();
151  const labelList& sc =
152  mesh.faceZones()[faceZoneID_.index()].slaveCells();
153 // Pout << "mc: " << mc << " sc: " << sc << endl;
154  const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
155  const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
156 
157  const labelList& own = mesh.faceOwner();
158  const labelList& nei = mesh.faceNeighbour();
159 
160  labelList addedCells(mf.size());
161 
162  forAll (mf, faceI)
163  {
164  addedCells[faceI] =
165  ref.setAction
166  (
167  polyAddCell
168  (
169  -1, // master point
170  -1, // master edge
171  mf[faceI], // master face
172  -1, // master cell
173  -1 // zone for cell
174  )
175  );
176  }
177 
178  // Create the new faces
179 
180  // Grab the local faces from the master zone
181  const faceList& zoneFaces = masterFaceLayer.localFaces();
182 
183  // Create the new faces by copying the master zone. All the added
184  // faces need to point into the newly created cells, which means
185  // all the faces from the master layer are flipped. The flux flip
186  // is determined from the original master layer face and the face
187  // owner: if the master cell is equal to the face owner the flux
188  // remains the same; otherwise it is flipped
189 
190  forAll (zoneFaces, faceI)
191  {
192  const face oldFace = zoneFaces[faceI].reverseFace();
193 
194  face newFace(oldFace.size());
195 
196  forAll (oldFace, pointI)
197  {
198  newFace[pointI] = addedPoints[oldFace[pointI]];
199  }
200 
201  bool flipFaceFlux = false;
202 
203  // Flip the face as necessary
204  if
205  (
206  mc[faceI] == nei[mf[faceI]]
207  || !mesh.isInternalFace(mf[faceI])
208  )
209  {
210  flipFaceFlux = true;
211  }
212 
213  // Add the face
214  ref.setAction
215  (
216  polyAddFace
217  (
218  newFace, // face
219  mc[faceI], // owner
220  addedCells[faceI], // neighbour
221  -1, // master point
222  -1, // master edge
223  mf[faceI], // master face for addition
224  flipFaceFlux, // flux flip
225  -1, // patch for face
226  -1, // zone for face
227  false // face zone flip
228  )
229  );
230 
231 // Pout << "adding face: " << newFace << " own: " << mc[faceI] << " nei: " << addedCells[faceI] << endl;
232  }
233 
234  // Modify the faces from the master zone for the new neighbour
235 
236  const faceList& faces = mesh.faces();
237 // Pout << "mfFlip: " << mfFlip << endl;
238  forAll (mf, faceI)
239  {
240  const label curfaceID = mf[faceI];
241 
242  // If the face is internal, modify its owner to be the newly
243  // created cell. No flip is necessary
244  if (!mesh.isInternalFace(curfaceID))
245  {
246  ref.setAction
247  (
248  polyModifyFace
249  (
250  faces[curfaceID], // modified face
251  curfaceID, // label of face being modified
252  addedCells[faceI], // owner
253  -1, // neighbour
254  false, // face flip
255  mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
256  false, // remove from zone
257  faceZoneID_.index(), // zone for face
258  mfFlip[faceI] // face flip in zone
259  )
260  );
261 // Pout << "Modifying a boundary face. Face: " << curfaceID << " flip: " << mfFlip[faceI] << endl;
262  }
263  // If slave cell is owner, the face remains the same (but with
264  // a new neighbour - the newly created cell). Otherwise, the
265  // face is flipped.
266  else if (sc[faceI] == own[curfaceID])
267  {
268  // Orientation is good, just change neighbour
269  ref.setAction
270  (
271  polyModifyFace
272  (
273  faces[curfaceID], // modified face
274  curfaceID, // label of face being modified
275  own[curfaceID], // owner
276  addedCells[faceI], // neighbour
277  false, // face flip
278  mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
279  false, // remove from zone
280  faceZoneID_.index(), // zone for face
281  mfFlip[faceI] // face flip in zone
282  )
283  );
284 
285 // Pout << "modify face, no flip " << curfaceID << " own: " << own[curfaceID] << " nei: " << addedCells[faceI] << endl;
286  }
287  else
288  {
289  // Change in orientation; flip face
290  ref.setAction
291  (
292  polyModifyFace
293  (
294  faces[curfaceID].reverseFace(), // modified face
295  curfaceID, // label of face being modified
296  nei[curfaceID], // owner
297  addedCells[faceI], // neighbour
298  true, // face flip
299  mesh.boundaryMesh().whichPatch(curfaceID), // patch for face
300  false, // remove from zone
301  faceZoneID_.index(), // zone for face
302  !mfFlip[faceI] // face flip in zone
303  )
304  );
305 // Pout << "modify face, with flip " << curfaceID << " own: " << own[curfaceID] << " nei: " << addedCells[faceI] << endl;
306  }
307  }
308 
309  // Create new faces from edges
310 
311  const edgeList& zoneLocalEdges = masterFaceLayer.edges();
312 
313  const labelListList& edgeFaces = masterFaceLayer.edgeFaces();
314 
315  label nInternalEdges = masterFaceLayer.nInternalEdges();
316 
317  const labelList& meshEdges =
318  mesh.faceZones()[faceZoneID_.index()].meshEdges();
319 
320  // Do all internal edges
321  for (label curEdgeID = 0; curEdgeID < nInternalEdges; curEdgeID++)
322  {
323  face newFace(4);
324 
325  newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
326  newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
327  newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
328  newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
329 
330  ref.setAction
331  (
332  polyAddFace
333  (
334  newFace, // face
335  addedCells[edgeFaces[curEdgeID][0]], // owner
336  addedCells[edgeFaces[curEdgeID][1]], // neighbour
337  -1, // master point
338  meshEdges[curEdgeID], // master edge
339  -1, // master face
340  false, // flip flux
341  -1, // patch for face
342  -1, // zone for face
343  false // face zone flip
344  )
345  );
346 
347 // Pout << "Add internal face off edge: " << newFace << " own: " << addedCells[edgeFaces[curEdgeID][0]] << " nei: " << addedCells[edgeFaces[curEdgeID][1]] << endl;
348  }
349 
350  // Prepare creation of faces from boundary edges.
351  // Note:
352  // A tricky part of the algorithm is finding the patch into which the
353  // newly created face will be added. For this purpose, take the edge
354  // and grab all the faces coming from it. From the set of faces
355  // eliminate the internal faces and find the boundary face from the rest.
356  // If there is more than one boundary face (excluding the ones in
357  // the master zone), the patch with the lower label is selected.
358 
359  const labelListList& meshEdgeFaces = mesh.edgeFaces();
360 
361  const faceZoneMesh& zoneMesh = mesh.faceZones();
362 
363  // Do all boundary edges
364 
365  for
366  (
367  label curEdgeID = nInternalEdges;
368  curEdgeID < zoneLocalEdges.size();
369  curEdgeID++
370  )
371  {
372  face newFace(4);
373  newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
374  newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
375  newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
376  newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
377 
378  // Determine the patch for insertion
379  const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
380 
381  label patchID = -1;
382  label zoneID = -1;
383 
384  forAll (curFaces, faceI)
385  {
386  const label cf = curFaces[faceI];
387 
388  if (!mesh.isInternalFace(cf))
389  {
390  // Face not internal. Check if it is in the zone
391  if (zoneMesh.whichZone(cf) != faceZoneID_.index())
392  {
393  // Found the face in a boundary patch which is not in zone
394  patchID = mesh.boundaryMesh().whichPatch(cf);
395  zoneID = mesh.faceZones().whichZone(cf);
396 
397  break;
398  }
399  }
400  }
401 
402  if (patchID < 0)
403  {
405  (
406  "void Foam::layerAdditionRemoval::setRefinement("
407  "polyTopoChange& ref)"
408  ) << "Cannot find patch for edge " << meshEdges[curEdgeID]
409  << ". Edge: " << mesh.edges()[meshEdges[curEdgeID]]
410  << abort(FatalError);
411  }
412 
413  ref.setAction
414  (
415  polyAddFace
416  (
417  newFace, // face
418  addedCells[edgeFaces[curEdgeID][0]], // owner
419  -1, // neighbour
420  -1, // master point
421  meshEdges[curEdgeID], // master edge
422  -1, // master face
423  false, // flip flux
424  patchID, // patch for face
425  zoneID, // zone
426  false // zone face flip
427  )
428  );
429 // Pout << "add boundary face: " << newFace << " into patch " << patchID << " own: " << addedCells[edgeFaces[curEdgeID][0]] << endl;
430  }
431 
432  // Modify the remaining faces of the master cells to reconnect to the new
433  // layer of faces.
434  // Algorithm: Go through all the cells of the master zone and make
435  // a map of faces to avoid duplicates. Do not insert the faces in
436  // the master patch (as they have already been dealt with). Make
437  // a master layer point renumbering map, which for every point in
438  // the master layer gives its new label. Loop through all faces in
439  // the map and attempt to renumber them using the master layer
440  // point renumbering map. Once the face is renumbered, compare it
441  // with the original face; if they are the same, the face has not
442  // changed; if not, modify the face but keep all of its old
443  // attributes (apart from the vertex numbers).
444 
445  // Create the map of faces in the master cell layer
446  labelHashSet masterCellFaceMap(primitiveMesh::facesPerCell_*mc.size());
447 
448  const cellList& cells = mesh.cells();
449 
450  forAll (mc, cellI)
451  {
452  const labelList& curFaces = cells[mc[cellI]];
453 
454  forAll (curFaces, faceI)
455  {
456  // Check if the face belongs to the master zone; if not add it
457  if (zoneMesh.whichZone(curFaces[faceI]) != faceZoneID_.index())
458  {
459  masterCellFaceMap.insert(curFaces[faceI]);
460  }
461  }
462  }
463 
464  // Create the master layer point map
465  Map<label> masterLayerPointMap(2*mp.size());
466 
467  forAll (mp, pointI)
468  {
469  masterLayerPointMap.insert
470  (
471  mp[pointI],
472  addedPoints[pointI]
473  );
474  }
475 
476  // Grab the list of faces of the master layer
477  const labelList masterCellFaces = masterCellFaceMap.toc();
478 
479  forAll (masterCellFaces, faceI)
480  {
481  // Attempt to renumber the face using the masterLayerPointMap.
482  // Missing point remain the same
483 
484  const label curFaceID = masterCellFaces[faceI];
485 
486  const face& oldFace = faces[curFaceID];
487 
488  face newFace(oldFace.size());
489 
490  bool changed = false;
491 
492  forAll (oldFace, pointI)
493  {
494  if (masterLayerPointMap.found(oldFace[pointI]))
495  {
496  changed = true;
497 
498  newFace[pointI] = masterLayerPointMap.find(oldFace[pointI])();
499  }
500  else
501  {
502  newFace[pointI] = oldFace[pointI];
503  }
504  }
505 
506  // If the face has changed, create a modification entry
507  if (changed)
508  {
509  // Get face zone and its flip
510  label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
511  bool modifiedFaceZoneFlip = false;
512 
513  if (modifiedFaceZone >= 0)
514  {
515  modifiedFaceZoneFlip =
516  mesh.faceZones()[modifiedFaceZone].flipMap()
517  [
518  mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
519  ];
520  }
521 
522  if (mesh.isInternalFace(curFaceID))
523  {
524  ref.setAction
525  (
526  polyModifyFace
527  (
528  newFace, // modified face
529  curFaceID, // label of face being modified
530  own[curFaceID], // owner
531  nei[curFaceID], // neighbour
532  false, // face flip
533  -1, // patch for face
534  false, // remove from zone
535  modifiedFaceZone, // zone for face
536  modifiedFaceZoneFlip // face flip in zone
537  )
538  );
539 // Pout << "modifying stick-out face. Internal Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " nei: " << nei[curFaceID] << endl;
540  }
541  else
542  {
543  ref.setAction
544  (
545  polyModifyFace
546  (
547  newFace, // modified face
548  curFaceID, // label of face being modified
549  own[curFaceID], // owner
550  -1, // neighbour
551  false, // face flip
552  mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
553  false, // remove from zone
554  modifiedFaceZone, // zone for face
555  modifiedFaceZoneFlip // face flip in zone
556  )
557  );
558 // Pout << "modifying stick-out face. Boundary Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " patch: " << mesh.boundaryMesh().whichPatch(curFaceID) << endl;
559  }
560  }
561  }
562 
563  if (debug)
564  {
565  Pout<< "void layerAdditionRemoval::addCellLayer("
566  << "polyTopoChange& ref) const "
567  << " for object " << name() << " : "
568  << "Finished adding cell layer" << endl;
569  }
570 }
571 
572 
573 // ************************ vim: set sw=4 sts=4 et: ************************ //