songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into molecularDynamics

Browse Source
This commit is contained in:
graham 2008-08-21 14:16:05 +01:00
commit d7981e2e99
43 changed files with 1314 additions and 310 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
View File

@ -23,7 +23,7 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
rhoSimpleFoam
rhoPimpleFoam
Description
Transient solver for turbulent flow of compressible fluids for

2
applications/solvers/incompressible/pimpleFoam/Make/files
View File

@ -1,3 +1,3 @@
pimpleFoam.C
EXE = $(FOAM_USER_APPBIN)/pimpleFoam
EXE = $(FOAM_APPBIN)/pimpleFoam

1
applications/test/speed/Make/options
View File

@ -0,0 +1 @@
EXE_INC = /* -ffast-math -mtune=core2 */

2
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
View File

@ -370,7 +370,7 @@ int main(int argc, char *argv[])
// Refinement parameters
refinementParameters refineParams(refineDict);
refineDriver.doRefine(refineDict, refineParams, wantSnap);
refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
writeMesh
(

4
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
View File

@ -491,6 +491,8 @@ void Foam::autoHexMeshDriver::doMesh()
if (wantRefine)
{
const dictionary& motionDict = dict_.subDict("motionDict");
autoRefineDriver refineDriver
(
meshRefinerPtr_(),
@ -502,7 +504,7 @@ void Foam::autoHexMeshDriver::doMesh()
// Get all the refinement specific params
refinementParameters refineParams(dict_, 1);
refineDriver.doRefine(dict_, refineParams, wantSnap);
refineDriver.doRefine(dict_, refineParams, wantSnap, motionDict);
// Write mesh
writeMesh("Refined mesh");

15
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
View File

@ -497,7 +497,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
void Foam::autoRefineDriver::baffleAndSplitMesh
(
const refinementParameters& refineParams,
const bool handleSnapProblems
const bool handleSnapProblems,
const dictionary& motionDict
)
{
Info<< nl
@ -514,6 +515,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
(
handleSnapProblems,
!handleSnapProblems, // merge free standing baffles?
motionDict,
const_cast<Time&>(mesh.time()),
globalToPatch_,
refineParams.keepPoints()[0]
@ -568,7 +570,8 @@ void Foam::autoRefineDriver::zonify
void Foam::autoRefineDriver::splitAndMergeBaffles
(
const refinementParameters& refineParams,
const bool handleSnapProblems
const bool handleSnapProblems,
const dictionary& motionDict
)
{
Info<< nl
@ -588,6 +591,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
(
handleSnapProblems,
false, // merge free standing baffles?
motionDict,
const_cast<Time&>(mesh.time()),
globalToPatch_,
refineParams.keepPoints()[0]
@ -685,7 +689,8 @@ void Foam::autoRefineDriver::doRefine
(
const dictionary& refineDict,
const refinementParameters& refineParams,
const bool prepareForSnapping
const bool prepareForSnapping,
const dictionary& motionDict
)
{
Info<< nl
@ -734,13 +739,13 @@ void Foam::autoRefineDriver::doRefine
);
// Introduce baffles at surface intersections
baffleAndSplitMesh(refineParams, prepareForSnapping);
baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
// Mesh is at its finest. Do optional zoning.
zonify(refineParams);
// Pull baffles apart
splitAndMergeBaffles(refineParams, prepareForSnapping);
splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
// Do something about cells with refined faces on the boundary
if (prepareForSnapping)

10
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
View File

@ -46,7 +46,6 @@ namespace Foam
class featureEdgeMesh;
class refinementParameters;
/*---------------------------------------------------------------------------*\
Class autoRefineDriver Declaration
\*---------------------------------------------------------------------------*/
@ -112,7 +111,8 @@ class autoRefineDriver
void baffleAndSplitMesh
(
const refinementParameters& refineParams,
const bool handleSnapProblems
const bool handleSnapProblems,
const dictionary& motionDict
);
//- Add zones
@ -121,7 +121,8 @@ class autoRefineDriver
void splitAndMergeBaffles
(
const refinementParameters& refineParams,
const bool handleSnapProblems
const bool handleSnapProblems,
const dictionary& motionDict
);
//- Merge refined boundary faces (from exposing coarser cell)
@ -163,7 +164,8 @@ public:
(
const dictionary& refineDict,
const refinementParameters& refineParams,
const bool prepareForSnapping
const bool prepareForSnapping,
const dictionary& motionDict
);
};

168
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -1294,6 +1294,169 @@ void Foam::autoSnapDriver::scaleMesh
}
// After snapping: correct patching according to nearest surface.
// Code is very similar to calcNearestSurface.
// - calculate face-wise snap distance as max of point-wise
// - calculate face-wise nearest surface point
// - repatch face according to patch for surface point.
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
(
const snapParameters& snapParams,
const labelList& adaptPatchIDs
)
{
const fvMesh& mesh = meshRefiner_.mesh();
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
Info<< "Repatching faces according to nearest surface ..." << endl;
// Get the labels of added patches.
autoPtr<indirectPrimitivePatch> ppPtr
(
meshRefinement::makePatch
(
mesh,
adaptPatchIDs
)
);
indirectPrimitivePatch& pp = ppPtr();
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces;
labelList unzonedSurfaces;
getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
// Faces that do not move
PackedList<1> isZonedFace(mesh.nFaces(), 0);
{
// 1. All faces on zoned surfaces
const wordList& faceZoneNames = surfaces.faceZoneNames();
const faceZoneMesh& fZones = mesh.faceZones();
forAll(zonedSurfaces, i)
{
label zoneSurfI = zonedSurfaces[i];
label zoneI = fZones.findZoneID(faceZoneNames[zoneSurfI]);
const faceZone& fZone = fZones[zoneI];
forAll(fZone, i)
{
isZonedFace.set(fZone[i], 1);
}
}
}
// Determine per pp face which patch it should be in
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Patch that face should be in
labelList closestPatch(pp.size(), -1);
{
// face snap distance as max of point snap distance
scalarField faceSnapDist(pp.size(), -GREAT);
{
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, pp));
const faceList& localFaces = pp.localFaces();
forAll(localFaces, faceI)
{
const face& f = localFaces[faceI];
forAll(f, fp)
{
faceSnapDist[faceI] = max
(
faceSnapDist[faceI],
snapDist[f[fp]]
);
}
}
}
pointField localFaceCentres(pp.size());
forAll(pp, i)
{
localFaceCentres[i] = mesh.faceCentres()[pp.addressing()[i]];
}
// Get nearest surface and region
labelList hitSurface;
labelList hitRegion;
surfaces.findNearestRegion
(
unzonedSurfaces,
localFaceCentres,
sqr(4*faceSnapDist), // sqr of attract distance
hitSurface,
hitRegion
);
// Get patch
forAll(pp, i)
{
label faceI = pp.addressing()[i];
if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
{
closestPatch[i] = globalToPatch_
[
surfaces.globalRegion
(
hitSurface[i],
hitRegion[i]
)
];
}
}
}
// Change those faces for which there is a different closest patch
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList ownPatch(mesh.nFaces(), -1);
labelList neiPatch(mesh.nFaces(), -1);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
ownPatch[pp.start()+i] = patchI;
neiPatch[pp.start()+i] = patchI;
}
}
label nChanged = 0;
forAll(closestPatch, i)
{
label faceI = pp.addressing()[i];
if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
{
ownPatch[faceI] = closestPatch[i];
neiPatch[faceI] = closestPatch[i];
nChanged++;
}
}
Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
<< " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
<< endl;
return meshRefiner_.createBaffles(ownPatch, neiPatch);
}
void Foam::autoSnapDriver::doSnap
(
const dictionary& snapDict,
@ -1329,7 +1492,7 @@ void Foam::autoSnapDriver::doSnap
);
indirectPrimitivePatch& pp = ppPtr();
// Distance to attact to nearest feature on surface
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, pp));
@ -1383,6 +1546,9 @@ void Foam::autoSnapDriver::doSnap
// Merge any introduced baffles.
mergeZoneBaffles(baffles);
// Repatch faces according to nearest.
repatchToSurface(snapParams, adaptPatchIDs);
}

