FreeFOAM The Cross-Platform CFD Toolkit
cellSplitter.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 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "cellSplitter.H"
29 #include <OpenFOAM/polyMesh.H>
35 #include <OpenFOAM/mapPolyMesh.H>
36 #include <meshTools/meshTools.H>
37 
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43 
44 defineTypeNameAndDebug(cellSplitter, 0);
45 
46 }
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void Foam::cellSplitter::getFaceInfo
51 (
52  const label faceI,
53  label& patchID,
54  label& zoneID,
55  label& zoneFlip
56 ) const
57 {
58  patchID = -1;
59 
60  if (!mesh_.isInternalFace(faceI))
61  {
62  patchID = mesh_.boundaryMesh().whichPatch(faceI);
63  }
64 
65  zoneID = mesh_.faceZones().whichZone(faceI);
66 
67  zoneFlip = false;
68 
69  if (zoneID >= 0)
70  {
71  const faceZone& fZone = mesh_.faceZones()[zoneID];
72 
73  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
74  }
75 }
76 
77 
78 // Find the new owner of faceI (since the original cell has been split into
79 // newCells
80 Foam::label Foam::cellSplitter::newOwner
81 (
82  const label faceI,
83  const Map<labelList>& cellToCells
84 ) const
85 {
86  label oldOwn = mesh_.faceOwner()[faceI];
87 
88  Map<labelList>::const_iterator fnd = cellToCells.find(oldOwn);
89 
90  if (fnd == cellToCells.end())
91  {
92  // Unsplit cell
93  return oldOwn;
94  }
95  else
96  {
97  // Look up index of face in the cells' faces.
98 
99  const labelList& newCells = fnd();
100 
101  const cell& cFaces = mesh_.cells()[oldOwn];
102 
103  return newCells[findIndex(cFaces, faceI)];
104  }
105 }
106 
107 
108 Foam::label Foam::cellSplitter::newNeighbour
109 (
110  const label faceI,
111  const Map<labelList>& cellToCells
112 ) const
113 {
114  label oldNbr = mesh_.faceNeighbour()[faceI];
115 
116  Map<labelList>::const_iterator fnd = cellToCells.find(oldNbr);
117 
118  if (fnd == cellToCells.end())
119  {
120  // Unsplit cell
121  return oldNbr;
122  }
123  else
124  {
125  // Look up index of face in the cells' faces.
126 
127  const labelList& newCells = fnd();
128 
129  const cell& cFaces = mesh_.cells()[oldNbr];
130 
131  return newCells[findIndex(cFaces, faceI)];
132  }
133 }
134 
135 
136 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
137 
138 // Construct from components
139 Foam::cellSplitter::cellSplitter(const polyMesh& mesh)
140 :
141  mesh_(mesh),
142  addedPoints_()
143 {}
144 
145 
146 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
147 
149 {}
150 
151 
152 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
153 
155 (
156  const Map<point>& cellToMidPoint,
157  polyTopoChange& meshMod
158 )
159 {
160  addedPoints_.clear();
161  addedPoints_.resize(cellToMidPoint.size());
162 
163 
164  //
165  // Introduce cellToMidPoints.
166  //
167 
168  forAllConstIter(Map<point>, cellToMidPoint, iter)
169  {
170  label cellI = iter.key();
171 
172  label anchorPoint = mesh_.cellPoints()[cellI][0];
173 
174  label addedPointI =
175  meshMod.setAction
176  (
177  polyAddPoint
178  (
179  iter(), // point
180  anchorPoint, // master point
181  -1, // zone for point
182  true // supports a cell
183  )
184  );
185  addedPoints_.insert(cellI, addedPointI);
186 
187  //Pout<< "Added point " << addedPointI
188  // << iter() << " in cell " << cellI << " with centre "
189  // << mesh_.cellCentres()[cellI] << endl;
190  }
191 
192 
193  //
194  // Add cells (first one is modified original cell)
195  //
196 
197  Map<labelList> cellToCells(cellToMidPoint.size());
198 
199  forAllConstIter(Map<point>, cellToMidPoint, iter)
200  {
201  label cellI = iter.key();
202 
203  const cell& cFaces = mesh_.cells()[cellI];
204 
205  // Cells created for this cell.
206  labelList newCells(cFaces.size());
207 
208  // First pyramid is the original cell
209  newCells[0] = cellI;
210 
211  // Add other pyramids
212  for (label i = 1; i < cFaces.size(); i++)
213  {
214  label addedCellI =
215  meshMod.setAction
216  (
217  polyAddCell
218  (
219  -1, // master point
220  -1, // master edge
221  -1, // master face
222  cellI, // master cell
223  -1 // zone
224  )
225  );
226 
227  newCells[i] = addedCellI;
228  }
229 
230  cellToCells.insert(cellI, newCells);
231 
232  //Pout<< "Split cell " << cellI
233  // << " with centre " << mesh_.cellCentres()[cellI] << nl
234  // << " faces:" << cFaces << nl
235  // << " into :" << newCells << endl;
236  }
237 
238 
239  //
240  // Introduce internal faces. These go from edges of the cell to the mid
241  // point.
242  //
243 
244  forAllConstIter(Map<point>, cellToMidPoint, iter)
245  {
246  label cellI = iter.key();
247 
248  label midPointI = addedPoints_[cellI];
249 
250  const cell& cFaces = mesh_.cells()[cellI];
251 
252  const labelList& cEdges = mesh_.cellEdges()[cellI];
253 
254  forAll(cEdges, i)
255  {
256  label edgeI = cEdges[i];
257  const edge& e = mesh_.edges()[edgeI];
258 
259  // Get the faces on the cell using the edge
260  label face0, face1;
261  meshTools::getEdgeFaces(mesh_, cellI, edgeI, face0, face1);
262 
263  // Get the cells on both sides of the face by indexing into cFaces.
264  // (since newly created cells are stored in cFaces order)
265  const labelList& newCells = cellToCells[cellI];
266 
267  label cell0 = newCells[findIndex(cFaces, face0)];
268  label cell1 = newCells[findIndex(cFaces, face1)];
269 
270  if (cell0 < cell1)
271  {
272  // Construct face to midpoint that is pointing away from
273  // (pyramid split off from) cellI
274 
275  const face& f0 = mesh_.faces()[face0];
276 
277  label index = findIndex(f0, e[0]);
278 
279  bool edgeInFaceOrder = (f0[f0.fcIndex(index)] == e[1]);
280 
281  // Check if cellI is the face owner
282 
283  face newF(3);
284  if (edgeInFaceOrder == (mesh_.faceOwner()[face0] == cellI))
285  {
286  // edge used in face order.
287  newF[0] = e[1];
288  newF[1] = e[0];
289  newF[2] = midPointI;
290  }
291  else
292  {
293  newF[0] = e[0];
294  newF[1] = e[1];
295  newF[2] = midPointI;
296  }
297 
298  // Now newF points away from cell0
299  meshMod.setAction
300  (
301  polyAddFace
302  (
303  newF, // face
304  cell0, // owner
305  cell1, // neighbour
306  -1, // master point
307  -1, // master edge
308  face0, // master face for addition
309  false, // flux flip
310  -1, // patch for face
311  -1, // zone for face
312  false // face zone flip
313  )
314  );
315  }
316  else
317  {
318  // Construct face to midpoint that is pointing away from
319  // (pyramid split off from) cellI
320 
321  const face& f1 = mesh_.faces()[face1];
322 
323  label index = findIndex(f1, e[0]);
324 
325  bool edgeInFaceOrder = (f1[f1.fcIndex(index)] == e[1]);
326 
327  // Check if cellI is the face owner
328 
329  face newF(3);
330  if (edgeInFaceOrder == (mesh_.faceOwner()[face1] == cellI))
331  {
332  // edge used in face order.
333  newF[0] = e[1];
334  newF[1] = e[0];
335  newF[2] = midPointI;
336  }
337  else
338  {
339  newF[0] = e[0];
340  newF[1] = e[1];
341  newF[2] = midPointI;
342  }
343 
344  // Now newF points away from cell1
345  meshMod.setAction
346  (
347  polyAddFace
348  (
349  newF, // face
350  cell1, // owner
351  cell0, // neighbour
352  -1, // master point
353  -1, // master edge
354  face0, // master face for addition
355  false, // flux flip
356  -1, // patch for face
357  -1, // zone for face
358  false // face zone flip
359  )
360  );
361  }
362  }
363  }
364 
365 
366  //
367  // Update all existing faces for split owner or neighbour.
368  //
369 
370 
371  // Mark off affected face.
372  boolList faceUpToDate(mesh_.nFaces(), true);
373 
374  forAllConstIter(Map<point>, cellToMidPoint, iter)
375  {
376  label cellI = iter.key();
377 
378  const cell& cFaces = mesh_.cells()[cellI];
379 
380  forAll(cFaces, i)
381  {
382  label faceI = cFaces[i];
383 
384  faceUpToDate[faceI] = false;
385  }
386  }
387 
388  forAll(faceUpToDate, faceI)
389  {
390  if (!faceUpToDate[faceI])
391  {
392  const face& f = mesh_.faces()[faceI];
393 
394  if (mesh_.isInternalFace(faceI))
395  {
396  label newOwn = newOwner(faceI, cellToCells);
397  label newNbr = newNeighbour(faceI, cellToCells);
398 
399  if (newOwn < newNbr)
400  {
401  meshMod.setAction
402  (
403  polyModifyFace
404  (
405  f,
406  faceI,
407  newOwn, // owner
408  newNbr, // neighbour
409  false, // flux flip
410  -1, // patch for face
411  false, // remove from zone
412  -1, // zone for face
413  false // face zone flip
414  )
415  );
416  }
417  else
418  {
419  meshMod.setAction
420  (
421  polyModifyFace
422  (
423  f.reverseFace(),
424  faceI,
425  newNbr, // owner
426  newOwn, // neighbour
427  false, // flux flip
428  -1, // patch for face
429  false, // remove from zone
430  -1, // zone for face
431  false // face zone flip
432  )
433  );
434  }
435 
436  }
437  else
438  {
439  label newOwn = newOwner(faceI, cellToCells);
440 
441  label patchID, zoneID, zoneFlip;
442  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
443 
444  meshMod.setAction
445  (
446  polyModifyFace
447  (
448  mesh_.faces()[faceI],
449  faceI,
450  newOwn, // owner
451  -1, // neighbour
452  false, // flux flip
453  patchID, // patch for face
454  false, // remove from zone
455  zoneID, // zone for face
456  zoneFlip // face zone flip
457  )
458  );
459  }
460 
461  faceUpToDate[faceI] = true;
462  }
463  }
464 }
465 
466 
467 void Foam::cellSplitter::updateMesh(const mapPolyMesh& morphMap)
468 {
469  // Create copy since we're deleting entries. Only if both cell and added
470  // point get mapped do they get inserted.
471  Map<label> newAddedPoints(addedPoints_.size());
472 
473  forAllConstIter(Map<label>, addedPoints_, iter)
474  {
475  label oldCellI = iter.key();
476 
477  label newCellI = morphMap.reverseCellMap()[oldCellI];
478 
479  label oldPointI = iter();
480 
481  label newPointI = morphMap.reversePointMap()[oldPointI];
482 
483  if (newCellI >= 0 && newPointI >= 0)
484  {
485  newAddedPoints.insert(newCellI, newPointI);
486  }
487  }
488 
489  // Copy
490  addedPoints_.transfer(newAddedPoints);
491 }
492 
493 
494 // ************************ vim: set sw=4 sts=4 et: ************************ //