Newer
Older
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "layerAdditionRemoval.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "polyTopoChange.H"
#include "polyTopoChanger.H"
#include "polyAddPoint.H"
#include "polyAddCell.H"
#include "polyAddFace.H"
#include "polyModifyFace.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
{
const polyMesh& mesh = topoChanger().mesh();
const primitiveFacePatch& masterFaceLayer =
mesh.faceZones()[faceZoneID_.index()]();
const pointField& points = mesh.points();
const labelList& mp = masterFaceLayer.meshPoints();
tmp<vectorField> textrusionDir(new vectorField(mp.size()));
vectorField& extrusionDir = textrusionDir.ref();
if (setLayerPairing())
{
if (debug)
{
Pout<< "void layerAdditionRemoval::extrusionDir() const "
<< " for object " << name() << " : "
<< "Using edges for point insertion" << endl;
}
// Detected a valid layer. Grab the point and face collapse mapping
const labelList& ptc = pointsPairing();
forAll(extrusionDir, mpI)
{
extrusionDir[mpI] = points[ptc[mpI]] - points[mp[mpI]];
}
}
else
{
if (debug)
{
Pout<< "void layerAdditionRemoval::extrusionDir() const "
<< " for object " << name() << " : "
<< "A valid layer could not be found in front of "
<< "the addition face layer. Using face-based "
<< "point normals for point addition"
<< endl;
}
extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
}
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
return textrusionDir;
}
void Foam::layerAdditionRemoval::addCellLayer
(
polyTopoChange& ref
) const
{
// Insert the layer addition instructions into the topological change
// Algorithm:
// 1. For every point in the master zone add a new point
// 2. For every face in the master zone add a cell
// 3. Go through all the edges of the master zone. For all internal edges,
// add a face with the correct orientation and make it internal.
// For all boundary edges, find the patch it belongs to and add the face
// to this patch
// 4. Create a set of new faces on the patch image and assign them to be
// between the old master cells and the newly created cells
// 5. Modify all the faces in the patch such that they are located
// between the old slave cells and newly created cells
if (debug)
{
Pout<< "void layerAdditionRemoval::addCellLayer("
<< "polyTopoChange& ref) const for object " << name() << " : "
<< "Adding cell layer" << endl;
}
// Create the points
const polyMesh& mesh = topoChanger().mesh();
const primitiveFacePatch& masterFaceLayer =
mesh.faceZones()[faceZoneID_.index()]();
const pointField& points = mesh.points();
const labelList& mp = masterFaceLayer.meshPoints();
// Get the extrusion direction for the added points
// Add the new points
labelList addedPoints(mp.size());
Henry Weller
committed
forAll(mp, pointi)
{
// Add the point nominal thickness away from the master zone point
// and grab the label
Henry Weller
committed
addedPoints[pointi] =
ref.setAction
(
polyAddPoint
(
Henry Weller
committed
points[mp[pointi]] // point
+ addDelta_*pointOffsets[pointi],
mp[pointi], // master point
-1, // zone for point
true // supports a cell
)
);
}
if (debug > 1)
{
Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
}
// Create the cells
const labelList& mc =
mesh.faceZones()[faceZoneID_.index()].masterCells();
const labelList& sc =
mesh.faceZones()[faceZoneID_.index()].slaveCells();
const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
labelList addedCells(mf.size());
label celli = mc[facei];
label zoneI = mesh.cellZones().whichZone(celli);
ref.setAction
(
polyAddCell
(
-1, // master point
-1, // master edge
mf[facei], // master face
zoneI // zone for cell
)
);
}
// Create the new faces
// Grab the local faces from the master zone
const faceList& zoneFaces = masterFaceLayer.localFaces();
// Create the new faces by copying the master zone. All the added
// faces need to point into the newly created cells, which means
// all the faces from the master layer are flipped. The flux flip
// is determined from the original master layer face and the face
// owner: if the master cell is equal to the face owner the flux
// remains the same; otherwise it is flipped
forAll(zoneFaces, facei)
const face oldFace = zoneFaces[facei].reverseFace();
face newFace(oldFace.size());
Henry Weller
committed
forAll(oldFace, pointi)
Henry Weller
committed
newFace[pointi] = addedPoints[oldFace[pointi]];
}
bool flipFaceFlux = false;
// Flip the face as necessary
if
(
!mesh.isInternalFace(mf[facei])
|| mc[facei] == nei[mf[facei]]
)
{
flipFaceFlux = true;
}
// Add the face
ref.setAction
(
polyAddFace
(
newFace, // face
mc[facei], // owner
addedCells[facei], // neighbour
-1, // master point
-1, // master edge
mf[facei], // master face for addition
flipFaceFlux, // flux flip
-1, // patch for face
-1, // zone for face
false // face zone flip
)
);
if (debug > 1)
{
Pout<< "adding face: " << newFace
<< " own: " << mc[facei]
<< " nei: " << addedCells[facei]
}
// Modify the faces from the master zone for the new neighbour
const faceList& faces = mesh.faces();
// Pout<< "mfFlip: " << mfFlip << endl;
const label curfaceID = mf[facei];
// If the face is internal, modify its owner to be the newly
// created cell. No flip is necessary
if (!mesh.isInternalFace(curfaceID))
{
ref.setAction
(
polyModifyFace
(
faces[curfaceID], // modified face
curfaceID, // label of face being modified
addedCells[facei], // owner
-1, // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
false, // remove from zone
faceZoneID_.index(), // zone for face
mfFlip[facei] // face flip in zone
if (debug > 1)
{
Pout<< "Modifying a boundary face. Face: " << curfaceID
<< " flip: " << mfFlip[facei]
// If slave cell is owner, the face remains the same (but with
// a new neighbour - the newly created cell). Otherwise, the
// face is flipped.
else if (sc[facei] == own[curfaceID])
{
// Orientation is good, just change neighbour
ref.setAction
(
polyModifyFace
(
faces[curfaceID], // modified face
curfaceID, // label of face being modified
own[curfaceID], // owner
addedCells[facei], // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curfaceID),// patch for face
false, // remove from zone
faceZoneID_.index(), // zone for face
mfFlip[facei] // face flip in zone
if (debug > 1)
{
Pout<< "modify face, no flip " << curfaceID
<< " own: " << own[curfaceID]
<< " nei: " << addedCells[facei]
}
else
{
// Change in orientation; flip face
ref.setAction
(
polyModifyFace
(
faces[curfaceID].reverseFace(), // modified face
curfaceID, // label of face being modified
nei[curfaceID], // owner
addedCells[facei], // neighbour
true, // face flip
mesh.boundaryMesh().whichPatch(curfaceID), // patch for face
false, // remove from zone
faceZoneID_.index(), // zone for face
!mfFlip[facei] // face flip in zone
if (debug > 1)
{
Pout<< "modify face, with flip " << curfaceID
<< " own: " << own[curfaceID]
<< " nei: " << addedCells[facei]
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
}
}
// Create new faces from edges
const edgeList& zoneLocalEdges = masterFaceLayer.edges();
const labelListList& edgeFaces = masterFaceLayer.edgeFaces();
label nInternalEdges = masterFaceLayer.nInternalEdges();
const labelList& meshEdges =
mesh.faceZones()[faceZoneID_.index()].meshEdges();
// Do all internal edges
for (label curEdgeID = 0; curEdgeID < nInternalEdges; curEdgeID++)
{
face newFace(4);
newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
ref.setAction
(
polyAddFace
(
newFace, // face
addedCells[edgeFaces[curEdgeID][0]], // owner
addedCells[edgeFaces[curEdgeID][1]], // neighbour
-1, // master point
meshEdges[curEdgeID], // master edge
-1, // master face
false, // flip flux
-1, // patch for face
-1, // zone for face
false // face zone flip
)
);
if (debug > 1)
{
Pout<< "Add internal face off edge: " << newFace
<< " own: " << addedCells[edgeFaces[curEdgeID][0]]
<< " nei: " << addedCells[edgeFaces[curEdgeID][1]]
<< endl;
}
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
}
// Prepare creation of faces from boundary edges.
// Note:
// A tricky part of the algorithm is finding the patch into which the
// newly created face will be added. For this purpose, take the edge
// and grab all the faces coming from it. From the set of faces
// eliminate the internal faces and find the boundary face from the rest.
// If there is more than one boundary face (excluding the ones in
// the master zone), the patch with the lower label is selected.
const labelListList& meshEdgeFaces = mesh.edgeFaces();
const faceZoneMesh& zoneMesh = mesh.faceZones();
// Do all boundary edges
for
(
label curEdgeID = nInternalEdges;
curEdgeID < zoneLocalEdges.size();
curEdgeID++
)
{
face newFace(4);
newFace[0] = mp[zoneLocalEdges[curEdgeID].start()];
newFace[1] = mp[zoneLocalEdges[curEdgeID].end()];
newFace[2] = addedPoints[zoneLocalEdges[curEdgeID].end()];
newFace[3] = addedPoints[zoneLocalEdges[curEdgeID].start()];
// Determine the patch for insertion
const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
label patchID = -1;
label zoneID = -1;
label extrudeFaceID = -1;
forAll(curFaces, curFacei)
extrudeFaceID = curFaces[curFacei];
if (!mesh.isInternalFace(extrudeFaceID))
// Face not internal. Check if it is not in the zone
const label cfZone = zoneMesh.whichZone(extrudeFaceID);
if (cfZone != faceZoneID_.index())
{
// Found the face in a boundary patch which is not in zone
patchID = mesh.boundaryMesh().whichPatch(extrudeFaceID);
zoneID = cfZone;
break;
}
}
}
if (patchID < 0)
{
FatalErrorInFunction
<< "Cannot find patch for edge " << meshEdges[curEdgeID]
<< ". Edge: " << mesh.edges()[meshEdges[curEdgeID]]
<< abort(FatalError);
}
ref.setAction
(
polyAddFace
(
newFace, // face
addedCells[edgeFaces[curEdgeID][0]], // owner
-1, // neighbour
-1, // master point
meshEdges[curEdgeID], // master edge
-1, // master face
false, // flip flux
patchID, // patch for face
zoneID, // zone
false // zone face flip
)
);
if (debug > 1)
{
Pout<< "add boundary face: " << newFace
<< " into patch " << patchID
<< " own: " << addedCells[edgeFaces[curEdgeID][0]]
<< endl;
}
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
// Handle duplicate boundary faces (for e.g. cyclicACMI)
if (patchID != -1)
{
for (const label cf : curFaces)
{
if (cf != extrudeFaceID && !mesh.isInternalFace(cf))
{
// Check if it is not in the zone and duplicate of
// extrudeFaceID
const label cfZone = zoneMesh.whichZone(cf);
if
(
cfZone != faceZoneID_.index()
&& face::compare(faces[cf], faces[extrudeFaceID]) == 1
)
{
// Found the face in a boundary patch which is not
// in zone (so would not be extruded above)
patchID = mesh.boundaryMesh().whichPatch(cf);
zoneID = cfZone;
ref.setAction
(
polyAddFace
(
newFace, // face
addedCells[edgeFaces[curEdgeID][0]],// owner
-1, // neighbour
-1, // point
meshEdges[curEdgeID], // edge
-1, // face
false, // flip flux
patchID, // patch
zoneID, // zone
false // zone flip
)
);
if (debug > 1)
{
Pout<< "add duplicate boundary face: " << newFace
<< " into patch " << patchID
<< " own: "
<< addedCells[edgeFaces[curEdgeID][0]] << endl;
}
}
}
}
}
}
// Modify the remaining faces of the master cells to reconnect to the new
// layer of faces.
// Algorithm: Go through all the cells of the master zone and make
// a map of faces to avoid duplicates. Do not insert the faces in
// the master patch (as they have already been dealt with). Make
// a master layer point renumbering map, which for every point in
// the master layer gives its new label. Loop through all faces in
// the map and attempt to renumber them using the master layer
// point renumbering map. Once the face is renumbered, compare it
// with the original face; if they are the same, the face has not
// changed; if not, modify the face but keep all of its old
// attributes (apart from the vertex numbers).
// Create the map of faces in the master cell layer
labelHashSet masterCellFaceMap(primitiveMesh::facesPerCell_*mc.size());
const cellList& cells = mesh.cells();
const labelList& curFaces = cells[mc[celli]];
forAll(curFaces, facei)
{
// Check if the face belongs to the master zone; if not add it
if (zoneMesh.whichZone(curFaces[facei]) != faceZoneID_.index())
masterCellFaceMap.insert(curFaces[facei]);
}
}
}
// Create the master layer point map
Map<label> masterLayerPointMap(2*mp.size());
Henry Weller
committed
forAll(mp, pointi)
{
masterLayerPointMap.insert
(
Henry Weller
committed
mp[pointi],
addedPoints[pointi]
);
}
// Grab the list of faces of the master layer
const labelList masterCellFaces = masterCellFaceMap.toc();
forAll(masterCellFaces, facei)
{
// Attempt to renumber the face using the masterLayerPointMap.
// Missing point remain the same
const label curFaceID = masterCellFaces[facei];
const face& oldFace = faces[curFaceID];
face newFace(oldFace.size());
bool changed = false;
Henry Weller
committed
forAll(oldFace, pointi)
Henry Weller
committed
if (masterLayerPointMap.found(oldFace[pointi]))
{
changed = true;
Henry Weller
committed
newFace[pointi] = masterLayerPointMap.find(oldFace[pointi])();
Henry Weller
committed
newFace[pointi] = oldFace[pointi];
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
}
}
// If the face has changed, create a modification entry
if (changed)
{
// Get face zone and its flip
label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
bool modifiedFaceZoneFlip = false;
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
mesh.faceZones()[modifiedFaceZone].flipMap()
[
mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
];
}
if (mesh.isInternalFace(curFaceID))
{
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
nei[curFaceID], // neighbour
false, // face flip
-1, // patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
if (debug > 1)
{
Pout<< "modifying stick-out face. Internal Old face: "
<< oldFace
<< " new face: " << newFace
<< " own: " << own[curFaceID]
<< " nei: " << nei[curFaceID]
<< endl;
}
}
else
{
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
-1, // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID),
// patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
if (debug > 1)
{
Pout<< "modifying stick-out face. Boundary Old face: "
<< oldFace
<< " new face: " << newFace
<< " own: " << own[curFaceID]
<< " patch: "
<< mesh.boundaryMesh().whichPatch(curFaceID)
<< endl;
}
}
}
}
if (debug)
{
Pout<< "void layerAdditionRemoval::addCellLayer(polyTopoChange&) const "
<< " for object " << name() << ": "
<< "Finished adding cell layer" << endl;
}
}
// ************************************************************************* //