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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "cellSplitter.H"
29#include "polyMesh.H"
30#include "polyTopoChange.H"
31#include "polyAddCell.H"
32#include "polyAddFace.H"
33#include "polyAddPoint.H"
34#include "polyModifyFace.H"
35#include "mapPolyMesh.H"
36#include "meshTools.H"
37
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
43defineTypeNameAndDebug(cellSplitter, 0);
44}
45
46
47// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48
49void Foam::cellSplitter::getFaceInfo
50(
51 const label facei,
52 label& patchID,
53 label& zoneID,
54 label& zoneFlip
55) const
56{
57 patchID = -1;
58
59 if (!mesh_.isInternalFace(facei))
60 {
61 patchID = mesh_.boundaryMesh().whichPatch(facei);
62 }
63
64 zoneID = mesh_.faceZones().whichZone(facei);
65
66 zoneFlip = false;
67
68 if (zoneID >= 0)
69 {
70 const faceZone& fZone = mesh_.faceZones()[zoneID];
71
72 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
73 }
74}
75
76
77// Find the new owner of facei (since the original cell has been split into
78// newCells
79Foam::label Foam::cellSplitter::newOwner
80(
81 const label facei,
82 const Map<labelList>& cellToCells
83) const
84{
85 const label old = mesh_.faceOwner()[facei];
86
87 const auto iter = cellToCells.cfind(old);
88
89 if (!iter.found())
90 {
91 // Unsplit cell
92 return old;
93 }
94
95
96 // Look up index of face in the cells' faces.
97
98 const labelList& newCells = *iter;
99
100 const cell& cFaces = mesh_.cells()[old];
101
102 return newCells[cFaces.find(facei)];
103}
104
105
106Foam::label Foam::cellSplitter::newNeighbour
107(
108 const label facei,
109 const Map<labelList>& cellToCells
110) const
111{
112 const label old = mesh_.faceNeighbour()[facei];
113
114 const auto iter = cellToCells.cfind(old);
115
116 if (!iter.found())
117 {
118 // Unsplit cell
119 return old;
120 }
121
122
123 // Look up index of face in the cells' faces.
124
125 const labelList& newCells = *iter;
126
127 const cell& cFaces = mesh_.cells()[old];
128
129 return newCells[cFaces.find(facei)];
130}
131
132
133// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
134
136:
137 mesh_(mesh),
138 addedPoints_()
139{}
140
141
142// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
143
145(
146 const Map<point>& cellToMidPoint,
147 polyTopoChange& meshMod
148)
149{
150 addedPoints_.clear();
151 addedPoints_.resize(cellToMidPoint.size());
152
153
154 //
155 // Introduce cellToMidPoints.
156 //
157
158 forAllConstIters(cellToMidPoint, iter)
159 {
160 const label celli = iter.key();
161
162 label anchorPoint = mesh_.cellPoints()[celli][0];
163
164 label addedPointi =
165 meshMod.setAction
166 (
167 polyAddPoint
168 (
169 iter.val(), // point
170 anchorPoint, // master point
171 -1, // zone for point
172 true // supports a cell
173 )
174 );
175 addedPoints_.insert(celli, addedPointi);
176
177 //Pout<< "Added point " << addedPointi
178 // << iter() << " in cell " << celli << " with centre "
179 // << mesh_.cellCentres()[celli] << endl;
180 }
181
182
183 //
184 // Add cells (first one is modified original cell)
185 //
186
187 Map<labelList> cellToCells(cellToMidPoint.size());
188
189 forAllConstIters(cellToMidPoint, iter)
190 {
191 const label celli = iter.key();
192
193 const cell& cFaces = mesh_.cells()[celli];
194
195 // Cells created for this cell.
196 labelList newCells(cFaces.size());
197
198 // First pyramid is the original cell
199 newCells[0] = celli;
200
201 // Add other pyramids
202 for (label i = 1; i < cFaces.size(); i++)
203 {
204 label addedCelli =
205 meshMod.setAction
206 (
207 polyAddCell
208 (
209 -1, // master point
210 -1, // master edge
211 -1, // master face
212 celli, // master cell
213 -1 // zone
214 )
215 );
216
217 newCells[i] = addedCelli;
218 }
219
220 cellToCells.insert(celli, newCells);
221
222 //Pout<< "Split cell " << celli
223 // << " with centre " << mesh_.cellCentres()[celli] << nl
224 // << " faces:" << cFaces << nl
225 // << " into :" << newCells << endl;
226 }
227
228
229 //
230 // Introduce internal faces. These go from edges of the cell to the mid
231 // point.
232 //
233
234 forAllConstIters(cellToMidPoint, iter)
235 {
236 const label celli = iter.key();
237
238 label midPointi = addedPoints_[celli];
239
240 const cell& cFaces = mesh_.cells()[celli];
241
242 const labelList& cEdges = mesh_.cellEdges()[celli];
243
244 forAll(cEdges, i)
245 {
246 label edgeI = cEdges[i];
247 const edge& e = mesh_.edges()[edgeI];
248
249 // Get the faces on the cell using the edge
250 label face0, face1;
251 meshTools::getEdgeFaces(mesh_, celli, edgeI, face0, face1);
252
253 // Get the cells on both sides of the face by indexing into cFaces.
254 // (since newly created cells are stored in cFaces order)
255 const labelList& newCells = cellToCells[celli];
256
257 label cell0 = newCells[cFaces.find(face0)];
258 label cell1 = newCells[cFaces.find(face1)];
259
260 if (cell0 < cell1)
261 {
262 // Construct face to midpoint that is pointing away from
263 // (pyramid split off from) celli
264
265 const face& f0 = mesh_.faces()[face0];
266
267 label index = f0.find(e[0]);
268
269 bool edgeInFaceOrder = (f0[f0.fcIndex(index)] == e[1]);
270
271 // Check if celli is the face owner
272
273 face newF(3);
274 if (edgeInFaceOrder == (mesh_.faceOwner()[face0] == celli))
275 {
276 // edge used in face order.
277 newF[0] = e[1];
278 newF[1] = e[0];
279 newF[2] = midPointi;
280 }
281 else
282 {
283 newF[0] = e[0];
284 newF[1] = e[1];
285 newF[2] = midPointi;
286 }
287
288 // Now newF points away from cell0
289 meshMod.setAction
290 (
291 polyAddFace
292 (
293 newF, // face
294 cell0, // owner
295 cell1, // neighbour
296 -1, // master point
297 -1, // master edge
298 face0, // master face for addition
299 false, // flux flip
300 -1, // patch for face
301 -1, // zone for face
302 false // face zone flip
303 )
304 );
305 }
306 else
307 {
308 // Construct face to midpoint that is pointing away from
309 // (pyramid split off from) celli
310
311 const face& f1 = mesh_.faces()[face1];
312
313 label index = f1.find(e[0]);
314
315 bool edgeInFaceOrder = (f1[f1.fcIndex(index)] == e[1]);
316
317 // Check if celli is the face owner
318
319 face newF(3);
320 if (edgeInFaceOrder == (mesh_.faceOwner()[face1] == celli))
321 {
322 // edge used in face order.
323 newF[0] = e[1];
324 newF[1] = e[0];
325 newF[2] = midPointi;
326 }
327 else
328 {
329 newF[0] = e[0];
330 newF[1] = e[1];
331 newF[2] = midPointi;
332 }
333
334 // Now newF points away from cell1
335 meshMod.setAction
336 (
337 polyAddFace
338 (
339 newF, // face
340 cell1, // owner
341 cell0, // neighbour
342 -1, // master point
343 -1, // master edge
344 face0, // master face for addition
345 false, // flux flip
346 -1, // patch for face
347 -1, // zone for face
348 false // face zone flip
349 )
350 );
351 }
352 }
353 }
354
355
356 //
357 // Update all existing faces for split owner or neighbour.
358 //
359
360
361 // Mark off affected face.
362 bitSet faceUpToDate(mesh_.nFaces(), true);
363
364 forAllConstIters(cellToMidPoint, iter)
365 {
366 const label celli = iter.key();
367
368 const cell& cFaces = mesh_.cells()[celli];
369
370 faceUpToDate.unset(cFaces);
371 }
372
373 forAll(faceUpToDate, facei)
374 {
375 if (!faceUpToDate.test(facei))
376 {
377 faceUpToDate.set(facei);
378
379 const face& f = mesh_.faces()[facei];
380
381 if (mesh_.isInternalFace(facei))
382 {
383 label newOwn = newOwner(facei, cellToCells);
384 label newNbr = newNeighbour(facei, cellToCells);
385
386 if (newOwn < newNbr)
387 {
388 meshMod.setAction
389 (
390 polyModifyFace
391 (
392 f,
393 facei,
394 newOwn, // owner
395 newNbr, // neighbour
396 false, // flux flip
397 -1, // patch for face
398 false, // remove from zone
399 -1, // zone for face
400 false // face zone flip
401 )
402 );
403 }
404 else
405 {
406 meshMod.setAction
407 (
408 polyModifyFace
409 (
410 f.reverseFace(),
411 facei,
412 newNbr, // owner
413 newOwn, // neighbour
414 false, // flux flip
415 -1, // patch for face
416 false, // remove from zone
417 -1, // zone for face
418 false // face zone flip
419 )
420 );
421 }
422
423 }
424 else
425 {
426 label newOwn = newOwner(facei, cellToCells);
427
428 label patchID, zoneID, zoneFlip;
429 getFaceInfo(facei, patchID, zoneID, zoneFlip);
430
431 meshMod.setAction
432 (
433 polyModifyFace
434 (
435 mesh_.faces()[facei],
436 facei,
437 newOwn, // owner
438 -1, // neighbour
439 false, // flux flip
440 patchID, // patch for face
441 false, // remove from zone
442 zoneID, // zone for face
443 zoneFlip // face zone flip
444 )
445 );
446 }
447 }
448 }
449}
450
451
452void Foam::cellSplitter::updateMesh(const mapPolyMesh& morphMap)
453{
454 // Create copy since we're deleting entries. Only if both cell and added
455 // point get mapped do they get inserted.
456 Map<label> newAddedPoints(addedPoints_.size());
457
458 forAllConstIters(addedPoints_, iter)
459 {
460 const label oldCelli = iter.key();
461 const label oldPointi = iter.val();
462
463 const label newCelli = morphMap.reverseCellMap()[oldCelli];
464 const label newPointi = morphMap.reversePointMap()[oldPointi];
465
466 if (newCelli >= 0 && newPointi >= 0)
467 {
468 newAddedPoints.insert(newCelli, newPointi);
469 }
470 }
471
472 // Copy
473 addedPoints_.transfer(newAddedPoints);
474}
475
476
477// ************************************************************************* //
void resize(const label sz)
Resize the hash table for efficiency.
Definition: HashTable.C:601
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
void clear()
Clear all entries from table.
Definition: HashTable.C:678
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type unset(const label i)
Definition: UList.H:539
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
Does pyramidal decomposition of selected cells. So all faces will become base of pyramid with as top ...
Definition: cellSplitter.H:60
void setRefinement(const Map< point > &cellToMidPoint, polyTopoChange &meshMod)
Insert mesh changes into meshMod.
void updateMesh()
Update for new mesh topology.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:498
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & cellEdges() const
const labelListList & cellPoints() const
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
const labelIOList & zoneID
void getEdgeFaces(const primitiveMesh &mesh, const label celli, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
Definition: meshTools.C:479
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278