zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature/foamyHexMesh-zoning'

Browse Source
This commit is contained in:
laurence
2013-09-04 12:28:16 +01:00
parent b8bb781e6d 040e6d9159
commit 0b0e0e69f4
10 changed files with 599 additions and 70 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/utilities/mesh/generation/foamyHexMesh/Make/options
View File

@ -21,6 +21,7 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/mesh/autoMesh/lnInclude \
-IvectorTools
EXE_LIBS = \

8
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
View File

@ -279,11 +279,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
const Triangulation& t
)
{
tmp<pointField> tpts(new pointField(t.number_of_vertices(), point::max));
tmp<pointField> tpts(new pointField(t.vertexCount(), point::max));
pointField& pts = tpts();
label nVert = 0;
for
(
typename Triangulation::Finite_vertices_iterator vit =
@ -292,9 +290,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
++vit
)
{
if (vit->internalOrBoundaryPoint())
if (vit->internalOrBoundaryPoint() && !vit->referred())
{
pts[nVert++] = topoint(vit->point());
pts[vit->index()] = topoint(vit->point());
}
}

4
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
View File

@ -22,6 +22,7 @@ EXE_INC = \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/mesh/autoMesh/lnInclude \
-IPrintTable \
-I../vectorTools
@ -32,4 +33,5 @@ LIB_LIBS = \
-ltriSurface \
-ldynamicMesh \
-lsurfMesh \
-lsampling
-lsampling \
-lautoMesh

26
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -600,6 +600,29 @@ private:
PackedBoolList& boundaryFacesToRemove
);
//- From meshRefinementBaffles.C. Use insidePoint for a surface to
// determine the cell zone.
void findCellZoneInsideWalk
(
const polyMesh& mesh,
const labelList& locationSurfaces,
const labelList& faceToSurface,
labelList& cellToSurface
) const;
//- Calculate the cell zones from cellCentres using all closed surfaces
labelList calcCellZones(const pointField& cellCentres) const;
//- Calculate the face zones
void calcFaceZones
(
const polyMesh& mesh,
const pointField& cellCentres,
const labelList& cellToSurface,
labelList& faceToSurface,
boolList& flipMap
) const;
//- Tet mesh calculation
void calcTetMesh
(
@ -712,9 +735,6 @@ private:
bool includeEmptyPatches = false
) const;
//- Create the cell centres to use for the mesh
void createCellCentres(pointField& cellCentres) const;
//- Sort the faces, owner and neighbour lists into
// upper-triangular order. For internal faces only, use
// before adding patch faces

28
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
View File