8
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -215,6 +215,14 @@ public:
motionSmoother&
);
//- Repatch faces according to surface nearest the face centre
autoPtr<mapPolyMesh> repatchToSurface
(
const snapParameters& snapParams,
const labelList& adaptPatchIDs
);
void doSnap
(
const dictionary& snapDict,

437
src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
View File

@ -50,6 +50,7 @@ License
#include "globalIndex.H"
#include "meshTools.H"
#include "OFstream.H"
#include "geomDecomp.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -426,15 +427,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
// Determine for multi-processor regions the lowest numbered cell on the lowest
// numbered processor.
void Foam::meshRefinement::getRegionMaster
void Foam::meshRefinement::getCoupledRegionMaster
(
const globalIndex& globalCells,
const boolList& blockedFace,
const regionSplit& globalRegion,
Map<label>& regionToMaster
) const
{
globalIndex globalCells(mesh_.nCells());
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
@ -483,6 +483,274 @@ void Foam::meshRefinement::getRegionMaster
}
void Foam::meshRefinement::calcLocalRegions
(
const globalIndex& globalCells,
const labelList& globalRegion,
const Map<label>& coupledRegionToMaster,
Map<label>& globalToLocalRegion,
pointField& localPoints
) const
{
globalToLocalRegion.resize(globalRegion.size());
DynamicList<point> localCc(globalRegion.size()/2);
forAll(globalRegion, cellI)
{
Map<label>::const_iterator fndMaster =
coupledRegionToMaster.find(globalRegion[cellI]);
if (fndMaster != coupledRegionToMaster.end())
{
// Multi-processor region.
if (globalCells.toGlobal(cellI) == fndMaster())
{
// I am master. Allocate region for me.
globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
localCc.append(mesh_.cellCentres()[cellI]);
}
}
else
{
// Single processor region.
if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
{
localCc.append(mesh_.cellCentres()[cellI]);
}
}
}
localCc.shrink();
localPoints.transfer(localCc);
if (localPoints.size() != globalToLocalRegion.size())
{
FatalErrorIn("calcLocalRegions(..)")
<< "localPoints:" << localPoints.size()
<< " globalToLocalRegion:" << globalToLocalRegion.size()
<< exit(FatalError);
}
}
Foam::label Foam::meshRefinement::getShiftedRegion
(
const globalIndex& indexer,
const Map<label>& globalToLocalRegion,
const Map<label>& coupledRegionToShifted,
const label globalRegion
)
{
Map<label>::const_iterator iter =
globalToLocalRegion.find(globalRegion);
if (iter != globalToLocalRegion.end())
{
// Region is 'owned' locally. Convert local region index into global.
return indexer.toGlobal(iter());
}
else
{
return coupledRegionToShifted[globalRegion];
}
}
// Add if not yet present
void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
{
if (findIndex(lst, elem) == -1)
{
label sz = lst.size();
lst.setSize(sz+1);
lst[sz] = elem;
}
}
void Foam::meshRefinement::calcRegionRegions
(
const labelList& globalRegion,
const Map<label>& globalToLocalRegion,
const Map<label>& coupledRegionToMaster,
labelListList& regionRegions
) const
{
// Global region indexing since we now know the shifted regions.
globalIndex shiftIndexer(globalToLocalRegion.size());
// Redo the coupledRegionToMaster to be in shifted region indexing.
Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
forAllConstIter(Map<label>, coupledRegionToMaster, iter)
{
label region = iter.key();
Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
if (fndRegion != globalToLocalRegion.end())
{
// A local cell is master of this region. Get its shifted region.
coupledRegionToShifted.insert
(
region,
shiftIndexer.toGlobal(fndRegion())
);
}
Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
Pstream::mapCombineScatter(coupledRegionToShifted);
}
// Determine region-region connectivity.
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This is for all my regions (so my local ones or the ones I am master
// of) the neighbouring region indices.
// Transfer lists.
PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
forAll(regionConnectivity, procI)
{
if (procI != Pstream::myProcNo())
{
regionConnectivity.set
(
procI,
new HashSet<edge, Hash<edge> >
(
coupledRegionToShifted.size()
/ Pstream::nProcs()
)
);
}
}
// Connectivity. For all my local regions the connected regions.
regionRegions.setSize(globalToLocalRegion.size());
// Add all local connectivity to regionRegions; add all non-local
// to the transferlists.
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
if (ownRegion != neiRegion)
{
label shiftOwnRegion = getShiftedRegion
(
shiftIndexer,
globalToLocalRegion,
coupledRegionToShifted,
ownRegion
);
label shiftNeiRegion = getShiftedRegion
(
shiftIndexer,
globalToLocalRegion,
coupledRegionToShifted,
neiRegion
);
// Connection between two regions. Send to owner of region.
// - is ownRegion 'owned' by me
// - is neiRegion 'owned' by me
if (shiftIndexer.isLocal(shiftOwnRegion))
{
label localI = shiftIndexer.toLocal(shiftOwnRegion);
addUnique(shiftNeiRegion, regionRegions[localI]);
}
else
{
label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
regionConnectivity[masterProc].insert
(
edge(shiftOwnRegion, shiftNeiRegion)
);
}
if (shiftIndexer.isLocal(shiftNeiRegion))
{
label localI = shiftIndexer.toLocal(shiftNeiRegion);
addUnique(shiftOwnRegion, regionRegions[localI]);
}
else
{
label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
regionConnectivity[masterProc].insert
(
edge(shiftOwnRegion, shiftNeiRegion)
);
}
}
}
// Send
forAll(regionConnectivity, procI)
{
if (procI != Pstream::myProcNo())
{
OPstream str(Pstream::blocking, procI);
str << regionConnectivity[procI];
}
}
// Receive
forAll(regionConnectivity, procI)
{
if (procI != Pstream::myProcNo())
{
IPstream str(Pstream::blocking, procI);
str >> regionConnectivity[procI];
}
}
// Add to addressing.
forAll(regionConnectivity, procI)
{
if (procI != Pstream::myProcNo())
{
for
(
HashSet<edge, Hash<edge> >::const_iterator iter =
regionConnectivity[procI].begin();
iter != regionConnectivity[procI].end();
++iter
)
{
const edge& e = iter.key();
bool someLocal = false;
if (shiftIndexer.isLocal(e[0]))
{
label localI = shiftIndexer.toLocal(e[0]);
addUnique(e[1], regionRegions[localI]);
someLocal = true;
}
if (shiftIndexer.isLocal(e[1]))
{
label localI = shiftIndexer.toLocal(e[1]);
addUnique(e[0], regionRegions[localI]);
someLocal = true;
}
if (!someLocal)
{
FatalErrorIn("calcRegionRegions(..)")
<< "Received from processor " << procI
<< " connection " << e
<< " where none of the elements is local to me."
<< abort(FatalError);
}
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -598,88 +866,88 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
// the region might span multiple domains so we want to get
// a "region master" per domain. Note that multi-processor
// regions can only occur on cells on coupled patches.
// Determine per coupled region the master cell (lowest numbered cell
// on lowest numbered processor)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Map<label> regionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
getRegionMaster(blockedFace, globalRegion, regionToMaster);
// Note: since the number of regions does not change by this the
// process can be seen as just a shift of a region onto a single
// processor.
// Global cell numbering engine
globalIndex globalCells(mesh_.nCells());
// Determine cell centre per region
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Now we can divide regions into
// - single-processor regions (almost all)
// - allocate local region + coordinate for it
// - multi-processor for which I am master
// - allocate local region + coordinate for it
// - multi-processor for which I am not master
// - do not allocate region.
// - but inherit new distribution from master.
// Determine per coupled region the master cell (lowest numbered cell
// on lowest numbered processor)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// (does not determine master for non-coupled=fully-local regions)
Map<label> globalToLocalRegion(mesh_.nCells());
DynamicList<point> localCc(mesh_.nCells()/10);
Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
getCoupledRegionMaster
(
globalCells,
blockedFace,
globalRegion,
coupledRegionToMaster
);
forAll(globalRegion, cellI)
{
Map<label>::const_iterator fndMaster =
regionToMaster.find(globalRegion[cellI]);
// Determine my regions
// ~~~~~~~~~~~~~~~~~~~~
// These are the fully local ones or the coupled ones of which I am master.
if (fndMaster != regionToMaster.end())
Map<label> globalToLocalRegion;
pointField localPoints;
calcLocalRegions
(
globalCells,
globalRegion,
coupledRegionToMaster,
globalToLocalRegion,
localPoints
);
// Find distribution for regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList regionDistribution;
if (isA<geomDecomp>(decomposer))
{
// Multi-processor region.
if (globalCells.toGlobal(cellI) == fndMaster())
{
// I am master. Allocate region for me.
globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
localCc.append(mesh_.cellCentres()[cellI]);
}
regionDistribution = decomposer.decompose(localPoints);
}
else
{
// Single processor region.
Map<label>::iterator iter =
globalToLocalRegion.find(globalRegion[cellI]);
if (iter == globalToLocalRegion.end())
{
globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
localCc.append(mesh_.cellCentres()[cellI]);
}
}
}
localCc.shrink();
pointField localPoints;
localPoints.transfer(localCc);
// Call decomposer on localCc
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList localDistribution = decomposer.decompose(localPoints);
// Distribute the destination processor for coupled regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Map<label> regionToDist(regionToMaster.size());
forAllConstIter(Map<label>, regionToMaster, iter)
{
if (globalCells.isLocal(iter()))
{
// Master cell is local.
regionToDist.insert
labelListList regionRegions;
calcRegionRegions
(
iter.key(),
localDistribution[globalToLocalRegion[iter.key()]]
globalRegion,
globalToLocalRegion,
coupledRegionToMaster,
regionRegions
);
regionDistribution = decomposer.decompose(regionRegions, localPoints);
}
// Convert region-based decomposition back to cell-based one
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Transfer destination processor back to all. Use global reduce for now.
Map<label> regionToDist(coupledRegionToMaster.size());
forAllConstIter(Map<label>, coupledRegionToMaster, iter)
{
label region = iter.key();
Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
if (regionFnd != globalToLocalRegion.end())
{
// Master cell is local. Store my distribution.
regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
}
else
{
@ -691,29 +959,27 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
Pstream::mapCombineScatter(regionToDist);
// Convert region-based decomposition back to cell-based one
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Determine destination for all cells
labelList distribution(mesh_.nCells());
forAll(globalRegion, cellI)
{
Map<label>::const_iterator fndMaster =
Map<label>::const_iterator fndRegion =
regionToDist.find(globalRegion[cellI]);
if (fndMaster != regionToDist.end())
if (fndRegion != regionToDist.end())
{
// Special handling
distribution[cellI] = fndMaster();
distribution[cellI] = fndRegion();
}
else
{
// region is local to the processor.
label localRegionI = globalToLocalRegion[globalRegion[cellI]];
distribution[cellI] = localDistribution[localRegionI];
distribution[cellI] = regionDistribution[localRegionI];
}
}
return distribution;
}
@ -839,9 +1105,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
if (debug)
{
Pout<< "Wanted distribution:"
<< distributor.countCells(distribution)
<< endl;
labelList nProcCells(distributor.countCells(distribution));
Pout<< "Wanted distribution:" << nProcCells << endl;
Pstream::listCombineGather(nProcCells, plusEqOp<label>());
Pstream::listCombineScatter(nProcCells);
Pout<< "Wanted resulting decomposition:" << endl;
forAll(nProcCells, procI)
{
Pout<< " " << procI << '\t' << nProcCells[procI] << endl;
}
Pout<< endl;
}
// Do actual sending/receiving of mesh
map = distributor.distribute(distribution);

57
src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -68,6 +68,7 @@ class featureEdgeMesh;
class fvMeshDistribute;
class searchableSurface;
class regionSplit;
class globalIndex;
/*---------------------------------------------------------------------------*\
Class meshRefinement Declaration
@ -157,15 +158,6 @@ private:
labelList& neiLevel,
pointField& neiCc
) const;
////- Calculate coupled boundary end vector and refinement level. Sort
//// so both sides of coupled face have data in same order.
//void calcCanonicalBoundaryData
//(
// labelList& leftLevel,
// pointField& leftCc,
// labelList& rightLevel,
// pointField& rightCc
//) const;
//- Find any intersection of surface. Store in surfaceIndex_.
void updateIntersections(const labelList& changedFaces);
@ -194,16 +186,52 @@ private:
const label exposedPatchI
);
// For decomposeCombineRegions
//- Used in decomposeCombineRegions. Given global region per cell
// determines master processor/cell for regions straddling
// procboundaries.
void getRegionMaster
void getCoupledRegionMaster
(
const globalIndex& globalCells,
const boolList& blockedFace,
const regionSplit& globalRegion,
Map<label>& regionToMaster
) const;
//- Determine regions that are local to me or coupled ones that
// are owned by me. Determine representative location.
void calcLocalRegions
(
const globalIndex& globalCells,
const labelList& globalRegion,
const Map<label>& coupledRegionToMaster,
Map<label>& globalToLocalRegion,
pointField& localPoints
) const;
//- Convert region into global index.
static label getShiftedRegion
(
const globalIndex& indexer,
const Map<label>& globalToLocalRegion,
const Map<label>& coupledRegionToShifted,
const label globalRegion
);
//- helper: add element if not in list. Linear search.
static void addUnique(const label, labelList&);
//- Calculate region connectivity. Major communication.
void calcRegionRegions
(
const labelList& globalRegion,
const Map<label>& globalToLocalRegion,
const Map<label>& coupledRegionToMaster,
labelListList& regionRegions
) const;
// Refinement candidate selection
@ -334,6 +362,14 @@ private:
const labelList& globalToPatch
) const;
//- Returns list with for every internal face -1 or the patch
// they should be baffled into.
labelList markFacesOnProblemCellsGeometric
(
const dictionary& motionDict,
const labelList& globalToPatch
) const;
// Baffle merging
@ -552,6 +588,7 @@ public:
(
const bool handleSnapProblems,
const bool mergeFreeStanding,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToPatch,
const point& keepPoint

334
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -44,9 +44,9 @@ License
#include "OFstream.H"
#include "regionSplit.H"
#include "removeCells.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
#include "motionSmoother.H"
#include "polyMeshGeometry.H"
#include "IOmanip.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -555,12 +555,13 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
boolList isBoundaryFace(mesh_.nFaces(), false);
// Fill boundary data. All elements on meshed patches get marked.
forAll(globalToPatch, i)
{
label patchI = globalToPatch[i];
// Get the labels of added patches.
labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
if (patchI != -1)
forAll(adaptPatchIDs, i)
{
label patchI = adaptPatchIDs[i];
const polyPatch& pp = patches[patchI];
label faceI = pp.start();
@ -578,7 +579,6 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
faceI++;
}
}
}
syncTools::syncPointList
(
@ -872,6 +872,264 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
return facePatch;
}
//XXXXXXXXXXXXXX
// Mark faces to be baffled to prevent snapping problems. Does
// test to find nearest surface and checks which faces would get squashed.
Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
(
const dictionary& motionDict,
const labelList& globalToPatch
) const
{
// Get the labels of added patches.
labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
// Construct addressing engine.
autoPtr<indirectPrimitivePatch> ppPtr
(
meshRefinement::makePatch
(
mesh_,
adaptPatchIDs
)
);
const indirectPrimitivePatch& pp = ppPtr();
const pointField& localPoints = pp.localPoints();
const labelList& meshPoints = pp.meshPoints();
// Find nearest (non-baffle) surface
pointField newPoints(mesh_.points());
{
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces_.findNearest
(
surfaces_.getUnnamedSurfaces(),
localPoints,
scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
hitSurface,
hitInfo
);
forAll(hitInfo, i)
{
if (hitInfo[i].hit())
{
//label pointI = meshPoints[i];
//Pout<< " " << pointI << " moved from "
// << mesh_.points()[pointI] << " by "
// << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
// << endl;
newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
}
}
}
// Per face (internal or coupled!) the patch that the
// baffle should get (or -1).
labelList facePatch(mesh_.nFaces(), -1);
// Count of baffled faces
label nBaffleFaces = 0;
// // Sync position? Or not since same face on both side so just sync
// // result of baffle.
//
// const scalar minArea(readScalar(motionDict.lookup("minArea")));
//
// Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
// << minArea << endl;
//
// pointField facePoints;
// for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
// {
// const face& f = mesh_.faces()[faceI];
//
// bool usesPatchPoint = false;
//
// facePoints.setSize(f.size());
// forAll(f, fp)
// {
// Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
//
// if (iter != pp.meshPointMap().end())
// {
// facePoints[fp] = newPosition[iter()];
// usesPatchPoint = true;
// }
// else
// {
// facePoints[fp] = mesh_.points()[f[fp]];
// }
// }
//
// if (usesPatchPoint)
// {
// // Check area of face wrt original area
// face identFace(identity(f.size()));
//
// if (identFace.mag(facePoints) < minArea)
// {
// facePatch[faceI] = getBafflePatch(facePatch, faceI);
// nBaffleFaces++;
// }
// }
// }
//
//
// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
//
// if (pp.coupled())
// {
// forAll(pp, i)
// {
// label faceI = pp.start()+i;
//
// const face& f = mesh_.faces()[faceI];
//
// bool usesPatchPoint = false;
//
// facePoints.setSize(f.size());
// forAll(f, fp)
// {
// Map<label>::const_iterator iter =
// pp.meshPointMap().find(f[fp]);
//
// if (iter != pp.meshPointMap().end())
// {
// facePoints[fp] = newPosition[iter()];
// usesPatchPoint = true;
// }
// else
// {
// facePoints[fp] = mesh_.points()[f[fp]];
// }
// }
//
// if (usesPatchPoint)
// {
// // Check area of face wrt original area
// face identFace(identity(f.size()));
//
// if (identFace.mag(facePoints) < minArea)
// {
// facePatch[faceI] = getBafflePatch(facePatch, faceI);
// nBaffleFaces++;
// }
// }
// }
// }
// }
{
pointField oldPoints(mesh_.points());
mesh_.movePoints(newPoints);
faceSet wrongFaces(mesh_, "wrongFaces", 100);
{
//motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
// Just check the errors from squashing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const labelList allFaces(identity(mesh_.nFaces()));
label nWrongFaces = 0;
scalar minArea(readScalar(motionDict.lookup("minArea")));
if (minArea > -SMALL)
{
polyMeshGeometry::checkFaceArea
(
false,
minArea,
mesh_,
mesh_.faceAreas(),
allFaces,
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces with area < "
<< setw(5) << minArea
<< " m^2 : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
scalar minDet = 0.01;
if (minDet > -1)
{
polyMeshGeometry::checkCellDeterminant
(
false,
minDet,
mesh_,
mesh_.faceAreas(),
allFaces,
polyMeshGeometry::affectedCells(mesh_, allFaces),
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces on cells with determinant < "
<< setw(5) << minDet << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
}
forAllConstIter(faceSet, wrongFaces, iter)
{
label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
{
facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
nBaffleFaces++;
//Pout<< " " << iter.key()
// //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
// << " is destined for patch " << facePatch[iter.key()]
// << endl;
}
}
// Restore points.
mesh_.movePoints(oldPoints);
}
Info<< "markFacesOnProblemCellsGeometric : marked "
<< returnReduce(nBaffleFaces, sumOp<label>())
<< " additional internal and coupled faces"
<< " to be converted into baffles." << endl;
syncTools::syncFaceList
(
mesh_,
facePatch,
maxEqOp<label>(),
false // no separation
);
return facePatch;
}
//XXXXXXXX
// Return a list of coupled face pairs, i.e. faces that use the same vertices.
@ -1554,6 +1812,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
(
const bool handleSnapProblems,
const bool mergeFreeStanding,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToPatch,
const point& keepPoint
@ -1584,6 +1843,12 @@ void Foam::meshRefinement::baffleAndSplitMesh
ownPatch,
neiPatch
);
if (debug)
{
runTime++;
}
createBaffles(ownPatch, neiPatch);
if (debug)
@ -1621,14 +1886,65 @@ void Foam::meshRefinement::baffleAndSplitMesh
labelList facePatch
(
markFacesOnProblemCells
//markFacesOnProblemCells
//(
// globalToPatch
//)
markFacesOnProblemCellsGeometric
(
motionDict,
globalToPatch
)
);
Info<< "Analyzed problem cells in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
if (debug)
{
// Dump all these faces to a faceSet.
faceSet problemGeom(mesh_, "problemFacesGeom", 100);
const labelList facePatchGeom
(
markFacesOnProblemCellsGeometric
(
motionDict,
globalToPatch
)
);
forAll(facePatchGeom, faceI)
{
if (facePatchGeom[faceI] != -1)
{
problemGeom.insert(faceI);
}
}
Pout<< "Dumping " << problemGeom.size()
<< " problem faces to " << problemGeom.objectPath() << endl;
problemGeom.write();
faceSet problemTopo(mesh_, "problemFacesTopo", 100);
const labelList facePatchTopo
(
markFacesOnProblemCells
(
globalToPatch
)
);
forAll(facePatchTopo, faceI)
{
if (facePatchTopo[faceI] != -1)
{
problemTopo.insert(faceI);
}
}
Pout<< "Dumping " << problemTopo.size()
<< " problem faces to " << problemTopo.objectPath() << endl;
problemTopo.write();
}
Info<< "Introducing baffles to delete problem cells." << nl << endl;
if (debug)

25
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -1277,22 +1277,13 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
if (Pstream::nProcs() > 1)
{
labelList distribution(decomposer.decompose(mesh_.cellCentres()));
// Get distribution such that baffle faces stay internal to the
// processor.
//labelList distribution(decomposePreserveBaffles(decomposer));
if (debug)
{
Pout<< "Wanted distribution:"
<< distributor.countCells(distribution)
<< endl;
}
// Do actual sending/receiving of mesh
distMap = distributor.distribute(distribution);
// Update numbering of meshRefiner
distribute(distMap);
distMap = balance
(
false, //keepZoneFaces
false, //keepBaffles
decomposer,
distributor
);
Info<< "Balanced mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
@ -1302,7 +1293,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
if (debug)
{
Pout<< "Writing " << msg
Pout<< "Writing balanced " << msg
<< " mesh to time " << mesh_.time().timeName() << endl;
write
(

90
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -377,15 +377,34 @@ Foam::refinementSurfaces::refinementSurfaces
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Get indices of named surfaces (surfaces with cellZoneName)
// Get indices of unnamed surfaces (surfaces without faceZoneName)
Foam::labelList Foam::refinementSurfaces::getUnnamedSurfaces() const
{
labelList anonymousSurfaces(faceZoneNames_.size());
label i = 0;
forAll(faceZoneNames_, surfI)
{
if (faceZoneNames_[surfI].size() == 0)
{
anonymousSurfaces[i++] = surfI;
}
}
anonymousSurfaces.setSize(i);
return anonymousSurfaces;
}
// Get indices of named surfaces (surfaces with faceZoneName)
Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
{
labelList namedSurfaces(cellZoneNames_.size());
labelList namedSurfaces(faceZoneNames_.size());
label namedI = 0;
forAll(cellZoneNames_, surfI)
forAll(faceZoneNames_, surfI)
{
if (cellZoneNames_[surfI].size() > 0)
if (faceZoneNames_[surfI].size() > 0)
{
namedSurfaces[namedI++] = surfI;
}
@ -846,6 +865,69 @@ void Foam::refinementSurfaces::findNearest
}
void Foam::refinementSurfaces::findNearestRegion
(
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& hitSurface,
labelList& hitRegion
) const
{
labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
// Do the tests. Note that findNearest returns index in geometries.
List<pointIndexHit> hitInfo;
searchableSurfacesQueries::findNearest
(
allGeometry_,
geometries,
samples,
nearestDistSqr,
hitSurface,
hitInfo
);
// Rework the hitSurface to be surface (i.e. index into surfaces_)
forAll(hitSurface, pointI)
{
if (hitSurface[pointI] != -1)
{
hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
}
}
// Collect the region
hitRegion.setSize(hitSurface.size());
hitRegion = -1;
forAll(surfacesToTest, i)
{
label surfI = surfacesToTest[i];
// Collect hits for surfI
const labelList localIndices(findIndices(hitSurface, surfI));
List<pointIndexHit> localHits
(
IndirectList<pointIndexHit>
(
hitInfo,
localIndices
)
);
labelList localRegion;
allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
forAll(localIndices, i)
{
hitRegion[localIndices[i]] = localRegion[i];
}
}
}
//// Find intersection with max of edge. Return -1 or the surface
//// with the highest maxLevel above currentLevel
//Foam::label Foam::refinementSurfaces::findHighestIntersection

18
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
View File

@ -154,7 +154,10 @@ public:
return cellZoneNames_;
}
//- Get indices of named surfaces (surfaces with cellZoneName)
//- Get indices of named surfaces (surfaces with faceZoneName)
labelList getUnnamedSurfaces() const;
//- Get indices of named surfaces (surfaces with faceZoneName)
labelList getNamedSurfaces() const;
//- Get indices of closed named surfaces
@ -249,6 +252,8 @@ public:
//- Find intersection nearest to the endpoints. surface1,2 are
// not indices into surfacesToTest but refinement surface indices.
// Returns surface, region on surface (so not global surface)
// and position on surface.
void findNearestIntersection
(
const labelList& surfacesToTest,
@ -282,6 +287,17 @@ public:
List<pointIndexHit>&
) const;
//- Find nearest point on surfaces. Return surface and region on
// surface (so not global surface)
void findNearestRegion
(
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& hitSurface,
labelList& hitRegion
) const;
//- Detect if a point is 'inside' (closed) surfaces.
// Returns -1 if not, returns first surface it is.
void findInside

2
src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
View File

@ -301,7 +301,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
// Is any of the faces a boundary face? If so, grab the patch
// A boundary-to-boundary collapse is checked for in validCollapse()
// and cannnot happen here.
// and cannot happen here.
if (!mesh.isInternalFace(mf[faceI]))
{

9
src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
View File

@ -50,7 +50,8 @@ Foam::Analytical<Type>::~Analytical()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::Analytical<Type>::integrate
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Analytical<Type>::integrate
(
const Type phi,
const scalar dt,
@ -58,7 +59,11 @@ Type Foam::Analytical<Type>::integrate
const scalar beta
) const
{
return alpha + (phi - alpha)*exp(-beta*dt);
typename IntegrationScheme<Type>::integrationResult retValue;
retValue.average() = alpha + (phi - alpha)*(1 - exp(-beta*dt))/(beta*dt);
retValue.value() = alpha + (phi - alpha)*exp(-beta*dt);
return retValue;
}

2
src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
View File

@ -74,7 +74,7 @@ public:
// Member Functions
//- Perform the integration
virtual Type integrate
virtual typename IntegrationScheme<Type>::integrationResult integrate
(
const Type phi,
const scalar dt,

9
src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
View File

@ -50,7 +50,8 @@ Foam::Euler<Type>::~Euler()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::Euler<Type>::integrate
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Euler<Type>::integrate
(
const Type phi,
const scalar dt,
@ -58,7 +59,11 @@ Type Foam::Euler<Type>::integrate
const scalar beta
) const
{
return (phi + dt*alpha)/(1.0 + dt/beta);
typename IntegrationScheme<Type>::integrationResult retValue;
retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
retValue.average() = 0.5*(phi + retValue.value());
return retValue;
}

2
src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
View File

@ -74,7 +74,7 @@ public:
// Member Functions
//- Perform the integration
virtual Type integrate
virtual typename IntegrationScheme<Type>::integrationResult integrate
(
const Type phi,
const scalar dt,

59
src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
View File

@ -53,6 +53,63 @@ namespace Foam
template<class Type>
class IntegrationScheme
{
public:
//- Helper class to supply results of integration
class integrationResult
{
//- Integration value
Type value_;
//- Average value across integration step
Type average_;
public:
//- Constructor
integrationResult()
:
value_(pTraits<Type>::zero),
average_(pTraits<Type>::zero)
{}
// Member functions
// Access
//- Return const access to the value
Type value() const
{
return value_;
}
//- Return const access to the average
Type average() const
{
return average_;
}
// Edit
//- Return access to the value for changing
Type& value()
{
return value_;
}
//- Return access to the average for changing
Type& average()
{
return average_;
}
};
private:
// Private data
//- Name of the Integration variable
@ -120,7 +177,7 @@ public:
// Member Functions
//- Perform the Integration
virtual Type integrate
virtual integrationResult integrate
(
const Type phi,
const scalar dt,

12
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -332,21 +332,11 @@ void Foam::KinematicCloud<ParcelType>::evolve()
inject(td);
move(td);
}
template<class ParcelType>
template<class TrackingData>
void Foam::KinematicCloud<ParcelType>::move
(
TrackingData& td
)
{
if (coupled_)
{
resetSourceTerms();
}
Cloud<ParcelType>::move(td);
}

6
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
View File

@ -257,10 +257,6 @@ protected:
ParcelType* p
);
//- Move the parcels
template<class TrackingData>
void move(TrackingData& td);
//- Post-injection checks
void postInjectCheck();
@ -436,7 +432,7 @@ public:
//- Reset the spray source terms
void resetSourceTerms();
//- Evolve the spray (move, inject)
//- Evolve the spray (inject, inject)
void evolve();
};

7
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
View File

@ -185,7 +185,12 @@ void Foam::ReactingCloud<ParcelType>::evolve()
inject(td);
this->move(td);
if (this->coupled())
{
resetSourceTerms();
}
Cloud<ParcelType>::move(td);
}

2
src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
View File

@ -204,7 +204,7 @@ public:
//- Reset the spray source terms
void resetSourceTerms();
//- Evolve the spray (move, inject)
//- Evolve the spray (inject, move)
void evolve();
};

7
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
View File

@ -176,7 +176,12 @@ void Foam::ThermoCloud<ParcelType>::evolve()
inject(td);
this->move(td);
if (this->coupled())
{
resetSourceTerms();
}
Cloud<ParcelType>::move(td);
}

6
src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
View File

@ -209,11 +209,11 @@ public:
// Cloud evolution functions
//- Evolve the spray (move, inject)
void evolve();
//- Reset the spray source terms
void resetSourceTerms();
//- Evolve the spray (inject, move)
void evolve();
};

13
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -155,15 +155,16 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const scalar bp = 1.0/(Cud + VSMALL);
const vector ap = Uc_/bp + rhoc_/rho_*td.g();
const vector ap = Uc_ + (1 - rhoc_/rho_)/(Cud + VSMALL)*td.g();
const scalar bp = Cud;
vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp);
vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp).value();
// Info<< "U_, Unew = " << U_ << ", " << Unew << endl;
// Calculate the momentum transfer to the continuous phase
dUTrans = -mass()*(Unew - U_);
// - do not include gravity impulse
dUTrans = -mass()*(Unew - U_ - dt*td.g());
// Make corrections for 2-D cases
if (meshInfo.caseIs2d())
@ -233,6 +234,9 @@ bool Foam::KinematicParcel<ParcelType>::move
// Update cell based properties
p.updateCellQuantities(td, dt, celli);
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
{
if (td.cloud().coupled())
{
p.calcCoupled(td, dt, celli);
@ -241,6 +245,7 @@ bool Foam::KinematicParcel<ParcelType>::move
{
p.calcUncoupled(td, dt, celli);
}
}
if (p.onBoundary() && td.keepParticle)
{

2
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
View File

@ -203,7 +203,7 @@ Foam::Ostream& Foam::operator<<
os << static_cast<const Particle<ParcelType>& >(p);
os.write
(
reinterpret_cast<const char*>(p.typeId()),
reinterpret_cast<const char*>(&p.typeId_),
sizeof(p.typeId())
+ sizeof(p.d())
+ sizeof(p.U())

149
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -59,7 +59,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
const scalar mass0 = this->mass();
const scalar np0 = this->nParticle_;
const scalar T0 = this->T_;
const scalar cp0 = this->cp_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
@ -79,9 +78,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
// - components do not necessarily exist in particle already
scalarList dMassSR(td.cloud().gases().size(), 0.0);
// Total mass lost from particle due to surface reactions
// - sub-model will adjust component mass fractions
scalar dMassMTSR = 0.0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -91,13 +93,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~
// Calculate mass transfer
// ~~~~~~~~~~~~~~~~~~~~~~~
@ -107,27 +102,36 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR);
// Initialise enthalpy retention to zero
scalar dhRet = 0.0;
calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
// New total mass
const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
// Ratio of mass devolatilised to the total volatile mass of the particle
const scalar fVol = 1 -
(YMixture_[0]*mass1)
/(td.cloud().composition().YMixture0()[0]*mass0_);
// Correct particle temperature to account for latent heat of
// devolatilisation
T1 -=
td.constProps().Ldevol()
*sum(dMassMT)
/(0.5*(mass0 + mass1)*this->cp_);
// Specific heat capacity of non-volatile components
const scalar cpNonVolatile =
(
YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
+ YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
)/(YMixture_[1] + YMixture_[2]);
// Add retained enthalpy from surface reaction to particle and remove
// from gas
T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_);
dhTrans -= dhRet;
// New specific heat capacity - linear variation until volatiles
// have evolved
const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
+ td.constProps().cp0();
// Correct dhTrans to account for enthalpy of evolved volatiles
dhTrans +=
sum(dMassMT)
*td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
// Correct dhTrans to account for enthalpy of consumed solids
dhTrans +=
sum(dMassSR)
*td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
// ~~~~~~~~~~~~~~~~~~~~~~~
@ -147,8 +151,8 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
td.cloud().UCoeff()[celli] += np0*mass0*Cud;
// Update enthalpy transfer
// td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans);
// - enthalpy of lost solids already accounted for
td.cloud().hTrans()[celli] += np0*dhTrans;
// Accumulate coefficient to be applied in carrier phase enthalpy coupling
td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@ -166,7 +170,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
{
td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
}
td.cloud().hTrans()[celli] += np0*mass1*cp1*T1;
td.cloud().hTrans()[celli] +=
np0*mass1
*(
YMixture_[0]*td.cloud().composition().HGas(YGas_, T1)
+ YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1)
);
td.cloud().UTrans()[celli] += np0*mass1*U1;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -176,7 +185,10 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = cp1;
this->cp_ =
YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
+ YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
+ YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
// Update particle density or diameter
if (td.cloud().massTransfer().changesVolume())
@ -205,7 +217,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
// const scalar cp0 = this->cp();
const scalar cp0 = this->cp_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Initialise transfer terms
@ -225,9 +237,12 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
// - components do not necessarily exist in particle already
scalarList dMassSR(td.cloud().gases().size(), 0.0);
// Total mass lost from particle due to surface reactions
// - sub-model will adjust component mass fractions
scalar dMassMTSR = 0.0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -236,16 +251,6 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
scalar htc = 0.0;
scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
// Limit new temp max by vapourisarion temperature
T1 = min(td.constProps().Tvap(), T1);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~
// Calculate mass transfer
@ -256,37 +261,23 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Calculate surface reactions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
calcSurfaceReactions
(
td,
dt,
celli,
T0,
T1,
dMassMTSR,
dMassSR
);
// Initialise enthalpy retention to zero
scalar dhRet = 0.0;
calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
// New total mass
const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
// Ratio of mass devolatilised to the total volatile mass of the particle
const scalar fVol = 1 -
(YMixture_[0]*mass1)
/(td.cloud().composition().YMixture0()[0]*mass0_);
// Specific heat capacity of non-volatile components
const scalar cpNonVolatile =
(
YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
+ YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
)/(YMixture_[1] + YMixture_[2]);
// New specific heat capacity - linear variation until volatiles
// have evolved
const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
+ td.constProps().cp0();
// New specific heat capacity
const scalar cp1 =
YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
+ YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
+ YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
// Add retained enthalpy to particle
T1 += dhRet/(mass0*0.5*(cp0 + cp1));
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Remove the particle when mass falls below minimum threshold
@ -389,8 +380,9 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
const label celli,
const scalar T0,
const scalar T1,
scalar& dMassMTSR,
scalarList& dMassMT
const scalarList& dMassMT,
scalarList& dMassSR,
scalar& dhRet
)
{
// Check that model is active
@ -409,21 +401,20 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
T0,
T1,
this->Tc_,
pc_,
this->rhoc_,
this->mass(),
dMassMT,
YGas_,
YLiquid_,
YSolid_,
YMixture_,
dMassMTSR,
dMassMT
dMassSR,
dhRet
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
#include "ReactingParcelIO.C"

11
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
View File

@ -89,6 +89,9 @@ public:
//- Boiling point [K]
const scalar Tbp_;
//- Latent heat of devolatilisation [J/kg]
const scalar Ldevol_;
public:
@ -102,6 +105,9 @@ public:
//- Return const access to the boiling point
inline scalar Tbp() const;
//- Return const access to the latent heat of devolatilisation
inline scalar Ldevol() const;
};
@ -210,8 +216,9 @@ protected:
const label celli,
const scalar T0,
const scalar T1,
scalar& dMassMTSR,
scalarList& dMassMT
const scalarList& dMassMT,
scalarList& dMassSR,
scalar& dhRet
);

11
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
View File

@ -34,7 +34,8 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
:
ThermoParcel<ParcelType>::constantProperties(dict),
Tvap_(dimensionedScalar(dict.lookup("Tvap")).value()),
Tbp_(dimensionedScalar(dict.lookup("Tbp")).value())
Tbp_(dimensionedScalar(dict.lookup("Tbp")).value()),
Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value())
{}
@ -127,6 +128,14 @@ Foam::ReactingParcel<ParcelType>::constantProperties::Tbp() const
}
template<class ParcelType>
inline Foam::scalar
Foam::ReactingParcel<ParcelType>::constantProperties::Ldevol() const
{
return Ldevol_;
}
// * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * //
template<class ParcelType>

5
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
View File

@ -297,10 +297,7 @@ Foam::Ostream& Foam::operator<<
os << static_cast<const ThermoParcel<ParcelType>& >(p);
os.write
(
reinterpret_cast<const char*>
(
&const_cast<ReactingParcel<ParcelType>&>(p).mass0()
),
reinterpret_cast<const char*>(&p.mass0_),
sizeof(p.mass0())
);
os << p.YMixture() << YGasLoc << YLiquidLoc << YSolidLoc;

12
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
View File

@ -203,17 +203,17 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
// Set new particle temperature
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Tnew = td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
// Integrate to find the new parcel temperature
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
dhTrans = -this->mass()*cp_*(Tnew - T_);
// Using average parcel temperature for enthalpy transfer calculation
dhTrans = dt*this->areaS()*htc*(Tres.average() - Tc_);
return Tnew;
return Tres.value();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
#include "ThermoParcelIO.C"

5
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
View File

@ -149,10 +149,7 @@ Foam::Ostream& Foam::operator<<
os << static_cast<const KinematicParcel<ParcelType>& >(p);
os.write
(
reinterpret_cast<const char*>
(
&const_cast<ThermoParcel<ParcelType>&>(p).T()
),
reinterpret_cast<const char*>(&p.T_),
sizeof(p.T()) + sizeof(p.cp())
);
}

9
src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
View File

@ -209,13 +209,20 @@ public:
//- Return the gas constant for the gas mixture
virtual scalar RGas(const scalarField& YGas) const = 0;
//- Return enthalpy for the gas mixture
//- Return enthalpy for the gas mixture [energy per unit mass]
virtual scalar HGas
(
const scalarField& YGas,
const scalar T
) const = 0;
//- Return enthalpy for the solid mixture [energy per unit mass]
virtual scalar HSolid
(
const scalarField& YSolid,
const scalar T
) const = 0;
//- Return specific heat caparcity for the gas mixture
virtual scalar cpGas
(

23
src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C
View File

@ -519,6 +519,29 @@ const
}
template<class CloudType>
Foam::scalar Foam::SingleMixtureFraction<CloudType>::HSolid
(
const scalarField& YSolid,
const scalar T
)
const
{
scalar HMixture = 0.0;
forAll(YSolid, i)
{
label id = solidGlobalIds_[i];
HMixture +=
YSolid[i]
*(
this->solids().properties()[id].Hf()
+ this->solids().properties()[id].cp()*T
);
}
return HMixture;
}
template<class CloudType>
Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpGas
(

5
src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H
View File

@ -188,9 +188,12 @@ public:
//- Return the gas constant for the gas mixture
scalar RGas(const scalarField& YGas) const;
//- Return enthalpy for the gas mixture
//- Return enthalpy for the gas mixture [energy per unit mass]
scalar HGas(const scalarField& YGas, const scalar T) const;
//- Return enthalpy for the solid mixture [energy per unit mass]
scalar HSolid(const scalarField& YSolid, const scalar T) const;
//- Return specific heat caparcity for the gas mixture
scalar cpGas(const scalarField& YGas, const scalar T) const;

10
src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
View File

@ -64,16 +64,20 @@ void Foam::NoSurfaceReaction<CloudType>::calculate
const scalar T0,
const scalar T1,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar massp,
const scalarList& dMassMT,
scalarField& YGas,
scalarField& YLiquid,
scalarField& YSolid,
scalarField& YMixture,
scalar& dMassMTSR,
scalarList& dMassSR
scalarList& dMassSR,
scalar& dhRet
) const
{}
{
// do nothing
}
// ************************************************************************* //

6
src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
View File

@ -84,14 +84,16 @@ public:
const scalar T0,
const scalar T1,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar massp,
const scalarList& dMassMT,
scalarField& YGas,
scalarField& YLiquid,
scalarField& YSolid,
scalarField& YMixture,
scalar& dMassMTSR,
scalarList& dMassSR
scalarList& dMassSR,
scalar& dhRet
) const;
};

6
src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
View File

@ -141,14 +141,16 @@ public:
const scalar T0,
const scalar T1,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar massp,
const scalarList& dMassMT,
scalarField& YGas,
scalarField& YLiquid,
scalarField& YSolid,
scalarField& YMixture,
scalar& dMassMTSR,
scalarList& dMassSR
scalarList& dMassSR,
scalar& dhRet
) const = 0;
};

8
src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.C
View File

@ -62,11 +62,6 @@ Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu
const scalar
) const
{
notImplemented
(
"Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu"
"(const scalar, const scalar)"
);
return 0.0;
}
@ -74,8 +69,7 @@ Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu
template <class CloudType>
Foam::scalar Foam::NoHeatTransfer<CloudType>::Pr() const
{
notImplemented("Foam::scalar Foam::NoHeatTransfer<CloudType>::Pr()");
return 0.0;
return 1.0;
}

1
src/meshTools/searchableSurface/searchableBox.C
View File

@ -436,7 +436,6 @@ void Foam::searchableBox::findLineAll
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
//XXX
// Tolerances:
// To find all intersections we add a small vector to the last intersection
// This is chosen such that