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

ENH: Refine cells with huge cellWeights. Store meshCutter in class to reuse

Browse Source 
later.
This commit is contained in:
graham
2011-06-03 16:18:29 +01:00
parent 1eb3cd5a52
commit fcfba3ade6
3 changed files with 346 additions and 135 deletions
Download Patch File Download Diff File Expand all files Collapse all files

245
src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "backgroundMeshDecomposition.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -58,13 +59,6 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
decompositionMethod& decomposer = decomposerPtr_();
hexRef8 meshCutter
(
mesh_,
labelList(mesh_.nCells(), 0),
labelList(mesh_.nPoints(), 0)
);
// For each cell in the mesh has it been determined if it is fully
// inside, outside, or overlaps the surface
labelList volumeStatus
@ -110,7 +104,6 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
{
labelList refCells = selectRefinementCells
(
meshCutter,
volumeStatus,
cellWeights
);
@ -118,7 +111,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
// Maintain 2:1 ratio
labelList newCellsToRefine
(
meshCutter.consistentRefinement
meshCutter_.consistentRefinement
(
refCells,
true // extend set
@ -141,12 +134,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
);
}
if
(
returnReduce(newCellsToRefine.size(), sumOp<label>())
== 0
)
if (returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
{
break;
}
@ -155,7 +143,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
polyTopoChange meshMod(mesh_);
// Play refinement commands into mesh changer.
meshCutter.setRefinement(newCellsToRefine, meshMod);
meshCutter_.setRefinement(newCellsToRefine, meshMod);
// Create mesh, return map from old to new mesh.
autoPtr<mapPolyMesh> map = meshMod.changeMesh
@ -163,7 +151,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
mesh_,
false, // inflate
true, // syncParallel
true, // orderCells (to reduce cell motion in scotch)
true, // orderCells (to reduce cell transfers)
false // orderPoints
);
@ -171,7 +159,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
mesh_.updateMesh(map);
// Update numbering of cells/vertices.
meshCutter.updateMesh(map);
meshCutter_.updateMesh(map);
{
// Map volumeStatus
@ -187,12 +175,11 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
if (oldCellI == -1)
{
newVolumeStatus[newCellI] =
searchableSurface::UNKNOWN;;
searchableSurface::UNKNOWN;
}
else
{
newVolumeStatus[newCellI] =
volumeStatus[oldCellI];
newVolumeStatus[newCellI] = volumeStatus[oldCellI];
}
}
@ -201,7 +188,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
if (debug)
{
Info<< "Refined from "
Info<< " Background mesh refined from "
<< returnReduce(map().nOldCells(), sumOp<label>())
<< " to " << mesh_.globalData().nTotalCells()
<< " cells." << endl;
@ -238,6 +225,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
}
}
// Hard code switch for this stage for testing
bool removeOutsideCells = false;
if (removeOutsideCells)
@ -277,7 +265,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
mesh_,
false, // inflate
true, // syncParallel
true, // orderCells (to reduce cell motion in scotch)
true, // orderCells (to reduce cell transfers)
false // orderPoints
);
@ -285,7 +273,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
mesh_.updateMesh(map);
// Update numbering of cells/vertices.
meshCutter.updateMesh(map);
meshCutter_.updateMesh(map);
cellRemover.updateMesh(map);
{
@ -323,7 +311,8 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
if (debug)
{
const_cast<Time&>(mesh_.time())++;
meshCutter.write();
Info<< "Time " << mesh_.time().timeName() << endl;
meshCutter_.write();
mesh_.write();
cellWeights.write();
}
@ -340,7 +329,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
autoPtr<mapDistributePolyMesh> mapDist =
distributor.distribute(newDecomp);
meshCutter.distribute(mapDist);
meshCutter_.distribute(mapDist);
mapDist().distributeCellData(volumeStatus);
@ -349,7 +338,8 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
printMeshData(mesh_);
const_cast<Time&>(mesh_.time())++;
meshCutter.write();
Info<< "Time " << mesh_.time().timeName() << endl;
meshCutter_.write();
mesh_.write();
cellWeights.write();
}
@ -359,6 +349,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
if (debug)
{
const_cast<Time&>(mesh_.time())++;
Info<< "Time " << mesh_.time().timeName() << endl;
cellWeights.write();
mesh_.write();
}
@ -424,7 +415,6 @@ void Foam::backgroundMeshDecomposition::printMeshData
bool Foam::backgroundMeshDecomposition::refineCell
(
const polyMesh& mesh,
label cellI,
label volType,
scalar& weightEstimate
@ -437,10 +427,10 @@ bool Foam::backgroundMeshDecomposition::refineCell
treeBoundBox cellBb
(
mesh.cells()[cellI].points
mesh_.cells()[cellI].points
(
mesh.faces(),
mesh.points()
mesh_.faces(),
mesh_.points()
)
);
@ -549,21 +539,18 @@ bool Foam::backgroundMeshDecomposition::refineCell
Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
(
const hexRef8& meshCutter,
labelList& volumeStatus,
volScalarField& cellWeights
) const
{
labelHashSet cellsToRefine;
const polyMesh& mesh = meshCutter.mesh();
// Determine/update the status of each cell
forAll(volumeStatus, cellI)
{
if (volumeStatus[cellI] == searchableSurface::MIXED)
{
if (meshCutter.cellLevel()[cellI] < minLevels_)
if (meshCutter_.cellLevel()[cellI] < minLevels_)
{
cellsToRefine.insert(cellI);
}
@ -575,7 +562,6 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
(
refineCell
(
mesh,
cellI,
volumeStatus[cellI],
cellWeights.internalField()[cellI]
@ -693,6 +679,12 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
IOobject::MUST_READ
)
),
meshCutter_
(
mesh_,
labelList(mesh_.nCells(), 0),
labelList(mesh_.nPoints(), 0)
),
boundaryFacesPtr_(),
bFTreePtr_(),
decomposeDict_
@ -714,7 +706,8 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
coeffsDict_.lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
),
minLevels_(readLabel(coeffsDict_.lookup("minLevels"))),
volRes_(readLabel(coeffsDict_.lookup("sampleResolution")))
volRes_(readLabel(coeffsDict_.lookup("sampleResolution"))),
maxCellWeightCoeff_(readScalar(coeffsDict_.lookup("maxCellWeightCoeff")))
{
if (!Pstream::parRun())
{
@ -761,16 +754,176 @@ Foam::backgroundMeshDecomposition::~backgroundMeshDecomposition()
Foam::autoPtr<Foam::mapDistributePolyMesh>
Foam::backgroundMeshDecomposition::distribute
(
volScalarField& cellWeights
volScalarField& cellWeights,
List<DynamicList<point> >& cellVertices
)
{
if (debug)
{
const_cast<Time&>(mesh_.time())++;
Info<< "Time " << mesh_.time().timeName() << endl;
cellWeights.write();
mesh_.write();
}
while (true)
{
// Refine large cells if necessary
scalar cellWeightLimit =
maxCellWeightCoeff_
*sum(cellWeights).value()
/mesh_.globalData().nTotalCells();
if (debug)
{
Info<< " cellWeightLimit " << cellWeightLimit << endl;
Pout<< " sum(cellWeights) "
<< sum(cellWeights.internalField())
<< " max(cellWeights) "
<< max(cellWeights.internalField())
<< endl;
}
labelHashSet cellsToRefine;
forAll(cellWeights, cWI)
{
if (cellWeights.internalField()[cWI] > cellWeightLimit)
{
cellsToRefine.insert(cWI);
}
}
if (returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
{
break;
}
// Maintain 2:1 ratio
labelList newCellsToRefine
(
meshCutter_.consistentRefinement
(
cellsToRefine.toc(),
true // extend set
)
);
if (debug && !cellsToRefine.empty())
{
Pout<< " cellWeights too large in " << cellsToRefine.size()
<< " cells" << endl;
}
forAll(newCellsToRefine, nCTRI)
{
label cellI = newCellsToRefine[nCTRI];
cellWeights.internalField()[cellI] /= 8.0;
}
// Mesh changing engine.
polyTopoChange meshMod(mesh_);
// Play refinement commands into mesh changer.
meshCutter_.setRefinement(newCellsToRefine, meshMod);
// Create mesh, return map from old to new mesh.
autoPtr<mapPolyMesh> map = meshMod.changeMesh
(
mesh_,
false, // inflate
true, // syncParallel
true, // orderCells (to reduce cell motion)
false // orderPoints
);
// Update fields
mesh_.updateMesh(map);
// Update numbering of cells/vertices.
meshCutter_.updateMesh(map);
{
// Map cellVertices
meshSearch cellSearch(mesh_);
const labelList& reverseCellMap = map().reverseCellMap();
List<DynamicList<point> > newCellVertices(mesh_.nCells());
forAll(cellVertices, oldCellI)
{
DynamicList<point>& oldCellVertices =
cellVertices[oldCellI];
if (findIndex(newCellsToRefine, oldCellI) >= 0)
{
// This old cell was refined so the cell for the vertices
// in the new mesh needs to be searched for.
forAll (oldCellVertices, oPI)
{
const point& v = oldCellVertices[oPI];
label newCellI = cellSearch.findCell(v);
if (newCellI == -1)
{
Pout<< "findCell backgroundMeshDecomposition "
<< v << " "
<< oldCellI
<< newCellI
<< endl;
}
newCellVertices[newCellI].append(v);
}
}
else
{
label newCellI = reverseCellMap[oldCellI];
forAll(oldCellVertices, oPI)
{
const point& v = oldCellVertices[oPI];
newCellVertices[newCellI].append(v);
}
}
}
cellVertices.transfer(newCellVertices);
}
if (debug)
{
Info<< " Background mesh refined from "
<< returnReduce(map().nOldCells(), sumOp<label>())
<< " to " << mesh_.globalData().nTotalCells()
<< " cells." << endl;
const_cast<Time&>(mesh_.time())++;
Info<< "Time " << mesh_.time().timeName() << endl;
cellWeights.write();
mesh_.write();
}
}
if (debug)
{
printMeshData(mesh_);
Pout<< " Pre distribute sum(cellWeights) "
<< sum(cellWeights.internalField())
<< " max(cellWeights) "
<< max(cellWeights.internalField())
<< endl;
}
labelList newDecomp = decomposerPtr_().decompose
(
mesh_,
@ -778,20 +931,32 @@ Foam::backgroundMeshDecomposition::distribute
cellWeights
);
Info<< " Redistributing background mesh cells" << endl;
fvMeshDistribute distributor(mesh_, mergeDist_);
autoPtr<mapDistributePolyMesh> mapDist =
distributor.distribute(newDecomp);
autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
meshCutter_.distribute(mapDist);
if (debug)
{
printMeshData(mesh_);
Pout<< " Post distribute sum(cellWeights) "
<< sum(cellWeights.internalField())
<< " max(cellWeights) "
<< max(cellWeights.internalField())
<< endl;
const_cast<Time&>(mesh_.time())++;
Info<< "Time " << mesh_.time().timeName() << endl;
mesh_.write();
cellWeights.write();
}
mapDist().distributeCellData(cellVertices);
buildPatchAndTree();
return mapDist;

12
src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
View File

@ -101,6 +101,9 @@ class backgroundMeshDecomposition
// is responsible for.
fvMesh mesh_;
//- Refinement object
hexRef8 meshCutter_;
//- Patch containing an independent representation of the surface of the
// mesh of this processor
autoPtr<bPatch> boundaryFacesPtr_;
@ -131,6 +134,10 @@ class backgroundMeshDecomposition
// investigate the local cell size
label volRes_;
//- Allowed factor above the average cell weight before a background
// cell needs to be split
scalar maxCellWeightCoeff_;
// Private Member Functions
@ -144,7 +151,6 @@ class backgroundMeshDecomposition
// it
bool refineCell
(
const polyMesh& mesh,
label cellI,
label volType,
scalar& weightEstimate
@ -154,7 +160,6 @@ class backgroundMeshDecomposition
// meshed
labelList selectRefinementCells
(
const hexRef8& meshCutter,
labelList& volumeStatus,
volScalarField& cellWeights
) const;
@ -195,7 +200,8 @@ public:
// returning a map to use to redistribute vertices.
autoPtr<mapDistributePolyMesh> distribute
(
volScalarField& cellWeights
volScalarField& cellWeights,
List<DynamicList<point> >& cellVertices
);
//- Is the given position inside the domain of this decomposition

224
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -859,111 +859,150 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
void Foam::conformalVoronoiMesh::distribute()
{
if (!Pstream::parRun())
{
return;
}
Info<< nl << "Redistributing points" << endl;
timeCheck("Before distribute");
const fvMesh& bMesh = decomposition_().mesh();
volScalarField cellWeights
(
IOobject
(
"cellWeights",
bMesh.time().timeName(),
bMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
bMesh,
dimensionedScalar("weight", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
meshSearch cellSearch(bMesh);
List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
while (true)
{
if (vit->real())
{
Foam::point v = topoint(vit->point());
label cellI = cellSearch.findCell(v);
if (cellI == -1)
{
Pout<< "findCell "
<< vit->type() << " "
<< vit->index() << " "
<< v << " "
<< cellI << endl;
}
else
{
cellVertices[cellI].append(topoint(vit->point()));
}
}
}
forAll(cellVertices, cI)
{
cellWeights.internalField()[cI] = min
scalar globalNumberOfVertices = returnReduce
(
max(cellVertices[cI].size(), 1e-3),
100
label(number_of_vertices()),
sumOp<label>()
);
}
autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
(
cellWeights
);
scalar unbalance = returnReduce
(
mag
(
1.0
- label(number_of_vertices())
/(globalNumberOfVertices/Pstream::nProcs())
),
maxOp<scalar>()
);
Pout<< "number_of_vertices() before distribution "
<< label(number_of_vertices()) << endl;
Info<< " Processor unbalance " << unbalance << endl;
mapDist().distributeCellData(cellVertices);
// Remove the entire tessellation
this->clear();
// Info<< "NEED TO MAP FEATURE POINTS" << endl;
// reinsertFeaturePoints();
startOfInternalPoints_ = number_of_vertices();
timeCheck("Distribution performed");
Info<< nl << "Inserting distributed tessellation" << endl;
std::vector<Point> pointsToInsert;
forAll(cellVertices, cI)
{
forAll(cellVertices[cI], cVPI)
if (unbalance <= cvMeshControls().maxLoadUnbalance())
{
pointsToInsert.push_back(toPoint(cellVertices[cI][cVPI]));
break;
}
Info<< " Total number of vertices before redistribution "
<< globalNumberOfVertices
<< endl;
// Pout<< " number_of_vertices() before distribution "
// << label(number_of_vertices()) << endl;
const fvMesh& bMesh = decomposition_().mesh();
volScalarField cellWeights
(
IOobject
(
"cellWeights",
bMesh.time().timeName(),
bMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
bMesh,
dimensionedScalar("weight", dimless, 1e-3),
zeroGradientFvPatchScalarField::typeName
);
meshSearch cellSearch(bMesh);
List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->real())
{
Foam::point v = topoint(vit->point());
label cellI = cellSearch.findCell(v);
if (cellI == -1)
{
Pout<< "findCell conformalVoronoiMesh::distribute findCell "
<< vit->type() << " "
<< vit->index() << " "
<< v << " "
<< cellI << endl;
}
else
{
cellVertices[cellI].append(topoint(vit->point()));
}
}
}
forAll(cellVertices, cI)
{
cellWeights.internalField()[cI] = max
(
cellVertices[cI].size(),
1e-2 // Small but finite weight for empty cells
);
}
autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
(
cellWeights,
cellVertices
);
// Remove the entire tessellation
this->clear();
Info<< "NEED TO MAP FEATURE POINTS" << endl;
// reinsertFeaturePoints();
startOfInternalPoints_ = number_of_vertices();
timeCheck("Distribution performed");
Info<< nl << " Inserting distributed tessellation" << endl;
std::vector<Point> pointsToInsert;
forAll(cellVertices, cI)
{
forAll(cellVertices[cI], cVPI)
{
pointsToInsert.push_back(toPoint(cellVertices[cI][cVPI]));
}
}
insertPoints(pointsToInsert);
Info<< " Total number of vertices after redistribution "
<< returnReduce(label(number_of_vertices()), sumOp<label>())
<< endl;
// Pout<< " number_of_vertices() after distribution "
// << label(number_of_vertices()) << endl;
if(cvMeshControls().objOutput())
{
writePoints("distributedPoints.obj", true);
}
timeCheck("After distribute");
}
insertPoints(pointsToInsert);
Pout<< "number_of_vertices() after distribution "
<< label(number_of_vertices()) << endl;
if(cvMeshControls().objOutput())
{
writePoints("distributedPoints.obj", true);
}
timeCheck("After distribute");
}
@ -1377,6 +1416,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
);
}
Info<< "CREATEFEATUREPOINTS MOVED" << endl;
// createFeaturePoints();
if (cvMeshControls().objOutput())