@ -2852,32 +2852,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
}
void Foam::conformalVoronoiMesh::createCellCentres
(
pointField& cellCentres
) const
{
cellCentres.setSize(number_of_vertices(), point::max);
label vertI = 0;
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (vit->internalOrBoundaryPoint())
{
cellCentres[vertI++] = topoint(vit->point());
}
}
cellCentres.setSize(vertI);
}
void Foam::conformalVoronoiMesh::sortFaces
(
faceList& faces,
@ -3095,7 +3069,7 @@ Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
{
Info<< nl << "Removing unused cells" << endl;
PackedBoolList cellUsed(number_of_vertices(), false);
PackedBoolList cellUsed(vertexCount(), false);
// Scan all faces to find all of the cells that are used

524
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -35,6 +35,11 @@ License
#include "pointMesh.H"
#include "indexedVertexOps.H"
#include "DelaunayMeshTools.H"
#include "surfaceZonesInfo.H"
#include "polyModifyCell.H"
#include "polyModifyFace.H"
#include "syncTools.H"
#include "regionSplit.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -402,6 +407,376 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
}
void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
(
const polyMesh& mesh,
const labelList& locationSurfaces, // indices of surfaces with inside point
const labelList& faceToSurface, // per face index of named surface
labelList& cellToSurface
) const
{
// Analyse regions. Reuse regionsplit
boolList blockedFace(mesh.nFaces());
//selectSeparatedCoupledFaces(blockedFace);
forAll(faceToSurface, faceI)
{
if (faceToSurface[faceI] == -1)
{
blockedFace[faceI] = false;
}
else
{
blockedFace[faceI] = true;
}
}
// No need to sync since namedSurfaceIndex already is synced
// Set region per cell based on walking
regionSplit cellRegion(mesh, blockedFace);
blockedFace.clear();
// Force calculation of face decomposition (used in findCell)
(void)mesh.tetBasePtIs();
const PtrList<surfaceZonesInfo>& surfZones =
geometryToConformTo().surfZones();
// For all locationSurface find the cell
forAll(locationSurfaces, i)
{
label surfI = locationSurfaces[i];
const Foam::point& insidePoint = surfZones[surfI].zoneInsidePoint();
const word& surfName = geometryToConformTo().geometry().names()[surfI];
Info<< " For surface " << surfName
<< " finding inside point " << insidePoint
<< endl;
// Find the region containing the insidePoint
label keepRegionI = -1;
label cellI = mesh.findCell(insidePoint);
if (cellI != -1)
{
keepRegionI = cellRegion[cellI];
}
reduce(keepRegionI, maxOp<label>());
Info<< " For surface " << surfName
<< " found point " << insidePoint << " in cell " << cellI
<< " in global region " << keepRegionI
<< " out of " << cellRegion.nRegions() << " regions." << endl;
if (keepRegionI == -1)
{
FatalErrorIn
(
"conformalVoronoiMesh::findCellZoneInsideWalk"
"(const polyMesh&, const labelList&"
", const labelList&, labelList&)"
) << "Point " << insidePoint
<< " is not inside the mesh." << nl
<< "Bounding box of the mesh:" << mesh.bounds()
<< exit(FatalError);
}
// Set all cells with this region
forAll(cellRegion, cellI)
{
if (cellRegion[cellI] == keepRegionI)
{
if (cellToSurface[cellI] == -2)
{
cellToSurface[cellI] = surfI;
}
else if (cellToSurface[cellI] != surfI)
{
WarningIn
(
"conformalVoronoiMesh::findCellZoneInsideWalk"
"(const labelList&, const labelList&"
", const labelList&, const labelList&)"
) << "Cell " << cellI
<< " at " << mesh.cellCentres()[cellI]
<< " is inside surface " << surfName
<< " but already marked as being in zone "
<< cellToSurface[cellI] << endl
<< "This can happen if your surfaces are not"
<< " (sufficiently) closed."
<< endl;
}
}
}
}
}
Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
(
const pointField& cellCentres
) const
{
labelList cellToSurface(cellCentres.size(), -1);
const PtrList<surfaceZonesInfo>& surfZones =
geometryToConformTo().surfZones();
// Get list of closed surfaces
labelList closedNamedSurfaces
(
surfaceZonesInfo::getAllClosedNamedSurfaces
(
surfZones,
geometryToConformTo().geometry(),
geometryToConformTo().surfaces()
)
);
forAll(closedNamedSurfaces, i)
{
label surfI = closedNamedSurfaces[i];
const searchableSurface& surface =
allGeometry()[geometryToConformTo().surfaces()[surfI]];
const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
surfZones[surfI].zoneInside();
if
(
selectionMethod != surfaceZonesInfo::INSIDE
&& selectionMethod != surfaceZonesInfo::OUTSIDE
&& selectionMethod != surfaceZonesInfo::INSIDEPOINT
)
{
FatalErrorIn("conformalVoronoiMesh::calcCellZones(..)")
<< "Trying to use surface "
<< surface.name()
<< " which has non-geometric inside selection method "
<< surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
<< exit(FatalError);
}
if (surface.hasVolumeType())
{
List<volumeType> volType;
surface.getVolumeType(cellCentres, volType);
bool selectInside = true;
if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
{
List<volumeType> volTypeInsidePoint;
surface.getVolumeType
(
pointField(1, surfZones[surfI].zoneInsidePoint()),
volTypeInsidePoint
);
if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
{
selectInside = false;
}
}
else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
{
selectInside = false;
}
forAll(volType, pointI)
{
if (cellToSurface[pointI] == -1)
{
if
(
(
volType[pointI] == volumeType::INSIDE
&& selectInside
)
|| (
volType[pointI] == volumeType::OUTSIDE
&& !selectInside
)
)
{
cellToSurface[pointI] = surfI;
}
}
}
}
}
return cellToSurface;
}
void Foam::conformalVoronoiMesh::calcFaceZones
(
const polyMesh& mesh,
const pointField& cellCentres,
const labelList& cellToSurface,
labelList& faceToSurface,
boolList& flipMap
) const
{
faceToSurface.setSize(mesh.nFaces(), -1);
flipMap.setSize(mesh.nFaces(), false);
const faceList& faces = mesh.faces();
const labelList& faceOwner = mesh.faceOwner();
const labelList& faceNeighbour = mesh.faceNeighbour();
forAll(faces, faceI)
{
const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
if (mesh.isInternalFace(faceI))
{
const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
flipMap[faceI] =
(
ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
? false
: true
);
if
(
(ownerSurfaceI >= 0 || neiSurfaceI >= 0)
&& ownerSurfaceI != neiSurfaceI
)
{
if (ownerSurfaceI > neiSurfaceI)
{
faceToSurface[faceI] = ownerSurfaceI;
}
else
{
faceToSurface[faceI] = neiSurfaceI;
}
}
}
else
{
if (ownerSurfaceI >= 0)
{
faceToSurface[faceI] = ownerSurfaceI;
}
}
}
const PtrList<surfaceZonesInfo>& surfZones =
geometryToConformTo().surfZones();
labelList insidePointNamedSurfaces
(
surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
);
// Use intersection of cellCentre connections
forAll(faces, faceI)
{
if
(
mesh.isInternalFace(faceI)
&& faceToSurface[faceI] < 0
)
{
const label own = faceOwner[faceI];
const label nei = faceNeighbour[faceI];
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo().findSurfaceAnyIntersection
(
cellCentres[own],
cellCentres[nei],
surfHit,
hitSurface
);
if (surfHit.hit())
{
if (findIndex(insidePointNamedSurfaces, hitSurface) != -1)
{
faceToSurface[faceI] = hitSurface;
vectorField norm;
geometryToConformTo().getNormal
(
hitSurface,
List<pointIndexHit>(1, surfHit),
norm
);
const vector fN = faces[faceI].normal(mesh.points());
if ((norm[0] & fN/(mag(fN) + SMALL)) < 0)
{
flipMap[faceI] = true;
}
else
{
flipMap[faceI] = false;
}
}
}
}
}
labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label ownSurface = cellToSurface[faceOwner[faceI]];
neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
}
}
}
syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label ownSurface = cellToSurface[faceOwner[faceI]];
label neiSurface = neiCellSurface[faceI-mesh.nInternalFaces()];
if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
{
// Give face the max cell zone
faceToSurface[faceI] = max(ownSurface, neiSurface);
}
}
}
}
// Sync
syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
}
Foam::autoPtr<Foam::fvMesh> Foam::conformalVoronoiMesh::createDummyMesh
(
const IOobject& io,
@ -883,7 +1258,11 @@ void Foam::conformalVoronoiMesh::writeMesh
// Check that the patch is not empty on every processor
reduce(totalPatchSize, sumOp<label>());
if (totalPatchSize > 0)
if
(
totalPatchSize > 0
// && !geometryToConformTo().surfZones().set(p)
)
{
patches[nValidPatches] = polyPatch::New
(
@ -898,6 +1277,145 @@ void Foam::conformalVoronoiMesh::writeMesh
}
}
patches.setSize(nValidPatches);
mesh.addFvPatches(patches);
// Add zones to mesh
{
Info<< " Adding zones to mesh" << endl;
const PtrList<surfaceZonesInfo>& surfZones =
geometryToConformTo().surfZones();
labelList cellToSurface(calcCellZones(cellCentres));
labelList faceToSurface;
boolList flipMap;
calcFaceZones
(
mesh,
cellCentres,
cellToSurface,
faceToSurface,
flipMap
);
labelList insidePointNamedSurfaces
(
surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
);
findCellZoneInsideWalk
(
mesh,
insidePointNamedSurfaces,
faceToSurface,
cellToSurface
);
labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
forAll(namedSurfaces, i)
{
label surfI = namedSurfaces[i];
Info<< incrIndent << indent << "Surface : "
<< geometryToConformTo().geometry().names()[surfI] << nl
<< indent << " faceZone : "
<< surfZones[surfI].faceZoneName() << nl
<< indent << " cellZone : "
<< surfZones[surfI].cellZoneName()
<< decrIndent << endl;
}
// Add zones to mesh
labelList surfaceToFaceZone =
surfaceZonesInfo::addFaceZonesToMesh
(
surfZones,
namedSurfaces,
mesh
);
labelList surfaceToCellZone =
surfaceZonesInfo::addCellZonesToMesh
(
surfZones,
namedSurfaces,
mesh
);
// Topochange container
polyTopoChange meshMod(mesh);
forAll(cellToSurface, cellI)
{
label surfaceI = cellToSurface[cellI];
if (surfaceI >= 0)
{
label zoneI = surfaceToCellZone[surfaceI];
if (zoneI >= 0)
{
meshMod.setAction
(
polyModifyCell
(
cellI,
false, // removeFromZone
zoneI
)
);
}
}
}
const labelList& faceOwner = mesh.faceOwner();
const labelList& faceNeighbour = mesh.faceNeighbour();
forAll(faceToSurface, faceI)
{
if (!mesh.isInternalFace(faceI))
{
continue;
}
label surfaceI = faceToSurface[faceI];
if (surfaceI >= 0)
{
label own = faceOwner[faceI];
label nei = faceNeighbour[faceI];
meshMod.setAction
(
polyModifyFace
(
mesh.faces()[faceI], // modified face
faceI, // label of face
own, // owner
nei, // neighbour
false, // face flip
-1, // patch for face
false, // remove from zone
surfaceToFaceZone[surfaceI], // zone for face
flipMap[faceI] // face flip in zone
)
);
}
}
// Change the mesh (no inflation, parallel sync)
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
}
// Add indirectPatchFaces to a face zone
{
labelList addr(boundaryFacesToRemove.count());
@ -928,10 +1446,6 @@ void Foam::conformalVoronoiMesh::writeMesh
);
}
patches.setSize(nValidPatches);
mesh.addFvPatches(patches);
timeCheck("Before fvMesh filtering");
autoPtr<polyMeshFilter> meshFilter;

