Skip to content
GitLab
Commit 3cf4d721 authored by henry's avatar henry
Browse files

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

parents 0dd8a2b8 d947aa9c
Hide whitespace changes
Inline Side-by-side
applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
View file @ 3cf4d721
......@@ -291,6 +291,22 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
if (allGeometry)
{
cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
if (mesh.checkConcaveCells(true, &cells))
{
noFailedChecks++;
label nCells = returnReduce(cells.size(), sumOp<label>());
Info<< " <<Writing " << nCells
<< " concave cells to set " << cells.name() << endl;
cells.instance() = mesh.pointsInstance();
cells.write();
}
}
return noFailedChecks;
}
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
View file @ 3cf4d721
......@@ -291,6 +291,9 @@ class primitiveMesh
//- Skewness warning threshold
static scalar skewThreshold_;
//- Threshold where faces are considered coplanar
static scalar planarCosAngle_;
protected:
......@@ -646,6 +649,13 @@ public:
labelHashSet* setPtr = NULL
) const;
//- Check for concave cells
bool checkConcaveCells
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check mesh topology for correctness.
// Returns false for no error.
......
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
View file @ 3cf4d721
......@@ -37,6 +37,7 @@ Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
......@@ -2093,6 +2094,129 @@ bool Foam::primitiveMesh::checkCellDeterminant
}
bool Foam::primitiveMesh::checkConcaveCells
(
const bool report,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool primitiveMesh::checkConcaveCells(const bool"
<< ", labelHashSet*) const: "
<< "checking for concave cells" << endl;
}
const cellList& c = cells();
const labelList& fOwner = faceOwner();
const vectorField& fAreas = faceAreas();
const pointField& fCentres = faceCentres();
const pointField& pts = points();
label nConcaveCells = 0;
forAll (c, cellI)
{
const cell& cFaces = c[cellI];
const labelList& cPoints = cellPoints()[cellI];
bool concave = false;
forAll(cFaces, i)
{
if (concave)
{
break;
}
label fI = cFaces[i];
const point& fC = fCentres[fI];
const face& f = faces()[fI];
vector fN = fAreas[fI];
fN /= max(mag(fN), VSMALL);
// Flip normal if required so that it is always pointing out of
// the cell
if (fOwner[fI] != cellI)
{
fN *= -1;
}
// Is any vertex of the cell on the wrong side of the
// plane of this face?
forAll(cPoints, cPtI)
{
label ptI = cPoints[cPtI];
// Skip points that are on this face
if (findIndex(f, ptI) > -1)
{
continue;
}
const point& pt = pts[ptI];
// If the cell is concave, the point will be on the
// positive normal side of the plane of f, defined by
// its centre and normal, and the angle between (pt -
// fC) and fN will be less than 90 degrees, so the dot
// product will be positive.
vector pC = (pt - fC);
pC /= max(mag(pC), VSMALL);
if ((pC & fN) > -planarCosAngle_)
{
// Concave or planar face
concave = true;
if (setPtr)
{
setPtr->insert(cellI);
}
nConcaveCells++;
break;
}
}
}
}
reduce(nConcaveCells, sumOp<label>());
if (nConcaveCells > 0)
{
if (debug || report)
{
Info<< " ***Concave cells found, number of cells: "
<< nConcaveCells << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Concave cell check OK." << endl;
}
return false;
}
return false;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::primitiveMesh::checkTopology(const bool report) const
......
src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
View file @ 3cf4d721
......@@ -202,12 +202,45 @@ Foam::label Foam::removePoints::countPointUsage
pointCanBeDeleted[pointI] = true;
nDeleted++;
}
}
edge0.clear();
edge1.clear();
// Protect any points on faces that would collapse down to nothing
// No particular intelligence so might protect too many points
forAll(mesh_.faces(), faceI)
{
const face& f = mesh_.faces()[faceI];
label nCollapse = 0;
forAll(f, fp)
{
if (pointCanBeDeleted[f[fp]])
{
nCollapse++;
}
}
if ((f.size() - nCollapse) < 3)
{
// Just unmark enough points
forAll(f, fp)
{
if (pointCanBeDeleted[f[fp]])
{
pointCanBeDeleted[f[fp]] = false;
--nCollapse;
if (nCollapse == 0)
{
break;
}
}
}
}
}
// Point can be deleted only if all processors want to delete it
syncTools::syncPointList
(
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View file @ 3cf4d721
......@@ -1375,7 +1375,7 @@ void Foam::autoLayerDriver::growNoExtrusion
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
)
) const
{
Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
......@@ -1406,7 +1406,7 @@ void Foam::autoLayerDriver::growNoExtrusion
{
if
(
extrudeStatus[f[fp]] == NOEXTRUDE
extrudeStatus[f[fp]] == EXTRUDE
&& grownExtrudeStatus[f[fp]] != NOEXTRUDE
)
{
......@@ -1419,6 +1419,31 @@ void Foam::autoLayerDriver::growNoExtrusion
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
false // no separation
);
forAll(status, i)
{
extrudeStatus[i] = extrudeMode(status[i]);
}
}
forAll(extrudeStatus, patchPointI)
{
if (extrudeStatus[patchPointI] == NOEXTRUDE)
......@@ -2711,8 +2736,6 @@ void Foam::autoLayerDriver::addLayers
const labelList& meshPoints = patches[patchI].meshPoints();
//scalar maxThickness = -VGREAT;
//scalar minThickness = VGREAT;
scalar sumThickness = 0;
scalar sumNearWallThickness = 0;
......@@ -2720,8 +2743,6 @@ void Foam::autoLayerDriver::addLayers
{
label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
//maxThickness = max(maxThickness, thickness[ppPointI]);
//minThickness = min(minThickness, thickness[ppPointI]);
sumThickness += thickness[ppPointI];
label nLay = patchNLayers[ppPointI];
......@@ -2749,8 +2770,6 @@ void Foam::autoLayerDriver::addLayers
if (totNPoints > 0)
{
//reduce(maxThickness, maxOp<scalar>());
//reduce(minThickness, minOp<scalar>());
avgThickness =
returnReduce(sumThickness, sumOp<scalar>())
/ totNPoints;
......@@ -2766,8 +2785,6 @@ void Foam::autoLayerDriver::addLayers
<< " " << setw(6) << layerParams.numLayers()[patchI]
<< " " << setw(8) << avgNearWallThickness
<< " " << setw(8) << avgThickness
//<< " " << setw(8) << minThickness
//<< " " << setw(8) << maxThickness
<< endl;
}
Info<< endl;
......@@ -3147,6 +3164,19 @@ void Foam::autoLayerDriver::addLayers
meshMover().movePoints(oldPoints);
meshMover().correct();
// Grow out region of non-extrusion
for (label i = 0; i < layerParams.nGrow(); i++)
{
growNoExtrusion
(
pp,
patchDisp,
patchNLayers,
extrudeStatus
);
}
Info<< endl;
}
......@@ -3293,7 +3323,6 @@ void Foam::autoLayerDriver::doLayers
// Balance
if (Pstream::parRun())
if (Pstream::parRun() && preBalance)
{
Info<< nl
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View file @ 3cf4d721
......@@ -232,13 +232,13 @@ class autoLayerDriver
) const;
//- Grow no-extrusion layer.
static void growNoExtrusion
void growNoExtrusion
(
const indirectPrimitivePatch& pp,
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
);
) const;
//- Calculate pointwise wanted and minimum thickness.
// thickness: wanted thickness
......
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View file @ 3cf4d721
......@@ -245,27 +245,8 @@ void Foam::meshRefinement::getBafflePatches
const pointField& cellCentres = mesh_.cellCentres();
// Build list of surfaces that are not to be baffled.
const wordList& cellZoneNames = surfaces_.cellZoneNames();
labelList surfacesToBaffle(cellZoneNames.size());
label baffleI = 0;
forAll(cellZoneNames, surfI)
{
if (cellZoneNames[surfI].size())
{
if (debug)
{
Pout<< "getBafflePatches : Not baffling surface "
<< surfaces_.names()[surfI] << endl;
}
}
else
{
surfacesToBaffle[baffleI++] = surfI;
}
}
surfacesToBaffle.setSize(baffleI);
// Surfaces that need to be baffled
const labelList surfacesToBaffle(surfaces_.getUnnamedSurfaces());
ownPatch.setSize(mesh_.nFaces());
ownPatch = -1;
......@@ -1253,7 +1234,7 @@ void Foam::meshRefinement::findCellZoneTopo
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
if (surfI != -1 && surfaceToCellZone[surfI] != -1)
{
// Calculate region to zone from cellRegions on either side
// of internal face.
......@@ -1305,7 +1286,7 @@ void Foam::meshRefinement::findCellZoneTopo
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
if (surfI != -1 && surfaceToCellZone[surfI] != -1)
{
bool changedCell = calcRegionToZone
(
......@@ -1439,6 +1420,15 @@ void Foam::meshRefinement::makeConsistentFaceIndex
}
}
}
else
{
// Unzonify boundary faces
forAll(pp, i)
{
label faceI = pp.start()+i;
namedSurfaceIndex[faceI] = -1;
}
}
}
}
......@@ -2042,14 +2032,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
labelList namedSurfaces(surfaces_.getNamedSurfaces());
boolList isNamedSurface(cellZoneNames.size(), false);
forAll(namedSurfaces, i)
{
label surfI = namedSurfaces[i];
isNamedSurface[surfI] = true;
Info<< "Surface : " << surfaces_.names()[surfI] << nl
<< " faceZone : " << faceZoneNames[surfI] << nl
<< " cellZone : " << cellZoneNames[surfI] << endl;
......@@ -2125,31 +2111,34 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
{
label surfI = namedSurfaces[i];
label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
if (zoneI == -1)
if (cellZoneNames[surfI] != word::null)
{
zoneI = cellZones.size();
cellZones.setSize(zoneI+1);
cellZones.set
(
zoneI,
new cellZone
label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
if (zoneI == -1)
{
zoneI = cellZones.size();
cellZones.setSize(zoneI+1);
cellZones.set
(
cellZoneNames[surfI], //name
labelList(0), //addressing
zoneI, //index
cellZones //cellZoneMesh
)
);
}
zoneI,
new cellZone
(
cellZoneNames[surfI], //name
labelList(0), //addressing
zoneI, //index
cellZones //cellZoneMesh
)
);
}
if (debug)
{
Pout<< "Cells inside " << surfaces_.names()[surfI]
<< " will go into cellZone " << zoneI << endl;
if (debug)
{
Pout<< "Cells inside " << surfaces_.names()[surfI]
<< " will go into cellZone " << zoneI << endl;
}
surfaceToCellZone[surfI] = zoneI;
}
surfaceToCellZone[surfI] = zoneI;
}
// Check they are synced
......@@ -2321,6 +2310,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (closedNamedSurfaces.size())
{
Info<< "Found " << closedNamedSurfaces.size()
<< " closed, named surfaces. Assigning cells in/outside"
<< " these surfaces to the corresponding cellZone."
<< nl << endl;
findCellZoneGeometric
(
closedNamedSurfaces, // indices of closed surfaces
......@@ -2333,8 +2327,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Set using walking
// ~~~~~~~~~~~~~~~~~
//if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
if (!allowFreeStandingZoneFaces)
{
Info<< "Walking to assign cellZones." << nl << endl;
// Topological walk
findCellZoneTopo
(
......@@ -2349,6 +2345,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
if (!allowFreeStandingZoneFaces)
{
Info<< "Only keeping zone faces inbetween different cellZones."
<< nl << endl;
makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
}
......
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View file @ 3cf4d721
......@@ -77,8 +77,18 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone"))
{
dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI];
dict.lookup("zoneInside") >> zoneInside_[surfI];
if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
{
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
else if (dict.found("zoneInside"))
{
IOWarningIn("refinementSurfaces::refinementSurfaces(..)", dict)
<< "Unused entry zoneInside for faceZone "
<< faceZoneNames_[surfI]
<< " since no cellZone specified."
<< endl;
}
}
// Global perpendicular angle
......@@ -314,8 +324,21 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone"))
{
dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI];
dict.lookup("zoneInside") >> zoneInside_[surfI];
if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
{
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
else if (dict.found("zoneInside"))
{
IOWarningIn
(
"refinementSurfaces::refinementSurfaces(..)",
dict
) << "Unused entry zoneInside for faceZone "
<< faceZoneNames_[surfI]
<< " since no cellZone specified."
<< endl;
}
}
// Global perpendicular angle
......@@ -475,18 +498,17 @@ Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
// Get indices of closed named surfaces
Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
{
labelList named(getNamedSurfaces());
labelList closed(cellZoneNames_.size());
labelList closed(named.size());
label closedI = 0;
forAll(named, i)
forAll(cellZoneNames_, surfI)
{
label surfI = named[i];
if (allGeometry_[surfaces_[surfI]].hasVolumeType())
if (cellZoneNames_[surfI].size())
{
closed[closedI++] = surfI;
if (allGeometry_[surfaces_[surfI]].hasVolumeType())
{
closed[closedI++] = surfI;
}
}
}