polyTopoChange: Removed used and deprecated functionality
This commit is contained in:
Henry Weller
@ -1,3 +0,0 @@
|
||||
collapseEdges.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/collapseEdges
|
||||
@ -1,10 +0,0 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-ldynamicMesh \
|
||||
-lmeshTools \
|
||||
-lgenericPatches \
|
||||
-lfiniteVolume
|
||||
@ -1,258 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Application
|
||||
collapseEdges
|
||||
|
||||
Description
|
||||
Collapses short edges and combines edges that are in line.
|
||||
|
||||
- collapse short edges. Length of edges to collapse provided as argument.
|
||||
- merge two edges if they are in line. Maximum angle provided as argument.
|
||||
- remove unused points.
|
||||
- collapse faces:
|
||||
- with small areas to a single point
|
||||
- that have a high aspect ratio (i.e. sliver face) to a single edge
|
||||
|
||||
Optionally checks the resulting mesh for bad faces and reduces the desired
|
||||
face length factor for those faces attached to the bad faces.
|
||||
|
||||
When collapsing an edge with one point on the boundary it will leave
|
||||
the boundary point intact. When both points inside it chooses random. When
|
||||
both points on boundary random again.
|
||||
|
||||
Usage
|
||||
- collapseEdges [OPTION]
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "argList.H"
|
||||
#include "Time.H"
|
||||
#include "timeSelector.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "fvMesh.H"
|
||||
#include "polyMeshFilter.H"
|
||||
#include "faceSet.H"
|
||||
#include "systemDict.H"
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
timeSelector::addOptions(true, false);
|
||||
argList::addNote
|
||||
(
|
||||
"Collapses small edges to a point.\n"
|
||||
"Optionally collapse small faces to a point and thin faces to an edge."
|
||||
);
|
||||
|
||||
argList::addBoolOption
|
||||
(
|
||||
"collapseFaces",
|
||||
"Collapse small and sliver faces as well as small edges"
|
||||
);
|
||||
|
||||
argList::addOption
|
||||
(
|
||||
"collapseFaceSet",
|
||||
"faceSet",
|
||||
"Collapse faces that are in the supplied face set"
|
||||
);
|
||||
|
||||
#include "addDictOption.H"
|
||||
#include "addOverwriteOption.H"
|
||||
#include "setRootCase.H"
|
||||
#include "createTimeNoFunctionObjects.H"
|
||||
const instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
|
||||
|
||||
#include "createMeshNoChangers.H"
|
||||
|
||||
const word oldInstance = mesh.pointsInstance();
|
||||
|
||||
const dictionary collapseDict(systemDict("collapseDict", args, mesh));
|
||||
|
||||
const bool overwrite = args.optionFound("overwrite");
|
||||
|
||||
const bool collapseFaces = args.optionFound("collapseFaces");
|
||||
const bool collapseFaceSet = args.optionFound("collapseFaceSet");
|
||||
|
||||
if (collapseFaces && collapseFaceSet)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Both face zone collapsing and face collapsing have been"
|
||||
<< "selected. Choose only one of:" << nl
|
||||
<< " -collapseFaces" << nl
|
||||
<< " -collapseFaceSet <faceSet>"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// maintain indirectPatchFaces if it is there (default) or force
|
||||
// (if collapseFaceSet option provided)
|
||||
word faceSetName("indirectPatchFaces");
|
||||
IOobject::readOption readFlag = IOobject::READ_IF_PRESENT;
|
||||
|
||||
if (args.optionReadIfPresent("collapseFaceSet", faceSetName))
|
||||
{
|
||||
readFlag = IOobject::MUST_READ;
|
||||
}
|
||||
|
||||
|
||||
|
||||
labelIOList pointPriority
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"pointPriority",
|
||||
runTime.name(),
|
||||
runTime,
|
||||
IOobject::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
labelList(mesh.nPoints(), labelMin)
|
||||
);
|
||||
forAll(timeDirs, timeI)
|
||||
{
|
||||
runTime.setTime(timeDirs[timeI], timeI);
|
||||
|
||||
Info<< "Time = " << runTime.userTimeName() << endl;
|
||||
|
||||
autoPtr<polyMeshFilter> meshFilterPtr;
|
||||
|
||||
label nBadFaces = 0;
|
||||
|
||||
faceSet indirectPatchFaces
|
||||
(
|
||||
mesh,
|
||||
faceSetName,
|
||||
readFlag,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
Info<< "Read faceSet " << indirectPatchFaces.name()
|
||||
<< " with "
|
||||
<< returnReduce(indirectPatchFaces.size(), sumOp<label>())
|
||||
<< " faces" << endl;
|
||||
|
||||
|
||||
{
|
||||
meshFilterPtr.set
|
||||
(
|
||||
new polyMeshFilter(mesh, pointPriority, collapseDict)
|
||||
);
|
||||
polyMeshFilter& meshFilter = meshFilterPtr();
|
||||
|
||||
// newMesh will be empty until it is filtered
|
||||
const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
|
||||
|
||||
// Filter small edges only. This reduces the number of faces so that
|
||||
// the face filtering is sped up.
|
||||
nBadFaces = meshFilter.filterEdges(0);
|
||||
{
|
||||
polyTopoChange meshMod(newMesh());
|
||||
|
||||
meshMod.changeMesh(mesh, false);
|
||||
|
||||
polyMeshFilter::copySets(newMesh(), mesh);
|
||||
}
|
||||
|
||||
pointPriority = meshFilter.pointPriority();
|
||||
}
|
||||
|
||||
if (collapseFaceSet)
|
||||
{
|
||||
meshFilterPtr.reset
|
||||
(
|
||||
new polyMeshFilter(mesh, pointPriority, collapseDict)
|
||||
);
|
||||
polyMeshFilter& meshFilter = meshFilterPtr();
|
||||
|
||||
const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
|
||||
|
||||
// Filter faces. Pass in the number of bad faces that are present
|
||||
// from the previous edge filtering to use as a stopping criterion.
|
||||
meshFilter.filter(indirectPatchFaces);
|
||||
{
|
||||
polyTopoChange meshMod(newMesh);
|
||||
|
||||
meshMod.changeMesh(mesh, false);
|
||||
|
||||
polyMeshFilter::copySets(newMesh(), mesh);
|
||||
}
|
||||
|
||||
pointPriority = meshFilter.pointPriority();
|
||||
}
|
||||
|
||||
if (collapseFaces)
|
||||
{
|
||||
meshFilterPtr.reset
|
||||
(
|
||||
new polyMeshFilter(mesh, pointPriority, collapseDict)
|
||||
);
|
||||
polyMeshFilter& meshFilter = meshFilterPtr();
|
||||
|
||||
const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
|
||||
|
||||
// Filter faces. Pass in the number of bad faces that are present
|
||||
// from the previous edge filtering to use as a stopping criterion.
|
||||
meshFilter.filter(nBadFaces);
|
||||
{
|
||||
polyTopoChange meshMod(newMesh);
|
||||
|
||||
meshMod.changeMesh(mesh, false);
|
||||
|
||||
polyMeshFilter::copySets(newMesh(), mesh);
|
||||
}
|
||||
|
||||
pointPriority = meshFilter.pointPriority();
|
||||
}
|
||||
|
||||
// Write resulting mesh
|
||||
if (!overwrite)
|
||||
{
|
||||
runTime++;
|
||||
}
|
||||
else
|
||||
{
|
||||
mesh.setInstance(oldInstance);
|
||||
}
|
||||
|
||||
Info<< nl << "Writing collapsed mesh to time "
|
||||
<< runTime.name() << nl << endl;
|
||||
|
||||
mesh.write();
|
||||
pointPriority.write();
|
||||
}
|
||||
|
||||
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,4 +0,0 @@
|
||||
cellSplitter.C
|
||||
modifyMesh.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/modifyMesh
|
||||
@ -1,8 +0,0 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lmeshTools \
|
||||
-lgenericPatches \
|
||||
-ldynamicMesh
|
||||
@ -1,491 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "cellSplitter.H"
|
||||
#include "polyMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "polyAddCell.H"
|
||||
#include "polyAddFace.H"
|
||||
#include "polyAddPoint.H"
|
||||
#include "polyModifyFace.H"
|
||||
#include "polyTopoChangeMap.H"
|
||||
#include "meshTools.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(cellSplitter, 0);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::cellSplitter::getFaceInfo
|
||||
(
|
||||
const label facei,
|
||||
label& patchID,
|
||||
label& zoneID,
|
||||
label& zoneFlip
|
||||
) const
|
||||
{
|
||||
patchID = -1;
|
||||
|
||||
if (!mesh_.isInternalFace(facei))
|
||||
{
|
||||
patchID = mesh_.boundaryMesh().whichPatch(facei);
|
||||
}
|
||||
|
||||
zoneID = mesh_.faceZones().whichZone(facei);
|
||||
|
||||
zoneFlip = false;
|
||||
|
||||
if (zoneID >= 0)
|
||||
{
|
||||
const faceZone& fZone = mesh_.faceZones()[zoneID];
|
||||
|
||||
zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Find the new owner of facei (since the original cell has been split into
|
||||
// newCells
|
||||
Foam::label Foam::cellSplitter::newOwner
|
||||
(
|
||||
const label facei,
|
||||
const Map<labelList>& cellToCells
|
||||
) const
|
||||
{
|
||||
label oldOwn = mesh_.faceOwner()[facei];
|
||||
|
||||
Map<labelList>::const_iterator fnd = cellToCells.find(oldOwn);
|
||||
|
||||
if (fnd == cellToCells.end())
|
||||
{
|
||||
// Unsplit cell
|
||||
return oldOwn;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Look up index of face in the cells' faces.
|
||||
|
||||
const labelList& newCells = fnd();
|
||||
|
||||
const cell& cFaces = mesh_.cells()[oldOwn];
|
||||
|
||||
return newCells[findIndex(cFaces, facei)];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::label Foam::cellSplitter::newNeighbour
|
||||
(
|
||||
const label facei,
|
||||
const Map<labelList>& cellToCells
|
||||
) const
|
||||
{
|
||||
label oldNbr = mesh_.faceNeighbour()[facei];
|
||||
|
||||
Map<labelList>::const_iterator fnd = cellToCells.find(oldNbr);
|
||||
|
||||
if (fnd == cellToCells.end())
|
||||
{
|
||||
// Unsplit cell
|
||||
return oldNbr;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Look up index of face in the cells' faces.
|
||||
|
||||
const labelList& newCells = fnd();
|
||||
|
||||
const cell& cFaces = mesh_.cells()[oldNbr];
|
||||
|
||||
return newCells[findIndex(cFaces, facei)];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
// Construct from components
|
||||
Foam::cellSplitter::cellSplitter(const polyMesh& mesh)
|
||||
:
|
||||
mesh_(mesh),
|
||||
addedPoints_()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::cellSplitter::~cellSplitter()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::cellSplitter::setRefinement
|
||||
(
|
||||
const Map<point>& cellToMidPoint,
|
||||
polyTopoChange& meshMod
|
||||
)
|
||||
{
|
||||
addedPoints_.clear();
|
||||
addedPoints_.resize(cellToMidPoint.size());
|
||||
|
||||
|
||||
//
|
||||
// Introduce cellToMidPoints.
|
||||
//
|
||||
|
||||
forAllConstIter(Map<point>, cellToMidPoint, iter)
|
||||
{
|
||||
label celli = iter.key();
|
||||
|
||||
label anchorPoint = mesh_.cellPoints()[celli][0];
|
||||
|
||||
label addedPointi =
|
||||
meshMod.setAction
|
||||
(
|
||||
polyAddPoint
|
||||
(
|
||||
iter(), // point
|
||||
anchorPoint, // master point
|
||||
-1, // zone for point
|
||||
true // supports a cell
|
||||
)
|
||||
);
|
||||
addedPoints_.insert(celli, addedPointi);
|
||||
|
||||
// Pout<< "Added point " << addedPointi
|
||||
// << iter() << " in cell " << celli << " with centre "
|
||||
// << mesh_.cellCentres()[celli] << endl;
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
// Add cells (first one is modified original cell)
|
||||
//
|
||||
|
||||
Map<labelList> cellToCells(cellToMidPoint.size());
|
||||
|
||||
forAllConstIter(Map<point>, cellToMidPoint, iter)
|
||||
{
|
||||
label celli = iter.key();
|
||||
|
||||
const cell& cFaces = mesh_.cells()[celli];
|
||||
|
||||
// Cells created for this cell.
|
||||
labelList newCells(cFaces.size());
|
||||
|
||||
// First pyramid is the original cell
|
||||
newCells[0] = celli;
|
||||
|
||||
// Add other pyramids
|
||||
for (label i = 1; i < cFaces.size(); i++)
|
||||
{
|
||||
label addedCelli =
|
||||
meshMod.setAction
|
||||
(
|
||||
polyAddCell
|
||||
(
|
||||
-1, // master point
|
||||
-1, // master edge
|
||||
-1, // master face
|
||||
celli, // master cell
|
||||
-1 // zone
|
||||
)
|
||||
);
|
||||
|
||||
newCells[i] = addedCelli;
|
||||
}
|
||||
|
||||
cellToCells.insert(celli, newCells);
|
||||
|
||||
// Pout<< "Split cell " << celli
|
||||
// << " with centre " << mesh_.cellCentres()[celli] << nl
|
||||
// << " faces:" << cFaces << nl
|
||||
// << " into :" << newCells << endl;
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
// Introduce internal faces. These go from edges of the cell to the mid
|
||||
// point.
|
||||
//
|
||||
|
||||
forAllConstIter(Map<point>, cellToMidPoint, iter)
|
||||
{
|
||||
label celli = iter.key();
|
||||
|
||||
label midPointi = addedPoints_[celli];
|
||||
|
||||
const cell& cFaces = mesh_.cells()[celli];
|
||||
|
||||
const labelList& cEdges = mesh_.cellEdges()[celli];
|
||||
|
||||
forAll(cEdges, i)
|
||||
{
|
||||
label edgeI = cEdges[i];
|
||||
const edge& e = mesh_.edges()[edgeI];
|
||||
|
||||
// Get the faces on the cell using the edge
|
||||
label face0, face1;
|
||||
meshTools::getEdgeFaces(mesh_, celli, edgeI, face0, face1);
|
||||
|
||||
// Get the cells on both sides of the face by indexing into cFaces.
|
||||
// (since newly created cells are stored in cFaces order)
|
||||
const labelList& newCells = cellToCells[celli];
|
||||
|
||||
label cell0 = newCells[findIndex(cFaces, face0)];
|
||||
label cell1 = newCells[findIndex(cFaces, face1)];
|
||||
|
||||
if (cell0 < cell1)
|
||||
{
|
||||
// Construct face to midpoint that is pointing away from
|
||||
// (pyramid split off from) celli
|
||||
|
||||
const face& f0 = mesh_.faces()[face0];
|
||||
|
||||
label index = findIndex(f0, e[0]);
|
||||
|
||||
bool edgeInFaceOrder = (f0[f0.fcIndex(index)] == e[1]);
|
||||
|
||||
// Check if celli is the face owner
|
||||
|
||||
face newF(3);
|
||||
if (edgeInFaceOrder == (mesh_.faceOwner()[face0] == celli))
|
||||
{
|
||||
// edge used in face order.
|
||||
newF[0] = e[1];
|
||||
newF[1] = e[0];
|
||||
newF[2] = midPointi;
|
||||
}
|
||||
else
|
||||
{
|
||||
newF[0] = e[0];
|
||||
newF[1] = e[1];
|
||||
newF[2] = midPointi;
|
||||
}
|
||||
|
||||
// Now newF points away from cell0
|
||||
meshMod.setAction
|
||||
(
|
||||
polyAddFace
|
||||
(
|
||||
newF, // face
|
||||
cell0, // owner
|
||||
cell1, // neighbour
|
||||
-1, // master point
|
||||
-1, // master edge
|
||||
face0, // master face for addition
|
||||
false, // flux flip
|
||||
-1, // patch for face
|
||||
-1, // zone for face
|
||||
false // face zone flip
|
||||
)
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Construct face to midpoint that is pointing away from
|
||||
// (pyramid split off from) celli
|
||||
|
||||
const face& f1 = mesh_.faces()[face1];
|
||||
|
||||
label index = findIndex(f1, e[0]);
|
||||
|
||||
bool edgeInFaceOrder = (f1[f1.fcIndex(index)] == e[1]);
|
||||
|
||||
// Check if celli is the face owner
|
||||
|
||||
face newF(3);
|
||||
if (edgeInFaceOrder == (mesh_.faceOwner()[face1] == celli))
|
||||
{
|
||||
// edge used in face order.
|
||||
newF[0] = e[1];
|
||||
newF[1] = e[0];
|
||||
newF[2] = midPointi;
|
||||
}
|
||||
else
|
||||
{
|
||||
newF[0] = e[0];
|
||||
newF[1] = e[1];
|
||||
newF[2] = midPointi;
|
||||
}
|
||||
|
||||
// Now newF points away from cell1
|
||||
meshMod.setAction
|
||||
(
|
||||
polyAddFace
|
||||
(
|
||||
newF, // face
|
||||
cell1, // owner
|
||||
cell0, // neighbour
|
||||
-1, // master point
|
||||
-1, // master edge
|
||||
face0, // master face for addition
|
||||
false, // flux flip
|
||||
-1, // patch for face
|
||||
-1, // zone for face
|
||||
false // face zone flip
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
// Update all existing faces for split owner or neighbour.
|
||||
//
|
||||
|
||||
|
||||
// Mark off affected face.
|
||||
boolList faceUpToDate(mesh_.nFaces(), true);
|
||||
|
||||
forAllConstIter(Map<point>, cellToMidPoint, iter)
|
||||
{
|
||||
label celli = iter.key();
|
||||
|
||||
const cell& cFaces = mesh_.cells()[celli];
|
||||
|
||||
forAll(cFaces, i)
|
||||
{
|
||||
label facei = cFaces[i];
|
||||
|
||||
faceUpToDate[facei] = false;
|
||||
}
|
||||
}
|
||||
|
||||
forAll(faceUpToDate, facei)
|
||||
{
|
||||
if (!faceUpToDate[facei])
|
||||
{
|
||||
const face& f = mesh_.faces()[facei];
|
||||
|
||||
if (mesh_.isInternalFace(facei))
|
||||
{
|
||||
label newOwn = newOwner(facei, cellToCells);
|
||||
label newNbr = newNeighbour(facei, cellToCells);
|
||||
|
||||
if (newOwn < newNbr)
|
||||
{
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
f,
|
||||
facei,
|
||||
newOwn, // owner
|
||||
newNbr, // neighbour
|
||||
false, // flux flip
|
||||
-1, // patch for face
|
||||
false, // remove from zone
|
||||
-1, // zone for face
|
||||
false // face zone flip
|
||||
)
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
f.reverseFace(),
|
||||
facei,
|
||||
newNbr, // owner
|
||||
newOwn, // neighbour
|
||||
false, // flux flip
|
||||
-1, // patch for face
|
||||
false, // remove from zone
|
||||
-1, // zone for face
|
||||
false // face zone flip
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
}
|
||||
else
|
||||
{
|
||||
label newOwn = newOwner(facei, cellToCells);
|
||||
|
||||
label patchID, zoneID, zoneFlip;
|
||||
getFaceInfo(facei, patchID, zoneID, zoneFlip);
|
||||
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
mesh_.faces()[facei],
|
||||
facei,
|
||||
newOwn, // owner
|
||||
-1, // neighbour
|
||||
false, // flux flip
|
||||
patchID, // patch for face
|
||||
false, // remove from zone
|
||||
zoneID, // zone for face
|
||||
zoneFlip // face zone flip
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
faceUpToDate[facei] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::cellSplitter::topoChange(const polyTopoChangeMap& map)
|
||||
{
|
||||
// Create copy since we're deleting entries. Only if both cell and added
|
||||
// point get mapped do they get inserted.
|
||||
Map<label> newAddedPoints(addedPoints_.size());
|
||||
|
||||
forAllConstIter(Map<label>, addedPoints_, iter)
|
||||
{
|
||||
label oldCelli = iter.key();
|
||||
|
||||
label newCelli = map.reverseCellMap()[oldCelli];
|
||||
|
||||
label oldPointi = iter();
|
||||
|
||||
label newPointi = map.reversePointMap()[oldPointi];
|
||||
|
||||
if (newCelli >= 0 && newPointi >= 0)
|
||||
{
|
||||
newAddedPoints.insert(newCelli, newPointi);
|
||||
}
|
||||
}
|
||||
|
||||
// Copy
|
||||
addedPoints_.transfer(newAddedPoints);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,153 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Class
|
||||
Foam::cellSplitter
|
||||
|
||||
Description
|
||||
Does pyramidal decomposition of selected cells. So all faces will become
|
||||
base of pyramid with as top a user-supplied point (usually the cell centre)
|
||||
|
||||
SourceFiles
|
||||
cellSplitter.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef cellSplitter_H
|
||||
#define cellSplitter_H
|
||||
|
||||
#include "Map.H"
|
||||
#include "edge.H"
|
||||
#include "typeInfo.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// Forward declaration of classes
|
||||
class polyTopoChange;
|
||||
class polyTopoChangeMap;
|
||||
class polyMesh;
|
||||
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class cellSplitter Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class cellSplitter
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Reference to mesh
|
||||
const polyMesh& mesh_;
|
||||
|
||||
//- Per cell the mid point added.
|
||||
Map<label> addedPoints_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Get patch and zone info for face
|
||||
void getFaceInfo
|
||||
(
|
||||
const label facei,
|
||||
label& patchID,
|
||||
label& zoneID,
|
||||
label& zoneFlip
|
||||
) const;
|
||||
|
||||
//- Find the new owner (if any) of the face.
|
||||
label newOwner
|
||||
(
|
||||
const label facei,
|
||||
const Map<labelList>& cellToCells
|
||||
) const;
|
||||
|
||||
//- Find the new neighbour (if any) of the face.
|
||||
label newNeighbour
|
||||
(
|
||||
const label facei,
|
||||
const Map<labelList>& cellToCells
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
ClassName("cellSplitter");
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from mesh
|
||||
cellSplitter(const polyMesh& mesh);
|
||||
|
||||
//- Disallow default bitwise copy construction
|
||||
cellSplitter(const cellSplitter&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
~cellSplitter();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Edit
|
||||
|
||||
//- Insert mesh changes into meshMod.
|
||||
// cellToMidPoint : cell to cut and position of its new midpoint
|
||||
void setRefinement
|
||||
(
|
||||
const Map<point>& cellToMidPoint,
|
||||
polyTopoChange& meshMod
|
||||
);
|
||||
|
||||
//- Force recalculation of locally stored data on topological change
|
||||
void topoChange(const polyTopoChangeMap&);
|
||||
|
||||
|
||||
// Access
|
||||
|
||||
//- Per cell the mid point added.
|
||||
const Map<label>& addedPoints() const
|
||||
{
|
||||
return addedPoints_;
|
||||
}
|
||||
|
||||
|
||||
// Member Operators
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const cellSplitter&) = delete;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,678 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Application
|
||||
modifyMesh
|
||||
|
||||
Description
|
||||
Manipulates mesh elements.
|
||||
|
||||
Actions are:
|
||||
(boundary)points:
|
||||
- move
|
||||
|
||||
(boundary)edges:
|
||||
- split and move introduced point
|
||||
|
||||
(boundary)faces:
|
||||
- split(triangulate) and move introduced point
|
||||
|
||||
edges:
|
||||
- collapse
|
||||
|
||||
cells:
|
||||
- split into polygonal base pyramids around newly introduced mid
|
||||
point
|
||||
|
||||
Is a bit of a loose collection of mesh change drivers.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "argList.H"
|
||||
#include "Time.H"
|
||||
#include "polyMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "polyTopoChangeMap.H"
|
||||
#include "boundaryCutter.H"
|
||||
#include "cellSplitter.H"
|
||||
#include "edgeCollapser.H"
|
||||
#include "meshTools.H"
|
||||
#include "Pair.H"
|
||||
#include "globalIndex.H"
|
||||
#include "systemDict.H"
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
// Locate point on patch. Returns (mesh) point label.
|
||||
label findPoint(const primitivePatch& pp, const point& nearPoint)
|
||||
{
|
||||
const pointField& points = pp.points();
|
||||
const labelList& meshPoints = pp.meshPoints();
|
||||
|
||||
// Find nearest and next nearest
|
||||
scalar minDistSqr = great;
|
||||
label minI = -1;
|
||||
|
||||
scalar almostMinDistSqr = great;
|
||||
label almostMinI = -1;
|
||||
|
||||
forAll(meshPoints, i)
|
||||
{
|
||||
label pointi = meshPoints[i];
|
||||
|
||||
scalar distSqr = magSqr(nearPoint - points[pointi]);
|
||||
|
||||
if (distSqr < minDistSqr)
|
||||
{
|
||||
almostMinDistSqr = minDistSqr;
|
||||
almostMinI = minI;
|
||||
|
||||
minDistSqr = distSqr;
|
||||
minI = pointi;
|
||||
}
|
||||
else if (distSqr < almostMinDistSqr)
|
||||
{
|
||||
almostMinDistSqr = distSqr;
|
||||
almostMinI = pointi;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Decide if nearPoint unique enough.
|
||||
Info<< "Found to point " << nearPoint << nl
|
||||
<< " nearest point : " << minI
|
||||
<< " distance " << Foam::sqrt(minDistSqr)
|
||||
<< " at " << points[minI] << nl
|
||||
<< " next nearest point : " << almostMinI
|
||||
<< " distance " << Foam::sqrt(almostMinDistSqr)
|
||||
<< " at " << points[almostMinI] << endl;
|
||||
|
||||
if (almostMinDistSqr < 4*minDistSqr)
|
||||
{
|
||||
Info<< "Next nearest too close to nearest. Aborting" << endl;
|
||||
|
||||
return -1;
|
||||
}
|
||||
else
|
||||
{
|
||||
return minI;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Locate edge on patch. Return mesh edge label.
|
||||
label findEdge
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const primitivePatch& pp,
|
||||
const point& nearPoint
|
||||
)
|
||||
{
|
||||
const pointField& localPoints = pp.localPoints();
|
||||
const pointField& points = pp.points();
|
||||
const labelList& meshPoints = pp.meshPoints();
|
||||
const edgeList& edges = pp.edges();
|
||||
|
||||
// Find nearest and next nearest
|
||||
scalar minDist = great;
|
||||
label minI = -1;
|
||||
|
||||
scalar almostMinDist = great;
|
||||
label almostMinI = -1;
|
||||
|
||||
forAll(edges, edgeI)
|
||||
{
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
pointHit pHit(e.line(localPoints).nearestDist(nearPoint));
|
||||
|
||||
if (pHit.hit())
|
||||
{
|
||||
if (pHit.distance() < minDist)
|
||||
{
|
||||
almostMinDist = minDist;
|
||||
almostMinI = minI;
|
||||
|
||||
minDist = pHit.distance();
|
||||
minI = meshTools::findEdge
|
||||
(
|
||||
mesh,
|
||||
meshPoints[e[0]],
|
||||
meshPoints[e[1]]
|
||||
);
|
||||
}
|
||||
else if (pHit.distance() < almostMinDist)
|
||||
{
|
||||
almostMinDist = pHit.distance();
|
||||
almostMinI = meshTools::findEdge
|
||||
(
|
||||
mesh,
|
||||
meshPoints[e[0]],
|
||||
meshPoints[e[1]]
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (minI == -1)
|
||||
{
|
||||
Info<< "Did not find edge close to point " << nearPoint << endl;
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
// Decide if nearPoint unique enough.
|
||||
Info<< "Found to point " << nearPoint << nl
|
||||
<< " nearest edge : " << minI
|
||||
<< " distance " << minDist << " endpoints "
|
||||
<< mesh.edges()[minI].line(points) << nl
|
||||
<< " next nearest edge : " << almostMinI
|
||||
<< " distance " << almostMinDist << " endpoints "
|
||||
<< mesh.edges()[almostMinI].line(points) << nl
|
||||
<< endl;
|
||||
|
||||
if (almostMinDist < 2*minDist)
|
||||
{
|
||||
Info<< "Next nearest too close to nearest. Aborting" << endl;
|
||||
|
||||
return -1;
|
||||
}
|
||||
else
|
||||
{
|
||||
return minI;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Find face on patch. Return mesh face label.
|
||||
label findFace
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const primitivePatch& pp,
|
||||
const point& nearPoint
|
||||
)
|
||||
{
|
||||
const pointField& points = pp.points();
|
||||
|
||||
// Find nearest and next nearest
|
||||
scalar minDist = great;
|
||||
label minI = -1;
|
||||
|
||||
scalar almostMinDist = great;
|
||||
label almostMinI = -1;
|
||||
|
||||
forAll(pp, patchFacei)
|
||||
{
|
||||
pointHit pHit(pp[patchFacei].nearestPoint(nearPoint, points));
|
||||
|
||||
if (pHit.hit())
|
||||
{
|
||||
if (pHit.distance() < minDist)
|
||||
{
|
||||
almostMinDist = minDist;
|
||||
almostMinI = minI;
|
||||
|
||||
minDist = pHit.distance();
|
||||
minI = patchFacei + mesh.nInternalFaces();
|
||||
}
|
||||
else if (pHit.distance() < almostMinDist)
|
||||
{
|
||||
almostMinDist = pHit.distance();
|
||||
almostMinI = patchFacei + mesh.nInternalFaces();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (minI == -1)
|
||||
{
|
||||
Info<< "Did not find face close to point " << nearPoint << endl;
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
// Decide if nearPoint unique enough.
|
||||
Info<< "Found to point " << nearPoint << nl
|
||||
<< " nearest face : " << minI
|
||||
<< " distance " << minDist
|
||||
<< " to face centre " << mesh.faceCentres()[minI] << nl
|
||||
<< " next nearest face : " << almostMinI
|
||||
<< " distance " << almostMinDist
|
||||
<< " to face centre " << mesh.faceCentres()[almostMinI] << nl
|
||||
<< endl;
|
||||
|
||||
if (almostMinDist < 2*minDist)
|
||||
{
|
||||
Info<< "Next nearest too close to nearest. Aborting" << endl;
|
||||
|
||||
return -1;
|
||||
}
|
||||
else
|
||||
{
|
||||
return minI;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Find cell with cell centre close to given point.
|
||||
label findCell(const primitiveMesh& mesh, const point& nearPoint)
|
||||
{
|
||||
label celli = mesh.findCell(nearPoint);
|
||||
|
||||
if (celli != -1)
|
||||
{
|
||||
scalar distToCcSqr = magSqr(nearPoint - mesh.cellCentres()[celli]);
|
||||
|
||||
const labelList& cPoints = mesh.cellPoints()[celli];
|
||||
|
||||
label minI = -1;
|
||||
scalar minDistSqr = great;
|
||||
|
||||
forAll(cPoints, i)
|
||||
{
|
||||
label pointi = cPoints[i];
|
||||
|
||||
scalar distSqr = magSqr(nearPoint - mesh.points()[pointi]);
|
||||
|
||||
if (distSqr < minDistSqr)
|
||||
{
|
||||
minDistSqr = distSqr;
|
||||
minI = pointi;
|
||||
}
|
||||
}
|
||||
|
||||
// Decide if nearPoint unique enough.
|
||||
Info<< "Found to point " << nearPoint << nl
|
||||
<< " nearest cell : " << celli
|
||||
<< " distance " << Foam::sqrt(distToCcSqr)
|
||||
<< " to cell centre " << mesh.cellCentres()[celli] << nl
|
||||
<< " nearest mesh point : " << minI
|
||||
<< " distance " << Foam::sqrt(minDistSqr)
|
||||
<< " to " << mesh.points()[minI] << nl
|
||||
<< endl;
|
||||
|
||||
if (minDistSqr < 4*distToCcSqr)
|
||||
{
|
||||
Info<< "Mesh point too close to nearest cell centre. Aborting"
|
||||
<< endl;
|
||||
|
||||
celli = -1;
|
||||
}
|
||||
}
|
||||
|
||||
return celli;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
#include "addOverwriteOption.H"
|
||||
|
||||
#include "setRootCase.H"
|
||||
#include "createTimeNoFunctionObjects.H"
|
||||
#include "createPolyMesh.H"
|
||||
const word oldInstance = mesh.pointsInstance();
|
||||
|
||||
const bool overwrite = args.optionFound("overwrite");
|
||||
|
||||
const dictionary dict(systemDict("modifyMeshDict", args, mesh));
|
||||
|
||||
// Read all from the dictionary.
|
||||
List<Pair<point>> pointsToMove(dict.lookup("pointsToMove"));
|
||||
List<Pair<point>> edgesToSplit(dict.lookup("edgesToSplit"));
|
||||
List<Pair<point>> facesToTriangulate
|
||||
(
|
||||
dict.lookup("facesToTriangulate")
|
||||
);
|
||||
|
||||
bool cutBoundary =
|
||||
(
|
||||
pointsToMove.size()
|
||||
|| edgesToSplit.size()
|
||||
|| facesToTriangulate.size()
|
||||
);
|
||||
|
||||
List<Pair<point>> edgesToCollapse(dict.lookup("edgesToCollapse"));
|
||||
|
||||
bool collapseEdge = edgesToCollapse.size();
|
||||
|
||||
List<Pair<point>> cellsToPyramidise(dict.lookup("cellsToSplit"));
|
||||
|
||||
bool cellsToSplit = cellsToPyramidise.size();
|
||||
|
||||
// List<Tuple2<pointField,point>>
|
||||
// cellsToCreate(dict.lookup("cellsToCreate"));
|
||||
|
||||
Info<< "Read from " << dict.name() << nl
|
||||
<< " Boundary cutting module:" << nl
|
||||
<< " points to move :" << pointsToMove.size() << nl
|
||||
<< " edges to split :" << edgesToSplit.size() << nl
|
||||
<< " faces to triangulate:" << facesToTriangulate.size() << nl
|
||||
<< " Cell splitting module:" << nl
|
||||
<< " cells to split :" << cellsToPyramidise.size() << nl
|
||||
<< " Edge collapsing module:" << nl
|
||||
<< " edges to collapse :" << edgesToCollapse.size() << nl
|
||||
//<< " cells to create :" << cellsToCreate.size() << nl
|
||||
<< endl;
|
||||
|
||||
if
|
||||
(
|
||||
(cutBoundary && collapseEdge)
|
||||
|| (cutBoundary && cellsToSplit)
|
||||
|| (collapseEdge && cellsToSplit)
|
||||
)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Used more than one mesh modifying module "
|
||||
<< "(boundary cutting, cell splitting, edge collapsing)" << nl
|
||||
<< "Please do them in separate passes." << exit(FatalError);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Get calculating engine for all of outside
|
||||
const SubList<face> outsideFaces
|
||||
(
|
||||
mesh.faces(),
|
||||
mesh.nFaces() - mesh.nInternalFaces(),
|
||||
mesh.nInternalFaces()
|
||||
);
|
||||
|
||||
primitivePatch allBoundary(outsideFaces, mesh.points());
|
||||
|
||||
|
||||
// Look up mesh labels and convert to input for boundaryCutter.
|
||||
|
||||
bool validInputs = true;
|
||||
|
||||
|
||||
Info<< nl << "Looking up points to move ..." << nl << endl;
|
||||
Map<point> pointToPos(pointsToMove.size());
|
||||
forAll(pointsToMove, i)
|
||||
{
|
||||
const Pair<point>& pts = pointsToMove[i];
|
||||
|
||||
label pointi = findPoint(allBoundary, pts.first());
|
||||
|
||||
if (pointi == -1 || !pointToPos.insert(pointi, pts.second()))
|
||||
{
|
||||
Info<< "Could not insert mesh point " << pointi
|
||||
<< " for input point " << pts.first() << nl
|
||||
<< "Perhaps the point is already marked for moving?" << endl;
|
||||
validInputs = false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Info<< nl << "Looking up edges to split ..." << nl << endl;
|
||||
Map<List<point>> edgeToCuts(edgesToSplit.size());
|
||||
forAll(edgesToSplit, i)
|
||||
{
|
||||
const Pair<point>& pts = edgesToSplit[i];
|
||||
|
||||
label edgeI = findEdge(mesh, allBoundary, pts.first());
|
||||
|
||||
if
|
||||
(
|
||||
edgeI == -1
|
||||
|| !edgeToCuts.insert(edgeI, List<point>(1, pts.second()))
|
||||
)
|
||||
{
|
||||
Info<< "Could not insert mesh edge " << edgeI
|
||||
<< " for input point " << pts.first() << nl
|
||||
<< "Perhaps the edge is already marked for cutting?" << endl;
|
||||
|
||||
validInputs = false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Info<< nl << "Looking up faces to triangulate ..." << nl << endl;
|
||||
Map<point> faceToDecompose(facesToTriangulate.size());
|
||||
forAll(facesToTriangulate, i)
|
||||
{
|
||||
const Pair<point>& pts = facesToTriangulate[i];
|
||||
|
||||
label facei = findFace(mesh, allBoundary, pts.first());
|
||||
|
||||
if (facei == -1 || !faceToDecompose.insert(facei, pts.second()))
|
||||
{
|
||||
Info<< "Could not insert mesh face " << facei
|
||||
<< " for input point " << pts.first() << nl
|
||||
<< "Perhaps the face is already marked for splitting?" << endl;
|
||||
|
||||
validInputs = false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
Info<< nl << "Looking up cells to convert to pyramids around"
|
||||
<< " cell centre ..." << nl << endl;
|
||||
Map<point> cellToPyrCentre(cellsToPyramidise.size());
|
||||
forAll(cellsToPyramidise, i)
|
||||
{
|
||||
const Pair<point>& pts = cellsToPyramidise[i];
|
||||
|
||||
label celli = findCell(mesh, pts.first());
|
||||
|
||||
if (celli == -1 || !cellToPyrCentre.insert(celli, pts.second()))
|
||||
{
|
||||
Info<< "Could not insert mesh cell " << celli
|
||||
<< " for input point " << pts.first() << nl
|
||||
<< "Perhaps the cell is already marked for splitting?" << endl;
|
||||
|
||||
validInputs = false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Info<< nl << "Looking up edges to collapse ..." << nl << endl;
|
||||
Map<point> edgeToPos(edgesToCollapse.size());
|
||||
forAll(edgesToCollapse, i)
|
||||
{
|
||||
const Pair<point>& pts = edgesToCollapse[i];
|
||||
|
||||
label edgeI = findEdge(mesh, allBoundary, pts.first());
|
||||
|
||||
if (edgeI == -1 || !edgeToPos.insert(edgeI, pts.second()))
|
||||
{
|
||||
Info<< "Could not insert mesh edge " << edgeI
|
||||
<< " for input point " << pts.first() << nl
|
||||
<< "Perhaps the edge is already marked for collapsing?" << endl;
|
||||
|
||||
validInputs = false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
if (!validInputs)
|
||||
{
|
||||
Info<< nl << "There was a problem in one of the inputs in the"
|
||||
<< " dictionary. Not modifying mesh." << endl;
|
||||
}
|
||||
else if (cellToPyrCentre.size())
|
||||
{
|
||||
Info<< nl << "All input cells located. Modifying mesh." << endl;
|
||||
|
||||
// Mesh change engine
|
||||
cellSplitter cutter(mesh);
|
||||
|
||||
// Topo change container
|
||||
polyTopoChange meshMod(mesh);
|
||||
|
||||
// Insert commands into meshMod
|
||||
cutter.setRefinement(cellToPyrCentre, meshMod);
|
||||
|
||||
// Do changes
|
||||
autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh, false);
|
||||
|
||||
if (map().hasMotionPoints())
|
||||
{
|
||||
mesh.setPoints(map().preMotionPoints());
|
||||
}
|
||||
|
||||
cutter.topoChange(map());
|
||||
|
||||
if (!overwrite)
|
||||
{
|
||||
runTime++;
|
||||
}
|
||||
else
|
||||
{
|
||||
mesh.setInstance(oldInstance);
|
||||
}
|
||||
|
||||
// Write resulting mesh
|
||||
Info<< "Writing modified mesh to time " << runTime.name() << endl;
|
||||
mesh.write();
|
||||
}
|
||||
else if (edgeToPos.size())
|
||||
{
|
||||
Info<< nl << "All input edges located. Modifying mesh." << endl;
|
||||
|
||||
// Mesh change engine
|
||||
edgeCollapser cutter(mesh);
|
||||
|
||||
const edgeList& edges = mesh.edges();
|
||||
const pointField& points = mesh.points();
|
||||
|
||||
pointField newPoints(points);
|
||||
|
||||
PackedBoolList collapseEdge(mesh.nEdges());
|
||||
Map<point> collapsePointToLocation(mesh.nPoints());
|
||||
|
||||
// Get new positions and construct collapse network
|
||||
forAllConstIter(Map<point>, edgeToPos, iter)
|
||||
{
|
||||
label edgeI = iter.key();
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
collapseEdge[edgeI] = true;
|
||||
collapsePointToLocation.set(e[1], points[e[0]]);
|
||||
|
||||
newPoints[e[0]] = iter();
|
||||
}
|
||||
|
||||
// Move master point to destination.
|
||||
mesh.setPoints(newPoints);
|
||||
|
||||
List<pointEdgeCollapse> allPointInfo;
|
||||
const globalIndex globalPoints(mesh.nPoints());
|
||||
labelList pointPriority(mesh.nPoints(), 0);
|
||||
|
||||
cutter.consistentCollapse
|
||||
(
|
||||
globalPoints,
|
||||
pointPriority,
|
||||
collapsePointToLocation,
|
||||
collapseEdge,
|
||||
allPointInfo
|
||||
);
|
||||
|
||||
// Topo change container
|
||||
polyTopoChange meshMod(mesh);
|
||||
|
||||
// Insert
|
||||
cutter.setRefinement(allPointInfo, meshMod);
|
||||
|
||||
// Do changes
|
||||
autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh, false);
|
||||
|
||||
if (map().hasMotionPoints())
|
||||
{
|
||||
mesh.setPoints(map().preMotionPoints());
|
||||
}
|
||||
|
||||
// Not implemented yet:
|
||||
// cutter.topoChange(map());
|
||||
|
||||
|
||||
if (!overwrite)
|
||||
{
|
||||
runTime++;
|
||||
}
|
||||
else
|
||||
{
|
||||
mesh.setInstance(oldInstance);
|
||||
}
|
||||
|
||||
// Write resulting mesh
|
||||
Info<< "Writing modified mesh to time " << runTime.name() << endl;
|
||||
mesh.write();
|
||||
}
|
||||
else
|
||||
{
|
||||
Info<< nl << "All input points located. Modifying mesh." << endl;
|
||||
|
||||
// Mesh change engine
|
||||
boundaryCutter cutter(mesh);
|
||||
|
||||
// Topo change container
|
||||
polyTopoChange meshMod(mesh);
|
||||
|
||||
// Insert commands into meshMod
|
||||
cutter.setRefinement
|
||||
(
|
||||
pointToPos,
|
||||
edgeToCuts,
|
||||
Map<labelPair>(0), // Faces to split diagonally
|
||||
faceToDecompose, // Faces to triangulate
|
||||
meshMod
|
||||
);
|
||||
|
||||
// Do changes
|
||||
autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh, false);
|
||||
|
||||
if (map().hasMotionPoints())
|
||||
{
|
||||
mesh.setPoints(map().preMotionPoints());
|
||||
}
|
||||
|
||||
cutter.topoChange(map());
|
||||
|
||||
if (!overwrite)
|
||||
{
|
||||
runTime++;
|
||||
}
|
||||
else
|
||||
{
|
||||
mesh.setInstance(oldInstance);
|
||||
}
|
||||
|
||||
// Write resulting mesh
|
||||
Info<< "Writing modified mesh to time " << runTime.name() << endl;
|
||||
mesh.write();
|
||||
}
|
||||
|
||||
|
||||
Info<< "\nEnd\n" << endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -4,7 +4,6 @@ polyTopoChange/polyTopoChange.C
|
||||
polyTopoChange/addPatchCellLayer.C
|
||||
polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
|
||||
polyTopoChange/edgeCollapser.C
|
||||
polyTopoChange/faceCollapser.C
|
||||
polyTopoChange/removeCells.C
|
||||
polyTopoChange/removeFaces.C
|
||||
polyTopoChange/refinementData.C
|
||||
@ -12,7 +11,6 @@ polyTopoChange/refinementDistanceData.C
|
||||
polyTopoChange/removePoints.C
|
||||
polyTopoChange/combineFaces.C
|
||||
polyTopoChange/duplicatePoints.C
|
||||
polyTopoChange/tetDecomposer.C
|
||||
|
||||
hexRef8 = polyTopoChange/hexRef8
|
||||
$(hexRef8)/hexRef8.C
|
||||
@ -31,11 +29,6 @@ polyTopoChanger/attachDetach/attachInterface.C
|
||||
polyTopoChanger/attachDetach/detachInterface.C
|
||||
polyTopoChanger/attachDetach/attachDetachPointMatchMap.C
|
||||
|
||||
polyTopoChanger/layerAdditionRemoval/layerAdditionRemoval.C
|
||||
polyTopoChanger/layerAdditionRemoval/setLayerPairing.C
|
||||
polyTopoChanger/layerAdditionRemoval/addCellLayer.C
|
||||
polyTopoChanger/layerAdditionRemoval/removeCellLayer.C
|
||||
|
||||
enrichedPatch = polyTopoChanger/slidingInterface/enrichedPatch
|
||||
$(enrichedPatch)/enrichedPatch.C
|
||||
$(enrichedPatch)/enrichedPatchPointMap.C
|
||||
|
||||
@ -1,505 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "faceCollapser.H"
|
||||
#include "polyMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "polyModifyPoint.H"
|
||||
#include "polyModifyFace.H"
|
||||
#include "polyRemoveFace.H"
|
||||
#include "SortableList.H"
|
||||
#include "meshTools.H"
|
||||
#include "OFstream.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
// Insert labelList into labelHashSet. Optional excluded element.
|
||||
void Foam::faceCollapser::insert
|
||||
(
|
||||
const labelList& elems,
|
||||
const label excludeElem,
|
||||
labelHashSet& set
|
||||
)
|
||||
{
|
||||
forAll(elems, i)
|
||||
{
|
||||
if (elems[i] != excludeElem)
|
||||
{
|
||||
set.insert(elems[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Find edge amongst candidate edges. FatalError if none.
|
||||
Foam::label Foam::faceCollapser::findEdge
|
||||
(
|
||||
const edgeList& edges,
|
||||
const labelList& edgeLabels,
|
||||
const label v0,
|
||||
const label v1
|
||||
)
|
||||
{
|
||||
forAll(edgeLabels, i)
|
||||
{
|
||||
label edgeI = edgeLabels[i];
|
||||
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
if
|
||||
(
|
||||
(e[0] == v0 && e[1] == v1)
|
||||
|| (e[0] == v1 && e[1] == v0)
|
||||
)
|
||||
{
|
||||
return edgeI;
|
||||
}
|
||||
}
|
||||
|
||||
FatalErrorInFunction
|
||||
<< " and " << v1 << " in edge labels " << edgeLabels
|
||||
<< abort(FatalError);
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
// Replace vertices in face
|
||||
void Foam::faceCollapser::filterFace
|
||||
(
|
||||
const Map<labelList>& splitEdges,
|
||||
const label facei,
|
||||
polyTopoChange& meshMod
|
||||
) const
|
||||
{
|
||||
const face& f = mesh_.faces()[facei];
|
||||
const labelList& fEdges = mesh_.faceEdges()[facei];
|
||||
|
||||
// Space for replaced vertices and split edges.
|
||||
DynamicList<label> newFace(10 * f.size());
|
||||
|
||||
forAll(f, fp)
|
||||
{
|
||||
label v0 = f[fp];
|
||||
|
||||
newFace.append(v0);
|
||||
|
||||
// Look ahead one to get edge.
|
||||
label fp1 = f.fcIndex(fp);
|
||||
|
||||
label v1 = f[fp1];
|
||||
|
||||
// Get split on edge if any.
|
||||
label edgeI = findEdge(mesh_.edges(), fEdges, v0, v1);
|
||||
|
||||
Map<labelList>::const_iterator edgeFnd =
|
||||
splitEdges.find(edgeI);
|
||||
|
||||
if (edgeFnd != splitEdges.end())
|
||||
{
|
||||
// edgeI has been split (by introducing new vertices).
|
||||
// Insert new vertices in face in correct order
|
||||
// (splitEdges was constructed to be from edge.start() to end())
|
||||
|
||||
const labelList& extraVerts = edgeFnd();
|
||||
|
||||
if (v0 == mesh_.edges()[edgeI].start())
|
||||
{
|
||||
forAll(extraVerts, i)
|
||||
{
|
||||
newFace.append(extraVerts[i]);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
forAllReverse(extraVerts, i)
|
||||
{
|
||||
newFace.append(extraVerts[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
face newF(newFace.shrink());
|
||||
|
||||
// Pout<< "Modifying face:" << facei << " from " << f << " to " << newFace
|
||||
// << endl;
|
||||
|
||||
if (newF != f)
|
||||
{
|
||||
label nei = -1;
|
||||
|
||||
label patchi = -1;
|
||||
|
||||
if (mesh_.isInternalFace(facei))
|
||||
{
|
||||
nei = mesh_.faceNeighbour()[facei];
|
||||
}
|
||||
else
|
||||
{
|
||||
patchi = mesh_.boundaryMesh().whichPatch(facei);
|
||||
}
|
||||
|
||||
// Get current zone info
|
||||
label zoneID = mesh_.faceZones().whichZone(facei);
|
||||
|
||||
bool zoneFlip = false;
|
||||
|
||||
if (zoneID >= 0)
|
||||
{
|
||||
const faceZone& fZone = mesh_.faceZones()[zoneID];
|
||||
|
||||
zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
|
||||
}
|
||||
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
newF, // modified face
|
||||
facei, // label of face being modified
|
||||
mesh_.faceOwner()[facei], // owner
|
||||
nei, // neighbour
|
||||
false, // face flip
|
||||
patchi, // patch for face
|
||||
false, // remove from zone
|
||||
zoneID, // zone for face
|
||||
zoneFlip // face flip in zone
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
// Construct from mesh
|
||||
Foam::faceCollapser::faceCollapser(const polyMesh& mesh)
|
||||
:
|
||||
mesh_(mesh)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::faceCollapser::setRefinement
|
||||
(
|
||||
const labelList& faceLabels,
|
||||
const labelList& fpStart,
|
||||
const labelList& fpEnd,
|
||||
polyTopoChange& meshMod
|
||||
) const
|
||||
{
|
||||
const pointField& points = mesh_.points();
|
||||
const edgeList& edges = mesh_.edges();
|
||||
const faceList& faces = mesh_.faces();
|
||||
const labelListList& edgeFaces = mesh_.edgeFaces();
|
||||
|
||||
|
||||
// From split edge to newly introduced point(s). Can be more than one per
|
||||
// edge!
|
||||
Map<labelList> splitEdges(faceLabels.size());
|
||||
|
||||
// Mark faces in any way affect by modifying points/edges. Used later on
|
||||
// to prevent having to redo all faces.
|
||||
labelHashSet affectedFaces(4*faceLabels.size());
|
||||
|
||||
|
||||
//
|
||||
// Add/remove vertices and construct mapping
|
||||
//
|
||||
|
||||
forAll(faceLabels, i)
|
||||
{
|
||||
const label facei = faceLabels[i];
|
||||
|
||||
const face& f = faces[facei];
|
||||
|
||||
const label fpA = fpStart[i];
|
||||
const label fpB = fpEnd[i];
|
||||
|
||||
const point& pA = points[f[fpA]];
|
||||
const point& pB = points[f[fpB]];
|
||||
|
||||
Pout<< "Face:" << f << " collapsed to fp:" << fpA << ' ' << fpB
|
||||
<< " with points:" << pA << ' ' << pB
|
||||
<< endl;
|
||||
|
||||
// Create line from fpA to fpB
|
||||
linePointRef lineAB(pA, pB);
|
||||
|
||||
// Get projections of all vertices onto line.
|
||||
|
||||
// Distance(squared) to pA for every point on face.
|
||||
SortableList<scalar> dist(f.size());
|
||||
|
||||
dist[fpA] = 0;
|
||||
dist[fpB] = magSqr(pB - pA);
|
||||
|
||||
// Step from fpA to fpB
|
||||
// ~~~~~~~~~~~~~~~~~~~~
|
||||
// (by incrementing)
|
||||
|
||||
label fpMin1 = fpA;
|
||||
label fp = f.fcIndex(fpMin1);
|
||||
|
||||
while (fp != fpB)
|
||||
{
|
||||
// See where fp sorts. Make sure it is above fpMin1!
|
||||
pointHit near = lineAB.nearestDist(points[f[fp]]);
|
||||
|
||||
scalar w = magSqr(near.rawPoint() - pA);
|
||||
|
||||
if (w <= dist[fpMin1])
|
||||
{
|
||||
// Offset.
|
||||
w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
|
||||
|
||||
point newPoint
|
||||
(
|
||||
pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
|
||||
);
|
||||
|
||||
Pout<< "Adapting position of vertex " << f[fp] << " on face "
|
||||
<< f << " from " << near.rawPoint() << " to " << newPoint
|
||||
<< endl;
|
||||
|
||||
near.setPoint(newPoint);
|
||||
}
|
||||
|
||||
// Responsibility of caller to make sure polyModifyPoint is only
|
||||
// called once per point. (so max only one collapse face per
|
||||
// edge)
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyPoint
|
||||
(
|
||||
f[fp],
|
||||
near.rawPoint(),
|
||||
false,
|
||||
-1,
|
||||
true
|
||||
)
|
||||
);
|
||||
|
||||
dist[fp] = w;
|
||||
|
||||
// Step to next vertex.
|
||||
fpMin1 = fp;
|
||||
fp = f.fcIndex(fpMin1);
|
||||
}
|
||||
|
||||
// Step from fpA to fpB
|
||||
// ~~~~~~~~~~~~~~~~~~~~
|
||||
// (by decrementing)
|
||||
|
||||
fpMin1 = fpA;
|
||||
fp = f.rcIndex(fpMin1);
|
||||
|
||||
while (fp != fpB)
|
||||
{
|
||||
// See where fp sorts. Make sure it is below fpMin1!
|
||||
pointHit near = lineAB.nearestDist(points[f[fp]]);
|
||||
|
||||
scalar w = magSqr(near.rawPoint() - pA);
|
||||
|
||||
if (w <= dist[fpMin1])
|
||||
{
|
||||
// Offset.
|
||||
w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
|
||||
|
||||
point newPoint
|
||||
(
|
||||
pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
|
||||
);
|
||||
|
||||
Pout<< "Adapting position of vertex " << f[fp] << " on face "
|
||||
<< f << " from " << near.rawPoint() << " to " << newPoint
|
||||
<< endl;
|
||||
|
||||
near.setPoint(newPoint);
|
||||
}
|
||||
|
||||
// Responsibility of caller to make sure polyModifyPoint is only
|
||||
// called once per point. (so max only one collapse face per
|
||||
// edge)
|
||||
meshMod.setAction
|
||||
(
|
||||
polyModifyPoint
|
||||
(
|
||||
f[fp],
|
||||
near.rawPoint(),
|
||||
false,
|
||||
-1,
|
||||
true
|
||||
)
|
||||
);
|
||||
|
||||
dist[fp] = w;
|
||||
|
||||
// Step to previous vertex.
|
||||
fpMin1 = fp;
|
||||
fp = f.rcIndex(fpMin1);
|
||||
}
|
||||
|
||||
dist.sort();
|
||||
|
||||
// Check that fpB sorts latest.
|
||||
if (dist.indices()[dist.size()-1] != fpB)
|
||||
{
|
||||
OFstream str("conflictingFace.obj");
|
||||
meshTools::writeOBJ(str, faceList(1, f), points);
|
||||
|
||||
FatalErrorInFunction
|
||||
<< "Trying to collapse face:" << facei << " vertices:" << f
|
||||
<< " to edges between vertices " << f[fpA] << " and "
|
||||
<< f[fpB] << " but " << f[fpB] << " does not seem to be the"
|
||||
<< " vertex furthest away from " << f[fpA] << endl
|
||||
<< "Dumped conflicting face to obj file conflictingFace.obj"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// From fp to index in sort:
|
||||
Pout<< "Face:" << f << " fpA:" << fpA << " fpB:" << fpB << nl;
|
||||
|
||||
labelList sortedFp(f.size());
|
||||
forAll(dist.indices(), i)
|
||||
{
|
||||
label fp = dist.indices()[i];
|
||||
|
||||
Pout<< " fp:" << fp << " distance:" << dist[i] << nl;
|
||||
|
||||
sortedFp[fp] = i;
|
||||
}
|
||||
|
||||
const labelList& fEdges = mesh_.faceEdges()[facei];
|
||||
|
||||
// Now look up all edges in the face and see if they get extra
|
||||
// vertices inserted and build an edge-to-intersected-points table.
|
||||
|
||||
// Order of inserted points is in edge order (from e.start to
|
||||
// e.end)
|
||||
|
||||
forAll(f, fp)
|
||||
{
|
||||
label fp1 = f.fcIndex(fp);
|
||||
|
||||
// Get index in sorted list
|
||||
label sorted0 = sortedFp[fp];
|
||||
label sorted1 = sortedFp[fp1];
|
||||
|
||||
// Get indices in increasing order
|
||||
DynamicList<label> edgePoints(f.size());
|
||||
|
||||
// Whether edgePoints are from fp to fp1
|
||||
bool fpToFp1;
|
||||
|
||||
if (sorted0 < sorted1)
|
||||
{
|
||||
fpToFp1 = true;
|
||||
|
||||
for (label j = sorted0+1; j < sorted1; j++)
|
||||
{
|
||||
edgePoints.append(f[dist.indices()[j]]);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
fpToFp1 = false;
|
||||
|
||||
for (label j = sorted1+1; j < sorted0; j++)
|
||||
{
|
||||
edgePoints.append(f[dist.indices()[j]]);
|
||||
}
|
||||
}
|
||||
|
||||
if (edgePoints.size())
|
||||
{
|
||||
edgePoints.shrink();
|
||||
|
||||
label edgeI = findEdge(edges, fEdges, f[fp], f[fp1]);
|
||||
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
if (fpToFp1 == (f[fp] == e.start()))
|
||||
{
|
||||
splitEdges.insert(edgeI, edgePoints);
|
||||
}
|
||||
else
|
||||
{
|
||||
reverse(edgePoints);
|
||||
splitEdges.insert(edgeI, edgePoints);
|
||||
}
|
||||
|
||||
// Mark all faces affected
|
||||
insert(edgeFaces[edgeI], facei, affectedFaces);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
forAllConstIter(Map<labelList>, splitEdges, iter)
|
||||
{
|
||||
Pout<< "Split edge:" << iter.key()
|
||||
<< " verts:" << mesh_.edges()[iter.key()]
|
||||
<< " in:" << nl;
|
||||
const labelList& edgePoints = iter();
|
||||
|
||||
forAll(edgePoints, i)
|
||||
{
|
||||
Pout<< " " << edgePoints[i] << nl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
// Remove faces.
|
||||
//
|
||||
|
||||
forAll(faceLabels, i)
|
||||
{
|
||||
const label facei = faceLabels[i];
|
||||
|
||||
meshMod.setAction(polyRemoveFace(facei));
|
||||
|
||||
// Update list of faces we still have to modify
|
||||
affectedFaces.erase(facei);
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
// Modify faces affected (but not removed)
|
||||
//
|
||||
|
||||
forAllConstIter(labelHashSet, affectedFaces, iter)
|
||||
{
|
||||
filterFace(splitEdges, iter.key(), meshMod);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,164 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Class
|
||||
Foam::faceCollapser
|
||||
|
||||
Description
|
||||
Collapses faces into edges. Used to remove sliver faces (faces with small
|
||||
area but non-zero span).
|
||||
|
||||
Passed in as
|
||||
- face label
|
||||
- the two indices in the face (fpA, fpB) which delimit the vertices to be
|
||||
kept.
|
||||
|
||||
Takes the vertices outside the range fpA..fpB and projects them onto the
|
||||
kept edges (edges using kept vertices only).
|
||||
|
||||
Note:
|
||||
- Use in combination with edge collapse to cleanup meshes.
|
||||
- Can not remove cells so will mess up trying to remove a face on a tet.
|
||||
- WIP. Should be combined with edge collapsing and cell collapsing into
|
||||
proper 'collapser'.
|
||||
- Caller is responsible for making sure kept vertices (fpA..fpB) for one
|
||||
face are not the vertices to be removed for another face.
|
||||
|
||||
SourceFiles
|
||||
faceCollapser.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef faceCollapser_H
|
||||
#define faceCollapser_H
|
||||
|
||||
#include "labelList.H"
|
||||
#include "point.H"
|
||||
#include "Map.H"
|
||||
#include "HashSet.H"
|
||||
#include "typeInfo.H"
|
||||
#include "edgeList.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// Forward declaration of classes
|
||||
class polyMesh;
|
||||
class polyTopoChange;
|
||||
class polyTopoChangeMap;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class faceCollapser Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class faceCollapser
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Reference to mesh
|
||||
const polyMesh& mesh_;
|
||||
|
||||
|
||||
// Static Functions
|
||||
|
||||
//- Insert labelList into labelHashSet. Optional excluded element.
|
||||
static void insert
|
||||
(
|
||||
const labelList& elems,
|
||||
const label excludeElem,
|
||||
labelHashSet& set
|
||||
);
|
||||
|
||||
//- Find edge amongst candidate edges.
|
||||
static label findEdge
|
||||
(
|
||||
const edgeList& edges,
|
||||
const labelList& edgeLabels,
|
||||
const label v0,
|
||||
const label v1
|
||||
);
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Replace vertices in face
|
||||
void filterFace
|
||||
(
|
||||
const Map<labelList>& splitEdges,
|
||||
const label facei,
|
||||
polyTopoChange& meshMod
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
ClassName("faceCollapser");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from mesh.
|
||||
faceCollapser(const polyMesh& mesh);
|
||||
|
||||
//- Disallow default bitwise copy construction
|
||||
faceCollapser(const faceCollapser&) = delete;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Edit
|
||||
|
||||
//- Collapse faces along endpoints. Play commands into
|
||||
// polyTopoChange to create mesh.
|
||||
void setRefinement
|
||||
(
|
||||
const labelList& faceLabels,
|
||||
const labelList& fpA,
|
||||
const labelList& fpB,
|
||||
polyTopoChange&
|
||||
) const;
|
||||
|
||||
//- Update stored quantities for new mesh labels.
|
||||
void topoChange(const polyTopoChangeMap&)
|
||||
{}
|
||||
|
||||
|
||||
// Member Operators
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const faceCollapser&) = delete;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,713 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "tetDecomposer.H"
|
||||
#include "meshTools.H"
|
||||
#include "polyMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "polyTopoChangeMap.H"
|
||||
#include "OFstream.H"
|
||||
#include "EdgeMap.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(tetDecomposer, 0);
|
||||
|
||||
template<>
|
||||
const char* NamedEnum<tetDecomposer::decompositionType, 2>::names[] =
|
||||
{
|
||||
"faceCentre",
|
||||
"faceDiagonal"
|
||||
};
|
||||
|
||||
const NamedEnum<tetDecomposer::decompositionType, 2>
|
||||
tetDecomposer::decompositionTypeNames;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::tetDecomposer::modifyFace
|
||||
(
|
||||
polyTopoChange& meshMod,
|
||||
const face& f,
|
||||
const label facei,
|
||||
const label own,
|
||||
const label nei,
|
||||
const label patchi,
|
||||
const label zoneI,
|
||||
const bool zoneFlip
|
||||
) const
|
||||
{
|
||||
// First usage of face. Modify.
|
||||
if (nei == -1 || own < nei)
|
||||
{
|
||||
meshMod.modifyFace
|
||||
(
|
||||
f, // modified face
|
||||
facei, // label of face
|
||||
own, // owner
|
||||
nei, // neighbour
|
||||
false, // face flip
|
||||
patchi, // patch for face
|
||||
zoneI, // zone for face
|
||||
zoneFlip // face flip in zone
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
meshMod.modifyFace
|
||||
(
|
||||
f.reverseFace(), // modified face
|
||||
facei, // label of face
|
||||
nei, // owner
|
||||
own, // neighbour
|
||||
true, // face flip
|
||||
patchi, // patch for face
|
||||
zoneI, // zone for face
|
||||
!zoneFlip // face flip in zone
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::tetDecomposer::addFace
|
||||
(
|
||||
polyTopoChange& meshMod,
|
||||
const face& f,
|
||||
const label own,
|
||||
const label nei,
|
||||
const label masterPointID,
|
||||
const label masterEdgeID,
|
||||
const label masterFaceID,
|
||||
const label patchi,
|
||||
const label zoneI,
|
||||
const bool zoneFlip
|
||||
) const
|
||||
{
|
||||
// Second or more usage of face. Add.
|
||||
if (nei == -1 || own < nei)
|
||||
{
|
||||
meshMod.addFace
|
||||
(
|
||||
f, // modified face
|
||||
own, // owner
|
||||
nei, // neighbour
|
||||
masterPointID, // master point
|
||||
masterEdgeID, // master edge
|
||||
masterFaceID, // master face
|
||||
false, // face flip
|
||||
patchi, // patch for face
|
||||
zoneI, // zone for face
|
||||
zoneFlip // face flip in zone
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
meshMod.addFace
|
||||
(
|
||||
f.reverseFace(), // modified face
|
||||
nei, // owner
|
||||
own, // neighbour
|
||||
masterPointID, // master point
|
||||
masterEdgeID, // master edge
|
||||
masterFaceID, // master face
|
||||
true, // face flip
|
||||
patchi, // patch for face
|
||||
zoneI, // zone for face
|
||||
!zoneFlip // face flip in zone
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Work out triangle index given the starting vertex in the face
|
||||
Foam::label Foam::tetDecomposer::triIndex(const label facei, const label fp)
|
||||
const
|
||||
{
|
||||
const face& f = mesh_.faces()[facei];
|
||||
const label fp0 = mesh_.tetBasePtIs()[facei];
|
||||
|
||||
// Work out triangle index on this face
|
||||
label thisTriI;
|
||||
if (fp == fp0)
|
||||
{
|
||||
thisTriI = 0;
|
||||
}
|
||||
else if (fp == f.rcIndex(fp0))
|
||||
{
|
||||
thisTriI = f.size()-3;
|
||||
}
|
||||
else
|
||||
{
|
||||
thisTriI = (fp-fp0-1) % (f.size()-2);
|
||||
}
|
||||
return thisTriI;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::tetDecomposer::tetDecomposer(const polyMesh& mesh)
|
||||
:
|
||||
mesh_(mesh)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::tetDecomposer::setRefinement
|
||||
(
|
||||
const decompositionType decomposeType,
|
||||
polyTopoChange& meshMod
|
||||
)
|
||||
{
|
||||
cellToPoint_.setSize(mesh_.nCells());
|
||||
forAll(mesh_.cellCentres(), celli)
|
||||
{
|
||||
// Any point on the cell
|
||||
label masterPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
|
||||
|
||||
cellToPoint_[celli] = meshMod.addPoint
|
||||
(
|
||||
mesh_.cellCentres()[celli],
|
||||
masterPointi,
|
||||
-1,
|
||||
true
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Add face centre points
|
||||
if (decomposeType == FACE_CENTRE_TRIS)
|
||||
{
|
||||
faceToPoint_.setSize(mesh_.nFaces());
|
||||
forAll(mesh_.faceCentres(), facei)
|
||||
{
|
||||
// Any point on the face
|
||||
const label masterPointi = mesh_.faces()[facei][0];
|
||||
|
||||
faceToPoint_[facei] = meshMod.addPoint
|
||||
(
|
||||
mesh_.faceCentres()[facei],
|
||||
masterPointi,
|
||||
-1,
|
||||
true
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Per face, per point (faceCentre) or triangle (faceDiag) the added cell
|
||||
faceOwnerCells_.setSize(mesh_.nFaces());
|
||||
faceNeighbourCells_.setSize(mesh_.nFaces());
|
||||
|
||||
if (decomposeType == FACE_CENTRE_TRIS)
|
||||
{
|
||||
forAll(faceOwnerCells_, facei)
|
||||
{
|
||||
const face& f = mesh_.faces()[facei];
|
||||
faceOwnerCells_[facei].setSize(f.size(), -1);
|
||||
faceNeighbourCells_[facei].setSize(f.size(), -1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Force construction of diagonal decomposition
|
||||
(void)mesh_.tetBasePtIs();
|
||||
|
||||
forAll(faceOwnerCells_, facei)
|
||||
{
|
||||
const face& f = mesh_.faces()[facei];
|
||||
faceOwnerCells_[facei].setSize(f.size()-2, -1);
|
||||
faceNeighbourCells_[facei].setSize(f.size()-2, -1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
forAll(mesh_.cells(), celli)
|
||||
{
|
||||
const cell& cFaces = mesh_.cells()[celli];
|
||||
|
||||
EdgeMap<label> edgeToFace(8*cFaces.size());
|
||||
|
||||
forAll(cFaces, cFacei)
|
||||
{
|
||||
label facei = cFaces[cFacei];
|
||||
const face& f = mesh_.faces()[facei];
|
||||
|
||||
// Get reference to either owner or neighbour
|
||||
labelList& added =
|
||||
(
|
||||
(mesh_.faceOwner()[facei] == celli)
|
||||
? faceOwnerCells_[facei]
|
||||
: faceNeighbourCells_[facei]
|
||||
);
|
||||
|
||||
if (decomposeType == FACE_CENTRE_TRIS)
|
||||
{
|
||||
forAll(f, fp)
|
||||
{
|
||||
if (cFacei == 0 && fp == 0)
|
||||
{
|
||||
// Reuse cell itself
|
||||
added[fp] = celli;
|
||||
}
|
||||
else
|
||||
{
|
||||
added[fp] = meshMod.addCell
|
||||
(
|
||||
-1, // masterPoint
|
||||
-1, // masterEdge
|
||||
-1, // masterFace
|
||||
celli, // masterCell
|
||||
mesh_.cellZones().whichZone(celli)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
for (label triI = 0; triI < f.size()-2; triI++)
|
||||
{
|
||||
if (cFacei == 0 && triI == 0)
|
||||
{
|
||||
// Reuse cell itself
|
||||
added[triI] = celli;
|
||||
}
|
||||
else
|
||||
{
|
||||
added[triI] = meshMod.addCell
|
||||
(
|
||||
-1, // masterPoint
|
||||
-1, // masterEdge
|
||||
-1, // masterFace
|
||||
celli, // masterCell
|
||||
mesh_.cellZones().whichZone(celli)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Add triangle faces
|
||||
face triangle(3);
|
||||
|
||||
forAll(mesh_.faces(), facei)
|
||||
{
|
||||
label own = mesh_.faceOwner()[facei];
|
||||
const labelList& addedOwn = faceOwnerCells_[facei];
|
||||
const labelList& addedNei = faceNeighbourCells_[facei];
|
||||
const face& f = mesh_.faces()[facei];
|
||||
|
||||
label patchi = -1;
|
||||
if (facei >= mesh_.nInternalFaces())
|
||||
{
|
||||
patchi = mesh_.boundaryMesh().whichPatch(facei);
|
||||
}
|
||||
|
||||
label zoneI = mesh_.faceZones().whichZone(facei);
|
||||
bool zoneFlip = false;
|
||||
if (zoneI != -1)
|
||||
{
|
||||
const faceZone& fz = mesh_.faceZones()[zoneI];
|
||||
zoneFlip = fz.flipMap()[fz.whichFace(facei)];
|
||||
}
|
||||
|
||||
|
||||
if (decomposeType == FACE_CENTRE_TRIS)
|
||||
{
|
||||
forAll(f, fp)
|
||||
{
|
||||
// 1. Front triangle (decomposition of face itself)
|
||||
// (between owner and neighbour cell)
|
||||
{
|
||||
triangle[0] = f[fp];
|
||||
triangle[1] = f[f.fcIndex(fp)];
|
||||
triangle[2] = faceToPoint_[facei];
|
||||
|
||||
if (fp == 0)
|
||||
{
|
||||
modifyFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
facei,
|
||||
addedOwn[fp],
|
||||
addedNei[fp],
|
||||
patchi,
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
addedOwn[fp],
|
||||
addedNei[fp],
|
||||
-1, // point
|
||||
-1, // edge
|
||||
facei, // face
|
||||
patchi,
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 2. Within owner cell - to cell centre
|
||||
{
|
||||
label newOwn = addedOwn[f.rcIndex(fp)];
|
||||
label newNei = addedOwn[fp];
|
||||
|
||||
triangle[0] = f[fp];
|
||||
triangle[1] = cellToPoint_[own];
|
||||
triangle[2] = faceToPoint_[facei];
|
||||
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
newOwn,
|
||||
newNei,
|
||||
f[fp], // point
|
||||
-1, // edge
|
||||
-1, // face
|
||||
-1, // patchi
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
// 2b. Within neighbour cell - to cell centre
|
||||
if (facei < mesh_.nInternalFaces())
|
||||
{
|
||||
label newOwn = addedNei[f.rcIndex(fp)];
|
||||
label newNei = addedNei[fp];
|
||||
|
||||
triangle[0] = f[fp];
|
||||
triangle[1] = faceToPoint_[facei];
|
||||
triangle[2] = cellToPoint_[mesh_.faceNeighbour()[facei]];
|
||||
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
newOwn,
|
||||
newNei,
|
||||
f[fp], // point
|
||||
-1, // edge
|
||||
-1, // face
|
||||
-1, // patchi
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
label fp0 = mesh_.tetBasePtIs()[facei];
|
||||
label fp = f.fcIndex(fp0);
|
||||
|
||||
for (label triI = 0; triI < f.size()-2; triI++)
|
||||
{
|
||||
label nextTri = triI+1;
|
||||
if (nextTri >= f.size()-2)
|
||||
{
|
||||
nextTri -= f.size()-2;
|
||||
}
|
||||
label nextFp = f.fcIndex(fp);
|
||||
|
||||
|
||||
// Triangle triI consisting of f[fp0], f[fp], f[nextFp]
|
||||
|
||||
|
||||
// 1. Front triangle (decomposition of face itself)
|
||||
// (between owner and neighbour cell)
|
||||
{
|
||||
triangle[0] = f[fp0];
|
||||
triangle[1] = f[fp];
|
||||
triangle[2] = f[nextFp];
|
||||
|
||||
if (triI == 0)
|
||||
{
|
||||
modifyFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
facei,
|
||||
addedOwn[triI],
|
||||
addedNei[triI],
|
||||
patchi,
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
addedOwn[triI],
|
||||
addedNei[triI],
|
||||
-1, // point
|
||||
-1, // edge
|
||||
facei, // face
|
||||
patchi,
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 2. Within owner cell - diagonal to cell centre
|
||||
if (triI < f.size()-3)
|
||||
{
|
||||
label newOwn = addedOwn[triI];
|
||||
label newNei = addedOwn[nextTri];
|
||||
|
||||
triangle[0] = f[fp0];
|
||||
triangle[1] = f[nextFp];
|
||||
triangle[2] = cellToPoint_[own];
|
||||
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
newOwn,
|
||||
newNei,
|
||||
f[fp], // point
|
||||
-1, // edge
|
||||
-1, // face
|
||||
-1, // patchi
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
|
||||
// 2b. Within neighbour cell - to cell centre
|
||||
if (facei < mesh_.nInternalFaces())
|
||||
{
|
||||
label newOwn = addedNei[triI];
|
||||
label newNei = addedNei[nextTri];
|
||||
|
||||
triangle[0] = f[nextFp];
|
||||
triangle[1] = f[fp0];
|
||||
triangle[2] =
|
||||
cellToPoint_[mesh_.faceNeighbour()[facei]];
|
||||
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
newOwn,
|
||||
newNei,
|
||||
f[fp], // point
|
||||
-1, // edge
|
||||
-1, // face
|
||||
-1, // patchi
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
fp = nextFp;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Add triangles for all edges.
|
||||
EdgeMap<label> edgeToFace;
|
||||
|
||||
forAll(mesh_.cells(), celli)
|
||||
{
|
||||
const cell& cFaces = mesh_.cells()[celli];
|
||||
|
||||
edgeToFace.clear();
|
||||
|
||||
forAll(cFaces, cFacei)
|
||||
{
|
||||
label facei = cFaces[cFacei];
|
||||
|
||||
label zoneI = mesh_.faceZones().whichZone(facei);
|
||||
bool zoneFlip = false;
|
||||
if (zoneI != -1)
|
||||
{
|
||||
const faceZone& fz = mesh_.faceZones()[zoneI];
|
||||
zoneFlip = fz.flipMap()[fz.whichFace(facei)];
|
||||
}
|
||||
|
||||
const face& f = mesh_.faces()[facei];
|
||||
// const labelList& fEdges = mesh_.faceEdges()[facei];
|
||||
forAll(f, fp)
|
||||
{
|
||||
label p0 = f[fp];
|
||||
label p1 = f[f.fcIndex(fp)];
|
||||
const edge e(p0, p1);
|
||||
|
||||
EdgeMap<label>::const_iterator edgeFnd = edgeToFace.find(e);
|
||||
if (edgeFnd == edgeToFace.end())
|
||||
{
|
||||
edgeToFace.insert(e, facei);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Found the other face on the edge.
|
||||
label otherFacei = edgeFnd();
|
||||
const face& otherF = mesh_.faces()[otherFacei];
|
||||
|
||||
// Found the other face on the edge. Note that since
|
||||
// we are looping in the same order the tets added for
|
||||
// otherFacei will be before those of facei
|
||||
|
||||
label otherFp = findIndex(otherF, p0);
|
||||
if (otherF.nextLabel(otherFp) == p1)
|
||||
{
|
||||
// ok. otherFp is first vertex of edge.
|
||||
}
|
||||
else if (otherF.prevLabel(otherFp) == p1)
|
||||
{
|
||||
otherFp = otherF.rcIndex(otherFp);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "problem." << abort(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// Triangle from edge to cell centre
|
||||
if (mesh_.faceOwner()[facei] == celli)
|
||||
{
|
||||
triangle[0] = p0;
|
||||
triangle[1] = p1;
|
||||
triangle[2] = cellToPoint_[celli];
|
||||
}
|
||||
else
|
||||
{
|
||||
triangle[0] = p1;
|
||||
triangle[1] = p0;
|
||||
triangle[2] = cellToPoint_[celli];
|
||||
}
|
||||
|
||||
// Determine tets on either side
|
||||
label thisTet, otherTet;
|
||||
|
||||
if (decomposeType == FACE_CENTRE_TRIS)
|
||||
{
|
||||
if (mesh_.faceOwner()[facei] == celli)
|
||||
{
|
||||
thisTet = faceOwnerCells_[facei][fp];
|
||||
}
|
||||
else
|
||||
{
|
||||
thisTet = faceNeighbourCells_[facei][fp];
|
||||
}
|
||||
|
||||
if (mesh_.faceOwner()[otherFacei] == celli)
|
||||
{
|
||||
otherTet = faceOwnerCells_[otherFacei][otherFp];
|
||||
}
|
||||
else
|
||||
{
|
||||
otherTet =
|
||||
faceNeighbourCells_[otherFacei][otherFp];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
label thisTriI = triIndex(facei, fp);
|
||||
if (mesh_.faceOwner()[facei] == celli)
|
||||
{
|
||||
thisTet = faceOwnerCells_[facei][thisTriI];
|
||||
}
|
||||
else
|
||||
{
|
||||
thisTet = faceNeighbourCells_[facei][thisTriI];
|
||||
}
|
||||
|
||||
label otherTriI = triIndex(otherFacei, otherFp);
|
||||
if (mesh_.faceOwner()[otherFacei] == celli)
|
||||
{
|
||||
otherTet = faceOwnerCells_[otherFacei][otherTriI];
|
||||
}
|
||||
else
|
||||
{
|
||||
otherTet =
|
||||
faceNeighbourCells_[otherFacei][otherTriI];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
addFace
|
||||
(
|
||||
meshMod,
|
||||
triangle,
|
||||
otherTet,
|
||||
thisTet,
|
||||
-1, // masterPoint
|
||||
-1, // fEdges[fp], // masterEdge
|
||||
facei, // masterFace
|
||||
-1, // patchi
|
||||
zoneI,
|
||||
zoneFlip
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::tetDecomposer::topoChange(const polyTopoChangeMap& map)
|
||||
{
|
||||
inplaceRenumber(map.reversePointMap(), cellToPoint_);
|
||||
inplaceRenumber(map.reversePointMap(), faceToPoint_);
|
||||
|
||||
forAll(faceOwnerCells_, facei)
|
||||
{
|
||||
inplaceRenumber(map.reverseCellMap(), faceOwnerCells_[facei]);
|
||||
}
|
||||
forAll(faceNeighbourCells_, facei)
|
||||
{
|
||||
inplaceRenumber(map.reverseCellMap(), faceNeighbourCells_[facei]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,206 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Class
|
||||
Foam::tetDecomposer
|
||||
|
||||
Description
|
||||
Decomposes polyMesh into tets.
|
||||
|
||||
SourceFiles
|
||||
tetDecomposer.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef tetDecomposer_H
|
||||
#define tetDecomposer_H
|
||||
|
||||
#include "DynamicList.H"
|
||||
#include "PackedBoolList.H"
|
||||
#include "boolList.H"
|
||||
#include "typeInfo.H"
|
||||
#include "NamedEnum.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
class polyMesh;
|
||||
class polyTopoChange;
|
||||
class face;
|
||||
class polyTopoChangeMap;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class tetDecomposer Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class tetDecomposer
|
||||
{
|
||||
public:
|
||||
|
||||
// Public data types
|
||||
|
||||
enum decompositionType
|
||||
{
|
||||
FACE_CENTRE_TRIS, //- Faces decomposed into triangles
|
||||
// using face-centre
|
||||
|
||||
FACE_DIAG_TRIS //- Faces decomposed into triangles diagonally
|
||||
};
|
||||
static const NamedEnum<decompositionType, 2> decompositionTypeNames;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
// Private Data
|
||||
|
||||
const polyMesh& mesh_;
|
||||
|
||||
//- From cell to tet point
|
||||
labelList cellToPoint_;
|
||||
|
||||
//- From face to tet point
|
||||
labelList faceToPoint_;
|
||||
|
||||
|
||||
// Per face, per point (faceCentre) or triangle (faceDiag)
|
||||
// the added tet on the owner side
|
||||
labelListList faceOwnerCells_;
|
||||
|
||||
// Per face, per point (faceCentre) or triangle (faceDiag)
|
||||
// the added tet on the neighbour side
|
||||
labelListList faceNeighbourCells_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Modify a face
|
||||
void modifyFace
|
||||
(
|
||||
polyTopoChange& meshMod,
|
||||
const face& f,
|
||||
const label facei,
|
||||
const label own,
|
||||
const label nei,
|
||||
const label patchi,
|
||||
const label zoneI,
|
||||
const bool zoneFlip
|
||||
) const;
|
||||
|
||||
//- Add a face
|
||||
void addFace
|
||||
(
|
||||
polyTopoChange& meshMod,
|
||||
const face& f,
|
||||
const label own,
|
||||
const label nei,
|
||||
const label masterPointID,
|
||||
const label masterEdgeID,
|
||||
const label masterFaceID,
|
||||
const label patchi,
|
||||
const label zoneI,
|
||||
const bool zoneFlip
|
||||
) const;
|
||||
|
||||
//- Work out triangle index given the starting vertex in the face
|
||||
label triIndex(const label facei, const label fp) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
ClassName("tetDecomposer");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from mesh
|
||||
tetDecomposer(const polyMesh&);
|
||||
|
||||
//- Disallow default bitwise copy construction
|
||||
tetDecomposer(const tetDecomposer&) = delete;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Access
|
||||
|
||||
//- From cell to tet point
|
||||
const labelList& cellToPoint() const
|
||||
{
|
||||
return cellToPoint_;
|
||||
}
|
||||
|
||||
//- From face to tet point
|
||||
const labelList& faceToPoint() const
|
||||
{
|
||||
return faceToPoint_;
|
||||
}
|
||||
|
||||
|
||||
//- Per face, per point (faceCentre) or triangle (faceDiag)
|
||||
// the added tet on the owner side
|
||||
const labelListList& faceOwnerCells() const
|
||||
{
|
||||
return faceOwnerCells_;
|
||||
}
|
||||
|
||||
//- Per face, per point (faceCentre) or triangle (faceDiag)
|
||||
// the added tet on the neighbour side
|
||||
const labelListList& faceNeighbourCells() const
|
||||
{
|
||||
return faceNeighbourCells_;
|
||||
}
|
||||
|
||||
|
||||
// Edit
|
||||
|
||||
//- Insert all changes into meshMod to convert the polyMesh into
|
||||
// tets.
|
||||
void setRefinement
|
||||
(
|
||||
const decompositionType decomposeType,
|
||||
polyTopoChange& meshMod
|
||||
);
|
||||
|
||||
//- Force recalculation of locally stored data on topological change
|
||||
void topoChange(const polyTopoChangeMap&);
|
||||
|
||||
|
||||
// Member Operators
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const tetDecomposer&) = delete;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,637 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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();
|
||||
}
|
||||
|
||||
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
|
||||
|
||||
const vectorField pointOffsets(extrusionDir());
|
||||
|
||||
// Add the new points
|
||||
labelList addedPoints(mp.size());
|
||||
|
||||
forAll(mp, pointi)
|
||||
{
|
||||
// Add the point nominal thickness away from the master zone point
|
||||
// and grab the label
|
||||
addedPoints[pointi] =
|
||||
ref.setAction
|
||||
(
|
||||
polyAddPoint
|
||||
(
|
||||
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());
|
||||
|
||||
forAll(mf, facei)
|
||||
{
|
||||
label celli = mc[facei];
|
||||
label zoneI = mesh.cellZones().whichZone(celli);
|
||||
|
||||
addedCells[facei] =
|
||||
ref.setAction
|
||||
(
|
||||
polyAddCell
|
||||
(
|
||||
-1, // master point
|
||||
-1, // master edge
|
||||
mf[facei], // master face
|
||||
-1, // master cell
|
||||
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());
|
||||
|
||||
forAll(oldFace, pointi)
|
||||
{
|
||||
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]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
// Modify the faces from the master zone for the new neighbour
|
||||
|
||||
const faceList& faces = mesh.faces();
|
||||
|
||||
// Pout<< "mfFlip: " << mfFlip << endl;
|
||||
|
||||
forAll(mf, facei)
|
||||
{
|
||||
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]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
// 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]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
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]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
// 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 meshFaceZones& meshZones = 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;
|
||||
|
||||
forAll(curFaces, facei)
|
||||
{
|
||||
const label cf = curFaces[facei];
|
||||
|
||||
if (!mesh.isInternalFace(cf))
|
||||
{
|
||||
// Face not internal. Check if it is in the zone
|
||||
if (meshZones.whichZone(cf) != faceZoneID_.index())
|
||||
{
|
||||
// Found the face in a boundary patch which is not in zone
|
||||
patchID = mesh.boundaryMesh().whichPatch(cf);
|
||||
zoneID = mesh.faceZones().whichZone(cf);
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
// 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();
|
||||
|
||||
forAll(mc, celli)
|
||||
{
|
||||
const labelList& curFaces = cells[mc[celli]];
|
||||
|
||||
forAll(curFaces, facei)
|
||||
{
|
||||
// Check if the face belongs to the master zone; if not add it
|
||||
if (meshZones.whichZone(curFaces[facei]) != faceZoneID_.index())
|
||||
{
|
||||
masterCellFaceMap.insert(curFaces[facei]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Create the master layer point map
|
||||
Map<label> masterLayerPointMap(2*mp.size());
|
||||
|
||||
forAll(mp, pointi)
|
||||
{
|
||||
masterLayerPointMap.insert
|
||||
(
|
||||
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;
|
||||
|
||||
forAll(oldFace, pointi)
|
||||
{
|
||||
if (masterLayerPointMap.found(oldFace[pointi]))
|
||||
{
|
||||
changed = true;
|
||||
|
||||
newFace[pointi] = masterLayerPointMap.find(oldFace[pointi])();
|
||||
}
|
||||
else
|
||||
{
|
||||
newFace[pointi] = oldFace[pointi];
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,518 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Description
|
||||
Cell layer addition/removal mesh modifier
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "layerAdditionRemoval.H"
|
||||
#include "polyTopoChanger.H"
|
||||
#include "polyMesh.H"
|
||||
#include "Time.H"
|
||||
#include "primitiveMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(layerAdditionRemoval, 0);
|
||||
addToRunTimeSelectionTable
|
||||
(
|
||||
polyMeshModifier,
|
||||
layerAdditionRemoval,
|
||||
dictionary
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
const Foam::scalar Foam::layerAdditionRemoval::addDelta_ = 0.3;
|
||||
const Foam::scalar Foam::layerAdditionRemoval::removeDelta_ = 0.1;
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::layerAdditionRemoval::checkDefinition()
|
||||
{
|
||||
if (!faceZoneID_.active())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Master face zone named " << faceZoneID_.name()
|
||||
<< " cannot be found."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
minLayerThickness_ < vSmall
|
||||
|| maxLayerThickness_ < minLayerThickness_
|
||||
)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Incorrect layer thickness definition."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
label nFaces = topoChanger().mesh().faceZones()[faceZoneID_.index()].size();
|
||||
|
||||
reduce(nFaces, sumOp<label>());
|
||||
|
||||
if (nFaces == 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Face extrusion zone contains no faces. "
|
||||
<< "Please check your mesh definition."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Cell layer addition/removal object " << name() << " :" << nl
|
||||
<< " faceZoneID: " << faceZoneID_ << endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::layerAdditionRemoval::readOldThickness
|
||||
(
|
||||
const dictionary& dict
|
||||
)
|
||||
{
|
||||
return dict.lookupOrDefault("oldLayerThickness", -1.0);
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::clearAddressing() const
|
||||
{
|
||||
if (pointsPairingPtr_)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "layerAdditionRemoval::clearAddressing()" << nl
|
||||
<< " clearing pointsPairingPtr_" << endl;
|
||||
}
|
||||
|
||||
deleteDemandDrivenData(pointsPairingPtr_);
|
||||
}
|
||||
|
||||
if (facesPairingPtr_)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "layerAdditionRemoval::clearAddressing()" << nl
|
||||
<< " clearing facesPairingPtr_" << endl;
|
||||
}
|
||||
|
||||
deleteDemandDrivenData(facesPairingPtr_);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::layerAdditionRemoval::layerAdditionRemoval
|
||||
(
|
||||
const word& name,
|
||||
const label index,
|
||||
const polyTopoChanger& ptc,
|
||||
const word& zoneName,
|
||||
const scalar minThickness,
|
||||
const scalar maxThickness,
|
||||
const Switch thicknessFromVolume
|
||||
)
|
||||
:
|
||||
polyMeshModifier(name, index, ptc, true),
|
||||
faceZoneID_(zoneName, ptc.mesh().faceZones()),
|
||||
minLayerThickness_(minThickness),
|
||||
maxLayerThickness_(maxThickness),
|
||||
thicknessFromVolume_(thicknessFromVolume),
|
||||
oldLayerThickness_(-1.0),
|
||||
pointsPairingPtr_(nullptr),
|
||||
facesPairingPtr_(nullptr),
|
||||
triggerRemoval_(-1),
|
||||
triggerAddition_(-1)
|
||||
{
|
||||
checkDefinition();
|
||||
}
|
||||
|
||||
|
||||
Foam::layerAdditionRemoval::layerAdditionRemoval
|
||||
(
|
||||
const word& name,
|
||||
const dictionary& dict,
|
||||
const label index,
|
||||
const polyTopoChanger& ptc
|
||||
)
|
||||
:
|
||||
polyMeshModifier(name, index, ptc, Switch(dict.lookup("active"))),
|
||||
faceZoneID_(dict.lookup("faceZoneName"), ptc.mesh().faceZones()),
|
||||
minLayerThickness_(dict.lookup<scalar>("minLayerThickness")),
|
||||
maxLayerThickness_(dict.lookup<scalar>("maxLayerThickness")),
|
||||
thicknessFromVolume_
|
||||
(
|
||||
dict.lookupOrDefault<Switch>("thicknessFromVolume", true)
|
||||
),
|
||||
oldLayerThickness_(readOldThickness(dict)),
|
||||
pointsPairingPtr_(nullptr),
|
||||
facesPairingPtr_(nullptr),
|
||||
triggerRemoval_(-1),
|
||||
triggerAddition_(-1)
|
||||
{
|
||||
checkDefinition();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::layerAdditionRemoval::~layerAdditionRemoval()
|
||||
{
|
||||
clearAddressing();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
bool Foam::layerAdditionRemoval::changeTopology() const
|
||||
{
|
||||
// Protect from multiple calculation in the same time-step
|
||||
if (triggerRemoval_ > -1 || triggerAddition_ > -1)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
// Go through all the cells in the master layer and calculate
|
||||
// approximate layer thickness as the ratio of the cell volume and
|
||||
// face area in the face zone.
|
||||
// Layer addition:
|
||||
// When the max thickness exceeds the threshold, trigger refinement.
|
||||
// Layer removal:
|
||||
// When the min thickness falls below the threshold, trigger removal.
|
||||
|
||||
const polyMesh& mesh = topoChanger().mesh();
|
||||
|
||||
const faceZone& fz = mesh.faceZones()[faceZoneID_.index()];
|
||||
const labelList& mc = fz.masterCells();
|
||||
|
||||
const scalarField& V = mesh.cellVolumes();
|
||||
const vectorField& S = mesh.faceAreas();
|
||||
|
||||
if (min(V) < -vSmall)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "negative cell volume. Error in mesh motion before "
|
||||
<< "topological change.\n V: " << V
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
scalar avgDelta = 0;
|
||||
scalar minDelta = great;
|
||||
scalar maxDelta = 0;
|
||||
label nDelta = 0;
|
||||
|
||||
if (thicknessFromVolume_)
|
||||
{
|
||||
// Thickness calculated from cell volume/face area
|
||||
forAll(fz, facei)
|
||||
{
|
||||
scalar curDelta = V[mc[facei]]/mag(S[fz[facei]]);
|
||||
avgDelta += curDelta;
|
||||
minDelta = min(minDelta, curDelta);
|
||||
maxDelta = max(maxDelta, curDelta);
|
||||
}
|
||||
|
||||
nDelta = fz.size();
|
||||
}
|
||||
else
|
||||
{
|
||||
// Thickness calculated from edges on layer
|
||||
const Map<label>& meshZonePointMap = fz().meshPointMap();
|
||||
|
||||
// Edges with only one point on zone
|
||||
forAll(mc, facei)
|
||||
{
|
||||
const cell& cFaces = mesh.cells()[mc[facei]];
|
||||
const edgeList cellEdges(cFaces.edges(mesh.faces()));
|
||||
|
||||
forAll(cellEdges, i)
|
||||
{
|
||||
const edge& e = cellEdges[i];
|
||||
|
||||
if (meshZonePointMap.found(e[0]))
|
||||
{
|
||||
if (!meshZonePointMap.found(e[1]))
|
||||
{
|
||||
scalar curDelta = e.mag(mesh.points());
|
||||
avgDelta += curDelta;
|
||||
nDelta++;
|
||||
minDelta = min(minDelta, curDelta);
|
||||
maxDelta = max(maxDelta, curDelta);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (meshZonePointMap.found(e[1]))
|
||||
{
|
||||
scalar curDelta = e.mag(mesh.points());
|
||||
avgDelta += curDelta;
|
||||
nDelta++;
|
||||
minDelta = min(minDelta, curDelta);
|
||||
maxDelta = max(maxDelta, curDelta);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
reduce(minDelta, minOp<scalar>());
|
||||
reduce(maxDelta, maxOp<scalar>());
|
||||
reduce(avgDelta, sumOp<scalar>());
|
||||
reduce(nDelta, sumOp<label>());
|
||||
|
||||
avgDelta /= nDelta;
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "bool layerAdditionRemoval::changeTopology() const "
|
||||
<< " for object " << name() << " : " << nl
|
||||
<< "Layer thickness: min: " << minDelta
|
||||
<< " max: " << maxDelta << " avg: " << avgDelta
|
||||
<< " old thickness: " << oldLayerThickness_ << nl
|
||||
<< "Removal threshold: " << minLayerThickness_
|
||||
<< " addition threshold: " << maxLayerThickness_ << endl;
|
||||
}
|
||||
|
||||
bool topologicalChange = false;
|
||||
|
||||
// If the thickness is decreasing and crosses the min thickness,
|
||||
// trigger removal
|
||||
if (oldLayerThickness_ < 0)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "First step. No addition/removal" << endl;
|
||||
}
|
||||
|
||||
// No topological changes allowed before first mesh motion
|
||||
oldLayerThickness_ = avgDelta;
|
||||
|
||||
topologicalChange = false;
|
||||
}
|
||||
else if (avgDelta < oldLayerThickness_)
|
||||
{
|
||||
// Layers moving towards removal
|
||||
if (minDelta < minLayerThickness_)
|
||||
{
|
||||
// Check layer pairing
|
||||
if (setLayerPairing())
|
||||
{
|
||||
// A mesh layer detected. Check that collapse is valid
|
||||
if (validCollapse())
|
||||
{
|
||||
// At this point, info about moving the old mesh
|
||||
// in a way to collapse the cells in the removed
|
||||
// layer is available. Not sure what to do with
|
||||
// it.
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "bool layerAdditionRemoval::changeTopology() "
|
||||
<< " const for object " << name() << " : "
|
||||
<< "Triggering layer removal" << endl;
|
||||
}
|
||||
|
||||
triggerRemoval_ = mesh.time().timeIndex();
|
||||
|
||||
// Old thickness looses meaning.
|
||||
// Set it up to indicate layer removal
|
||||
oldLayerThickness_ = great;
|
||||
|
||||
topologicalChange = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
// No removal, clear addressing
|
||||
clearAddressing();
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
oldLayerThickness_ = avgDelta;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Layers moving towards addition
|
||||
if (maxDelta > maxLayerThickness_)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "bool layerAdditionRemoval::changeTopology() const "
|
||||
<< " for object " << name() << " : "
|
||||
<< "Triggering layer addition" << endl;
|
||||
}
|
||||
|
||||
triggerAddition_ = mesh.time().timeIndex();
|
||||
|
||||
// Old thickness looses meaning.
|
||||
// Set it up to indicate layer removal
|
||||
oldLayerThickness_ = 0;
|
||||
|
||||
topologicalChange = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
oldLayerThickness_ = avgDelta;
|
||||
}
|
||||
}
|
||||
|
||||
return topologicalChange;
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
|
||||
{
|
||||
// Insert the layer addition/removal instructions
|
||||
// into the topological change
|
||||
|
||||
if (triggerRemoval_ == topoChanger().mesh().time().timeIndex())
|
||||
{
|
||||
removeCellLayer(ref);
|
||||
|
||||
// Clear addressing. This also resets the addition/removal data
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
|
||||
<< "for object " << name() << " : "
|
||||
<< "Clearing addressing after layer removal" << endl;
|
||||
}
|
||||
|
||||
triggerRemoval_ = -1;
|
||||
clearAddressing();
|
||||
}
|
||||
|
||||
if (triggerAddition_ == topoChanger().mesh().time().timeIndex())
|
||||
{
|
||||
addCellLayer(ref);
|
||||
|
||||
// Clear addressing. This also resets the addition/removal data
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
|
||||
<< "for object " << name() << " : "
|
||||
<< "Clearing addressing after layer addition" << endl;
|
||||
}
|
||||
|
||||
triggerAddition_ = -1;
|
||||
clearAddressing();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::topoChange(const polyTopoChangeMap&)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "layerAdditionRemoval::topoChange(const polyTopoChangeMap&) "
|
||||
<< "for object " << name() << " : "
|
||||
<< "Clearing addressing on external request";
|
||||
|
||||
if (pointsPairingPtr_ || facesPairingPtr_)
|
||||
{
|
||||
Pout<< "Pointers set." << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
Pout<< "Pointers not set." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
// Mesh has changed topologically. Update local topological data
|
||||
faceZoneID_.update(topoChanger().mesh().faceZones());
|
||||
|
||||
clearAddressing();
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
|
||||
{
|
||||
if (t < vSmall || maxLayerThickness_ < t)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Incorrect layer thickness definition."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
minLayerThickness_ = t;
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
|
||||
{
|
||||
if (t < minLayerThickness_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Incorrect layer thickness definition."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
maxLayerThickness_ = t;
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::write(Ostream& os) const
|
||||
{
|
||||
os << nl << type() << nl
|
||||
<< name()<< nl
|
||||
<< faceZoneID_ << nl
|
||||
<< minLayerThickness_ << nl
|
||||
<< oldLayerThickness_ << nl
|
||||
<< maxLayerThickness_ << nl
|
||||
<< thicknessFromVolume_ << endl;
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
|
||||
{
|
||||
os << nl << name() << nl << token::BEGIN_BLOCK << nl
|
||||
<< " type " << type()
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " faceZoneName " << faceZoneID_.name()
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " minLayerThickness " << minLayerThickness_
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " maxLayerThickness " << maxLayerThickness_
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " thicknessFromVolume " << thicknessFromVolume_
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " oldLayerThickness " << oldLayerThickness_
|
||||
<< token::END_STATEMENT << nl
|
||||
<< " active " << active()
|
||||
<< token::END_STATEMENT << nl
|
||||
<< token::END_BLOCK << endl;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,234 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Class
|
||||
Foam::layerAdditionRemoval
|
||||
|
||||
Description
|
||||
Cell layer addition mesh modifier
|
||||
|
||||
SourceFiles
|
||||
layerAdditionRemoval.C
|
||||
addCellLayer.C
|
||||
removeCellLayer.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef layerAdditionRemoval_H
|
||||
#define layerAdditionRemoval_H
|
||||
|
||||
#include "polyMeshModifier.H"
|
||||
#include "primitiveFacePatch.H"
|
||||
#include "ZoneIDs.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class layerAdditionRemoval Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class layerAdditionRemoval
|
||||
:
|
||||
public polyMeshModifier
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Master face zone ID
|
||||
faceZoneID faceZoneID_;
|
||||
|
||||
//- Min thickness of extrusion layer. Triggers layer removal
|
||||
mutable scalar minLayerThickness_;
|
||||
|
||||
//- Max thickness of extrusion layer. Triggers layer addition
|
||||
mutable scalar maxLayerThickness_;
|
||||
|
||||
//- Switch to calculate thickness as volume/area
|
||||
// If false, thickness calculated from edges
|
||||
const bool thicknessFromVolume_;
|
||||
|
||||
//- Layer thickness from previous step
|
||||
// Used to decide whether to add or remove layers
|
||||
mutable scalar oldLayerThickness_;
|
||||
|
||||
//- Point pairing
|
||||
mutable labelList* pointsPairingPtr_;
|
||||
|
||||
//- Face pairing
|
||||
mutable labelList* facesPairingPtr_;
|
||||
|
||||
//- Layer removal trigger time index
|
||||
mutable label triggerRemoval_;
|
||||
|
||||
//- Layer addition trigger time index
|
||||
mutable label triggerAddition_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Check validity of construction data
|
||||
void checkDefinition();
|
||||
|
||||
|
||||
// Topological changes
|
||||
|
||||
//- Check for valid layer
|
||||
bool validCollapse() const;
|
||||
|
||||
//- Set layer pairing. Return true if a valid layer exists
|
||||
bool setLayerPairing() const;
|
||||
|
||||
//- Return points pairing in a layer (not automatic!)
|
||||
const labelList& pointsPairing() const;
|
||||
|
||||
//- Return faces pairing in a layer (not automatic!)
|
||||
const labelList& facesPairing() const;
|
||||
|
||||
//- Calculate the offset to the next layer
|
||||
tmp<vectorField> extrusionDir() const;
|
||||
|
||||
//- Add a layer of cells
|
||||
void addCellLayer(polyTopoChange&) const;
|
||||
|
||||
//- Remove a layer of cells
|
||||
void removeCellLayer(polyTopoChange&) const;
|
||||
|
||||
//- Clear addressing
|
||||
void clearAddressing() const;
|
||||
|
||||
|
||||
// Helpers
|
||||
|
||||
//- Optionally read old thickness
|
||||
static scalar readOldThickness(const dictionary&);
|
||||
|
||||
|
||||
// Static Data Members
|
||||
|
||||
//- Thickness insertion fraction for the pre-motion
|
||||
static const scalar addDelta_;
|
||||
|
||||
//- Thickness removal fraction for the cell collapse
|
||||
// Note: the cell will be collapsed to this relative
|
||||
// thickness before the layer is removed.
|
||||
static const scalar removeDelta_;
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("layerAdditionRemoval");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
layerAdditionRemoval
|
||||
(
|
||||
const word& name,
|
||||
const label index,
|
||||
const polyTopoChanger& ptc,
|
||||
const word& zoneName,
|
||||
const scalar minThickness,
|
||||
const scalar maxThickness,
|
||||
const Switch thicknessFromVolume = true
|
||||
);
|
||||
|
||||
//- Construct from dictionary
|
||||
layerAdditionRemoval
|
||||
(
|
||||
const word& name,
|
||||
const dictionary& dict,
|
||||
const label index,
|
||||
const polyTopoChanger& ptc
|
||||
);
|
||||
|
||||
//- Disallow default bitwise copy construction
|
||||
layerAdditionRemoval(const layerAdditionRemoval&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~layerAdditionRemoval();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Check for topology change
|
||||
virtual bool changeTopology() const;
|
||||
|
||||
//- Insert the layer addition/removal instructions
|
||||
// into the topological change
|
||||
virtual void setRefinement(polyTopoChange&) const;
|
||||
|
||||
//- Modify motion points to comply with the topological change
|
||||
virtual void modifyMotionPoints(pointField& motionPoints) const;
|
||||
|
||||
//- Force recalculation of locally stored data on topological change
|
||||
virtual void topoChange(const polyTopoChangeMap&);
|
||||
|
||||
|
||||
// Edit
|
||||
|
||||
//- Return min layer thickness which triggers removal
|
||||
scalar minLayerThickness() const
|
||||
{
|
||||
return minLayerThickness_;
|
||||
}
|
||||
|
||||
//- Set min layer thickness which triggers removal
|
||||
void setMinLayerThickness(const scalar t) const;
|
||||
|
||||
//- Return max layer thickness which triggers removal
|
||||
scalar maxLayerThickness() const
|
||||
{
|
||||
return maxLayerThickness_;
|
||||
}
|
||||
|
||||
//- Set max layer thickness which triggers removal
|
||||
void setMaxLayerThickness(const scalar t) const;
|
||||
|
||||
|
||||
//- Write
|
||||
virtual void write(Ostream&) const;
|
||||
|
||||
//- Write dictionary
|
||||
virtual void writeDict(Ostream&) const;
|
||||
|
||||
|
||||
// Member Operators
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const layerAdditionRemoval&) = delete;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,420 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "oppositeFace.H"
|
||||
#include "polyTopoChanger.H"
|
||||
#include "polyRemoveCell.H"
|
||||
#include "polyRemoveFace.H"
|
||||
#include "polyRemovePoint.H"
|
||||
#include "polyModifyFace.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
bool Foam::layerAdditionRemoval::validCollapse() const
|
||||
{
|
||||
// Check for valid layer collapse
|
||||
// - no boundary-to-boundary collapse
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Checking layer collapse for object " << name() << endl;
|
||||
}
|
||||
|
||||
// Grab the face collapse mapping
|
||||
const polyMesh& mesh = topoChanger().mesh();
|
||||
|
||||
const labelList& ftc = facesPairing();
|
||||
const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
|
||||
|
||||
label nBoundaryHits = 0;
|
||||
|
||||
forAll(mf, facei)
|
||||
{
|
||||
if
|
||||
(
|
||||
!mesh.isInternalFace(mf[facei])
|
||||
&& !mesh.isInternalFace(ftc[facei])
|
||||
)
|
||||
{
|
||||
nBoundaryHits++;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Finished checking layer collapse for object "
|
||||
<< name() <<". Number of boundary-on-boundary hits: "
|
||||
<< nBoundaryHits << endl;
|
||||
}
|
||||
|
||||
if (returnReduce(nBoundaryHits, sumOp<label>()) > 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::layerAdditionRemoval::removeCellLayer
|
||||
(
|
||||
polyTopoChange& ref
|
||||
) const
|
||||
{
|
||||
// Algorithm for layer removal. Second phase: topological change
|
||||
// 2) Go through all the faces of the master cell layer and remove
|
||||
// the ones that are not in the master face zone.
|
||||
//
|
||||
// 3) Grab all the faces coming out of points that are collapsed
|
||||
// and select the ones that are not marked for removal. Go
|
||||
// through the remaining faces and replace the point to remove by
|
||||
// the equivalent point in the master face zone.
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Removing the cell layer for object " << name() << endl;
|
||||
}
|
||||
|
||||
const polyMesh& mesh = topoChanger().mesh();
|
||||
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
|
||||
const labelList& ptc = pointsPairing();
|
||||
const labelList& ftc = facesPairing();
|
||||
|
||||
// Remove all the cells from the master layer
|
||||
const labelList& mc =
|
||||
mesh.faceZones()[faceZoneID_.index()].masterCells();
|
||||
|
||||
forAll(mc, facei)
|
||||
{
|
||||
label slaveSideCell = own[ftc[facei]];
|
||||
|
||||
if (mesh.isInternalFace(ftc[facei]) && slaveSideCell == mc[facei])
|
||||
{
|
||||
// Owner cell of the face is being removed.
|
||||
// Grab the neighbour instead
|
||||
slaveSideCell = nei[ftc[facei]];
|
||||
}
|
||||
|
||||
ref.setAction(polyRemoveCell(mc[facei], slaveSideCell));
|
||||
}
|
||||
|
||||
// Remove all the faces from the master layer cells which are not in
|
||||
// the master face layer
|
||||
labelHashSet facesToRemoveMap(mc.size()*primitiveMesh::facesPerCell_);
|
||||
|
||||
const cellList& cells = mesh.cells();
|
||||
|
||||
forAll(mc, celli)
|
||||
{
|
||||
const cell& curCell = cells[mc[celli]];
|
||||
|
||||
forAll(curCell, facei)
|
||||
{
|
||||
// Check if the face is in the master zone. If not, remove it
|
||||
if
|
||||
(
|
||||
mesh.faceZones().whichZone(curCell[facei])
|
||||
!= faceZoneID_.index()
|
||||
)
|
||||
{
|
||||
facesToRemoveMap.insert(curCell[facei]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
forAllConstIter(labelHashSet, facesToRemoveMap, iter)
|
||||
{
|
||||
ref.setAction(polyRemoveFace(iter.key()));
|
||||
}
|
||||
|
||||
// Remove all points that will be collapsed
|
||||
forAll(ptc, pointi)
|
||||
{
|
||||
ref.setAction(polyRemovePoint(ptc[pointi]));
|
||||
}
|
||||
|
||||
// Grab all faces coming off points to be deleted. If the face
|
||||
// has not been deleted, replace the removed point with the
|
||||
// equivalent from the master layer.
|
||||
|
||||
// Make a map of all point to be removed, giving the master point label
|
||||
// for each of them
|
||||
|
||||
Map<label> removedPointMap(2*ptc.size());
|
||||
|
||||
const labelList& meshPoints =
|
||||
mesh.faceZones()[faceZoneID_.index()]().meshPoints();
|
||||
|
||||
forAll(ptc, pointi)
|
||||
{
|
||||
removedPointMap.insert(ptc[pointi], meshPoints[pointi]);
|
||||
}
|
||||
|
||||
const labelListList& pf = mesh.pointFaces();
|
||||
|
||||
const faceList& faces = mesh.faces();
|
||||
|
||||
// Make a list of faces to be modified using the map to avoid duplicates
|
||||
labelHashSet facesToModify(ptc.size()*primitiveMesh::facesPerPoint_);
|
||||
|
||||
forAll(ptc, pointi)
|
||||
{
|
||||
const labelList& curFaces = pf[ptc[pointi]];
|
||||
|
||||
forAll(curFaces, facei)
|
||||
{
|
||||
if (!facesToRemoveMap.found(curFaces[facei]))
|
||||
{
|
||||
facesToModify.insert(curFaces[facei]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
labelList ftm = facesToModify.toc();
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "faces to modify: " << ftm << endl;
|
||||
}
|
||||
|
||||
forAll(ftm, facei)
|
||||
{
|
||||
// For every face to modify, copy the face and re-map the vertices.
|
||||
// It is known all the faces will be changed since they hang off
|
||||
// re-mapped vertices
|
||||
label curFaceID = ftm[facei];
|
||||
|
||||
face newFace(faces[curFaceID]);
|
||||
|
||||
forAll(newFace, pointi)
|
||||
{
|
||||
Map<label>::iterator rpmIter =
|
||||
removedPointMap.find(newFace[pointi]);
|
||||
|
||||
if (rpmIter != removedPointMap.end())
|
||||
{
|
||||
// Point mapped. Replace it
|
||||
newFace[pointi] = rpmIter();
|
||||
}
|
||||
}
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "face label: " << curFaceID
|
||||
<< " old face: " << faces[curFaceID]
|
||||
<< " new face: " << newFace << endl;
|
||||
}
|
||||
|
||||
// Get face zone and its flip
|
||||
label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
|
||||
bool modifiedFaceZoneFlip = false;
|
||||
|
||||
if (modifiedFaceZone >= 0)
|
||||
{
|
||||
const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
|
||||
modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
|
||||
}
|
||||
|
||||
label newNeighbour = -1;
|
||||
|
||||
if (curFaceID < mesh.nInternalFaces())
|
||||
{
|
||||
newNeighbour = nei[curFaceID];
|
||||
}
|
||||
|
||||
// Modify the face
|
||||
ref.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
newFace, // modified face
|
||||
curFaceID, // label of face being modified
|
||||
own[curFaceID], // owner
|
||||
newNeighbour, // neighbour
|
||||
false, // face flip
|
||||
mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
|
||||
false, // remove from zone
|
||||
modifiedFaceZone, // zone for face
|
||||
modifiedFaceZoneFlip // face flip in zone
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
// Modify the faces in the master layer to point past the removed cells
|
||||
|
||||
const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
|
||||
const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
|
||||
|
||||
forAll(mf, facei)
|
||||
{
|
||||
// Grab the owner and neighbour of the faces to be collapsed and get rid
|
||||
// of the cell to be removed
|
||||
label masterSideCell = own[mf[facei]];
|
||||
|
||||
if (masterSideCell == mc[facei])
|
||||
{
|
||||
if (mesh.isInternalFace(mf[facei]))
|
||||
{
|
||||
// Owner cell of the face is being removed.
|
||||
// Grab the neighbour instead
|
||||
masterSideCell = nei[mf[facei]];
|
||||
}
|
||||
else
|
||||
{
|
||||
masterSideCell = -1;
|
||||
}
|
||||
}
|
||||
|
||||
label slaveSideCell = own[ftc[facei]];
|
||||
|
||||
if (slaveSideCell == mc[facei])
|
||||
{
|
||||
if (mesh.isInternalFace(ftc[facei]))
|
||||
{
|
||||
// Owner cell of the face is being removed.
|
||||
// Grab the neighbour instead
|
||||
slaveSideCell = nei[ftc[facei]];
|
||||
}
|
||||
else
|
||||
{
|
||||
slaveSideCell = -1;
|
||||
}
|
||||
}
|
||||
|
||||
// Find out if the face needs to be flipped
|
||||
label newOwner = -1;
|
||||
label newNeighbour = -1;
|
||||
bool flipFace = false;
|
||||
label newPatchID = -1;
|
||||
label newZoneID = -1;
|
||||
|
||||
// Is any of the faces a boundary face? If so, grab the patch
|
||||
// A boundary-to-boundary collapse is checked for in validCollapse()
|
||||
// and cannot happen here.
|
||||
|
||||
if (!mesh.isInternalFace(mf[facei]))
|
||||
{
|
||||
// Master is the boundary face: it gets a new owner but no flip
|
||||
newOwner = slaveSideCell;
|
||||
newNeighbour = -1;
|
||||
flipFace = false;
|
||||
newPatchID = mesh.boundaryMesh().whichPatch(mf[facei]);
|
||||
newZoneID = mesh.faceZones().whichZone(mf[facei]);
|
||||
}
|
||||
else if (!mesh.isInternalFace(ftc[facei]))
|
||||
{
|
||||
// Slave is the boundary face: grab its patch
|
||||
newOwner = slaveSideCell;
|
||||
newNeighbour = -1;
|
||||
|
||||
// Find out if the face flip is necessary
|
||||
if (own[mf[facei]] == slaveSideCell)
|
||||
{
|
||||
flipFace = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
flipFace = true;
|
||||
}
|
||||
|
||||
newPatchID = mesh.boundaryMesh().whichPatch(ftc[facei]);
|
||||
|
||||
// The zone of the master face is preserved
|
||||
newZoneID = mesh.faceZones().whichZone(mf[facei]);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Both faces are internal: flip is decided based on which of the
|
||||
// new cells around it has a lower label
|
||||
newOwner = min(masterSideCell, slaveSideCell);
|
||||
newNeighbour = max(masterSideCell, slaveSideCell);
|
||||
|
||||
if (newOwner == own[mf[facei]] || newNeighbour == nei[mf[facei]])
|
||||
{
|
||||
flipFace = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
flipFace = true;
|
||||
}
|
||||
|
||||
newPatchID = -1;
|
||||
|
||||
// The zone of the master face is preserved
|
||||
newZoneID = mesh.faceZones().whichZone(mf[facei]);
|
||||
}
|
||||
|
||||
// Modify the face and flip if necessary
|
||||
face newFace = faces[mf[facei]];
|
||||
bool zoneFlip = mfFlip[facei];
|
||||
|
||||
if (flipFace)
|
||||
{
|
||||
newFace.flip();
|
||||
zoneFlip = !zoneFlip;
|
||||
}
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "Modifying face " << mf[facei]
|
||||
<< " newFace: " << newFace << nl
|
||||
<< " newOwner: " << newOwner
|
||||
<< " newNeighbour: " << newNeighbour
|
||||
<< " flipFace: " << flipFace
|
||||
<< " newPatchID: " << newPatchID
|
||||
<< " newZoneID: " << newZoneID << nl
|
||||
<< " oldOwn: " << own[mf[facei]]
|
||||
<< " oldNei: " << nei[mf[facei]] << endl;
|
||||
}
|
||||
|
||||
ref.setAction
|
||||
(
|
||||
polyModifyFace
|
||||
(
|
||||
newFace, // modified face
|
||||
mf[facei], // label of face being modified
|
||||
newOwner, // owner
|
||||
newNeighbour, // neighbour
|
||||
flipFace, // flip
|
||||
newPatchID, // patch for face
|
||||
false, // remove from zone
|
||||
newZoneID, // new zone
|
||||
zoneFlip // face flip in zone
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,232 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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/>.
|
||||
|
||||
Description
|
||||
Remove a layer of cells and prepare addressing data
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "layerAdditionRemoval.H"
|
||||
#include "polyMesh.H"
|
||||
#include "primitiveMesh.H"
|
||||
#include "polyTopoChange.H"
|
||||
#include "oppositeFace.H"
|
||||
#include "polyTopoChanger.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
bool Foam::layerAdditionRemoval::setLayerPairing() const
|
||||
{
|
||||
// Note:
|
||||
// This is also the most complex part of the topological change.
|
||||
// Therefore it will be calculated here and stored as temporary
|
||||
// data until the actual topological change, after which it will
|
||||
// be cleared.
|
||||
|
||||
// Algorithm for point collapse
|
||||
// 1) Go through the master cell layer and for every face of
|
||||
// the face zone find the opposite face in the master cell.
|
||||
// Check the direction of the opposite face and adjust as
|
||||
// necessary. Check other faces to find an edge defining
|
||||
// relative orientation of the two faces and adjust the face
|
||||
// as necessary. Once the face is adjusted, record the
|
||||
// addressing between the master and slave vertex layer.
|
||||
|
||||
const polyMesh& mesh = topoChanger().mesh();
|
||||
|
||||
const labelList& mc =
|
||||
mesh.faceZones()[faceZoneID_.index()].masterCells();
|
||||
|
||||
const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
|
||||
|
||||
const boolList& mfFlip =
|
||||
mesh.faceZones()[faceZoneID_.index()].flipMap();
|
||||
|
||||
const faceList& faces = mesh.faces();
|
||||
const cellList& cells = mesh.cells();
|
||||
|
||||
// Grab the local faces from the master zone
|
||||
const faceList& mlf =
|
||||
mesh.faceZones()[faceZoneID_.index()]().localFaces();
|
||||
|
||||
const labelList& meshPoints =
|
||||
mesh.faceZones()[faceZoneID_.index()]().meshPoints();
|
||||
|
||||
// Create a list of points to collapse for every point of
|
||||
// the master patch
|
||||
if (pointsPairingPtr_ || facesPairingPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Problem with layer pairing data"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
pointsPairingPtr_ = new labelList(meshPoints.size(), -1);
|
||||
labelList& ptc = *pointsPairingPtr_;
|
||||
|
||||
facesPairingPtr_ = new labelList(mf.size(), -1);
|
||||
labelList& ftc = *facesPairingPtr_;
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "meshPoints: " << meshPoints << nl
|
||||
<< "localPoints: "
|
||||
<< mesh.faceZones()[faceZoneID_.index()]().localPoints()
|
||||
<< endl;
|
||||
}
|
||||
|
||||
// For all faces, create the mapping
|
||||
label nPointErrors = 0;
|
||||
label nFaceErrors = 0;
|
||||
|
||||
forAll(mf, facei)
|
||||
{
|
||||
// Get the local master face
|
||||
face curLocalFace = mlf[facei];
|
||||
|
||||
// Flip face based on flip index to recover original orientation
|
||||
if (mfFlip[facei])
|
||||
{
|
||||
curLocalFace.flip();
|
||||
}
|
||||
|
||||
// Get the opposing face from the master cell
|
||||
oppositeFace lidFace =
|
||||
cells[mc[facei]].opposingFace(mf[facei], faces);
|
||||
|
||||
if (!lidFace.found())
|
||||
{
|
||||
// This is not a valid layer; cannot continue
|
||||
nFaceErrors++;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "curMasterFace: " << faces[mf[facei]] << nl
|
||||
<< "cell shape: " << mesh.cellShapes()[mc[facei]] << nl
|
||||
<< "curLocalFace: " << curLocalFace << nl
|
||||
<< "lidFace: " << lidFace
|
||||
<< " master index: " << lidFace.masterIndex()
|
||||
<< " oppositeIndex: " << lidFace.oppositeIndex() << endl;
|
||||
}
|
||||
|
||||
// Grab the opposite face for face collapse addressing
|
||||
ftc[facei] = lidFace.oppositeIndex();
|
||||
|
||||
// Using the local face insert the points into the lid list
|
||||
forAll(curLocalFace, pointi)
|
||||
{
|
||||
const label clp = curLocalFace[pointi];
|
||||
|
||||
if (ptc[clp] == -1)
|
||||
{
|
||||
// Point not mapped yet. Insert the label
|
||||
ptc[clp] = lidFace[pointi];
|
||||
}
|
||||
else
|
||||
{
|
||||
// Point mapped from some other face. Check the label
|
||||
if (ptc[clp] != lidFace[pointi])
|
||||
{
|
||||
nPointErrors++;
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
Pout<< "Topological error in cell layer pairing. "
|
||||
<< "This mesh is either topologically incorrect "
|
||||
<< "or the master face layer is not defined "
|
||||
<< "consistently. Please check the "
|
||||
<< "face zone flip map." << nl
|
||||
<< "First index: " << ptc[clp]
|
||||
<< " new index: " << lidFace[pointi] << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
reduce(nPointErrors, sumOp<label>());
|
||||
reduce(nFaceErrors, sumOp<label>());
|
||||
|
||||
if (nPointErrors > 0 || nFaceErrors > 0)
|
||||
{
|
||||
clearAddressing();
|
||||
|
||||
return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Valid layer
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
const Foam::labelList& Foam::layerAdditionRemoval::pointsPairing() const
|
||||
{
|
||||
if (!pointsPairingPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Problem with layer pairing data for object " << name()
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
return *pointsPairingPtr_;
|
||||
}
|
||||
|
||||
const Foam::labelList& Foam::layerAdditionRemoval::facesPairing() const
|
||||
{
|
||||
if (!facesPairingPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Problem with layer pairing data for object " << name()
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
return *facesPairingPtr_;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::layerAdditionRemoval::modifyMotionPoints
|
||||
(
|
||||
pointField& motionPoints
|
||||
) const
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "void layerAdditionRemoval::modifyMotionPoints("
|
||||
<< "pointField& motionPoints) const for object "
|
||||
<< name() << " : ";
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "No motion point adjustment" << endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user