24
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
View File

@ -579,30 +579,6 @@ Foam::conformationSurfaces::conformationSurfaces
Info<< features_[fI].name() << endl;
}
}
Info<< "ZONES" << endl;
forAll(surfZones_, surfI)
{
if (surfZones_.set(surfI))
{
const surfaceZonesInfo& sInfo = surfZones_[surfI];
Info<< " " << surfI << nl
<< " faceZoneName = " << sInfo.faceZoneName() << nl
<< " cellZoneName = " << sInfo.cellZoneName() << nl
<< " zoneInside = "
<< surfaceZonesInfo::areaSelectionAlgoNames[sInfo.zoneInside()]
<< nl
<< " zoneInsidePoint = " << sInfo.zoneInsidePoint() << nl
<< " faceType = "
<< surfaceZonesInfo::faceZoneTypeNames[sInfo.faceType()]
<< endl;
}
else
{
Info<< " " << surfI << " EMPTY" << endl;
}
}
}

1
applications/utilities/mesh/generation/foamyQuadMesh/Make/options
View File

@ -23,6 +23,7 @@ EXE_INC = \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/mesh/autoMesh/lnInclude
EXE_LIBS = \
$(CGAL_LIBS) \

45
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
View File

@ -124,7 +124,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
{
IOWarningIn
(
"refinementSurfaces::refinementSurfaces(..)",
"surfaceZonesInfo::surfaceZonesInfo(..)",
surfacesDict
) << "Illegal entry zoneInside "
<< areaSelectionAlgoNames[zoneInside_]
@ -137,7 +137,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
{
IOWarningIn
(
"refinementSurfaces::refinementSurfaces(..)",
"surfaceZonesInfo::surfaceZonesInfo(..)",
surfacesDict
) << "Unused entry zoneInside for faceZone "
<< faceZoneName_
@ -215,7 +215,11 @@ Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
label namedI = 0;
forAll(surfList, surfI)
{
if (surfList[surfI].faceZoneName().size())
if
(
surfList.set(surfI)
&& surfList[surfI].faceZoneName().size()
)
{
namedSurfaces[namedI++] = surfI;
}
@ -241,7 +245,8 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
{
if
(
surfList[surfI].cellZoneName().size()
surfList.set(surfI)
&& surfList[surfI].cellZoneName().size()
&& (
surfList[surfI].zoneInside() == surfaceZonesInfo::INSIDE
|| surfList[surfI].zoneInside() == surfaceZonesInfo::OUTSIDE
@ -258,6 +263,35 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
}
// Get indices of closed named surfaces
Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
(
const PtrList<surfaceZonesInfo>& surfList,
const searchableSurfaces& allGeometry,
const labelList& surfaces
)
{
labelList closed(surfList.size());
label closedI = 0;
forAll(surfList, surfI)
{
if
(
surfList.set(surfI)
&& surfList[surfI].cellZoneName().size()
&& allGeometry[surfaces[surfI]].hasVolumeType()
)
{
closed[closedI++] = surfI;
}
}
closed.setSize(closedI);
return closed;
}
// Get indices of named surfaces with a
Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
(
@ -271,7 +305,8 @@ Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
{
if
(
surfList[surfI].cellZoneName().size()
surfList.set(surfI)
&& surfList[surfI].cellZoneName().size()
&& surfList[surfI].zoneInside() == surfaceZonesInfo::INSIDEPOINT
)
{

8
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
View File

@ -199,6 +199,14 @@ public:
const labelList& surfaces
);
//- Get indices of surfaces with a cellZone that are closed.
static labelList getAllClosedNamedSurfaces
(
const PtrList<surfaceZonesInfo>& surfList,
const searchableSurfaces& allGeometry,
const labelList& surfaces
);
//- Get indices of surfaces with a cellZone that have 'insidePoint'
// section.
static labelList getInsidePointNamedSurfaces