c9333a5ac8
- provide a plain stream() method on messageStream to reduce reliance
on casting operators and slightly opaque operator()() calls etc
- support alternative stream for messageStream serial output.
This can be used to support local redirection of output.
For example,
refPtr<OFstream> logging; // or autoPtr, unique_ptr etc
// Later...
Info.stream(logging.get())
<< "Detailed output ..." << endl;
This will use the stdout semantics in the normal case, or allow
redirection to an output file if a target output stream is defined,
but still effectively use /dev/null on non-master processes.
This is mostly the same as this ternary
(logging ? *logging : Info())
except that the ternary could be incorrect on sub-processes,
requires more typing etc.
ENH: use case-relative names of dictionary, IOstream for FatalIOError
- normally yields more easily understandable information
5049 lines
146 KiB
C
5049 lines
146 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2015 OpenFOAM Foundation
|
|
Copyright (C) 2015-2021 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Description
|
|
All to do with adding cell layers
|
|
|
|
\*----------------------------------------------------------------------------*/
|
|
|
|
#include "snappyLayerDriver.H"
|
|
#include "fvMesh.H"
|
|
#include "Time.H"
|
|
#include "meshRefinement.H"
|
|
#include "removePoints.H"
|
|
#include "pointFields.H"
|
|
#include "motionSmoother.H"
|
|
#include "unitConversion.H"
|
|
#include "pointSet.H"
|
|
#include "faceSet.H"
|
|
#include "cellSet.H"
|
|
#include "polyTopoChange.H"
|
|
#include "mapPolyMesh.H"
|
|
#include "addPatchCellLayer.H"
|
|
#include "mapDistributePolyMesh.H"
|
|
#include "OBJstream.H"
|
|
#include "layerParameters.H"
|
|
#include "combineFaces.H"
|
|
#include "IOmanip.H"
|
|
#include "globalIndex.H"
|
|
#include "DynamicField.H"
|
|
#include "PatchTools.H"
|
|
#include "slipPointPatchFields.H"
|
|
#include "fixedValuePointPatchFields.H"
|
|
#include "zeroFixedValuePointPatchFields.H"
|
|
#include "calculatedPointPatchFields.H"
|
|
#include "cyclicSlipPointPatchFields.H"
|
|
#include "fixedValueFvPatchFields.H"
|
|
#include "localPointRegion.H"
|
|
#include "externalDisplacementMeshMover.H"
|
|
#include "scalarIOField.H"
|
|
#include "profiling.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
|
|
defineTypeNameAndDebug(snappyLayerDriver, 0);
|
|
|
|
} // End namespace Foam
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
// For debugging: Dump displacement to .obj files
|
|
void Foam::snappyLayerDriver::dumpDisplacement
|
|
(
|
|
const fileName& prefix,
|
|
const indirectPrimitivePatch& pp,
|
|
const vectorField& patchDisp,
|
|
const List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
OBJstream dispStr(prefix + "_disp.obj");
|
|
Info<< "Writing all displacements to " << dispStr.name() << endl;
|
|
|
|
forAll(patchDisp, patchPointi)
|
|
{
|
|
const point& pt = pp.localPoints()[patchPointi];
|
|
dispStr.write(linePointRef(pt, pt + patchDisp[patchPointi]));
|
|
}
|
|
|
|
|
|
OBJstream illStr(prefix + "_illegal.obj");
|
|
Info<< "Writing invalid displacements to " << illStr.name() << endl;
|
|
|
|
forAll(patchDisp, patchPointi)
|
|
{
|
|
if (extrudeStatus[patchPointi] != EXTRUDE)
|
|
{
|
|
const point& pt = pp.localPoints()[patchPointi];
|
|
illStr.write(linePointRef(pt, pt + patchDisp[patchPointi]));
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::tmp<Foam::scalarField> Foam::snappyLayerDriver::avgPointData
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const scalarField& pointFld
|
|
)
|
|
{
|
|
tmp<scalarField> tfaceFld(new scalarField(pp.size(), Zero));
|
|
scalarField& faceFld = tfaceFld.ref();
|
|
|
|
forAll(pp.localFaces(), facei)
|
|
{
|
|
const face& f = pp.localFaces()[facei];
|
|
if (f.size())
|
|
{
|
|
forAll(f, fp)
|
|
{
|
|
faceFld[facei] += pointFld[f[fp]];
|
|
}
|
|
faceFld[facei] /= f.size();
|
|
}
|
|
}
|
|
return tfaceFld;
|
|
}
|
|
|
|
|
|
// Check that primitivePatch is not multiply connected. Collect non-manifold
|
|
// points in pointSet.
|
|
void Foam::snappyLayerDriver::checkManifold
|
|
(
|
|
const indirectPrimitivePatch& fp,
|
|
pointSet& nonManifoldPoints
|
|
)
|
|
{
|
|
// Check for non-manifold points (surface pinched at point)
|
|
fp.checkPointManifold(false, &nonManifoldPoints);
|
|
|
|
// Check for edge-faces (surface pinched at edge)
|
|
const labelListList& edgeFaces = fp.edgeFaces();
|
|
|
|
forAll(edgeFaces, edgei)
|
|
{
|
|
const labelList& eFaces = edgeFaces[edgei];
|
|
|
|
if (eFaces.size() > 2)
|
|
{
|
|
const edge& e = fp.edges()[edgei];
|
|
|
|
nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
|
|
nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::checkMeshManifold() const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
Info<< nl << "Checking mesh manifoldness ..." << endl;
|
|
|
|
pointSet nonManifoldPoints
|
|
(
|
|
mesh,
|
|
"nonManifoldPoints",
|
|
mesh.nPoints() / 100
|
|
);
|
|
|
|
// Build primitivePatch out of faces and check it for problems.
|
|
checkManifold
|
|
(
|
|
indirectPrimitivePatch
|
|
(
|
|
IndirectList<face>
|
|
(
|
|
mesh.faces(),
|
|
identity(mesh.boundaryMesh().range()) // All outside faces
|
|
),
|
|
mesh.points()
|
|
),
|
|
nonManifoldPoints
|
|
);
|
|
|
|
label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
|
|
|
|
if (nNonManif > 0)
|
|
{
|
|
Info<< "Outside of mesh is multiply connected across edges or"
|
|
<< " points." << nl
|
|
<< "This is not a fatal error but might cause some unexpected"
|
|
<< " behaviour." << nl
|
|
//<< "Writing " << nNonManif
|
|
//<< " points where this happens to pointSet "
|
|
//<< nonManifoldPoints.name()
|
|
<< endl;
|
|
|
|
//nonManifoldPoints.instance() = meshRefiner_.timeName();
|
|
//nonManifoldPoints.write();
|
|
}
|
|
Info<< endl;
|
|
}
|
|
|
|
|
|
|
|
// Unset extrusion on point. Returns true if anything unset.
|
|
bool Foam::snappyLayerDriver::unmarkExtrusion
|
|
(
|
|
const label patchPointi,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
if (extrudeStatus[patchPointi] == EXTRUDE)
|
|
{
|
|
extrudeStatus[patchPointi] = NOEXTRUDE;
|
|
patchNLayers[patchPointi] = 0;
|
|
patchDisp[patchPointi] = Zero;
|
|
return true;
|
|
}
|
|
else if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
|
|
{
|
|
extrudeStatus[patchPointi] = NOEXTRUDE;
|
|
patchNLayers[patchPointi] = 0;
|
|
patchDisp[patchPointi] = Zero;
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
// Unset extrusion on face. Returns true if anything unset.
|
|
bool Foam::snappyLayerDriver::unmarkExtrusion
|
|
(
|
|
const face& localFace,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
bool unextruded = false;
|
|
|
|
forAll(localFace, fp)
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
localFace[fp],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
unextruded = true;
|
|
}
|
|
}
|
|
return unextruded;
|
|
}
|
|
|
|
|
|
Foam::label Foam::snappyLayerDriver::constrainFp(const label sz, const label fp)
|
|
{
|
|
if (fp >= sz)
|
|
{
|
|
return 0;
|
|
}
|
|
else if (fp < 0)
|
|
{
|
|
return sz-1;
|
|
}
|
|
else
|
|
{
|
|
return fp;
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::countCommonPoints
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const label facei,
|
|
|
|
Map<label>& nCommonPoints
|
|
) const
|
|
{
|
|
const faceList& localFaces = pp.localFaces();
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
|
|
const face& f = localFaces[facei];
|
|
|
|
nCommonPoints.clear();
|
|
|
|
forAll(f, fp)
|
|
{
|
|
label pointi = f[fp];
|
|
const labelList& pFaces = pointFaces[pointi];
|
|
|
|
forAll(pFaces, pFacei)
|
|
{
|
|
label nbFacei = pFaces[pFacei];
|
|
|
|
if (facei < nbFacei)
|
|
{
|
|
// Only check once for each combination of two faces.
|
|
++(nCommonPoints(nbFacei, 0));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
bool Foam::snappyLayerDriver::checkCommonOrder
|
|
(
|
|
const label nCommon,
|
|
const face& curFace,
|
|
const face& nbFace
|
|
) const
|
|
{
|
|
forAll(curFace, fp)
|
|
{
|
|
// Get the index in the neighbouring face shared with curFace
|
|
const label nb = nbFace.find(curFace[fp]);
|
|
|
|
if (nb != -1)
|
|
{
|
|
|
|
// Check the whole face from nb onwards for shared vertices
|
|
// with neighbouring face. Rule is that any shared vertices
|
|
// should be consecutive on both faces i.e. if they are
|
|
// vertices fp,fp+1,fp+2 on one face they should be
|
|
// vertices nb, nb+1, nb+2 (or nb+2, nb+1, nb) on the
|
|
// other face.
|
|
|
|
|
|
// Vertices before and after on curFace
|
|
label fpPlus1 = curFace.fcIndex(fp);
|
|
label fpMin1 = curFace.rcIndex(fp);
|
|
|
|
// Vertices before and after on nbFace
|
|
label nbPlus1 = nbFace.fcIndex(nb);
|
|
label nbMin1 = nbFace.rcIndex(nb);
|
|
|
|
// Find order of walking by comparing next points on both
|
|
// faces.
|
|
label curInc = labelMax;
|
|
label nbInc = labelMax;
|
|
|
|
if (nbFace[nbPlus1] == curFace[fpPlus1])
|
|
{
|
|
curInc = 1;
|
|
nbInc = 1;
|
|
}
|
|
else if (nbFace[nbPlus1] == curFace[fpMin1])
|
|
{
|
|
curInc = -1;
|
|
nbInc = 1;
|
|
}
|
|
else if (nbFace[nbMin1] == curFace[fpMin1])
|
|
{
|
|
curInc = -1;
|
|
nbInc = -1;
|
|
}
|
|
else
|
|
{
|
|
curInc = 1;
|
|
nbInc = -1;
|
|
}
|
|
|
|
|
|
// Pass1: loop until start of common vertices found.
|
|
label curNb = nb;
|
|
label curFp = fp;
|
|
|
|
do
|
|
{
|
|
curFp = constrainFp(curFace.size(), curFp+curInc);
|
|
curNb = constrainFp(nbFace.size(), curNb+nbInc);
|
|
} while (curFace[curFp] == nbFace[curNb]);
|
|
|
|
// Pass2: check equality walking from curFp, curNb
|
|
// in opposite order.
|
|
|
|
curInc = -curInc;
|
|
nbInc = -nbInc;
|
|
|
|
for (label commonI = 0; commonI < nCommon; commonI++)
|
|
{
|
|
curFp = constrainFp(curFace.size(), curFp+curInc);
|
|
curNb = constrainFp(nbFace.size(), curNb+nbInc);
|
|
|
|
if (curFace[curFp] != nbFace[curNb])
|
|
{
|
|
// Error: gap in string of connected vertices
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// Done the curFace - nbFace combination.
|
|
break;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::checkCommonOrder
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const label facei,
|
|
const Map<label>& nCommonPoints,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
forAllConstIters(nCommonPoints, iter)
|
|
{
|
|
const label nbFacei = iter.key();
|
|
const label nCommon = iter.val();
|
|
|
|
const face& curFace = pp[facei];
|
|
const face& nbFace = pp[nbFacei];
|
|
|
|
if
|
|
(
|
|
nCommon >= 2
|
|
&& nCommon != nbFace.size()
|
|
&& nCommon != curFace.size()
|
|
)
|
|
{
|
|
bool stringOk = checkCommonOrder(nCommon, curFace, nbFace);
|
|
|
|
if (!stringOk)
|
|
{
|
|
// Note: unmark whole face or just the common points?
|
|
// For now unmark the whole face
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[facei],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[nbFacei],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::handleNonStringConnected
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
// Detect faces which are connected on non-consecutive vertices.
|
|
// This is the "<<Number of faces with non-consecutive shared points"
|
|
// warning from checkMesh. These faces cannot be extruded so
|
|
// there is no need to even attempt it.
|
|
|
|
List<extrudeMode> oldExtrudeStatus;
|
|
autoPtr<OBJstream> str;
|
|
if (debug&meshRefinement::LAYERINFO)
|
|
{
|
|
oldExtrudeStatus = extrudeStatus;
|
|
str.reset
|
|
(
|
|
new OBJstream
|
|
(
|
|
meshRefiner_.mesh().time().path()
|
|
/"nonStringConnected.obj"
|
|
)
|
|
);
|
|
Pout<< "Dumping string edges to " << str().name();
|
|
}
|
|
|
|
|
|
// 1) Local
|
|
Map<label> nCommonPoints(128);
|
|
|
|
forAll(pp, facei)
|
|
{
|
|
countCommonPoints(pp, facei, nCommonPoints);
|
|
|
|
// Faces share pointi. Find any more shared points
|
|
// and if not in single string unmark all. See
|
|
// primitiveMesh::checkCommonOrder
|
|
checkCommonOrder
|
|
(
|
|
pp,
|
|
facei,
|
|
nCommonPoints,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
|
|
// 2) TDB. Other face remote
|
|
|
|
|
|
|
|
if (debug&meshRefinement::LAYERINFO)
|
|
{
|
|
forAll(extrudeStatus, pointi)
|
|
{
|
|
if (extrudeStatus[pointi] != oldExtrudeStatus[pointi])
|
|
{
|
|
str().write
|
|
(
|
|
meshRefiner_.mesh().points()[pp.meshPoints()[pointi]]
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// No extrusion at non-manifold points.
|
|
void Foam::snappyLayerDriver::handleNonManifolds
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& meshEdges,
|
|
const labelListList& edgeGlobalFaces,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
Info<< nl << "Handling non-manifold points ..." << endl;
|
|
|
|
// Detect non-manifold points
|
|
Info<< nl << "Checking patch manifoldness ..." << endl;
|
|
|
|
pointSet nonManifoldPoints(mesh, "nonManifoldPoints", pp.nPoints());
|
|
|
|
// 1. Local check. Note that we do not check for e.g. two patch faces
|
|
// being connected via a point since their connection might be ok
|
|
// through a coupled patch. The ultimate is to do a proper point-face
|
|
// walk which is done when actually duplicating the points. Here we just
|
|
// do the obvious problems.
|
|
{
|
|
// Check for edge-faces (surface pinched at edge)
|
|
const labelListList& edgeFaces = pp.edgeFaces();
|
|
|
|
forAll(edgeFaces, edgei)
|
|
{
|
|
const labelList& eFaces = edgeFaces[edgei];
|
|
if (eFaces.size() > 2)
|
|
{
|
|
const edge& e = pp.edges()[edgei];
|
|
nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
|
|
nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
|
|
}
|
|
}
|
|
}
|
|
|
|
// 2. Remote check for boundary edges on coupled boundaries
|
|
forAll(edgeGlobalFaces, edgei)
|
|
{
|
|
if (edgeGlobalFaces[edgei].size() > 2)
|
|
{
|
|
// So boundary edges that are connected to more than 2 processors
|
|
// i.e. a non-manifold edge which is exactly on a processor
|
|
// boundary.
|
|
const edge& e = pp.edges()[edgei];
|
|
nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
|
|
nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
|
|
}
|
|
}
|
|
|
|
|
|
label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
|
|
|
|
Info<< "Outside of local patch is multiply connected across edges or"
|
|
<< " points at " << nNonManif << " points." << endl;
|
|
|
|
if (nNonManif > 0)
|
|
{
|
|
// Make sure all processors use the same information. The edge might
|
|
// not exist locally but remotely there might be a problem with this
|
|
// edge.
|
|
nonManifoldPoints.sync(mesh);
|
|
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
|
|
forAll(meshPoints, patchPointi)
|
|
{
|
|
if (nonManifoldPoints.found(meshPoints[patchPointi]))
|
|
{
|
|
unmarkExtrusion
|
|
(
|
|
patchPointi,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero for all " << nNonManif
|
|
<< " non-manifold points" << endl;
|
|
|
|
|
|
|
|
// 4. Check for extrusion of baffles i.e. all edges of a face having the
|
|
// same two neighbouring faces (one of which is the current face).
|
|
// Note: this is detected locally already before - this test is for the
|
|
// extremely rare occurrence where the baffle faces are on
|
|
// different processors.
|
|
{
|
|
label nBaffleFaces = 0;
|
|
|
|
const labelListList& faceEdges = pp.faceEdges();
|
|
forAll(pp, facei)
|
|
{
|
|
const labelList& fEdges = faceEdges[facei];
|
|
|
|
const labelList& globFaces0 = edgeGlobalFaces[fEdges[0]];
|
|
if (globFaces0.size() == 2)
|
|
{
|
|
const edge e0(globFaces0[0], globFaces0[1]);
|
|
bool isBaffle = true;
|
|
for (label fp = 1; fp < fEdges.size(); fp++)
|
|
{
|
|
const labelList& globFaces = edgeGlobalFaces[fEdges[fp]];
|
|
if
|
|
(
|
|
(globFaces.size() != 2)
|
|
|| (edge(globFaces[0], globFaces[1]) != e0)
|
|
)
|
|
{
|
|
isBaffle = false;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (isBaffle)
|
|
{
|
|
bool unextrude = unmarkExtrusion
|
|
(
|
|
pp.localFaces()[facei],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
if (unextrude)
|
|
{
|
|
//Pout<< "Detected extrusion of baffle face "
|
|
// << pp.faceCentres()[facei]
|
|
// << " since all edges have the same neighbours "
|
|
// << e0 << endl;
|
|
|
|
nBaffleFaces++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
reduce(nBaffleFaces, sumOp<label>());
|
|
|
|
if (nBaffleFaces)
|
|
{
|
|
Info<< "Set displacement to zero for all points on " << nBaffleFaces
|
|
<< " baffle faces" << endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Parallel feature edge detection. Assumes non-manifold edges already handled.
|
|
void Foam::snappyLayerDriver::handleFeatureAngle
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& meshEdges,
|
|
const scalar minAngle,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
const scalar minCos = Foam::cos(degToRad(minAngle));
|
|
|
|
Info<< nl << "Handling feature edges (angle < " << minAngle
|
|
<< ") ..." << endl;
|
|
|
|
if (minCos < 1-SMALL && minCos > -1+SMALL)
|
|
{
|
|
// Normal component of normals of connected faces.
|
|
vectorField edgeNormal(mesh.nEdges(), point::max);
|
|
|
|
const labelListList& edgeFaces = pp.edgeFaces();
|
|
|
|
forAll(edgeFaces, edgei)
|
|
{
|
|
const labelList& eFaces = pp.edgeFaces()[edgei];
|
|
|
|
label meshEdgei = meshEdges[edgei];
|
|
|
|
forAll(eFaces, i)
|
|
{
|
|
nomalsCombine()
|
|
(
|
|
edgeNormal[meshEdgei],
|
|
pp.faceNormals()[eFaces[i]]
|
|
);
|
|
}
|
|
}
|
|
|
|
syncTools::syncEdgeList
|
|
(
|
|
mesh,
|
|
edgeNormal,
|
|
nomalsCombine(),
|
|
point::max // null value
|
|
);
|
|
|
|
autoPtr<OBJstream> str;
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
str.reset
|
|
(
|
|
new OBJstream
|
|
(
|
|
mesh.time().path()
|
|
/ "featureEdges_"
|
|
+ meshRefiner_.timeName()
|
|
+ ".obj"
|
|
)
|
|
);
|
|
Info<< "Writing feature edges to " << str().name() << endl;
|
|
}
|
|
|
|
label nFeats = 0;
|
|
|
|
// Now on coupled edges the edgeNormal will have been truncated and
|
|
// only be still be the old value where two faces have the same normal
|
|
forAll(edgeFaces, edgei)
|
|
{
|
|
const labelList& eFaces = pp.edgeFaces()[edgei];
|
|
|
|
label meshEdgei = meshEdges[edgei];
|
|
|
|
const vector& n = edgeNormal[meshEdgei];
|
|
|
|
if (n != point::max)
|
|
{
|
|
scalar cos = n & pp.faceNormals()[eFaces[0]];
|
|
|
|
if (cos < minCos)
|
|
{
|
|
const edge& e = pp.edges()[edgei];
|
|
|
|
unmarkExtrusion
|
|
(
|
|
e[0],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
unmarkExtrusion
|
|
(
|
|
e[1],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
nFeats++;
|
|
|
|
if (str)
|
|
{
|
|
const point& p0 = pp.localPoints()[e[0]];
|
|
const point& p1 = pp.localPoints()[e[1]];
|
|
str().write(linePointRef(p0, p1));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero for points on "
|
|
<< returnReduce(nFeats, sumOp<label>())
|
|
<< " feature edges" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
// No extrusion on cells with warped faces. Calculates the thickness of the
|
|
// layer and compares it to the space the warped face takes up. Disables
|
|
// extrusion if layer thickness is more than faceRatio of the thickness of
|
|
// the face.
|
|
void Foam::snappyLayerDriver::handleWarpedFaces
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const scalar faceRatio,
|
|
const boolList& relativeSizes,
|
|
const scalar edge0Len,
|
|
const labelList& cellLevel,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
Info<< nl << "Handling cells with warped patch faces ..." << nl;
|
|
|
|
const pointField& points = mesh.points();
|
|
|
|
// Local reference to face centres also used to trigger consistent
|
|
// [re-]building of demand-driven face centres and areas
|
|
const vectorField& faceCentres = mesh.faceCentres();
|
|
|
|
label nWarpedFaces = 0;
|
|
|
|
forAll(pp, i)
|
|
{
|
|
const face& f = pp[i];
|
|
label faceI = pp.addressing()[i];
|
|
label patchI = patches.patchID()[faceI-mesh.nInternalFaces()];
|
|
|
|
// It is hard to calculate some length scale if not in relative
|
|
// mode so disable this check.
|
|
|
|
if (relativeSizes[patchI] && f.size() > 3)
|
|
{
|
|
label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
|
|
scalar edgeLen = edge0Len/(1<<ownLevel);
|
|
|
|
// Normal distance to face centre plane
|
|
const point& fc = faceCentres[faceI];
|
|
const vector& fn = pp.faceNormals()[i];
|
|
|
|
scalarField vProj(f.size());
|
|
|
|
forAll(f, fp)
|
|
{
|
|
vector n = points[f[fp]] - fc;
|
|
vProj[fp] = (n & fn);
|
|
}
|
|
|
|
// Get normal 'span' of face
|
|
scalar minVal = min(vProj);
|
|
scalar maxVal = max(vProj);
|
|
|
|
if ((maxVal - minVal) > faceRatio * edgeLen)
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[i],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nWarpedFaces++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero on "
|
|
<< returnReduce(nWarpedFaces, sumOp<label>())
|
|
<< " warped faces since layer would be > " << faceRatio
|
|
<< " of the size of the bounding box." << endl;
|
|
}
|
|
|
|
|
|
//// No extrusion on cells with multiple patch faces. There usually is a reason
|
|
//// why combinePatchFaces hasn't succeeded.
|
|
//void Foam::snappyLayerDriver::handleMultiplePatchFaces
|
|
//(
|
|
// const indirectPrimitivePatch& pp,
|
|
// pointField& patchDisp,
|
|
// labelList& patchNLayers,
|
|
// List<extrudeMode>& extrudeStatus
|
|
//) const
|
|
//{
|
|
// const fvMesh& mesh = meshRefiner_.mesh();
|
|
//
|
|
// Info<< nl << "Handling cells with multiple patch faces ..." << nl;
|
|
//
|
|
// const labelListList& pointFaces = pp.pointFaces();
|
|
//
|
|
// // Cells that should not get an extrusion layer
|
|
// cellSet multiPatchCells(mesh, "multiPatchCells", pp.size());
|
|
//
|
|
// // Detect points that use multiple faces on same cell.
|
|
// forAll(pointFaces, patchPointi)
|
|
// {
|
|
// const labelList& pFaces = pointFaces[patchPointi];
|
|
//
|
|
// labelHashSet pointCells(pFaces.size());
|
|
//
|
|
// forAll(pFaces, i)
|
|
// {
|
|
// label celli = mesh.faceOwner()[pp.addressing()[pFaces[i]]];
|
|
//
|
|
// if (!pointCells.insert(celli))
|
|
// {
|
|
// // Second or more occurrence of cell so cell has two or more
|
|
// // pp faces connected to this point.
|
|
// multiPatchCells.insert(celli);
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// label nMultiPatchCells = returnReduce
|
|
// (
|
|
// multiPatchCells.size(),
|
|
// sumOp<label>()
|
|
// );
|
|
//
|
|
// Info<< "Detected " << nMultiPatchCells
|
|
// << " cells with multiple (connected) patch faces." << endl;
|
|
//
|
|
// label nChanged = 0;
|
|
//
|
|
// if (nMultiPatchCells > 0)
|
|
// {
|
|
// multiPatchCells.instance() = meshRefiner_.timeName();
|
|
// Info<< "Writing " << nMultiPatchCells
|
|
// << " cells with multiple (connected) patch faces to cellSet "
|
|
// << multiPatchCells.objectPath() << endl;
|
|
// multiPatchCells.write();
|
|
//
|
|
//
|
|
// // Go through all points and remove extrusion on any cell in
|
|
// // multiPatchCells
|
|
// // (has to be done in separate loop since having one point on
|
|
// // multipatches has to reset extrusion on all points of cell)
|
|
//
|
|
// forAll(pointFaces, patchPointi)
|
|
// {
|
|
// if (extrudeStatus[patchPointi] != NOEXTRUDE)
|
|
// {
|
|
// const labelList& pFaces = pointFaces[patchPointi];
|
|
//
|
|
// forAll(pFaces, i)
|
|
// {
|
|
// label celli =
|
|
// mesh.faceOwner()[pp.addressing()[pFaces[i]]];
|
|
//
|
|
// if (multiPatchCells.found(celli))
|
|
// {
|
|
// if
|
|
// (
|
|
// unmarkExtrusion
|
|
// (
|
|
// patchPointi,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
// )
|
|
// {
|
|
// nChanged++;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// reduce(nChanged, sumOp<label>());
|
|
// }
|
|
//
|
|
// Info<< "Prevented extrusion on " << nChanged
|
|
// << " points due to multiple patch faces." << nl << endl;
|
|
//}
|
|
|
|
|
|
void Foam::snappyLayerDriver::setNumLayers
|
|
(
|
|
const labelList& patchToNLayers,
|
|
const labelList& patchIDs,
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus,
|
|
label& nAddedCells
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
Info<< nl << "Handling points with inconsistent layer specification ..."
|
|
<< endl;
|
|
|
|
// Get for every point (really only necessary on patch external points)
|
|
// the max and min of any patch faces using it.
|
|
labelList maxLayers(patchNLayers.size(), labelMin);
|
|
labelList minLayers(patchNLayers.size(), labelMax);
|
|
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
|
|
const labelList& meshPoints = mesh.boundaryMesh()[patchi].meshPoints();
|
|
|
|
label wantedLayers = patchToNLayers[patchi];
|
|
|
|
forAll(meshPoints, patchPointi)
|
|
{
|
|
label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
|
|
|
|
maxLayers[ppPointi] = max(wantedLayers, maxLayers[ppPointi]);
|
|
minLayers[ppPointi] = min(wantedLayers, minLayers[ppPointi]);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
maxLayers,
|
|
maxEqOp<label>(),
|
|
labelMin // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
minLayers,
|
|
minEqOp<label>(),
|
|
labelMax // null value
|
|
);
|
|
|
|
// Unmark any point with different min and max
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
//label nConflicts = 0;
|
|
|
|
forAll(maxLayers, i)
|
|
{
|
|
if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
|
|
<< " maxLayers:" << maxLayers
|
|
<< " minLayers:" << minLayers
|
|
<< abort(FatalError);
|
|
}
|
|
else if (maxLayers[i] == minLayers[i])
|
|
{
|
|
// Ok setting.
|
|
patchNLayers[i] = maxLayers[i];
|
|
}
|
|
else
|
|
{
|
|
// Inconsistent num layers between patch faces using point
|
|
//if
|
|
//(
|
|
// unmarkExtrusion
|
|
// (
|
|
// i,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
//)
|
|
//{
|
|
// nConflicts++;
|
|
//}
|
|
patchNLayers[i] = maxLayers[i];
|
|
}
|
|
}
|
|
|
|
|
|
// Calculate number of cells to create
|
|
nAddedCells = 0;
|
|
forAll(pp.localFaces(), facei)
|
|
{
|
|
const face& f = pp.localFaces()[facei];
|
|
|
|
// Get max of extrusion per point
|
|
label nCells = 0;
|
|
forAll(f, fp)
|
|
{
|
|
nCells = max(nCells, patchNLayers[f[fp]]);
|
|
}
|
|
|
|
nAddedCells += nCells;
|
|
}
|
|
reduce(nAddedCells, sumOp<label>());
|
|
|
|
//reduce(nConflicts, sumOp<label>());
|
|
//
|
|
//Info<< "Set displacement to zero for " << nConflicts
|
|
// << " points due to points being on multiple regions"
|
|
// << " with inconsistent nLayers specification." << endl;
|
|
}
|
|
|
|
|
|
// Construct pointVectorField with correct boundary conditions for adding
|
|
// layers
|
|
Foam::tmp<Foam::pointVectorField>
|
|
Foam::snappyLayerDriver::makeLayerDisplacementField
|
|
(
|
|
const pointMesh& pMesh,
|
|
const labelList& numLayers
|
|
)
|
|
{
|
|
// Construct displacement field.
|
|
const pointBoundaryMesh& pointPatches = pMesh.boundary();
|
|
|
|
wordList patchFieldTypes
|
|
(
|
|
pointPatches.size(),
|
|
slipPointPatchVectorField::typeName
|
|
);
|
|
wordList actualPatchTypes(patchFieldTypes.size());
|
|
forAll(pointPatches, patchi)
|
|
{
|
|
actualPatchTypes[patchi] = pointPatches[patchi].type();
|
|
}
|
|
|
|
forAll(numLayers, patchi)
|
|
{
|
|
// 0 layers: do not allow slip so fixedValue 0
|
|
// >0 layers: fixedValue which gets adapted
|
|
if (numLayers[patchi] == 0)
|
|
{
|
|
patchFieldTypes[patchi] =
|
|
zeroFixedValuePointPatchVectorField::typeName;
|
|
}
|
|
else if (numLayers[patchi] > 0)
|
|
{
|
|
patchFieldTypes[patchi] = fixedValuePointPatchVectorField::typeName;
|
|
}
|
|
}
|
|
|
|
forAll(pointPatches, patchi)
|
|
{
|
|
if (isA<processorPointPatch>(pointPatches[patchi]))
|
|
{
|
|
patchFieldTypes[patchi] = calculatedPointPatchVectorField::typeName;
|
|
}
|
|
else if (isA<cyclicPointPatch>(pointPatches[patchi]))
|
|
{
|
|
patchFieldTypes[patchi] = cyclicSlipPointPatchVectorField::typeName;
|
|
}
|
|
}
|
|
|
|
|
|
const polyMesh& mesh = pMesh();
|
|
|
|
// Note: time().timeName() instead of meshRefinement::timeName() since
|
|
// postprocessable field.
|
|
|
|
return tmp<pointVectorField>::New
|
|
(
|
|
IOobject
|
|
(
|
|
"pointDisplacement",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
pMesh,
|
|
dimensionedVector(dimLength, Zero),
|
|
patchFieldTypes,
|
|
actualPatchTypes
|
|
);
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::growNoExtrusion
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
|
|
|
|
List<extrudeMode> grownExtrudeStatus(extrudeStatus);
|
|
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
label nGrown = 0;
|
|
|
|
forAll(localFaces, facei)
|
|
{
|
|
const face& f = localFaces[facei];
|
|
|
|
bool hasSqueeze = false;
|
|
forAll(f, fp)
|
|
{
|
|
if (extrudeStatus[f[fp]] == NOEXTRUDE)
|
|
{
|
|
hasSqueeze = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (hasSqueeze)
|
|
{
|
|
// Squeeze all points of face
|
|
forAll(f, fp)
|
|
{
|
|
if
|
|
(
|
|
extrudeStatus[f[fp]] == EXTRUDE
|
|
&& grownExtrudeStatus[f[fp]] != NOEXTRUDE
|
|
)
|
|
{
|
|
grownExtrudeStatus[f[fp]] = NOEXTRUDE;
|
|
nGrown++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
extrudeStatus.transfer(grownExtrudeStatus);
|
|
|
|
|
|
// Synchronise since might get called multiple times.
|
|
// Use the fact that NOEXTRUDE is the minimum value.
|
|
{
|
|
labelList status(extrudeStatus.size());
|
|
forAll(status, i)
|
|
{
|
|
status[i] = extrudeStatus[i];
|
|
}
|
|
syncTools::syncPointList
|
|
(
|
|
meshRefiner_.mesh(),
|
|
pp.meshPoints(),
|
|
status,
|
|
minEqOp<label>(),
|
|
labelMax // null value
|
|
);
|
|
forAll(status, i)
|
|
{
|
|
extrudeStatus[i] = extrudeMode(status[i]);
|
|
}
|
|
}
|
|
|
|
|
|
forAll(extrudeStatus, patchPointi)
|
|
{
|
|
if (extrudeStatus[patchPointi] == NOEXTRUDE)
|
|
{
|
|
patchDisp[patchPointi] = Zero;
|
|
patchNLayers[patchPointi] = 0;
|
|
}
|
|
}
|
|
|
|
reduce(nGrown, sumOp<label>());
|
|
|
|
Info<< "Set displacement to zero for an additional " << nGrown
|
|
<< " points." << endl;
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::determineSidePatches
|
|
(
|
|
const globalIndex& globalFaces,
|
|
const labelListList& edgeGlobalFaces,
|
|
const indirectPrimitivePatch& pp,
|
|
|
|
labelList& edgePatchID,
|
|
labelList& edgeZoneID,
|
|
boolList& edgeFlip,
|
|
labelList& inflateFaceID
|
|
)
|
|
{
|
|
// Sometimes edges-to-be-extruded are on more than 2 processors.
|
|
// Work out which 2 hold the faces to be extruded and thus which procpatch
|
|
// the edge-face should be in. As an additional complication this might
|
|
// mean that 2 procesors that were only edge-connected now suddenly need
|
|
// to become face-connected i.e. have a processor patch between them.
|
|
|
|
fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
// Determine edgePatchID. Any additional processor boundary gets added to
|
|
// patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
|
|
// of patches.
|
|
// Note: faceZones are at this point split into baffles so any zone
|
|
// information might also come from boundary faces (hence
|
|
// zoneFromAnyFace set in call to calcExtrudeInfo)
|
|
label nPatches;
|
|
Map<label> nbrProcToPatch;
|
|
Map<label> patchToNbrProc;
|
|
addPatchCellLayer::calcExtrudeInfo
|
|
(
|
|
true, // zoneFromAnyFace
|
|
|
|
mesh,
|
|
globalFaces,
|
|
edgeGlobalFaces,
|
|
pp,
|
|
|
|
edgePatchID,
|
|
nPatches,
|
|
nbrProcToPatch,
|
|
patchToNbrProc,
|
|
edgeZoneID,
|
|
edgeFlip,
|
|
inflateFaceID
|
|
);
|
|
|
|
label nOldPatches = mesh.boundaryMesh().size();
|
|
label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
|
|
Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
|
|
<< " handle extrusion of non-manifold processor boundaries."
|
|
<< endl;
|
|
|
|
if (nAdded > 0)
|
|
{
|
|
// We might not add patches in same order as in patchToNbrProc
|
|
// so prepare to renumber edgePatchID
|
|
Map<label> wantedToAddedPatch;
|
|
|
|
for (label patchi = nOldPatches; patchi < nPatches; patchi++)
|
|
{
|
|
label nbrProci = patchToNbrProc[patchi];
|
|
word name
|
|
(
|
|
processorPolyPatch::newName(Pstream::myProcNo(), nbrProci)
|
|
);
|
|
|
|
dictionary patchDict;
|
|
patchDict.add("type", processorPolyPatch::typeName);
|
|
patchDict.add("myProcNo", Pstream::myProcNo());
|
|
patchDict.add("neighbProcNo", nbrProci);
|
|
patchDict.add("nFaces", 0);
|
|
patchDict.add("startFace", mesh.nFaces());
|
|
|
|
//Pout<< "Adding patch " << patchi
|
|
// << " name:" << name
|
|
// << " between " << Pstream::myProcNo()
|
|
// << " and " << nbrProci << endl;
|
|
|
|
label procPatchi = meshRefiner_.appendPatch
|
|
(
|
|
mesh,
|
|
mesh.boundaryMesh().size(), // new patch index
|
|
name,
|
|
patchDict
|
|
);
|
|
wantedToAddedPatch.insert(patchi, procPatchi);
|
|
}
|
|
|
|
// Renumber edgePatchID
|
|
forAll(edgePatchID, i)
|
|
{
|
|
label patchi = edgePatchID[i];
|
|
const auto fnd = wantedToAddedPatch.cfind(patchi);
|
|
if (fnd.found())
|
|
{
|
|
edgePatchID[i] = fnd.val();
|
|
}
|
|
}
|
|
|
|
mesh.clearOut();
|
|
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::calculateLayerThickness
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& patchIDs,
|
|
const layerParameters& layerParams,
|
|
const labelList& cellLevel,
|
|
const labelList& patchNLayers,
|
|
const scalar edge0Len,
|
|
|
|
scalarField& thickness,
|
|
scalarField& minThickness,
|
|
scalarField& expansionRatio
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
|
|
// Rework patch-wise layer parameters into minimum per point
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
// Note: only layer parameters consistent with layer specification
|
|
// method (see layerParameters) will be correct.
|
|
scalarField firstLayerThickness(pp.nPoints(), GREAT);
|
|
scalarField finalLayerThickness(pp.nPoints(), GREAT);
|
|
scalarField totalThickness(pp.nPoints(), GREAT);
|
|
scalarField expRatio(pp.nPoints(), GREAT);
|
|
|
|
minThickness.setSize(pp.nPoints());
|
|
minThickness = GREAT;
|
|
|
|
thickness.setSize(pp.nPoints());
|
|
thickness = GREAT;
|
|
|
|
expansionRatio.setSize(pp.nPoints());
|
|
expansionRatio = GREAT;
|
|
|
|
for (const label patchi : patchIDs)
|
|
{
|
|
const labelList& meshPoints = patches[patchi].meshPoints();
|
|
|
|
forAll(meshPoints, patchPointi)
|
|
{
|
|
label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
|
|
|
|
firstLayerThickness[ppPointi] = min
|
|
(
|
|
firstLayerThickness[ppPointi],
|
|
layerParams.firstLayerThickness()[patchi]
|
|
);
|
|
finalLayerThickness[ppPointi] = min
|
|
(
|
|
finalLayerThickness[ppPointi],
|
|
layerParams.finalLayerThickness()[patchi]
|
|
);
|
|
totalThickness[ppPointi] = min
|
|
(
|
|
totalThickness[ppPointi],
|
|
layerParams.thickness()[patchi]
|
|
);
|
|
expRatio[ppPointi] = min
|
|
(
|
|
expRatio[ppPointi],
|
|
layerParams.expansionRatio()[patchi]
|
|
);
|
|
minThickness[ppPointi] = min
|
|
(
|
|
minThickness[ppPointi],
|
|
layerParams.minThickness()[patchi]
|
|
);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
firstLayerThickness,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
finalLayerThickness,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
totalThickness,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
expRatio,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
minThickness,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
|
|
|
|
// Now the thicknesses are set according to the minimum of connected
|
|
// patches.
|
|
|
|
// Determine per point the max cell level of connected cells
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
labelList maxPointLevel(pp.nPoints(), labelMin);
|
|
{
|
|
forAll(pp, i)
|
|
{
|
|
label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
|
|
|
|
const face& f = pp.localFaces()[i];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
maxPointLevel,
|
|
maxEqOp<label>(),
|
|
labelMin // null value
|
|
);
|
|
}
|
|
|
|
|
|
// Rework relative thickness into absolute
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
// by multiplying with the internal cell size.
|
|
// Note that we cannot loop over the patches since then points on
|
|
// multiple patches would get multiplied with edgeLen twice ..
|
|
{
|
|
// Multiplication factor for relative sizes
|
|
scalarField edgeLen(pp.nPoints(), GREAT);
|
|
|
|
labelList spec(pp.nPoints(), layerParameters::FIRST_AND_TOTAL);
|
|
|
|
bitSet isRelativePoint(mesh.nPoints());
|
|
|
|
for (const label patchi : patchIDs)
|
|
{
|
|
const labelList& meshPoints = patches[patchi].meshPoints();
|
|
const layerParameters::thicknessModelType patchSpec =
|
|
layerParams.layerModels()[patchi];
|
|
const bool relSize = layerParams.relativeSizes()[patchi];
|
|
|
|
for (const label meshPointi : meshPoints)
|
|
{
|
|
const label ppPointi = pp.meshPointMap()[meshPointi];
|
|
|
|
// Note: who wins if different specs?
|
|
|
|
// Calculate undistorted edge size for this level.
|
|
edgeLen[ppPointi] = min
|
|
(
|
|
edgeLen[ppPointi],
|
|
edge0Len/(1<<maxPointLevel[ppPointi])
|
|
);
|
|
spec[ppPointi] = max(spec[ppPointi], patchSpec);
|
|
isRelativePoint[meshPointi] =
|
|
isRelativePoint[meshPointi]
|
|
|| relSize;
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
edgeLen,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
spec,
|
|
maxEqOp<label>(),
|
|
label(layerParameters::FIRST_AND_TOTAL) // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
isRelativePoint,
|
|
orEqOp<unsigned int>(),
|
|
0
|
|
);
|
|
|
|
|
|
|
|
|
|
forAll(pp.meshPoints(), pointi)
|
|
{
|
|
const label meshPointi = pp.meshPoints()[pointi];
|
|
const layerParameters::thicknessModelType pointSpec =
|
|
static_cast<layerParameters::thicknessModelType>(spec[pointi]);
|
|
|
|
if (pointSpec == layerParameters::FIRST_AND_RELATIVE_FINAL)
|
|
{
|
|
// This overrules the relative sizes flag for
|
|
// first (always absolute) and final (always relative)
|
|
finalLayerThickness[pointi] *= edgeLen[pointi];
|
|
if (isRelativePoint[meshPointi])
|
|
{
|
|
totalThickness[pointi] *= edgeLen[pointi];
|
|
minThickness[pointi] *= edgeLen[pointi];
|
|
}
|
|
}
|
|
else if (isRelativePoint[meshPointi])
|
|
{
|
|
firstLayerThickness[pointi] *= edgeLen[pointi];
|
|
finalLayerThickness[pointi] *= edgeLen[pointi];
|
|
totalThickness[pointi] *= edgeLen[pointi];
|
|
minThickness[pointi] *= edgeLen[pointi];
|
|
}
|
|
|
|
thickness[pointi] = min
|
|
(
|
|
thickness[pointi],
|
|
layerParameters::layerThickness
|
|
(
|
|
pointSpec,
|
|
patchNLayers[pointi],
|
|
firstLayerThickness[pointi],
|
|
finalLayerThickness[pointi],
|
|
totalThickness[pointi],
|
|
expRatio[pointi]
|
|
)
|
|
);
|
|
expansionRatio[pointi] = min
|
|
(
|
|
expansionRatio[pointi],
|
|
layerParameters::layerExpansionRatio
|
|
(
|
|
pointSpec,
|
|
patchNLayers[pointi],
|
|
firstLayerThickness[pointi],
|
|
finalLayerThickness[pointi],
|
|
totalThickness[pointi],
|
|
expRatio[pointi]
|
|
)
|
|
);
|
|
}
|
|
}
|
|
|
|
// Synchronise the determined thicknes. Note that this should not be
|
|
// necessary since the inputs to the calls to layerThickness,
|
|
// layerExpansionRatio above are already parallel consistent
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
thickness,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
expansionRatio,
|
|
minEqOp<scalar>(),
|
|
GREAT // null value
|
|
);
|
|
|
|
//Info<< "calculateLayerThickness : min:" << gMin(thickness)
|
|
// << " max:" << gMax(thickness) << endl;
|
|
|
|
// Print a bit
|
|
{
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
const int oldPrecision = Info.stream().precision();
|
|
|
|
// Find maximum length of a patch name, for a nicer output
|
|
label maxPatchNameLen = 0;
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
word patchName = patches[patchi].name();
|
|
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
|
|
}
|
|
|
|
Info<< nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
|
|
<< setw(0) << " faces layers avg thickness[m]" << nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << " "
|
|
<< setw(0) << " near-wall overall" << nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
|
|
<< setw(0) << " ----- ------ --------- -------" << endl;
|
|
|
|
|
|
const bitSet isMasterPoint(syncTools::getMasterPoints(mesh));
|
|
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
|
|
const labelList& meshPoints = patches[patchi].meshPoints();
|
|
const layerParameters::thicknessModelType spec =
|
|
layerParams.layerModels()[patchi];
|
|
|
|
scalar sumThickness = 0;
|
|
scalar sumNearWallThickness = 0;
|
|
label nMasterPoints = 0;
|
|
|
|
forAll(meshPoints, patchPointi)
|
|
{
|
|
label meshPointi = meshPoints[patchPointi];
|
|
if (isMasterPoint[meshPointi])
|
|
{
|
|
label ppPointi = pp.meshPointMap()[meshPointi];
|
|
|
|
sumThickness += thickness[ppPointi];
|
|
sumNearWallThickness += layerParams.firstLayerThickness
|
|
(
|
|
spec,
|
|
patchNLayers[ppPointi],
|
|
firstLayerThickness[ppPointi],
|
|
finalLayerThickness[ppPointi],
|
|
thickness[ppPointi],
|
|
expansionRatio[ppPointi]
|
|
);
|
|
nMasterPoints++;
|
|
}
|
|
}
|
|
|
|
label totNPoints = returnReduce(nMasterPoints, sumOp<label>());
|
|
|
|
// For empty patches, totNPoints is 0.
|
|
scalar avgThickness = 0;
|
|
scalar avgNearWallThickness = 0;
|
|
|
|
if (totNPoints > 0)
|
|
{
|
|
avgThickness =
|
|
returnReduce(sumThickness, sumOp<scalar>())
|
|
/ totNPoints;
|
|
avgNearWallThickness =
|
|
returnReduce(sumNearWallThickness, sumOp<scalar>())
|
|
/ totNPoints;
|
|
}
|
|
|
|
Info<< setf(ios_base::left) << setw(maxPatchNameLen)
|
|
<< patches[patchi].name() << setprecision(3)
|
|
<< " " << setw(8)
|
|
<< returnReduce(patches[patchi].size(), sumOp<scalar>())
|
|
<< " " << setw(6) << layerParams.numLayers()[patchi]
|
|
<< " " << setw(8) << avgNearWallThickness
|
|
<< " " << setw(8) << avgThickness
|
|
<< endl;
|
|
}
|
|
Info<< setprecision(oldPrecision) << endl;
|
|
}
|
|
}
|
|
|
|
|
|
// Synchronize displacement among coupled patches.
|
|
void Foam::snappyLayerDriver::syncPatchDisplacement
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const scalarField& minThickness,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
|
|
label nChangedTotal = 0;
|
|
|
|
while (true)
|
|
{
|
|
label nChanged = 0;
|
|
|
|
// Sync displacement (by taking min)
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
meshPoints,
|
|
patchDisp,
|
|
minMagSqrEqOp<vector>(),
|
|
point::rootMax // null value
|
|
);
|
|
|
|
// Unmark if displacement too small
|
|
forAll(patchDisp, i)
|
|
{
|
|
if (mag(patchDisp[i]) < minThickness[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
labelList syncPatchNLayers(patchNLayers);
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
meshPoints,
|
|
syncPatchNLayers,
|
|
minEqOp<label>(),
|
|
labelMax // null value
|
|
);
|
|
|
|
// Reset if differs
|
|
// 1. take max
|
|
forAll(syncPatchNLayers, i)
|
|
{
|
|
if (syncPatchNLayers[i] != patchNLayers[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
meshPoints,
|
|
syncPatchNLayers,
|
|
maxEqOp<label>(),
|
|
labelMin // null value
|
|
);
|
|
|
|
// Reset if differs
|
|
// 2. take min
|
|
forAll(syncPatchNLayers, i)
|
|
{
|
|
if (syncPatchNLayers[i] != patchNLayers[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
nChangedTotal += nChanged;
|
|
|
|
if (!returnReduce(nChanged, sumOp<label>()))
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
//Info<< "Prevented extrusion on "
|
|
// << returnReduce(nChangedTotal, sumOp<label>())
|
|
// << " coupled patch points during syncPatchDisplacement." << endl;
|
|
}
|
|
|
|
|
|
// Calculate displacement vector for all patch points. Uses pointNormal.
|
|
// Checks that displaced patch point would be visible from all centres
|
|
// of the faces using it.
|
|
// extrudeStatus is both input and output and gives the status of each
|
|
// patch point.
|
|
void Foam::snappyLayerDriver::getPatchDisplacement
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const scalarField& thickness,
|
|
const scalarField& minThickness,
|
|
const scalarField& expansionRatio,
|
|
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Determining displacement for added points"
|
|
<< " according to pointNormal ..." << endl;
|
|
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const vectorField& faceNormals = pp.faceNormals();
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
const pointField& localPoints = pp.localPoints();
|
|
|
|
// Determine pointNormal
|
|
// ~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
pointField pointNormals(PatchTools::pointNormals(mesh, pp));
|
|
|
|
|
|
// Determine local length scale on patch
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Start off from same thickness everywhere (except where no extrusion)
|
|
patchDisp = thickness*pointNormals;
|
|
|
|
|
|
label nNoVisNormal = 0;
|
|
label nExtrudeRemove = 0;
|
|
|
|
|
|
////XXXXXXXX
|
|
// {
|
|
// OBJstream twoStr
|
|
// (
|
|
// mesh.time().path()
|
|
// / "twoFacePoints_"
|
|
// + meshRefiner_.timeName()
|
|
// + ".obj"
|
|
// );
|
|
// OBJstream multiStr
|
|
// (
|
|
// mesh.time().path()
|
|
// / "multiFacePoints_"
|
|
// + meshRefiner_.timeName()
|
|
// + ".obj"
|
|
// );
|
|
// Pout<< "Writing points inbetween two faces on same cell to "
|
|
// << twoStr.name() << endl;
|
|
// Pout<< "Writing points inbetween three or more faces on same cell to "
|
|
// << multiStr.name() << endl;
|
|
// // Check whether inbetween patch faces on same cell
|
|
// Map<labelList> cellToFaces;
|
|
// forAll(pointNormals, patchPointi)
|
|
// {
|
|
// const labelList& pFaces = pointFaces[patchPointi];
|
|
//
|
|
// cellToFaces.clear();
|
|
// forAll(pFaces, pFacei)
|
|
// {
|
|
// const label patchFacei = pFaces[pFacei];
|
|
// const label meshFacei = pp.addressing()[patchFacei];
|
|
// const label celli = mesh.faceOwner()[meshFacei];
|
|
// Map<labelList>::iterator faceFnd = cellToFaces.find(celli);
|
|
// if (faceFnd.found())
|
|
// {
|
|
// labelList& faces = faceFnd();
|
|
// faces.appendUniq(patchFacei);
|
|
// }
|
|
// else
|
|
// {
|
|
// cellToFaces.insert(celli, labelList(one{}, patchFacei));
|
|
// }
|
|
// }
|
|
//
|
|
// forAllConstIter(Map<labelList>, cellToFaces, iter)
|
|
// {
|
|
// if (iter().size() == 2)
|
|
// {
|
|
// twoStr.write(pp.localPoints()[patchPointi]);
|
|
// }
|
|
// else if (iter().size() > 2)
|
|
// {
|
|
// multiStr.write(pp.localPoints()[patchPointi]);
|
|
//
|
|
// const scalar ratio =
|
|
// layerParameters::finalLayerThicknessRatio
|
|
// (
|
|
// patchNLayers[patchPointi],
|
|
// expansionRatio[patchPointi]
|
|
// );
|
|
// // Get thickness of cell next to bulk
|
|
// const vector finalDisp
|
|
// (
|
|
// ratio*patchDisp[patchPointi]
|
|
// );
|
|
//
|
|
// //Pout<< "** point:" << pp.localPoints()[patchPointi]
|
|
// // << " on cell:" << iter.key()
|
|
// // << " faces:" << iter()
|
|
// // << " displacement was:" << patchDisp[patchPointi]
|
|
// // << " ratio:" << ratio
|
|
// // << " finalDispl:" << finalDisp;
|
|
//
|
|
// // Half this thickness
|
|
// patchDisp[patchPointi] -= 0.8*finalDisp;
|
|
//
|
|
// //Pout<< " new displacement:"
|
|
// // << patchDisp[patchPointi] << endl;
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// Pout<< "Written " << multiStr.nVertices()
|
|
// << " points inbetween three or more faces on same cell to "
|
|
// << multiStr.name() << endl;
|
|
// }
|
|
////XXXXXXXX
|
|
|
|
|
|
// Check if no extrude possible.
|
|
forAll(pointNormals, patchPointi)
|
|
{
|
|
label meshPointi = pp.meshPoints()[patchPointi];
|
|
|
|
if (extrudeStatus[patchPointi] == NOEXTRUDE)
|
|
{
|
|
// Do not use unmarkExtrusion; forcibly set to zero extrusion.
|
|
patchNLayers[patchPointi] = 0;
|
|
patchDisp[patchPointi] = Zero;
|
|
}
|
|
else
|
|
{
|
|
// Get normal
|
|
const vector& n = pointNormals[patchPointi];
|
|
|
|
if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointi]))
|
|
{
|
|
if (debug&meshRefinement::ATTRACTION)
|
|
{
|
|
Pout<< "No valid normal for point " << meshPointi
|
|
<< ' ' << pp.points()[meshPointi]
|
|
<< "; setting displacement to "
|
|
<< patchDisp[patchPointi]
|
|
<< endl;
|
|
}
|
|
|
|
extrudeStatus[patchPointi] = EXTRUDEREMOVE;
|
|
nNoVisNormal++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// At illegal points make displacement average of new neighbour positions
|
|
forAll(extrudeStatus, patchPointi)
|
|
{
|
|
if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
|
|
{
|
|
point avg(Zero);
|
|
label nPoints = 0;
|
|
|
|
const labelList& pEdges = pp.pointEdges()[patchPointi];
|
|
|
|
forAll(pEdges, i)
|
|
{
|
|
label edgei = pEdges[i];
|
|
|
|
label otherPointi = pp.edges()[edgei].otherVertex(patchPointi);
|
|
|
|
if (extrudeStatus[otherPointi] != NOEXTRUDE)
|
|
{
|
|
avg += localPoints[otherPointi] + patchDisp[otherPointi];
|
|
nPoints++;
|
|
}
|
|
}
|
|
|
|
if (nPoints > 0)
|
|
{
|
|
if (debug&meshRefinement::ATTRACTION)
|
|
{
|
|
Pout<< "Displacement at illegal point "
|
|
<< localPoints[patchPointi]
|
|
<< " set to "
|
|
<< (avg / nPoints - localPoints[patchPointi])
|
|
<< endl;
|
|
}
|
|
|
|
patchDisp[patchPointi] =
|
|
avg / nPoints
|
|
- localPoints[patchPointi];
|
|
|
|
nExtrudeRemove++;
|
|
}
|
|
else
|
|
{
|
|
// All surrounding points are not extruded. Leave patchDisp
|
|
// intact.
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Detected " << returnReduce(nNoVisNormal, sumOp<label>())
|
|
<< " points with point normal pointing through faces." << nl
|
|
<< "Reset displacement at "
|
|
<< returnReduce(nExtrudeRemove, sumOp<label>())
|
|
<< " points to average of surrounding points." << endl;
|
|
|
|
// Make sure displacement is equal on both sides of coupled patches.
|
|
syncPatchDisplacement
|
|
(
|
|
pp,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
Info<< endl;
|
|
}
|
|
|
|
|
|
bool Foam::snappyLayerDriver::sameEdgeNeighbour
|
|
(
|
|
const labelListList& globalEdgeFaces,
|
|
const label myGlobalFacei,
|
|
const label nbrGlobFacei,
|
|
const label edgei
|
|
) const
|
|
{
|
|
const labelList& eFaces = globalEdgeFaces[edgei];
|
|
if (eFaces.size() == 2)
|
|
{
|
|
return edge(myGlobalFacei, nbrGlobFacei) == edge(eFaces[0], eFaces[1]);
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::getVertexString
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelListList& globalEdgeFaces,
|
|
const label facei,
|
|
const label edgei,
|
|
const label myGlobFacei,
|
|
const label nbrGlobFacei,
|
|
DynamicList<label>& vertices
|
|
) const
|
|
{
|
|
const labelList& fEdges = pp.faceEdges()[facei];
|
|
label fp = fEdges.find(edgei);
|
|
|
|
if (fp == -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "problem." << abort(FatalError);
|
|
}
|
|
|
|
// Search back
|
|
label startFp = fp;
|
|
|
|
forAll(fEdges, i)
|
|
{
|
|
label prevFp = fEdges.rcIndex(startFp);
|
|
if
|
|
(
|
|
!sameEdgeNeighbour
|
|
(
|
|
globalEdgeFaces,
|
|
myGlobFacei,
|
|
nbrGlobFacei,
|
|
fEdges[prevFp]
|
|
)
|
|
)
|
|
{
|
|
break;
|
|
}
|
|
startFp = prevFp;
|
|
}
|
|
|
|
label endFp = fp;
|
|
forAll(fEdges, i)
|
|
{
|
|
label nextFp = fEdges.fcIndex(endFp);
|
|
if
|
|
(
|
|
!sameEdgeNeighbour
|
|
(
|
|
globalEdgeFaces,
|
|
myGlobFacei,
|
|
nbrGlobFacei,
|
|
fEdges[nextFp]
|
|
)
|
|
)
|
|
{
|
|
break;
|
|
}
|
|
endFp = nextFp;
|
|
}
|
|
|
|
const face& f = pp.localFaces()[facei];
|
|
vertices.clear();
|
|
fp = startFp;
|
|
while (fp != endFp)
|
|
{
|
|
vertices.append(f[fp]);
|
|
fp = f.fcIndex(fp);
|
|
}
|
|
vertices.append(f[fp]);
|
|
fp = f.fcIndex(fp);
|
|
vertices.append(f[fp]);
|
|
}
|
|
|
|
|
|
// Truncates displacement
|
|
// - for all patchFaces in the faceset displacement gets set to zero
|
|
// - all displacement < minThickness gets set to zero
|
|
Foam::label Foam::snappyLayerDriver::truncateDisplacement
|
|
(
|
|
const globalIndex& globalFaces,
|
|
const labelListList& edgeGlobalFaces,
|
|
const indirectPrimitivePatch& pp,
|
|
const scalarField& minThickness,
|
|
const faceSet& illegalPatchFaces,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
label nChanged = 0;
|
|
|
|
const Map<label>& meshPointMap = pp.meshPointMap();
|
|
|
|
for (const label facei : illegalPatchFaces)
|
|
{
|
|
if (mesh.isInternalFace(facei))
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Faceset " << illegalPatchFaces.name()
|
|
<< " contains internal face " << facei << nl
|
|
<< "It should only contain patch faces" << abort(FatalError);
|
|
}
|
|
|
|
const face& f = mesh.faces()[facei];
|
|
|
|
|
|
forAll(f, fp)
|
|
{
|
|
const auto fnd = meshPointMap.cfind(f[fp]);
|
|
if (fnd.found())
|
|
{
|
|
const label patchPointi = fnd.val();
|
|
|
|
if (extrudeStatus[patchPointi] != NOEXTRUDE)
|
|
{
|
|
unmarkExtrusion
|
|
(
|
|
patchPointi,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
forAll(patchDisp, patchPointi)
|
|
{
|
|
if (mag(patchDisp[patchPointi]) < minThickness[patchPointi])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
patchPointi,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
else if (extrudeStatus[patchPointi] == NOEXTRUDE)
|
|
{
|
|
// Make sure displacement is 0. Should already be so but ...
|
|
patchDisp[patchPointi] = Zero;
|
|
patchNLayers[patchPointi] = 0;
|
|
}
|
|
}
|
|
|
|
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
while (true)
|
|
{
|
|
syncPatchDisplacement
|
|
(
|
|
pp,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
|
|
// Pinch
|
|
// ~~~~~
|
|
|
|
// Make sure that a face doesn't have two non-consecutive areas
|
|
// not extruded (e.g. quad where vertex 0 and 2 are not extruded
|
|
// but 1 and 3 are) since this gives topological errors.
|
|
|
|
label nPinched = 0;
|
|
|
|
forAll(localFaces, i)
|
|
{
|
|
const face& localF = localFaces[i];
|
|
|
|
// Count number of transitions from unsnapped to snapped.
|
|
label nTrans = 0;
|
|
|
|
extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
|
|
|
|
forAll(localF, fp)
|
|
{
|
|
extrudeMode fpMode = extrudeStatus[localF[fp]];
|
|
|
|
if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
|
|
{
|
|
nTrans++;
|
|
}
|
|
prevMode = fpMode;
|
|
}
|
|
|
|
if (nTrans > 1)
|
|
{
|
|
// Multiple pinches. Reset whole face as unextruded.
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
localF,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nPinched++;
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
reduce(nPinched, sumOp<label>());
|
|
|
|
Info<< "truncateDisplacement : Unextruded " << nPinched
|
|
<< " faces due to non-consecutive vertices being extruded." << endl;
|
|
|
|
|
|
// Butterfly
|
|
// ~~~~~~~~~
|
|
|
|
// Make sure that a string of edges becomes a single face so
|
|
// not a butterfly. Occasionally an 'edge' will have a single dangling
|
|
// vertex due to face combining. These get extruded as a single face
|
|
// (with a dangling vertex) so make sure this extrusion forms a single
|
|
// shape.
|
|
// - continuous i.e. no butterfly:
|
|
// + +
|
|
// |\ /|
|
|
// | \ / |
|
|
// +--+--+
|
|
// - extrudes from all but the endpoints i.e. no partial
|
|
// extrude
|
|
// +
|
|
// /|
|
|
// / |
|
|
// +--+--+
|
|
// The common error topology is a pinch somewhere in the middle
|
|
label nButterFly = 0;
|
|
{
|
|
DynamicList<label> stringedVerts;
|
|
forAll(pp.edges(), edgei)
|
|
{
|
|
const labelList& globFaces = edgeGlobalFaces[edgei];
|
|
|
|
if (globFaces.size() == 2)
|
|
{
|
|
label myFacei = pp.edgeFaces()[edgei][0];
|
|
label myGlobalFacei = globalFaces.toGlobal
|
|
(
|
|
pp.addressing()[myFacei]
|
|
);
|
|
label nbrGlobalFacei =
|
|
(
|
|
globFaces[0] != myGlobalFacei
|
|
? globFaces[0]
|
|
: globFaces[1]
|
|
);
|
|
getVertexString
|
|
(
|
|
pp,
|
|
edgeGlobalFaces,
|
|
myFacei,
|
|
edgei,
|
|
myGlobalFacei,
|
|
nbrGlobalFacei,
|
|
stringedVerts
|
|
);
|
|
|
|
if
|
|
(
|
|
extrudeStatus[stringedVerts[0]] != NOEXTRUDE
|
|
|| extrudeStatus[stringedVerts.last()] != NOEXTRUDE
|
|
)
|
|
{
|
|
// Any pinch in the middle
|
|
bool pinch = false;
|
|
for (label i = 1; i < stringedVerts.size()-1; i++)
|
|
{
|
|
if (extrudeStatus[stringedVerts[i]] == NOEXTRUDE)
|
|
{
|
|
pinch = true;
|
|
break;
|
|
}
|
|
}
|
|
if (pinch)
|
|
{
|
|
forAll(stringedVerts, i)
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
stringedVerts[i],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nButterFly++;
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
reduce(nButterFly, sumOp<label>());
|
|
|
|
Info<< "truncateDisplacement : Unextruded " << nButterFly
|
|
<< " faces due to stringed edges with inconsistent extrusion."
|
|
<< endl;
|
|
|
|
|
|
|
|
// Consistent number of layers
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Make sure that a face has consistent number of layers for all
|
|
// its vertices.
|
|
|
|
label nDiffering = 0;
|
|
|
|
//forAll(localFaces, i)
|
|
//{
|
|
// const face& localF = localFaces[i];
|
|
//
|
|
// label numLayers = -1;
|
|
//
|
|
// forAll(localF, fp)
|
|
// {
|
|
// if (patchNLayers[localF[fp]] > 0)
|
|
// {
|
|
// if (numLayers == -1)
|
|
// {
|
|
// numLayers = patchNLayers[localF[fp]];
|
|
// }
|
|
// else if (numLayers != patchNLayers[localF[fp]])
|
|
// {
|
|
// // Differing number of layers
|
|
// if
|
|
// (
|
|
// unmarkExtrusion
|
|
// (
|
|
// localF,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
// )
|
|
// {
|
|
// nDiffering++;
|
|
// nChanged++;
|
|
// }
|
|
// break;
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
//
|
|
//reduce(nDiffering, sumOp<label>());
|
|
//
|
|
//Info<< "truncateDisplacement : Unextruded " << nDiffering
|
|
// << " faces due to having differing number of layers." << endl;
|
|
|
|
if (nPinched+nButterFly+nDiffering == 0)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
return nChanged;
|
|
}
|
|
|
|
|
|
// Setup layer information (at points and faces) to modify mesh topology in
|
|
// regions where layer mesh terminates.
|
|
void Foam::snappyLayerDriver::setupLayerInfoTruncation
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& patchNLayers,
|
|
const List<extrudeMode>& extrudeStatus,
|
|
const label nBufferCellsNoExtrude,
|
|
labelList& nPatchPointLayers,
|
|
labelList& nPatchFaceLayers
|
|
) const
|
|
{
|
|
Info<< nl << "Setting up information for layer truncation ..." << endl;
|
|
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
if (nBufferCellsNoExtrude < 0)
|
|
{
|
|
Info<< nl << "Performing no layer truncation."
|
|
<< " nBufferCellsNoExtrude set to less than 0 ..." << endl;
|
|
|
|
// Face layers if any point gets extruded
|
|
forAll(pp.localFaces(), patchFacei)
|
|
{
|
|
const face& f = pp.localFaces()[patchFacei];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
const label nPointLayers = patchNLayers[f[fp]];
|
|
if (nPointLayers > 0)
|
|
{
|
|
if (nPatchFaceLayers[patchFacei] == -1)
|
|
{
|
|
nPatchFaceLayers[patchFacei] = nPointLayers;
|
|
}
|
|
else
|
|
{
|
|
nPatchFaceLayers[patchFacei] = min
|
|
(
|
|
nPatchFaceLayers[patchFacei],
|
|
nPointLayers
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
nPatchPointLayers = patchNLayers;
|
|
|
|
// Set any unset patch face layers
|
|
forAll(nPatchFaceLayers, patchFacei)
|
|
{
|
|
if (nPatchFaceLayers[patchFacei] == -1)
|
|
{
|
|
nPatchFaceLayers[patchFacei] = 0;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Determine max point layers per face.
|
|
labelList maxLevel(pp.size(), Zero);
|
|
|
|
forAll(pp.localFaces(), patchFacei)
|
|
{
|
|
const face& f = pp.localFaces()[patchFacei];
|
|
|
|
// find patch faces where layer terminates (i.e contains extrude
|
|
// and noextrude points).
|
|
|
|
bool noExtrude = false;
|
|
label mLevel = 0;
|
|
|
|
forAll(f, fp)
|
|
{
|
|
if (extrudeStatus[f[fp]] == NOEXTRUDE)
|
|
{
|
|
noExtrude = true;
|
|
}
|
|
mLevel = max(mLevel, patchNLayers[f[fp]]);
|
|
}
|
|
|
|
if (mLevel > 0)
|
|
{
|
|
// So one of the points is extruded. Check if all are extruded
|
|
// or is a mix.
|
|
|
|
if (noExtrude)
|
|
{
|
|
nPatchFaceLayers[patchFacei] = 1;
|
|
maxLevel[patchFacei] = mLevel;
|
|
}
|
|
else
|
|
{
|
|
maxLevel[patchFacei] = mLevel;
|
|
}
|
|
}
|
|
}
|
|
|
|
// We have the seed faces (faces with nPatchFaceLayers != maxLevel)
|
|
// Now do a meshwave across the patch where we pick up neighbours
|
|
// of seed faces.
|
|
// Note: quite inefficient. Could probably be coded better.
|
|
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
|
|
label nLevels = gMax(patchNLayers);
|
|
|
|
// flag neighbouring patch faces with number of layers to grow
|
|
for (label ilevel = 1; ilevel < nLevels; ilevel++)
|
|
{
|
|
label nBuffer;
|
|
|
|
if (ilevel == 1)
|
|
{
|
|
nBuffer = nBufferCellsNoExtrude - 1;
|
|
}
|
|
else
|
|
{
|
|
nBuffer = nBufferCellsNoExtrude;
|
|
}
|
|
|
|
for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
|
|
{
|
|
labelList tempCounter(nPatchFaceLayers);
|
|
|
|
boolList foundNeighbour(pp.nPoints(), false);
|
|
|
|
forAll(pp.meshPoints(), patchPointi)
|
|
{
|
|
forAll(pointFaces[patchPointi], pointFacei)
|
|
{
|
|
label facei = pointFaces[patchPointi][pointFacei];
|
|
|
|
if
|
|
(
|
|
nPatchFaceLayers[facei] != -1
|
|
&& maxLevel[facei] > 0
|
|
)
|
|
{
|
|
foundNeighbour[patchPointi] = true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
foundNeighbour,
|
|
orEqOp<bool>(),
|
|
false // null value
|
|
);
|
|
|
|
forAll(pp.meshPoints(), patchPointi)
|
|
{
|
|
if (foundNeighbour[patchPointi])
|
|
{
|
|
forAll(pointFaces[patchPointi], pointFacei)
|
|
{
|
|
label facei = pointFaces[patchPointi][pointFacei];
|
|
if
|
|
(
|
|
nPatchFaceLayers[facei] == -1
|
|
&& maxLevel[facei] > 0
|
|
&& ilevel < maxLevel[facei]
|
|
)
|
|
{
|
|
tempCounter[facei] = ilevel;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
nPatchFaceLayers = tempCounter;
|
|
}
|
|
}
|
|
|
|
forAll(pp.localFaces(), patchFacei)
|
|
{
|
|
if (nPatchFaceLayers[patchFacei] == -1)
|
|
{
|
|
nPatchFaceLayers[patchFacei] = maxLevel[patchFacei];
|
|
}
|
|
}
|
|
|
|
forAll(pp.meshPoints(), patchPointi)
|
|
{
|
|
if (extrudeStatus[patchPointi] != NOEXTRUDE)
|
|
{
|
|
forAll(pointFaces[patchPointi], pointFacei)
|
|
{
|
|
label face = pointFaces[patchPointi][pointFacei];
|
|
nPatchPointLayers[patchPointi] = max
|
|
(
|
|
nPatchPointLayers[patchPointi],
|
|
nPatchFaceLayers[face]
|
|
);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
nPatchPointLayers[patchPointi] = 0;
|
|
}
|
|
}
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
nPatchPointLayers,
|
|
maxEqOp<label>(),
|
|
label(0) // null value
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Does any of the cells use a face from faces?
|
|
bool Foam::snappyLayerDriver::cellsUseFace
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& cellLabels,
|
|
const labelHashSet& faces
|
|
)
|
|
{
|
|
forAll(cellLabels, i)
|
|
{
|
|
const cell& cFaces = mesh.cells()[cellLabels[i]];
|
|
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
if (faces.found(cFaces[cFacei]))
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
// Checks the newly added cells and locally unmarks points so they
|
|
// will not get extruded next time round. Returns global number of unmarked
|
|
// points (0 if all was fine)
|
|
Foam::label Foam::snappyLayerDriver::checkAndUnmark
|
|
(
|
|
const addPatchCellLayer& addLayer,
|
|
const dictionary& meshQualityDict,
|
|
const bool additionalReporting,
|
|
const List<labelPair>& baffles,
|
|
const indirectPrimitivePatch& pp,
|
|
const fvMesh& newMesh,
|
|
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
// Check the resulting mesh for errors
|
|
Info<< nl << "Checking mesh with layer ..." << endl;
|
|
faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
|
|
motionSmoother::checkMesh
|
|
(
|
|
false,
|
|
newMesh,
|
|
meshQualityDict,
|
|
identity(newMesh.nFaces()),
|
|
baffles,
|
|
wrongFaces,
|
|
false // dryRun_
|
|
);
|
|
Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
|
|
<< " illegal faces"
|
|
<< " (concave, zero area or negative cell pyramid volume)"
|
|
<< endl;
|
|
|
|
// Undo local extrusion if
|
|
// - any of the added cells in error
|
|
|
|
label nChanged = 0;
|
|
|
|
// Get all cells in the layer.
|
|
labelListList addedCells
|
|
(
|
|
addPatchCellLayer::addedCells
|
|
(
|
|
newMesh,
|
|
addLayer.layerFaces()
|
|
)
|
|
);
|
|
|
|
// Check if any of the faces in error uses any face of an added cell
|
|
// - if additionalReporting print the few remaining areas for ease of
|
|
// finding out where the problems are.
|
|
|
|
const label nReportMax = 10;
|
|
DynamicField<point> disabledFaceCentres(nReportMax);
|
|
|
|
forAll(addedCells, oldPatchFacei)
|
|
{
|
|
// Get the cells (in newMesh labels) per old patch face (in mesh
|
|
// labels)
|
|
const labelList& fCells = addedCells[oldPatchFacei];
|
|
|
|
if (cellsUseFace(newMesh, fCells, wrongFaces))
|
|
{
|
|
// Unmark points on old mesh
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[oldPatchFacei],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
if (additionalReporting && (nChanged < nReportMax))
|
|
{
|
|
disabledFaceCentres.append
|
|
(
|
|
pp.faceCentres()[oldPatchFacei]
|
|
);
|
|
}
|
|
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
label nChangedTotal = returnReduce(nChanged, sumOp<label>());
|
|
|
|
if (additionalReporting)
|
|
{
|
|
// Limit the number of points to be printed so that
|
|
// not too many points are reported when running in parallel
|
|
// Not accurate, i.e. not always nReportMax points are written,
|
|
// but this estimation avoid some communication here.
|
|
// The important thing, however, is that when only a few faces
|
|
// are disabled, their coordinates are printed, and this should be
|
|
// the case
|
|
label nReportLocal = nChanged;
|
|
if (nChangedTotal > nReportMax)
|
|
{
|
|
nReportLocal = min
|
|
(
|
|
max(nChangedTotal / Pstream::nProcs(), 1),
|
|
min
|
|
(
|
|
nChanged,
|
|
max(nReportMax / Pstream::nProcs(), 1)
|
|
)
|
|
);
|
|
}
|
|
|
|
if (nReportLocal)
|
|
{
|
|
Pout<< "Checked mesh with layers. Disabled extrusion at " << endl;
|
|
for (label i=0; i < nReportLocal; i++)
|
|
{
|
|
Pout<< " " << disabledFaceCentres[i] << endl;
|
|
}
|
|
}
|
|
|
|
label nReportTotal = returnReduce(nReportLocal, sumOp<label>());
|
|
|
|
if (nReportTotal < nChangedTotal)
|
|
{
|
|
Info<< "Suppressed disabled extrusion message for other "
|
|
<< nChangedTotal - nReportTotal << " faces." << endl;
|
|
}
|
|
}
|
|
|
|
return nChangedTotal;
|
|
}
|
|
|
|
|
|
//- Count global number of extruded faces
|
|
Foam::label Foam::snappyLayerDriver::countExtrusion
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
// Count number of extruded patch faces
|
|
label nExtruded = 0;
|
|
{
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
forAll(localFaces, i)
|
|
{
|
|
const face& localFace = localFaces[i];
|
|
|
|
forAll(localFace, fp)
|
|
{
|
|
if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
|
|
{
|
|
nExtruded++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return returnReduce(nExtruded, sumOp<label>());
|
|
}
|
|
|
|
|
|
Foam::List<Foam::labelPair> Foam::snappyLayerDriver::getBafflesOnAddedMesh
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& newToOldFaces,
|
|
const List<labelPair>& baffles
|
|
)
|
|
{
|
|
// The problem is that the baffle faces are now inside the
|
|
// mesh (addPatchCellLayer modifies original boundary faces and
|
|
// adds new ones. So 2 pass:
|
|
// - find the boundary face for all faces originating from baffle
|
|
// - use the boundary face for the new baffles
|
|
|
|
Map<label> baffleSet(4*baffles.size());
|
|
forAll(baffles, bafflei)
|
|
{
|
|
baffleSet.insert(baffles[bafflei][0], bafflei);
|
|
baffleSet.insert(baffles[bafflei][1], bafflei);
|
|
}
|
|
|
|
|
|
List<labelPair> newBaffles(baffles.size(), labelPair(-1, -1));
|
|
for
|
|
(
|
|
label facei = mesh.nInternalFaces();
|
|
facei < mesh.nFaces();
|
|
facei++
|
|
)
|
|
{
|
|
label oldFacei = newToOldFaces[facei];
|
|
|
|
const auto faceFnd = baffleSet.find(oldFacei);
|
|
if (faceFnd.found())
|
|
{
|
|
label bafflei = faceFnd();
|
|
labelPair& p = newBaffles[bafflei];
|
|
if (p[0] == -1)
|
|
{
|
|
p[0] = facei;
|
|
}
|
|
else if (p[1] == -1)
|
|
{
|
|
p[1] = facei;
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Problem:" << facei << " at:"
|
|
<< mesh.faceCentres()[facei]
|
|
<< " is on same baffle as " << p[0]
|
|
<< " at:" << mesh.faceCentres()[p[0]]
|
|
<< " and " << p[1]
|
|
<< " at:" << mesh.faceCentres()[p[1]]
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
return newBaffles;
|
|
}
|
|
|
|
|
|
// Collect layer faces and layer cells into mesh fields for ease of handling
|
|
void Foam::snappyLayerDriver::getLayerCellsFaces
|
|
(
|
|
const polyMesh& mesh,
|
|
const addPatchCellLayer& addLayer,
|
|
const scalarField& oldRealThickness,
|
|
|
|
labelList& cellNLayers,
|
|
scalarField& faceRealThickness
|
|
)
|
|
{
|
|
cellNLayers.setSize(mesh.nCells());
|
|
cellNLayers = 0;
|
|
faceRealThickness.setSize(mesh.nFaces());
|
|
faceRealThickness = 0;
|
|
|
|
// Mark all faces in the layer
|
|
const labelListList& layerFaces = addLayer.layerFaces();
|
|
|
|
// Mark all cells in the layer.
|
|
labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
|
|
|
|
forAll(addedCells, oldPatchFacei)
|
|
{
|
|
const labelList& added = addedCells[oldPatchFacei];
|
|
|
|
const labelList& layer = layerFaces[oldPatchFacei];
|
|
|
|
if (layer.size())
|
|
{
|
|
// Leave out original internal face
|
|
forAll(added, i)
|
|
{
|
|
cellNLayers[added[i]] = layer.size()-1;
|
|
}
|
|
}
|
|
}
|
|
|
|
forAll(layerFaces, oldPatchFacei)
|
|
{
|
|
const labelList& layer = layerFaces[oldPatchFacei];
|
|
const scalar realThickness = oldRealThickness[oldPatchFacei];
|
|
|
|
if (layer.size())
|
|
{
|
|
// Layer contains both original boundary face and new boundary
|
|
// face so is nLayers+1. Leave out old internal face.
|
|
for (label i = 1; i < layer.size(); i++)
|
|
{
|
|
faceRealThickness[layer[i]] = realThickness;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::printLayerData
|
|
(
|
|
const fvMesh& mesh,
|
|
const labelList& patchIDs,
|
|
const labelList& cellNLayers,
|
|
const scalarField& faceWantedThickness,
|
|
const scalarField& faceRealThickness
|
|
) const
|
|
{
|
|
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
|
|
|
const int oldPrecision = Info.stream().precision();
|
|
|
|
// Find maximum length of a patch name, for a nicer output
|
|
label maxPatchNameLen = 0;
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
word patchName = pbm[patchi].name();
|
|
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
|
|
}
|
|
|
|
Info<< nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
|
|
<< setw(0) << " faces layers overall thickness" << nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << " "
|
|
<< setw(0) << " [m] [%]" << nl
|
|
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
|
|
<< setw(0) << " ----- ------ --- ---" << endl;
|
|
|
|
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
const polyPatch& pp = pbm[patchi];
|
|
|
|
label sumSize = pp.size();
|
|
|
|
// Number of layers
|
|
const labelList& faceCells = pp.faceCells();
|
|
label sumNLayers = 0;
|
|
forAll(faceCells, i)
|
|
{
|
|
sumNLayers += cellNLayers[faceCells[i]];
|
|
}
|
|
|
|
// Thickness
|
|
scalarField::subField patchWanted = pbm[patchi].patchSlice
|
|
(
|
|
faceWantedThickness
|
|
);
|
|
scalarField::subField patchReal = pbm[patchi].patchSlice
|
|
(
|
|
faceRealThickness
|
|
);
|
|
|
|
scalar sumRealThickness = sum(patchReal);
|
|
scalar sumFraction = 0;
|
|
forAll(patchReal, i)
|
|
{
|
|
if (patchWanted[i] > VSMALL)
|
|
{
|
|
sumFraction += (patchReal[i]/patchWanted[i]);
|
|
}
|
|
}
|
|
|
|
|
|
reduce(sumSize, sumOp<label>());
|
|
reduce(sumNLayers, sumOp<label>());
|
|
reduce(sumRealThickness, sumOp<scalar>());
|
|
reduce(sumFraction, sumOp<scalar>());
|
|
|
|
|
|
scalar avgLayers = 0;
|
|
scalar avgReal = 0;
|
|
scalar avgFraction = 0;
|
|
if (sumSize > 0)
|
|
{
|
|
avgLayers = scalar(sumNLayers)/sumSize;
|
|
avgReal = sumRealThickness/sumSize;
|
|
avgFraction = sumFraction/sumSize;
|
|
}
|
|
|
|
Info<< setf(ios_base::left) << setw(maxPatchNameLen)
|
|
<< pbm[patchi].name() << setprecision(3)
|
|
<< " " << setw(8) << sumSize
|
|
<< " " << setw(8) << avgLayers
|
|
<< " " << setw(8) << avgReal
|
|
<< " " << setw(8) << 100*avgFraction
|
|
<< endl;
|
|
}
|
|
Info<< setprecision(oldPrecision) << endl;
|
|
}
|
|
|
|
|
|
bool Foam::snappyLayerDriver::writeLayerSets
|
|
(
|
|
const fvMesh& mesh,
|
|
const labelList& cellNLayers,
|
|
const scalarField& faceRealThickness
|
|
) const
|
|
{
|
|
bool allOk = true;
|
|
{
|
|
label nAdded = 0;
|
|
forAll(cellNLayers, celli)
|
|
{
|
|
if (cellNLayers[celli] > 0)
|
|
{
|
|
nAdded++;
|
|
}
|
|
}
|
|
cellSet addedCellSet(mesh, "addedCells", nAdded);
|
|
forAll(cellNLayers, celli)
|
|
{
|
|
if (cellNLayers[celli] > 0)
|
|
{
|
|
addedCellSet.insert(celli);
|
|
}
|
|
}
|
|
addedCellSet.instance() = meshRefiner_.timeName();
|
|
Info<< "Writing "
|
|
<< returnReduce(addedCellSet.size(), sumOp<label>())
|
|
<< " added cells to cellSet "
|
|
<< addedCellSet.name() << endl;
|
|
bool ok = addedCellSet.write();
|
|
allOk = allOk && ok;
|
|
}
|
|
{
|
|
label nAdded = 0;
|
|
for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
|
|
{
|
|
if (faceRealThickness[facei] > 0)
|
|
{
|
|
nAdded++;
|
|
}
|
|
}
|
|
|
|
faceSet layerFacesSet(mesh, "layerFaces", nAdded);
|
|
for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
|
|
{
|
|
if (faceRealThickness[facei] > 0)
|
|
{
|
|
layerFacesSet.insert(facei);
|
|
}
|
|
}
|
|
layerFacesSet.instance() = meshRefiner_.timeName();
|
|
Info<< "Writing "
|
|
<< returnReduce(layerFacesSet.size(), sumOp<label>())
|
|
<< " faces inside added layer to faceSet "
|
|
<< layerFacesSet.name() << endl;
|
|
bool ok = layerFacesSet.write();
|
|
allOk = allOk && ok;
|
|
}
|
|
return allOk;
|
|
}
|
|
|
|
|
|
bool Foam::snappyLayerDriver::writeLayerData
|
|
(
|
|
const fvMesh& mesh,
|
|
const labelList& patchIDs,
|
|
const labelList& cellNLayers,
|
|
const scalarField& faceWantedThickness,
|
|
const scalarField& faceRealThickness
|
|
) const
|
|
{
|
|
bool allOk = true;
|
|
|
|
if (meshRefinement::writeLevel() & meshRefinement::WRITELAYERSETS)
|
|
{
|
|
bool ok = writeLayerSets(mesh, cellNLayers, faceRealThickness);
|
|
allOk = allOk && ok;
|
|
}
|
|
|
|
if (meshRefinement::writeLevel() & meshRefinement::WRITELAYERFIELDS)
|
|
{
|
|
Info<< nl << "Writing fields with layer information:" << incrIndent
|
|
<< endl;
|
|
{
|
|
volScalarField fld
|
|
(
|
|
IOobject
|
|
(
|
|
"nSurfaceLayers",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE,
|
|
false
|
|
),
|
|
mesh,
|
|
dimensionedScalar(dimless, Zero),
|
|
fixedValueFvPatchScalarField::typeName
|
|
);
|
|
forAll(fld, celli)
|
|
{
|
|
fld[celli] = cellNLayers[celli];
|
|
}
|
|
volScalarField::Boundary& fldBf = fld.boundaryFieldRef();
|
|
|
|
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
const polyPatch& pp = pbm[patchi];
|
|
const labelList& faceCells = pp.faceCells();
|
|
scalarField pfld(faceCells.size());
|
|
forAll(faceCells, i)
|
|
{
|
|
pfld[i] = cellNLayers[faceCells[i]];
|
|
}
|
|
fldBf[patchi] == pfld;
|
|
}
|
|
Info<< indent << fld.name() << " : actual number of layers"
|
|
<< endl;
|
|
bool ok = fld.write();
|
|
allOk = allOk && ok;
|
|
}
|
|
{
|
|
volScalarField fld
|
|
(
|
|
IOobject
|
|
(
|
|
"thickness",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE,
|
|
false
|
|
),
|
|
mesh,
|
|
dimensionedScalar(dimless, Zero),
|
|
fixedValueFvPatchScalarField::typeName
|
|
);
|
|
volScalarField::Boundary& fldBf = fld.boundaryFieldRef();
|
|
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
fldBf[patchi] == pbm[patchi].patchSlice(faceRealThickness);
|
|
}
|
|
Info<< indent << fld.name() << " : overall layer thickness"
|
|
<< endl;
|
|
bool ok = fld.write();
|
|
allOk = allOk && ok;
|
|
}
|
|
{
|
|
volScalarField fld
|
|
(
|
|
IOobject
|
|
(
|
|
"thicknessFraction",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE,
|
|
false
|
|
),
|
|
mesh,
|
|
dimensionedScalar(dimless, Zero),
|
|
fixedValueFvPatchScalarField::typeName
|
|
);
|
|
volScalarField::Boundary& fldBf = fld.boundaryFieldRef();
|
|
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchi = patchIDs[i];
|
|
|
|
scalarField::subField patchWanted = pbm[patchi].patchSlice
|
|
(
|
|
faceWantedThickness
|
|
);
|
|
scalarField::subField patchReal = pbm[patchi].patchSlice
|
|
(
|
|
faceRealThickness
|
|
);
|
|
|
|
// Convert patchReal to relative thickness
|
|
scalarField pfld(patchReal.size(), Zero);
|
|
forAll(patchReal, i)
|
|
{
|
|
if (patchWanted[i] > VSMALL)
|
|
{
|
|
pfld[i] = patchReal[i]/patchWanted[i];
|
|
}
|
|
}
|
|
|
|
fldBf[patchi] == pfld;
|
|
}
|
|
Info<< indent << fld.name()
|
|
<< " : overall layer thickness (fraction"
|
|
<< " of desired thickness)" << endl;
|
|
bool ok = fld.write();
|
|
allOk = allOk && ok;
|
|
}
|
|
Info<< decrIndent<< endl;
|
|
}
|
|
|
|
return allOk;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::snappyLayerDriver::snappyLayerDriver
|
|
(
|
|
meshRefinement& meshRefiner,
|
|
const labelList& globalToMasterPatch,
|
|
const labelList& globalToSlavePatch,
|
|
const bool dryRun
|
|
)
|
|
:
|
|
meshRefiner_(meshRefiner),
|
|
globalToMasterPatch_(globalToMasterPatch),
|
|
globalToSlavePatch_(globalToSlavePatch),
|
|
dryRun_(dryRun)
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::snappyLayerDriver::mergePatchFacesUndo
|
|
(
|
|
const layerParameters& layerParams,
|
|
const dictionary& motionDict,
|
|
const meshRefinement::FaceMergeType mergeType
|
|
)
|
|
{
|
|
// Clip to 30 degrees. Not helpful!
|
|
//scalar planarAngle = min(30.0, layerParams.featureAngle());
|
|
scalar planarAngle = layerParams.mergePatchFacesAngle();
|
|
scalar minCos = Foam::cos(degToRad(planarAngle));
|
|
|
|
scalar concaveCos = Foam::cos(degToRad(layerParams.concaveAngle()));
|
|
|
|
Info<< nl
|
|
<< "Merging all faces of a cell" << nl
|
|
<< "---------------------------" << nl
|
|
<< " - which are on the same patch" << nl
|
|
<< " - which make an angle < " << planarAngle
|
|
<< " degrees"
|
|
<< " (cos:" << minCos << ')' << nl
|
|
<< " - as long as the resulting face doesn't become concave"
|
|
<< " by more than "
|
|
<< layerParams.concaveAngle() << " degrees" << nl
|
|
<< " (0=straight, 180=fully concave)" << nl
|
|
<< endl;
|
|
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
|
|
|
|
labelList duplicateFace(mesh.nFaces(), -1);
|
|
forAll(couples, i)
|
|
{
|
|
const labelPair& cpl = couples[i];
|
|
duplicateFace[cpl[0]] = cpl[1];
|
|
duplicateFace[cpl[1]] = cpl[0];
|
|
}
|
|
|
|
label nChanged = meshRefiner_.mergePatchFacesUndo
|
|
(
|
|
minCos,
|
|
concaveCos,
|
|
meshRefiner_.meshedPatches(),
|
|
motionDict,
|
|
duplicateFace,
|
|
mergeType // How to merge co-planar patch faces
|
|
);
|
|
|
|
nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::addLayers
|
|
(
|
|
const layerParameters& layerParams,
|
|
const dictionary& motionDict,
|
|
const labelList& patchIDs,
|
|
const label nAllowableErrors,
|
|
decompositionMethod& decomposer,
|
|
fvMeshDistribute& distributor
|
|
)
|
|
{
|
|
fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
|
|
// faceZones of type internal or baffle (for merging points across)
|
|
labelList internalOrBaffleFaceZones;
|
|
{
|
|
List<surfaceZonesInfo::faceZoneType> fzTypes(2);
|
|
fzTypes[0] = surfaceZonesInfo::INTERNAL;
|
|
fzTypes[1] = surfaceZonesInfo::BAFFLE;
|
|
internalOrBaffleFaceZones = meshRefiner_.getZones(fzTypes);
|
|
}
|
|
|
|
// faceZones of type internal (for checking mesh quality across and
|
|
// merging baffles)
|
|
const labelList internalFaceZones
|
|
(
|
|
meshRefiner_.getZones
|
|
(
|
|
List<surfaceZonesInfo::faceZoneType>
|
|
(
|
|
1,
|
|
surfaceZonesInfo::INTERNAL
|
|
)
|
|
)
|
|
);
|
|
|
|
// Create baffles (pairs of faces that share the same points)
|
|
// Baffles stored as owner and neighbour face that have been created.
|
|
List<labelPair> baffles;
|
|
{
|
|
labelList originatingFaceZone;
|
|
meshRefiner_.createZoneBaffles
|
|
(
|
|
identity(mesh.faceZones().size()),
|
|
baffles,
|
|
originatingFaceZone
|
|
);
|
|
|
|
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing baffled mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Duplicate points on faceZones of type boundary. Should normally already
|
|
// be done by snapping phase
|
|
{
|
|
autoPtr<mapPolyMesh> map = meshRefiner_.dupNonManifoldBoundaryPoints();
|
|
if (map)
|
|
{
|
|
const labelList& reverseFaceMap = map->reverseFaceMap();
|
|
forAll(baffles, i)
|
|
{
|
|
label f0 = reverseFaceMap[baffles[i].first()];
|
|
label f1 = reverseFaceMap[baffles[i].second()];
|
|
baffles[i] = labelPair(f0, f1);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Now we have all patches determine the number of layer per patch
|
|
// This will be layerParams.numLayers() except for faceZones where one
|
|
// side does get layers and the other not in which case we want to
|
|
// suppress movement by explicitly setting numLayers 0
|
|
labelList numLayers(layerParams.numLayers());
|
|
{
|
|
labelHashSet layerIDs(patchIDs);
|
|
forAll(mesh.faceZones(), zonei)
|
|
{
|
|
label mpi, spi;
|
|
surfaceZonesInfo::faceZoneType fzType;
|
|
bool hasInfo = meshRefiner_.getFaceZoneInfo
|
|
(
|
|
mesh.faceZones()[zonei].name(),
|
|
mpi,
|
|
spi,
|
|
fzType
|
|
);
|
|
if (hasInfo)
|
|
{
|
|
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
|
if (layerIDs.found(mpi) && !layerIDs.found(spi))
|
|
{
|
|
// Only layers on master side. Fix points on slave side
|
|
Info<< "On faceZone " << mesh.faceZones()[zonei].name()
|
|
<< " adding layers to master patch " << pbm[mpi].name()
|
|
<< " only. Freezing points on slave patch "
|
|
<< pbm[spi].name() << endl;
|
|
numLayers[spi] = 0;
|
|
}
|
|
else if (!layerIDs.found(mpi) && layerIDs.found(spi))
|
|
{
|
|
// Only layers on slave side. Fix points on master side
|
|
Info<< "On faceZone " << mesh.faceZones()[zonei].name()
|
|
<< " adding layers to slave patch " << pbm[spi].name()
|
|
<< " only. Freezing points on master patch "
|
|
<< pbm[mpi].name() << endl;
|
|
numLayers[mpi] = 0;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Duplicate points on faceZones that layers are added to
|
|
labelList pointToMaster;
|
|
|
|
{
|
|
// Check outside of baffles for non-manifoldness
|
|
|
|
// Points that are candidates for duplication
|
|
labelList candidatePoints;
|
|
{
|
|
// Do full analysis to see if we need to extrude points
|
|
// so have to duplicate them
|
|
autoPtr<indirectPrimitivePatch> pp
|
|
(
|
|
meshRefinement::makePatch
|
|
(
|
|
mesh,
|
|
patchIDs
|
|
)
|
|
);
|
|
|
|
// Displacement for all pp.localPoints. Set to large value to
|
|
// avoid truncation in syncPatchDisplacement because of
|
|
// minThickness.
|
|
vectorField patchDisp(pp().nPoints(), vector(GREAT, GREAT, GREAT));
|
|
labelList patchNLayers(pp().nPoints(), Zero);
|
|
label nIdealTotAddedCells = 0;
|
|
List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
|
|
// Get number of layers per point from number of layers per patch
|
|
setNumLayers
|
|
(
|
|
numLayers, // per patch the num layers
|
|
patchIDs, // patches that are being moved
|
|
*pp, // indirectpatch for all faces moving
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus,
|
|
nIdealTotAddedCells
|
|
);
|
|
// Make sure displacement is equal on both sides of coupled patches.
|
|
// Note that we explicitly disable the minThickness truncation
|
|
// of the patchDisp here.
|
|
syncPatchDisplacement
|
|
(
|
|
*pp,
|
|
scalarField(patchDisp.size(), Zero), //minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
|
|
// Do duplication only if all patch points decide to extrude. Ignore
|
|
// contribution from non-patch points. Note that we need to
|
|
// apply this to all mesh points
|
|
labelList minPatchState(mesh.nPoints(), labelMax);
|
|
forAll(extrudeStatus, patchPointi)
|
|
{
|
|
label pointi = pp().meshPoints()[patchPointi];
|
|
minPatchState[pointi] = extrudeStatus[patchPointi];
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
minPatchState,
|
|
minEqOp<label>(), // combine op
|
|
labelMax // null value
|
|
);
|
|
|
|
// So now minPatchState:
|
|
// - labelMax on non-patch points
|
|
// - NOEXTRUDE if any patch point was not extruded
|
|
// - EXTRUDE or EXTRUDEREMOVE if all patch points are extruded/
|
|
// extrudeRemove.
|
|
|
|
label n = 0;
|
|
forAll(minPatchState, pointi)
|
|
{
|
|
label state = minPatchState[pointi];
|
|
if (state == EXTRUDE || state == EXTRUDEREMOVE)
|
|
{
|
|
n++;
|
|
}
|
|
}
|
|
candidatePoints.setSize(n);
|
|
n = 0;
|
|
forAll(minPatchState, pointi)
|
|
{
|
|
label state = minPatchState[pointi];
|
|
if (state == EXTRUDE || state == EXTRUDEREMOVE)
|
|
{
|
|
candidatePoints[n++] = pointi;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Not duplicate points on either side of baffles that don't get any
|
|
// layers
|
|
labelPairList nonDupBaffles;
|
|
|
|
{
|
|
// faceZones that are not being duplicated
|
|
DynamicList<label> nonDupZones(mesh.faceZones().size());
|
|
|
|
labelHashSet layerIDs(patchIDs);
|
|
forAll(mesh.faceZones(), zonei)
|
|
{
|
|
label mpi, spi;
|
|
surfaceZonesInfo::faceZoneType fzType;
|
|
bool hasInfo = meshRefiner_.getFaceZoneInfo
|
|
(
|
|
mesh.faceZones()[zonei].name(),
|
|
mpi,
|
|
spi,
|
|
fzType
|
|
);
|
|
if (hasInfo && !layerIDs.found(mpi) && !layerIDs.found(spi))
|
|
{
|
|
nonDupZones.append(zonei);
|
|
}
|
|
}
|
|
nonDupBaffles = meshRefinement::subsetBaffles
|
|
(
|
|
mesh,
|
|
nonDupZones,
|
|
localPointRegion::findDuplicateFacePairs(mesh)
|
|
);
|
|
}
|
|
|
|
|
|
const localPointRegion regionSide(mesh, nonDupBaffles, candidatePoints);
|
|
|
|
autoPtr<mapPolyMesh> map = meshRefiner_.dupNonManifoldPoints
|
|
(
|
|
regionSide
|
|
);
|
|
|
|
if (map)
|
|
{
|
|
// Store point duplication
|
|
pointToMaster.setSize(mesh.nPoints(), -1);
|
|
|
|
const labelList& pointMap = map().pointMap();
|
|
const labelList& reversePointMap = map().reversePointMap();
|
|
|
|
forAll(pointMap, pointi)
|
|
{
|
|
label oldPointi = pointMap[pointi];
|
|
label newMasterPointi = reversePointMap[oldPointi];
|
|
|
|
if (newMasterPointi != pointi)
|
|
{
|
|
// Found slave. Mark both master and slave
|
|
pointToMaster[pointi] = newMasterPointi;
|
|
pointToMaster[newMasterPointi] = newMasterPointi;
|
|
}
|
|
}
|
|
|
|
// Update baffle numbering
|
|
{
|
|
const labelList& reverseFaceMap = map().reverseFaceMap();
|
|
forAll(baffles, i)
|
|
{
|
|
label f0 = reverseFaceMap[baffles[i].first()];
|
|
label f1 = reverseFaceMap[baffles[i].second()];
|
|
baffles[i] = labelPair(f0, f1);
|
|
}
|
|
}
|
|
|
|
|
|
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing point-duplicate mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
|
|
OBJstream str
|
|
(
|
|
mesh.time().path()
|
|
/ "duplicatePoints_"
|
|
+ meshRefiner_.timeName()
|
|
+ ".obj"
|
|
);
|
|
Info<< "Writing point-duplicates to " << str.name() << endl;
|
|
const pointField& p = mesh.points();
|
|
forAll(pointMap, pointi)
|
|
{
|
|
label newMasteri = reversePointMap[pointMap[pointi]];
|
|
|
|
if (newMasteri != pointi)
|
|
{
|
|
str.write(linePointRef(p[pointi], p[newMasteri]));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Add layers to patches
|
|
// ~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Now we have
|
|
// - mesh with optional baffles and duplicated points for faceZones
|
|
// where layers are to be added
|
|
// - pointToMaster : correspondence for duplicated points
|
|
// - baffles : list of pairs of faces
|
|
|
|
|
|
autoPtr<indirectPrimitivePatch> pp
|
|
(
|
|
meshRefinement::makePatch
|
|
(
|
|
mesh,
|
|
patchIDs
|
|
)
|
|
);
|
|
|
|
|
|
// Global face indices engine
|
|
const globalIndex globalFaces(mesh.nFaces());
|
|
|
|
// Determine extrudePatch.edgeFaces in global numbering (so across
|
|
// coupled patches). This is used only to string up edges between coupled
|
|
// faces (all edges between same (global)face indices get extruded).
|
|
labelListList edgeGlobalFaces
|
|
(
|
|
addPatchCellLayer::globalEdgeFaces
|
|
(
|
|
mesh,
|
|
globalFaces,
|
|
*pp
|
|
)
|
|
);
|
|
|
|
// Determine patches for extruded boundary edges. Calculates if any
|
|
// additional processor patches need to be constructed.
|
|
|
|
labelList edgePatchID;
|
|
labelList edgeZoneID;
|
|
boolList edgeFlip;
|
|
labelList inflateFaceID;
|
|
determineSidePatches
|
|
(
|
|
globalFaces,
|
|
edgeGlobalFaces,
|
|
*pp,
|
|
|
|
edgePatchID,
|
|
edgeZoneID,
|
|
edgeFlip,
|
|
inflateFaceID
|
|
);
|
|
|
|
|
|
// Point-wise extrusion data
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Displacement for all pp.localPoints. Set to large value so does
|
|
// not trigger the minThickness truncation (see syncPatchDisplacement below)
|
|
vectorField patchDisp(pp().nPoints(), vector(GREAT, GREAT, GREAT));
|
|
|
|
// Number of layers for all pp.localPoints. Note: only valid if
|
|
// extrudeStatus = EXTRUDE.
|
|
labelList patchNLayers(pp().nPoints(), Zero);
|
|
|
|
// Ideal number of cells added
|
|
label nIdealTotAddedCells = 0;
|
|
|
|
// Whether to add edge for all pp.localPoints.
|
|
List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
|
|
|
|
|
|
{
|
|
// Get number of layers per point from number of layers per patch
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
setNumLayers
|
|
(
|
|
numLayers, // per patch the num layers
|
|
patchIDs, // patches that are being moved
|
|
*pp, // indirectpatch for all faces moving
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus,
|
|
nIdealTotAddedCells
|
|
);
|
|
|
|
// Precalculate mesh edge labels for patch edges
|
|
labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
|
|
|
|
|
|
// Disable extrusion on split strings of common points
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
handleNonStringConnected
|
|
(
|
|
*pp,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
|
|
// Disable extrusion on non-manifold points
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
handleNonManifolds
|
|
(
|
|
*pp,
|
|
meshEdges,
|
|
edgeGlobalFaces,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Disable extrusion on feature angles
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
handleFeatureAngle
|
|
(
|
|
*pp,
|
|
meshEdges,
|
|
layerParams.featureAngle(),
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Disable extrusion on warped faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
// It is hard to calculate some length scale if not in relative
|
|
// mode so disable this check.
|
|
if (!layerParams.relativeSizes().found(false))
|
|
{
|
|
// Undistorted edge length
|
|
const scalar edge0Len =
|
|
meshRefiner_.meshCutter().level0EdgeLength();
|
|
const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
|
|
|
|
handleWarpedFaces
|
|
(
|
|
*pp,
|
|
layerParams.maxFaceThicknessRatio(),
|
|
layerParams.relativeSizes(),
|
|
edge0Len,
|
|
cellLevel,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
|
|
//// Disable extrusion on cells with multiple patch faces
|
|
//// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
//
|
|
//handleMultiplePatchFaces
|
|
//(
|
|
// *pp,
|
|
//
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
//);
|
|
|
|
addProfiling(grow, "snappyHexMesh::layers::grow");
|
|
|
|
// Grow out region of non-extrusion
|
|
for (label i = 0; i < layerParams.nGrow(); i++)
|
|
{
|
|
growNoExtrusion
|
|
(
|
|
*pp,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Undistorted edge length
|
|
const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
|
|
const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
|
|
|
|
// Determine (wanted) point-wise overall layer thickness and expansion
|
|
// ratio
|
|
scalarField thickness(pp().nPoints());
|
|
scalarIOField minThickness
|
|
(
|
|
IOobject
|
|
(
|
|
"minThickness",
|
|
meshRefiner_.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ
|
|
),
|
|
pp().nPoints()
|
|
);
|
|
scalarField expansionRatio(pp().nPoints());
|
|
calculateLayerThickness
|
|
(
|
|
*pp,
|
|
patchIDs,
|
|
layerParams,
|
|
cellLevel,
|
|
patchNLayers,
|
|
edge0Len,
|
|
|
|
thickness,
|
|
minThickness,
|
|
expansionRatio
|
|
);
|
|
|
|
|
|
|
|
// Current set of topology changes. (changing mesh clears out
|
|
// polyTopoChange)
|
|
polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
|
|
// Per cell 0 or number of layers in the cell column it is part of
|
|
labelList cellNLayers;
|
|
// Per face actual overall layer thickness
|
|
scalarField faceRealThickness;
|
|
// Per face wanted overall layer thickness
|
|
scalarField faceWantedThickness(mesh.nFaces(), Zero);
|
|
{
|
|
UIndirectList<scalar>(faceWantedThickness, pp->addressing()) =
|
|
avgPointData(*pp, thickness);
|
|
}
|
|
|
|
|
|
{
|
|
// Overall displacement field
|
|
pointVectorField displacement
|
|
(
|
|
makeLayerDisplacementField
|
|
(
|
|
pointMesh::New(mesh),
|
|
numLayers
|
|
)
|
|
);
|
|
|
|
// Allocate run-time selectable mesh mover
|
|
autoPtr<externalDisplacementMeshMover> medialAxisMoverPtr;
|
|
{
|
|
// Set up controls for meshMover
|
|
dictionary combinedDict(layerParams.dict());
|
|
// Add mesh quality constraints
|
|
combinedDict.merge(motionDict);
|
|
// Where to get minThickness from
|
|
combinedDict.add("minThicknessName", minThickness.name());
|
|
|
|
const List<labelPair> internalBaffles
|
|
(
|
|
meshRefinement::subsetBaffles
|
|
(
|
|
mesh,
|
|
internalFaceZones,
|
|
localPointRegion::findDuplicateFacePairs(mesh)
|
|
)
|
|
);
|
|
|
|
// Take over patchDisp as boundary conditions on displacement
|
|
// pointVectorField
|
|
medialAxisMoverPtr = externalDisplacementMeshMover::New
|
|
(
|
|
layerParams.meshShrinker(),
|
|
combinedDict,
|
|
internalBaffles,
|
|
displacement
|
|
);
|
|
|
|
|
|
if (dryRun_)
|
|
{
|
|
string errorMsg(FatalError.message());
|
|
string IOerrorMsg(FatalIOError.message());
|
|
|
|
if (errorMsg.size() || IOerrorMsg.size())
|
|
{
|
|
//errorMsg = "[dryRun] " + errorMsg;
|
|
//errorMsg.replaceAll("\n", "\n[dryRun] ");
|
|
//IOerrorMsg = "[dryRun] " + IOerrorMsg;
|
|
//IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
|
|
|
|
IOWarningInFunction(combinedDict)
|
|
<< nl
|
|
<< "Missing/incorrect required dictionary entries:"
|
|
<< nl << nl
|
|
<< IOerrorMsg.c_str() << nl
|
|
<< errorMsg.c_str() << nl << nl
|
|
<< "Exiting dry-run" << nl << endl;
|
|
|
|
if (Pstream::parRun())
|
|
{
|
|
Perr<< "\nFOAM parallel run exiting\n" << endl;
|
|
Pstream::exit(0);
|
|
}
|
|
else
|
|
{
|
|
Perr<< "\nFOAM exiting\n" << endl;
|
|
std::exit(0);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Saved old points
|
|
const pointField oldPoints(mesh.points());
|
|
|
|
for
|
|
(
|
|
label iteration = 0;
|
|
iteration < layerParams.nLayerIter();
|
|
iteration++
|
|
)
|
|
{
|
|
Info<< nl
|
|
<< "Layer addition iteration " << iteration << nl
|
|
<< "--------------------------" << endl;
|
|
|
|
|
|
// Unset the extrusion at the pp.
|
|
const dictionary& meshQualityDict =
|
|
(
|
|
iteration < layerParams.nRelaxedIter()
|
|
? motionDict
|
|
: motionDict.subDict("relaxed")
|
|
);
|
|
|
|
if (iteration >= layerParams.nRelaxedIter())
|
|
{
|
|
Info<< "Switched to relaxed meshQuality constraints." << endl;
|
|
}
|
|
|
|
|
|
|
|
// Make sure displacement is equal on both sides of coupled patches.
|
|
// Note that this also does the patchDisp < minThickness truncation
|
|
// so for the first pass make sure the patchDisp is larger than
|
|
// that.
|
|
syncPatchDisplacement
|
|
(
|
|
*pp,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Displacement acc. to pointnormals
|
|
getPatchDisplacement
|
|
(
|
|
*pp,
|
|
thickness,
|
|
minThickness,
|
|
expansionRatio,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Shrink mesh by displacement value first.
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
{
|
|
const pointField oldPatchPos(pp().localPoints());
|
|
|
|
// We have patchDisp which is the outwards pointing
|
|
// extrusion distance. Convert into an inwards pointing
|
|
// shrink distance
|
|
patchDisp = -patchDisp;
|
|
|
|
// Take over patchDisp into pointDisplacement field and
|
|
// adjust both for multi-patch constraints
|
|
motionSmootherAlgo::setDisplacement
|
|
(
|
|
patchIDs,
|
|
*pp,
|
|
patchDisp,
|
|
displacement
|
|
);
|
|
|
|
|
|
// Move mesh
|
|
// ~~~~~~~~~
|
|
|
|
// Set up controls for meshMover
|
|
dictionary combinedDict(layerParams.dict());
|
|
// Add standard quality constraints
|
|
combinedDict.merge(motionDict);
|
|
// Add relaxed constraints (overrides standard ones)
|
|
combinedDict.merge(meshQualityDict);
|
|
// Where to get minThickness from
|
|
combinedDict.add("minThicknessName", minThickness.name());
|
|
|
|
labelList checkFaces(identity(mesh.nFaces()));
|
|
medialAxisMoverPtr().move
|
|
(
|
|
combinedDict,
|
|
nAllowableErrors,
|
|
checkFaces
|
|
);
|
|
|
|
pp().movePoints(mesh.points());
|
|
|
|
// Update patchDisp (since not all might have been honoured)
|
|
patchDisp = oldPatchPos - pp().localPoints();
|
|
}
|
|
|
|
// Truncate displacements that are too small (this will do internal
|
|
// ones, coupled ones have already been truncated by
|
|
// syncPatchDisplacement)
|
|
faceSet dummySet(mesh, "wrongPatchFaces", 0);
|
|
truncateDisplacement
|
|
(
|
|
globalFaces,
|
|
edgeGlobalFaces,
|
|
*pp,
|
|
minThickness,
|
|
dummySet,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
|
|
// Dump to .obj file for debugging.
|
|
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
|
|
{
|
|
dumpDisplacement
|
|
(
|
|
mesh.time().path()/"layer_" + meshRefiner_.timeName(),
|
|
pp(),
|
|
patchDisp,
|
|
extrudeStatus
|
|
);
|
|
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing shrunk mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
|
|
// See comment in snappySnapDriver why we should not remove
|
|
// meshPhi using mesh.clearOut().
|
|
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
}
|
|
|
|
|
|
// Mesh topo change engine. Insert current mesh.
|
|
polyTopoChange meshMod(mesh);
|
|
|
|
// Grow layer of cells on to patch. Handles zero sized displacement.
|
|
addPatchCellLayer addLayer(mesh);
|
|
|
|
// Determine per point/per face number of layers to extrude. Also
|
|
// handles the slow termination of layers when going switching
|
|
// layers
|
|
|
|
labelList nPatchPointLayers(pp().nPoints(), -1);
|
|
labelList nPatchFaceLayers(pp().size(), -1);
|
|
setupLayerInfoTruncation
|
|
(
|
|
*pp,
|
|
patchNLayers,
|
|
extrudeStatus,
|
|
layerParams.nBufferCellsNoExtrude(),
|
|
nPatchPointLayers,
|
|
nPatchFaceLayers
|
|
);
|
|
|
|
// Calculate displacement for final layer for addPatchLayer.
|
|
// (layer of cells next to the original mesh)
|
|
vectorField finalDisp(patchNLayers.size(), Zero);
|
|
|
|
forAll(nPatchPointLayers, i)
|
|
{
|
|
scalar ratio = layerParameters::finalLayerThicknessRatio
|
|
(
|
|
nPatchPointLayers[i],
|
|
expansionRatio[i]
|
|
);
|
|
finalDisp[i] = ratio*patchDisp[i];
|
|
}
|
|
|
|
|
|
const scalarField invExpansionRatio(1.0/expansionRatio);
|
|
|
|
// Add topo regardless of whether extrudeStatus is extruderemove.
|
|
// Not add layer if patchDisp is zero.
|
|
addLayer.setRefinement
|
|
(
|
|
globalFaces,
|
|
edgeGlobalFaces,
|
|
|
|
invExpansionRatio,
|
|
pp(),
|
|
|
|
edgePatchID, // boundary patch for extruded boundary edges
|
|
edgeZoneID, // zone for extruded edges
|
|
edgeFlip,
|
|
inflateFaceID,
|
|
|
|
|
|
labelList(0), // exposed patchIDs, not used for adding layers
|
|
nPatchFaceLayers, // layers per face
|
|
nPatchPointLayers, // layers per point
|
|
finalDisp, // thickness of layer nearest internal mesh
|
|
meshMod
|
|
);
|
|
|
|
if (debug)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
}
|
|
|
|
// Store mesh changes for if mesh is correct.
|
|
savedMeshMod = meshMod;
|
|
|
|
|
|
// With the stored topo changes we create a new mesh so we can
|
|
// undo if necessary.
|
|
|
|
autoPtr<fvMesh> newMeshPtr;
|
|
autoPtr<mapPolyMesh> mapPtr = meshMod.makeMesh
|
|
(
|
|
newMeshPtr,
|
|
IOobject
|
|
(
|
|
//mesh.name()+"_layer",
|
|
mesh.name(),
|
|
static_cast<polyMesh&>(mesh).instance(),
|
|
mesh.time(), // register with runTime
|
|
IOobject::READ_IF_PRESENT, // read fv* if present
|
|
static_cast<polyMesh&>(mesh).writeOpt()
|
|
), // io params from original mesh but new name
|
|
mesh, // original mesh
|
|
true // parallel sync
|
|
);
|
|
fvMesh& newMesh = *newMeshPtr;
|
|
mapPolyMesh& map = *mapPtr;
|
|
|
|
// Get timing, but more importantly get memory information
|
|
addProfiling(grow, "snappyHexMesh::layers::updateMesh");
|
|
|
|
//?necessary? Update fields
|
|
newMesh.updateMesh(map);
|
|
|
|
newMesh.setInstance(meshRefiner_.timeName());
|
|
|
|
// Update numbering on addLayer:
|
|
// - cell/point labels to be newMesh.
|
|
// - patchFaces to remain in oldMesh order.
|
|
addLayer.updateMesh
|
|
(
|
|
map,
|
|
identity(pp().size()),
|
|
identity(pp().nPoints())
|
|
);
|
|
|
|
// Collect layer faces and cells for outside loop.
|
|
getLayerCellsFaces
|
|
(
|
|
newMesh,
|
|
addLayer,
|
|
avgPointData(*pp, mag(patchDisp))(), // current thickness
|
|
|
|
cellNLayers,
|
|
faceRealThickness
|
|
);
|
|
|
|
|
|
// Count number of added cells
|
|
label nAddedCells = 0;
|
|
forAll(cellNLayers, celli)
|
|
{
|
|
if (cellNLayers[celli] > 0)
|
|
{
|
|
nAddedCells++;
|
|
}
|
|
}
|
|
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
|
|
<< endl;
|
|
newMesh.write();
|
|
writeLayerSets(newMesh, cellNLayers, faceRealThickness);
|
|
|
|
// Reset the instance of the original mesh so next iteration
|
|
// it dumps a complete mesh. This is just so that the inbetween
|
|
// newMesh does not upset e.g. paraFoam cycling through the
|
|
// times.
|
|
mesh.setInstance(meshRefiner_.timeName());
|
|
}
|
|
|
|
|
|
//- Get baffles in newMesh numbering. Note that we cannot detect
|
|
// baffles here since the points are duplicated
|
|
List<labelPair> internalBaffles;
|
|
{
|
|
// From old mesh face to corresponding newMesh boundary face
|
|
labelList meshToNewMesh(mesh.nFaces(), -1);
|
|
for
|
|
(
|
|
label facei = newMesh.nInternalFaces();
|
|
facei < newMesh.nFaces();
|
|
facei++
|
|
)
|
|
{
|
|
label newMeshFacei = map.faceMap()[facei];
|
|
if (newMeshFacei != -1)
|
|
{
|
|
meshToNewMesh[newMeshFacei] = facei;
|
|
}
|
|
}
|
|
|
|
List<labelPair> newMeshBaffles(baffles.size());
|
|
label newi = 0;
|
|
forAll(baffles, i)
|
|
{
|
|
const labelPair& p = baffles[i];
|
|
labelPair newMeshBaffle
|
|
(
|
|
meshToNewMesh[p[0]],
|
|
meshToNewMesh[p[1]]
|
|
);
|
|
if (newMeshBaffle[0] != -1 && newMeshBaffle[1] != -1)
|
|
{
|
|
newMeshBaffles[newi++] = newMeshBaffle;
|
|
}
|
|
}
|
|
newMeshBaffles.setSize(newi);
|
|
|
|
internalBaffles = meshRefinement::subsetBaffles
|
|
(
|
|
newMesh,
|
|
internalFaceZones,
|
|
newMeshBaffles
|
|
);
|
|
|
|
Info<< "Detected "
|
|
<< returnReduce(internalBaffles.size(), sumOp<label>())
|
|
<< " baffles across faceZones of type internal" << nl
|
|
<< endl;
|
|
}
|
|
|
|
label nTotChanged = checkAndUnmark
|
|
(
|
|
addLayer,
|
|
meshQualityDict,
|
|
layerParams.additionalReporting(),
|
|
internalBaffles,
|
|
pp(),
|
|
newMesh,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
label nTotExtruded = countExtrusion(*pp, extrudeStatus);
|
|
label nTotFaces = returnReduce(pp().size(), sumOp<label>());
|
|
label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
|
|
|
|
Info<< "Extruding " << nTotExtruded
|
|
<< " out of " << nTotFaces
|
|
<< " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
|
|
<< " Removed extrusion at " << nTotChanged << " faces."
|
|
<< endl
|
|
<< "Added " << nTotAddedCells << " out of "
|
|
<< nIdealTotAddedCells
|
|
<< " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells
|
|
<< "%)." << endl;
|
|
|
|
if (nTotChanged == 0)
|
|
{
|
|
break;
|
|
}
|
|
|
|
// Reset mesh points and start again
|
|
mesh.movePoints(oldPoints);
|
|
pp().movePoints(mesh.points());
|
|
medialAxisMoverPtr().movePoints(mesh.points());
|
|
|
|
// Grow out region of non-extrusion
|
|
for (label i = 0; i < layerParams.nGrow(); i++)
|
|
{
|
|
growNoExtrusion
|
|
(
|
|
*pp,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
|
|
Info<< endl;
|
|
}
|
|
}
|
|
|
|
|
|
// At this point we have a (shrunk) mesh and a set of topology changes
|
|
// which will make a valid mesh with layer. Apply these changes to the
|
|
// current mesh.
|
|
|
|
{
|
|
// Apply the stored topo changes to the current mesh.
|
|
autoPtr<mapPolyMesh> mapPtr = savedMeshMod.changeMesh(mesh, false);
|
|
mapPolyMesh& map = *mapPtr;
|
|
|
|
// Hack to remove meshPhi - mapped incorrectly. TBD.
|
|
mesh.clearOut();
|
|
|
|
// Update fields
|
|
mesh.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map.hasMotionPoints())
|
|
{
|
|
mesh.movePoints(map.preMotionPoints());
|
|
}
|
|
else
|
|
{
|
|
// Delete mesh volumes.
|
|
mesh.clearOut();
|
|
}
|
|
|
|
// Reset the instance for if in overwrite mode
|
|
mesh.setInstance(meshRefiner_.timeName());
|
|
|
|
meshRefiner_.updateMesh(map, labelList(0));
|
|
|
|
// Update numbering of faceWantedThickness
|
|
meshRefinement::updateList
|
|
(
|
|
map.faceMap(),
|
|
scalar(0),
|
|
faceWantedThickness
|
|
);
|
|
|
|
// Print data now that we still have patches for the zones
|
|
//if (meshRefinement::outputLevel() & meshRefinement::OUTPUTLAYERINFO)
|
|
printLayerData
|
|
(
|
|
mesh,
|
|
patchIDs,
|
|
cellNLayers,
|
|
faceWantedThickness,
|
|
faceRealThickness
|
|
);
|
|
|
|
|
|
// Dump for debugging
|
|
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing mesh with layers but disconnected to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
}
|
|
|
|
|
|
// Use geometric detection of points-to-be-merged
|
|
// - detect any boundary face created from a duplicated face (=baffle)
|
|
// - on these mark any point created from a duplicated point
|
|
if (returnReduce(pointToMaster.size(), sumOp<label>()))
|
|
{
|
|
// Estimate number of points-to-be-merged
|
|
DynamicList<label> candidates(baffles.size()*4);
|
|
|
|
// Mark whether old face was on baffle
|
|
bitSet oldBaffleFace(map.nOldFaces());
|
|
forAll(baffles, i)
|
|
{
|
|
const labelPair& baffle = baffles[i];
|
|
oldBaffleFace.set(baffle[0]);
|
|
oldBaffleFace.set(baffle[1]);
|
|
}
|
|
|
|
// Collect candidate if
|
|
// - point on boundary face originating from baffle
|
|
// - and point originating from duplicate
|
|
for
|
|
(
|
|
label facei = mesh.nInternalFaces();
|
|
facei < mesh.nFaces();
|
|
facei++
|
|
)
|
|
{
|
|
label oldFacei = map.faceMap()[facei];
|
|
if (oldFacei != -1 && oldBaffleFace.test(oldFacei))
|
|
{
|
|
const face& f = mesh.faces()[facei];
|
|
forAll(f, fp)
|
|
{
|
|
label pointi = f[fp];
|
|
label oldPointi = map.pointMap()[pointi];
|
|
|
|
if (pointToMaster[oldPointi] != -1)
|
|
{
|
|
candidates.append(pointi);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Do geometric merge. Ideally we'd like to use a topological
|
|
// merge but we've thrown away all layer-wise addressing when
|
|
// throwing away addPatchCellLayer engine. Also the addressing
|
|
// is extremely complicated. There is no problem with merging
|
|
// too many points; the problem would be if merging baffles.
|
|
// Trust mergeZoneBaffles to do sufficient checks.
|
|
labelList oldToNew;
|
|
label nNew = mergePoints
|
|
(
|
|
pointField(mesh.points(), candidates),
|
|
meshRefiner_.mergeDistance(),
|
|
false,
|
|
oldToNew
|
|
);
|
|
|
|
// Extract points to be merged (i.e. multiple points originating
|
|
// from a single one)
|
|
|
|
labelListList newToOld(invertOneToMany(nNew, oldToNew));
|
|
|
|
// Extract points with more than one old one
|
|
pointToMaster.setSize(mesh.nPoints());
|
|
pointToMaster = -1;
|
|
|
|
forAll(newToOld, newi)
|
|
{
|
|
const labelList& oldPoints = newToOld[newi];
|
|
if (oldPoints.size() > 1)
|
|
{
|
|
labelList meshPoints
|
|
(
|
|
labelUIndList(candidates, oldPoints)
|
|
);
|
|
label masteri = min(meshPoints);
|
|
forAll(meshPoints, i)
|
|
{
|
|
pointToMaster[meshPoints[i]] = masteri;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// Count duplicate points
|
|
label nPointPairs = 0;
|
|
forAll(pointToMaster, pointi)
|
|
{
|
|
label otherPointi = pointToMaster[pointi];
|
|
if (otherPointi != -1)
|
|
{
|
|
nPointPairs++;
|
|
}
|
|
}
|
|
reduce(nPointPairs, sumOp<label>());
|
|
if (nPointPairs > 0)
|
|
{
|
|
// Merge any duplicated points
|
|
Info<< "Merging " << nPointPairs << " duplicated points ..." << endl;
|
|
|
|
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
|
|
{
|
|
OBJstream str
|
|
(
|
|
mesh.time().path()
|
|
/ "mergePoints_"
|
|
+ meshRefiner_.timeName()
|
|
+ ".obj"
|
|
);
|
|
Info<< "Points to be merged to " << str.name() << endl;
|
|
forAll(pointToMaster, pointi)
|
|
{
|
|
label otherPointi = pointToMaster[pointi];
|
|
if (otherPointi != -1)
|
|
{
|
|
const point& pt = mesh.points()[pointi];
|
|
const point& otherPt = mesh.points()[otherPointi];
|
|
str.write(linePointRef(pt, otherPt));
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
autoPtr<mapPolyMesh> map = meshRefiner_.mergePoints(pointToMaster);
|
|
if (map)
|
|
{
|
|
inplaceReorder(map().reverseCellMap(), cellNLayers);
|
|
|
|
const labelList& reverseFaceMap = map().reverseFaceMap();
|
|
inplaceReorder(reverseFaceMap, faceWantedThickness);
|
|
inplaceReorder(reverseFaceMap, faceRealThickness);
|
|
|
|
Info<< "Merged points in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
}
|
|
}
|
|
|
|
if (mesh.faceZones().size() > 0)
|
|
{
|
|
// Merge any baffles
|
|
Info<< "Converting baffles back into zoned faces ..."
|
|
<< endl;
|
|
|
|
autoPtr<mapPolyMesh> map = meshRefiner_.mergeZoneBaffles
|
|
(
|
|
true, // internal zones
|
|
false // baffle zones
|
|
);
|
|
if (map)
|
|
{
|
|
inplaceReorder(map().reverseCellMap(), cellNLayers);
|
|
|
|
const labelList& faceMap = map().faceMap();
|
|
|
|
// Make sure to keep the max since on two patches only one has
|
|
// layers.
|
|
scalarField newFaceRealThickness(mesh.nFaces(), Zero);
|
|
scalarField newFaceWantedThickness(mesh.nFaces(), Zero);
|
|
forAll(newFaceRealThickness, facei)
|
|
{
|
|
label oldFacei = faceMap[facei];
|
|
if (oldFacei >= 0)
|
|
{
|
|
scalar& realThick = newFaceRealThickness[facei];
|
|
realThick = max(realThick, faceRealThickness[oldFacei]);
|
|
scalar& wanted = newFaceWantedThickness[facei];
|
|
wanted = max(wanted, faceWantedThickness[oldFacei]);
|
|
}
|
|
}
|
|
faceRealThickness.transfer(newFaceRealThickness);
|
|
faceWantedThickness.transfer(newFaceWantedThickness);
|
|
}
|
|
|
|
Info<< "Converted baffles in = "
|
|
<< meshRefiner_.mesh().time().cpuTimeIncrement()
|
|
<< " s\n" << nl << endl;
|
|
}
|
|
|
|
// Do final balancing
|
|
// ~~~~~~~~~~~~~~~~~~
|
|
|
|
if (Pstream::parRun())
|
|
{
|
|
Info<< nl
|
|
<< "Doing final balancing" << nl
|
|
<< "---------------------" << nl
|
|
<< endl;
|
|
|
|
if (debug)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
}
|
|
|
|
// Balance. No restriction on face zones and baffles.
|
|
autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
|
|
(
|
|
false,
|
|
false,
|
|
scalarField(mesh.nCells(), 1.0),
|
|
decomposer,
|
|
distributor
|
|
);
|
|
|
|
// Re-distribute flag of layer faces and cells
|
|
map().distributeCellData(cellNLayers);
|
|
map().distributeFaceData(faceWantedThickness);
|
|
map().distributeFaceData(faceRealThickness);
|
|
}
|
|
|
|
|
|
// Write mesh data
|
|
// ~~~~~~~~~~~~~~~
|
|
|
|
if (!dryRun_)
|
|
{
|
|
writeLayerData
|
|
(
|
|
mesh,
|
|
patchIDs,
|
|
cellNLayers,
|
|
faceWantedThickness,
|
|
faceRealThickness
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::snappyLayerDriver::doLayers
|
|
(
|
|
const dictionary& shrinkDict,
|
|
const dictionary& motionDict,
|
|
const layerParameters& layerParams,
|
|
const meshRefinement::FaceMergeType mergeType,
|
|
const bool preBalance,
|
|
decompositionMethod& decomposer,
|
|
fvMeshDistribute& distributor
|
|
)
|
|
{
|
|
addProfiling(layers, "snappyHexMesh::layers");
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
Info<< nl
|
|
<< "Shrinking and layer addition phase" << nl
|
|
<< "----------------------------------" << nl
|
|
<< endl;
|
|
|
|
|
|
Info<< "Using mesh parameters " << motionDict << nl << endl;
|
|
|
|
// Merge coplanar boundary faces
|
|
if
|
|
(
|
|
mergeType == meshRefinement::FaceMergeType::GEOMETRIC
|
|
|| mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
|
|
)
|
|
{
|
|
mergePatchFacesUndo(layerParams, motionDict, mergeType);
|
|
}
|
|
|
|
|
|
// Per patch the number of layers (-1 or 0 if no layer)
|
|
const labelList& numLayers = layerParams.numLayers();
|
|
|
|
// Patches that need to get a layer
|
|
DynamicList<label> patchIDs(numLayers.size());
|
|
label nFacesWithLayers = 0;
|
|
forAll(numLayers, patchi)
|
|
{
|
|
if (numLayers[patchi] > 0)
|
|
{
|
|
const polyPatch& pp = mesh.boundaryMesh()[patchi];
|
|
|
|
if (!pp.coupled())
|
|
{
|
|
patchIDs.append(patchi);
|
|
nFacesWithLayers += mesh.boundaryMesh()[patchi].size();
|
|
}
|
|
else
|
|
{
|
|
WarningInFunction
|
|
<< "Ignoring layers on coupled patch " << pp.name()
|
|
<< endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Add contributions from faceZones that get layers
|
|
const faceZoneMesh& fZones = mesh.faceZones();
|
|
forAll(fZones, zonei)
|
|
{
|
|
label mpi, spi;
|
|
surfaceZonesInfo::faceZoneType fzType;
|
|
meshRefiner_.getFaceZoneInfo(fZones[zonei].name(), mpi, spi, fzType);
|
|
|
|
if (numLayers[mpi] > 0)
|
|
{
|
|
nFacesWithLayers += fZones[zonei].size();
|
|
}
|
|
if (numLayers[spi] > 0)
|
|
{
|
|
nFacesWithLayers += fZones[zonei].size();
|
|
}
|
|
}
|
|
|
|
|
|
patchIDs.shrink();
|
|
|
|
if (returnReduce(nFacesWithLayers, sumOp<label>()) == 0)
|
|
{
|
|
Info<< nl << "No layers to generate ..." << endl;
|
|
}
|
|
else
|
|
{
|
|
// Check that outside of mesh is not multiply connected.
|
|
checkMeshManifold();
|
|
|
|
// Check initial mesh
|
|
Info<< "Checking initial mesh ..." << endl;
|
|
labelHashSet wrongFaces(mesh.nFaces()/100);
|
|
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun_);
|
|
const label nInitErrors = returnReduce
|
|
(
|
|
wrongFaces.size(),
|
|
sumOp<label>()
|
|
);
|
|
|
|
Info<< "Detected " << nInitErrors << " illegal faces"
|
|
<< " (concave, zero area or negative cell pyramid volume)"
|
|
<< endl;
|
|
|
|
|
|
bool faceZoneOnCoupledFace = false;
|
|
|
|
if (!preBalance)
|
|
{
|
|
// Check if there are faceZones on processor boundaries. This
|
|
// requires balancing to move them off the processor boundaries.
|
|
|
|
// Is face on a faceZone
|
|
bitSet isExtrudedZoneFace(mesh.nFaces());
|
|
{
|
|
// Add contributions from faceZones that get layers
|
|
const faceZoneMesh& fZones = mesh.faceZones();
|
|
forAll(fZones, zonei)
|
|
{
|
|
const faceZone& fZone = fZones[zonei];
|
|
const word& fzName = fZone.name();
|
|
|
|
label mpi, spi;
|
|
surfaceZonesInfo::faceZoneType fzType;
|
|
meshRefiner_.getFaceZoneInfo(fzName, mpi, spi, fzType);
|
|
|
|
if (numLayers[mpi] > 0 || numLayers[spi])
|
|
{
|
|
isExtrudedZoneFace.set(fZone);
|
|
}
|
|
}
|
|
}
|
|
|
|
bitSet intOrCoupled
|
|
(
|
|
syncTools::getInternalOrCoupledFaces(mesh)
|
|
);
|
|
|
|
for
|
|
(
|
|
label facei = mesh.nInternalFaces();
|
|
facei < mesh.nFaces();
|
|
facei++
|
|
)
|
|
{
|
|
if (intOrCoupled[facei] && isExtrudedZoneFace.test(facei))
|
|
{
|
|
faceZoneOnCoupledFace = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
reduce(faceZoneOnCoupledFace, orOp<bool>());
|
|
}
|
|
|
|
|
|
|
|
|
|
// Balance
|
|
if (Pstream::parRun() && (preBalance || faceZoneOnCoupledFace))
|
|
{
|
|
Info<< nl
|
|
<< "Doing initial balancing" << nl
|
|
<< "-----------------------" << nl
|
|
<< endl;
|
|
|
|
scalarField cellWeights(mesh.nCells(), 1);
|
|
forAll(numLayers, patchi)
|
|
{
|
|
if (numLayers[patchi] > 0)
|
|
{
|
|
const polyPatch& pp = mesh.boundaryMesh()[patchi];
|
|
forAll(pp.faceCells(), i)
|
|
{
|
|
cellWeights[pp.faceCells()[i]] += numLayers[patchi];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Add contributions from faceZones that get layers
|
|
const faceZoneMesh& fZones = mesh.faceZones();
|
|
forAll(fZones, zonei)
|
|
{
|
|
const faceZone& fZone = fZones[zonei];
|
|
const word& fzName = fZone.name();
|
|
|
|
label mpi, spi;
|
|
surfaceZonesInfo::faceZoneType fzType;
|
|
meshRefiner_.getFaceZoneInfo(fzName, mpi, spi, fzType);
|
|
|
|
if (numLayers[mpi] > 0)
|
|
{
|
|
// Get the owner side for unflipped faces, neighbour side
|
|
// for flipped ones
|
|
const labelList& cellIDs = fZone.slaveCells();
|
|
forAll(cellIDs, i)
|
|
{
|
|
if (cellIDs[i] >= 0)
|
|
{
|
|
cellWeights[cellIDs[i]] += numLayers[mpi];
|
|
}
|
|
}
|
|
}
|
|
if (numLayers[spi] > 0)
|
|
{
|
|
const labelList& cellIDs = fZone.masterCells();
|
|
forAll(cellIDs, i)
|
|
{
|
|
if (cellIDs[i] >= 0)
|
|
{
|
|
cellWeights[cellIDs[i]] += numLayers[mpi];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Balance mesh (and meshRefinement). Restrict faceZones to
|
|
// be on internal faces only since they will be converted into
|
|
// baffles.
|
|
autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
|
|
(
|
|
true, // keepZoneFaces
|
|
false,
|
|
cellWeights,
|
|
decomposer,
|
|
distributor
|
|
);
|
|
}
|
|
|
|
|
|
// Do all topo changes
|
|
addLayers
|
|
(
|
|
layerParams,
|
|
motionDict,
|
|
patchIDs,
|
|
nInitErrors,
|
|
decomposer,
|
|
distributor
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|