zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Patches contributed by Mattijs Janssens:

Browse Source 
splitMeshRegions: handle flipping of faces for surface fields

subsetMesh: subset dimensionedFields

decomposePar: use run-time selection of decomposition constraints. Used to
    keep cells on particular processors. See the decomposeParDict in

$FOAM_UTILITIES/parallel/decomposePar:
  - preserveBaffles: keep baffle faces on same processor
  - preserveFaceZones: keep faceZones owner and neighbour on same processor
  - preservePatches: keep owner and neighbour on same processor. Note: not
    suitable for cyclicAMI since these are not coupled on the patch level
  - singleProcessorFaceSets: keep complete faceSet on a single processor
  - refinementHistory: keep cells originating from a single cell on the
    same processor.

decomposePar: clean up decomposition of refinement data from snappyHexMesh

reconstructPar: reconstruct refinement data (refineHexMesh, snappyHexMesh)

reconstructParMesh: reconstruct refinement data (refineHexMesh, snappyHexMesh)

redistributePar:
  - corrected mapping surfaceFields
  - adding processor patches in order consistent with decomposePar

argList: check that slaves are running same version as master

fvMeshSubset: move to dynamicMesh library

fvMeshDistribute:
  - support for mapping dimensionedFields
  - corrected mapping of surfaceFields

parallel routines: allow parallel running on single processor

Field: support for
  - distributed mapping
  - mapping with flipping

mapDistribute: support for flipping

AMIInterpolation: avoid constructing localPoints
This commit is contained in:
Henry Weller
2016-05-15 16:36:48 +01:00
parent 26658647fa
commit 1441f8cab0
93 changed files with 9230 additions and 2650 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/dynamicMesh/Make/files
View File

@ -28,17 +28,22 @@ polyTopoChange/polyTopoChange/addPatchCellLayer.C
polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
polyTopoChange/polyTopoChange/edgeCollapser.C
polyTopoChange/polyTopoChange/faceCollapser.C
polyTopoChange/polyTopoChange/hexRef8.C
polyTopoChange/polyTopoChange/removeCells.C
polyTopoChange/polyTopoChange/removeFaces.C
polyTopoChange/polyTopoChange/refinementData.C
polyTopoChange/polyTopoChange/refinementDistanceData.C
polyTopoChange/polyTopoChange/refinementHistory.C
polyTopoChange/polyTopoChange/removePoints.C
polyTopoChange/polyTopoChange/combineFaces.C
polyTopoChange/polyTopoChange/duplicatePoints.C
polyTopoChange/polyTopoChange/tetDecomposer.C
hexRef8 = polyTopoChange/polyTopoChange/hexRef8
$(hexRef8)/hexRef8.C
$(hexRef8)/hexRef8Data.C
$(hexRef8)/refinementHistory.C
slidingInterface/slidingInterface.C
slidingInterface/slidingInterfaceProjectPoints.C
slidingInterface/coupleSlidingInterface.C
@ -83,6 +88,8 @@ polyMeshAdder/polyMeshAdder.C
fvMeshTools/fvMeshTools.C
fvMeshSubset/fvMeshSubset.C
motionSmoother/motionSmoother.C
motionSmoother/motionSmootherAlgo.C
motionSmoother/motionSmootherAlgoCheck.C

14
src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
View File

@ -28,6 +28,14 @@ License
#include "faceCoupleInfo.H"
#include "fvMesh.H"
/* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */
namespace Foam
{
defineTypeNameAndDebug(fvMeshAdder, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::labelList Foam::fvMeshAdder::calcPatchMap
@ -105,6 +113,12 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::fvMeshAdder::add
fvMeshAdder::MapSurfaceFields<symmTensor>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapSurfaceFields<tensor>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapDimFields<scalar>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapDimFields<vector>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapDimFields<sphericalTensor>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapDimFields<symmTensor>(mapPtr, mesh0, mesh1);
fvMeshAdder::MapDimFields<tensor>(mapPtr, mesh0, mesh1);
return mapPtr;
}

26
src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
View File

@ -26,7 +26,7 @@ Class
Description
Adds two fvMeshes without using any polyMesh morphing.
Uses fvMeshAdder.
Uses polyMeshAdder.
SourceFiles
fvMeshAdder.C
@ -42,6 +42,7 @@ SourceFiles
#include "fvPatchFieldsFwd.H"
#include "fvsPatchFieldsFwd.H"
#include "fvPatchFieldMapper.H"
#include "DimensionedField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -101,8 +102,22 @@ private:
const GeometricField<Type, fvsPatchField, surfaceMesh>& fldToAdd
);
//- Update single dimensionedField.
template<class Type>
static void MapDimField
(
const mapAddedPolyMesh& meshMap,
DimensionedField<Type, volMesh>& fld,
const DimensionedField<Type, volMesh>& fldToAdd
);
public:
// Declare name of the class and its debug switch
ClassName("fvMeshAdder");
// Member Functions
//- Inplace add mesh to fvMesh. Maps all stored fields. Returns map.
@ -131,6 +146,15 @@ public:
const fvMesh& mesh,
const fvMesh& meshToAdd
);
//- Map all DimensionedFields of Type
template<class Type>
static void MapDimFields
(
const mapAddedPolyMesh&,
const fvMesh& mesh,
const fvMesh& meshToAdd
);
};

103
src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
View File

@ -280,6 +280,12 @@ void Foam::fvMeshAdder::MapVolFields
++fieldIter
)
{
if (debug)
{
Pout<< "MapVolFields : Storing old time for " << fieldIter()->name()
<< endl;
}
const_cast<GeometricField<Type, fvPatchField, volMesh>*>(fieldIter())
->storeOldTimes();
}
@ -304,6 +310,12 @@ void Foam::fvMeshAdder::MapVolFields
const GeometricField<Type, fvPatchField, volMesh>& fldToAdd =
*fieldsToAdd[fld.name()];
if (debug)
{
Pout<< "MapVolFields : mapping " << fld.name()
<< " and " << fldToAdd.name() << endl;
}
MapVolField<Type>(meshMap, fld, fldToAdd);
}
else
@ -585,8 +597,13 @@ void Foam::fvMeshAdder::MapSurfaceFields
++fieldIter
)
{
const_cast<fldType*>(fieldIter())
->storeOldTimes();
if (debug)
{
Pout<< "MapSurfaceFields : Storing old time for "
<< fieldIter()->name() << endl;
}
const_cast<fldType*>(fieldIter())->storeOldTimes();
}
@ -604,6 +621,12 @@ void Foam::fvMeshAdder::MapSurfaceFields
{
const fldType& fldToAdd = *fieldsToAdd[fld.name()];
if (debug)
{
Pout<< "MapSurfaceFields : mapping " << fld.name()
<< " and " << fldToAdd.name() << endl;
}
MapSurfaceField<Type>(meshMap, fld, fldToAdd);
}
else
@ -617,4 +640,80 @@ void Foam::fvMeshAdder::MapSurfaceFields
}
template<class Type>
void Foam::fvMeshAdder::MapDimField
(
const mapAddedPolyMesh& meshMap,
DimensionedField<Type, volMesh>& fld,
const DimensionedField<Type, volMesh>& fldToAdd
)
{
const fvMesh& mesh = fld.mesh();
// Store old field
Field<Type> oldField(fld);
fld.setSize(mesh.nCells());
fld.rmap(oldField, meshMap.oldCellMap());
fld.rmap(fldToAdd, meshMap.addedCellMap());
}
template<class Type>
void Foam::fvMeshAdder::MapDimFields
(
const mapAddedPolyMesh& meshMap,
const fvMesh& mesh,
const fvMesh& meshToAdd
)
{
typedef DimensionedField<Type, volMesh> fldType;
// Note: use strict flag on lookupClass to avoid picking up
// volFields
HashTable<const fldType*> fields
(
mesh.objectRegistry::lookupClass<fldType>(true)
);
HashTable<const fldType*> fieldsToAdd
(
meshToAdd.objectRegistry::lookupClass<fldType>(true)
);
for
(
typename HashTable<const fldType*>::
iterator fieldIter = fields.begin();
fieldIter != fields.end();
++fieldIter
)
{
fldType& fld = const_cast<fldType&>(*fieldIter());
if (fieldsToAdd.found(fld.name()))
{
const fldType& fldToAdd = *fieldsToAdd[fld.name()];
if (debug)
{
Pout<< "MapDimFields : mapping " << fld.name()
<< " and " << fldToAdd.name() << endl;
}
MapDimField<Type>(meshMap, fld, fldToAdd);
}
else
{
WarningIn("fvMeshAdder::MapDimFields(..)")
<< "Not mapping field " << fld.name()
<< " since not present on mesh to add"
<< endl;
}
}
}
// ************************************************************************* //

414
src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
View File

@ -40,17 +40,125 @@ License
#include "syncTools.H"
#include "CompactListList.H"
#include "fvMeshTools.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(fvMeshDistribute, 0);
//- Less function class that can be used for sorting processor patches
class lessProcPatches
{
const labelList& nbrProc_;
const labelList& referPatchID_;
public:
lessProcPatches( const labelList& nbrProc, const labelList& referPatchID)
:
nbrProc_(nbrProc),
referPatchID_(referPatchID)
{}
bool operator()(const label a, const label b)
{
if (nbrProc_[a] < nbrProc_[b])
{
return true;
}
else if (nbrProc_[a] > nbrProc_[b])
{
return false;
}
else
{
// Equal neighbour processor
return referPatchID_[a] < referPatchID_[b];
}
}
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::fvMeshDistribute::inplaceRenumberWithFlip
(
const labelUList& oldToNew,
const bool oldToNewHasFlip,
const bool lstHasFlip,
labelUList& lst
)
{
if (!lstHasFlip && !oldToNewHasFlip)
{
Foam::inplaceRenumber(oldToNew, lst);
}
else
{
// Either input data or map encodes sign so result encodes sign
forAll(lst, elemI)
{
// Extract old value and sign
label val = lst[elemI];
label sign = 1;
if (lstHasFlip)
{
if (val > 0)
{
val = val-1;
}
else if (val < 0)
{
val = -val-1;
sign = -1;
}
else
{
FatalErrorInFunction
<< "Problem : zero value " << val
<< " at index " << elemI << " out of " << lst.size()
<< " list with flip bit" << exit(FatalError);
}
}
// Lookup new value and possibly change sign
label newVal = oldToNew[val];
if (oldToNewHasFlip)
{
if (newVal > 0)
{
newVal = newVal-1;
}
else if (newVal < 0)
{
newVal = -newVal-1;
sign = -sign;
}
else
{
FatalErrorInFunction
<< "Problem : zero value " << newVal
<< " at index " << elemI << " out of "
<< oldToNew.size()
<< " list with flip bit" << exit(FatalError);
}
}
// Encode new value and sign
lst[elemI] = sign*(newVal+1);
}
}
}
Foam::labelList Foam::fvMeshDistribute::select
(
const bool selectEqual,
@ -453,7 +561,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
forAll(constructFaceMap, proci)
{
inplaceRenumber(map().reverseFaceMap(), constructFaceMap[proci]);
inplaceRenumberWithFlip
(
map().reverseFaceMap(),
false,
true,
constructFaceMap[proci]
);
}
@ -882,12 +996,45 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
meshMod
);
//// Generate test field
//tmp<surfaceScalarField> sfld(generateTestField(mesh_));
// Save internal fields (note: not as DimensionedFields since would
// get mapped)
PtrList<Field<scalar>> sFlds;
saveInternalFields(sFlds);
PtrList<Field<vector>> vFlds;
saveInternalFields(vFlds);
PtrList<Field<sphericalTensor>> sptFlds;
saveInternalFields(sptFlds);
PtrList<Field<symmTensor>> sytFlds;
saveInternalFields(sytFlds);
PtrList<Field<tensor>> tFlds;
saveInternalFields(tFlds);
// Change the mesh. No inflation. Note: no parallel comms allowed.
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
// Update fields
mesh_.updateMesh(map);
// Any exposed faces in a surfaceField will not be mapped. Map the value
// of these separately (until there is support in all PatchFields for
// mapping from internal faces ...)
mapExposedFaces(map(), sFlds);
mapExposedFaces(map(), vFlds);
mapExposedFaces(map(), sptFlds);
mapExposedFaces(map(), sytFlds);
mapExposedFaces(map(), tFlds);
//// Test test field
//testField(sfld);
// Move mesh (since morphing does not do this)
if (map().hasMotionPoints())
{
@ -911,10 +1058,18 @@ void Foam::fvMeshDistribute::addProcPatches
// contain for all current boundary faces the global patchID (for non-proc
// patch) or the processor.
// Determine a visit order such that the processor patches get added
// in order of increasing neighbour processor (and for same neighbour
// processor (in case of processor cyclics) in order of increasing
// 'refer' patch)
labelList indices;
sortedOrder(nbrProc, indices, lessProcPatches(nbrProc, referPatchID));
procPatchID.setSize(Pstream::nProcs());
forAll(nbrProc, bFacei)
forAll(indices, i)
{
label bFacei = indices[i];
label proci = nbrProc[bFacei];
if (proci != -1 && proci != Pstream::myProcNo())
@ -927,6 +1082,7 @@ void Foam::fvMeshDistribute::addProcPatches
if (referPatchID[bFacei] == -1)
{
// Ordinary processor boundary
processorPolyPatch pp
(
0, // size
@ -934,7 +1090,7 @@ void Foam::fvMeshDistribute::addProcPatches
mesh_.boundaryMesh().size(),
mesh_.boundaryMesh(),
Pstream::myProcNo(),
nbrProc[bFacei]
proci
);
procPatchID[proci].insert
@ -957,7 +1113,6 @@ void Foam::fvMeshDistribute::addProcPatches
(
mesh_.boundaryMesh()[referPatchID[bFacei]]
);
processorCyclicPolyPatch pp
(
0, // size
@ -965,7 +1120,7 @@ void Foam::fvMeshDistribute::addProcPatches
mesh_.boundaryMesh().size(),
mesh_.boundaryMesh(),
Pstream::myProcNo(),
nbrProc[bFacei],
proci,
pcPatch.name(),
pcPatch.transform()
);
@ -1500,6 +1655,33 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
checkEqualWordList("surfaceTensorFields", surfTensors);
typedef volScalarField::Internal dimScalType;
const wordList dimScalars(mesh_.names(dimScalType::typeName));
checkEqualWordList("volScalarField::Internal", dimScalars);
typedef volVectorField::Internal dimVecType;
const wordList dimVectors(mesh_.names(dimVecType::typeName));
checkEqualWordList("volVectorField::Internal", dimVectors);
typedef volSphericalTensorField::Internal dimSphereType;
const wordList dimSphereTensors(mesh_.names(dimSphereType::typeName));
checkEqualWordList
(
"volSphericalTensorField::Internal",
dimSphereTensors
);
typedef volSymmTensorField::Internal dimSymmTensorType;
const wordList dimSymmTensors(mesh_.names(dimSymmTensorType::typeName));
checkEqualWordList
(
"volSymmTensorField::Internal",
dimSymmTensors
);
typedef volTensorField::Internal dimTensorType;
const wordList dimTensors(mesh_.names(dimTensorType::typeName));
checkEqualWordList("volTensorField::Internal", dimTensors);
@ -1626,10 +1808,13 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
);
subCellMap[recvProc] = subsetter.cellMap();
subFaceMap[recvProc] = renumber
subFaceMap[recvProc] = subsetter.faceFlipMap();
inplaceRenumberWithFlip
(
repatchFaceMap,
subsetter.faceMap()
false, // oldToNew has flip
true, // subFaceMap has flip
subFaceMap[recvProc]
);
subPointMap[recvProc] = subsetter.pointMap();
subPatchMap[recvProc] = subsetter.patchMap();
@ -1681,6 +1866,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
procSourceNewNbrProc,
str
);
// volFields
sendFields<volScalarField>(recvProc, volScalars, subsetter, str);
sendFields<volVectorField>(recvProc, volVectors, subsetter, str);
sendFields<volSphericalTensorField>
@ -1699,6 +1886,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
);
sendFields<volTensorField>(recvProc, volTensors, subsetter, str);
// surfaceFields
sendFields<surfaceScalarField>
(
recvProc,
@ -1734,6 +1922,43 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
subsetter,
str
);
// dimensionedFields
sendFields<volScalarField::Internal>
(
recvProc,
dimScalars,
subsetter,
str
);
sendFields<volVectorField::Internal>
(
recvProc,
dimVectors,
subsetter,
str
);
sendFields<volSphericalTensorField::Internal>
(
recvProc,
dimSphereTensors,
subsetter,
str
);
sendFields<volSymmTensorField::Internal>
(
recvProc,
dimSymmTensors,
subsetter,
str
);
sendFields<volTensorField::Internal>
(
recvProc,
dimTensors,
subsetter,
str
);
}
}
@ -1771,12 +1996,24 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
repatchFaceMap,
subMap().faceMap()
);
// Insert the sign bit from face flipping
labelList& faceMap = subFaceMap[Pstream::myProcNo()];
forAll(faceMap, faceI)
{
faceMap[faceI] += 1;
}
const labelHashSet& flip = subMap().flipFaceFlux();
forAllConstIter(labelHashSet, flip, iter)
{
label faceI = iter.key();
faceMap[faceI] = -faceMap[faceI];
}
subPointMap[Pstream::myProcNo()] = subMap().pointMap();
subPatchMap[Pstream::myProcNo()] = identity(patches.size());
// Initialize all addressing into current mesh
constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces());
constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces()) + 1;
constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
@ -1872,17 +2109,26 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
labelList domainSourceNewNbrProc;
autoPtr<fvMesh> domainMeshPtr;
PtrList<volScalarField> vsf;
PtrList<volVectorField> vvf;
PtrList<volSphericalTensorField> vsptf;
PtrList<volSymmTensorField> vsytf;
PtrList<volTensorField> vtf;
PtrList<surfaceScalarField> ssf;
PtrList<surfaceVectorField> svf;
PtrList<surfaceSphericalTensorField> ssptf;
PtrList<surfaceSymmTensorField> ssytf;
PtrList<surfaceTensorField> stf;
PtrList<volScalarField::Internal> dsf;
PtrList<volVectorField::Internal> dvf;
PtrList<volSphericalTensorField::Internal> dstf;
PtrList<volSymmTensorField::Internal> dsytf;
PtrList<volTensorField::Internal> dtf;
// Opposite of sendMesh
{
domainMeshPtr = receiveMesh
@ -1908,6 +2154,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
// of problems reading consecutive fields from single stream.
dictionary fieldDicts(str);
// Vol fields
receiveFields<volScalarField>
(
sendProc,
@ -1949,6 +2196,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
fieldDicts.subDict(volTensorField::typeName)
);
// Surface fields
receiveFields<surfaceScalarField>
(
sendProc,
@ -1989,12 +2237,70 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
stf,
fieldDicts.subDict(surfaceTensorField::typeName)
);
// Dimensioned fields
receiveFields<volScalarField::Internal>
(
sendProc,
dimScalars,
domainMesh,
dsf,
fieldDicts.subDict
(
volScalarField::Internal::typeName
)
);
receiveFields<volVectorField::Internal>
(
sendProc,
dimVectors,
domainMesh,
dvf,
fieldDicts.subDict
(
volVectorField::Internal::typeName
)
);
receiveFields<volSphericalTensorField::Internal>
(
sendProc,
dimSphereTensors,
domainMesh,
dstf,
fieldDicts.subDict
(
volSphericalTensorField::Internal::
typeName
)
);
receiveFields<volSymmTensorField::Internal>
(
sendProc,
dimSymmTensors,
domainMesh,
dsytf,
fieldDicts.subDict
(
volSymmTensorField::Internal::typeName
)
);
receiveFields<volTensorField::Internal>
(
sendProc,
dimTensors,
domainMesh,
dtf,
fieldDicts.subDict
(
volTensorField::Internal::typeName
)
);
}
const fvMesh& domainMesh = domainMeshPtr();
constructCellMap[sendProc] = identity(domainMesh.nCells());
constructFaceMap[sendProc] = identity(domainMesh.nFaces());
constructFaceMap[sendProc] = identity(domainMesh.nFaces()) + 1;
constructPointMap[sendProc] = identity(domainMesh.nPoints());
constructPatchMap[sendProc] =
identity(domainMesh.boundaryMesh().size());
@ -2105,28 +2411,76 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
domainSourceNewNbrProc
);
// Update all addressing so xxProcAddressing points to correct item
// in masterMesh.
// Update all addressing so xxProcAddressing points to correct
// item in masterMesh.
const labelList& oldCellMap = map().oldCellMap();
const labelList& oldFaceMap = map().oldFaceMap();
const labelList& oldPointMap = map().oldPointMap();
const labelList& oldPatchMap = map().oldPatchMap();
//Note: old mesh faces never flipped!
forAll(constructPatchMap, proci)
{
if (proci != sendProc && constructPatchMap[proci].size())
{
// Processor already in mesh (either myProcNo or received)
inplaceRenumber(oldCellMap, constructCellMap[proci]);
inplaceRenumber(oldFaceMap, constructFaceMap[proci]);
inplaceRenumberWithFlip
(
oldFaceMap,
false,
true,
constructFaceMap[proci]
);
inplaceRenumber(oldPointMap, constructPointMap[proci]);
inplaceRenumber(oldPatchMap, constructPatchMap[proci]);
}
}
labelHashSet flippedAddedFaces;
{
// Find out if any faces of domain mesh were flipped (boundary
// faces becoming internal)
label nBnd = domainMesh.nFaces()-domainMesh.nInternalFaces();
flippedAddedFaces.resize(nBnd/4);
for
(
label domainFaceI = domainMesh.nInternalFaces();
domainFaceI < domainMesh.nFaces();
domainFaceI++
)
{
label newFaceI = map().addedFaceMap()[domainFaceI];
label newCellI = mesh_.faceOwner()[newFaceI];
label domainCellI = domainMesh.faceOwner()[domainFaceI];
if (newCellI != map().addedCellMap()[domainCellI])
{
flippedAddedFaces.insert(domainFaceI);
}
}
}
// Added processor
inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
inplaceRenumber(map().addedFaceMap(), constructFaceMap[sendProc]);
// Add flip
forAllConstIter(labelHashSet, flippedAddedFaces, iter)
{
label domainFaceI = iter.key();
label& val = constructFaceMap[sendProc][domainFaceI];
val = -val;
}
inplaceRenumberWithFlip
(
map().addedFaceMap(),
false,
true, // constructFaceMap has flip sign
constructFaceMap[sendProc]
);
inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
@ -2236,35 +2590,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
Zero
);
initPatchFields<surfaceScalarField, processorFvsPatchField<scalar>>
(
Zero
);
initPatchFields<surfaceVectorField, processorFvsPatchField<vector>>
(
Zero
);
initPatchFields
<
surfaceSphericalTensorField,
processorFvsPatchField<sphericalTensor>
>
(
Zero
);
initPatchFields
<
surfaceSymmTensorField,
processorFvsPatchField<symmTensor>
>
(
Zero
);
initPatchFields<surfaceTensorField, processorFvsPatchField<tensor>>
(
Zero
);
mesh_.setInstance(mesh_.time().timeName());
@ -2308,7 +2633,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
constructPointMap.xfer(),
constructFaceMap.xfer(),
constructCellMap.xfer(),
constructPatchMap.xfer()
constructPatchMap.xfer(),
true, // subFaceMap has flip
true // constructFaceMap has flip
)
);
}

23
src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
View File

@ -82,6 +82,14 @@ class fvMeshDistribute
// Private Member Functions
static void inplaceRenumberWithFlip
(
const labelUList& oldToNew,
const bool oldToNewHasFlip,
const bool lstHasFlip,
labelUList& lst
);
//- Find indices with value
static labelList select
(
@ -117,6 +125,18 @@ class fvMeshDistribute
const PtrList<FieldField<fvsPatchField, T>>& oldBflds
);
//- Save internal fields of surfaceFields
template<class T>
void saveInternalFields(PtrList<Field<T>>& iflds) const;
//- Set value of patch faces resulting from internal faces
template<class T>
void mapExposedFaces
(
const mapPolyMesh& map,
const PtrList<Field<T>>& oldFlds
);
//- Init patch fields of certain type
template<class GeoField, class PatchFieldType>
void initPatchFields
@ -151,6 +171,7 @@ class fvMeshDistribute
labelListList& constructPointMap
);
// Coupling information
//- Construct the local environment of all boundary faces.
@ -240,7 +261,7 @@ class fvMeshDistribute
(
const labelList& neighbourNewProc, // new processor per b. face
const labelList& referPatchID, // -1 or original patch
const List<Map<label>>& procPatchID// patchID
const List<Map<label>>& procPatchID // patchID
);
//- Send mesh and coupling data.

142
src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C
View File

@ -55,13 +55,14 @@ void Foam::fvMeshDistribute::printFieldInfo(const fvMesh& mesh)
}
// Save whole boundary field
template<class T, class Mesh>
void Foam::fvMeshDistribute::saveBoundaryFields
(
PtrList<FieldField<fvsPatchField, T>>& bflds
) const
{
// Save whole boundary field
typedef GeometricField<T, fvsPatchField, Mesh> fldType;
HashTable<const fldType*> flds
@ -84,7 +85,6 @@ void Foam::fvMeshDistribute::saveBoundaryFields
}
// Map boundary field
template<class T, class Mesh>
void Foam::fvMeshDistribute::mapBoundaryFields
(
@ -92,6 +92,8 @@ void Foam::fvMeshDistribute::mapBoundaryFields
const PtrList<FieldField<fvsPatchField, T>>& oldBflds
)
{
// Map boundary field
const labelList& oldPatchStarts = map.oldPatchStarts();
const labelList& faceMap = map.faceMap();
@ -145,13 +147,103 @@ void Foam::fvMeshDistribute::mapBoundaryFields
}
// Init patch fields of certain type
template<class T>
void Foam::fvMeshDistribute::saveInternalFields
(
PtrList<Field<T> >& iflds
) const
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> fldType;
HashTable<const fldType*> flds
(
static_cast<const fvMesh&>(mesh_).objectRegistry::lookupClass<fldType>()
);
iflds.setSize(flds.size());
label i = 0;
forAllConstIter(typename HashTable<const fldType*>, flds, iter)
{
const fldType& fld = *iter();
iflds.set(i, fld.primitiveField().clone());
i++;
}
}
template<class T>
void Foam::fvMeshDistribute::mapExposedFaces
(
const mapPolyMesh& map,
const PtrList<Field<T> >& oldFlds
)
{
// Set boundary values of exposed internal faces
const labelList& faceMap = map.faceMap();
typedef GeometricField<T, fvsPatchField, surfaceMesh> fldType;
HashTable<fldType*> flds
(
mesh_.objectRegistry::lookupClass<fldType>()
);
if (flds.size() != oldFlds.size())
{
FatalErrorIn("fvMeshDistribute::mapExposedFaces(..)") << "problem"
<< abort(FatalError);
}
label fieldI = 0;
forAllIter(typename HashTable<fldType*>, flds, iter)
{
fldType& fld = *iter();
typename fldType::Boundary& bfld = fld.boundaryFieldRef();
const Field<T>& oldInternal = oldFlds[fieldI++];
// Pull from old internal field into bfld.
forAll(bfld, patchI)
{
fvsPatchField<T>& patchFld = bfld[patchI];
forAll(patchFld, i)
{
const label faceI = patchFld.patch().start()+i;
label oldFaceI = faceMap[faceI];
if (oldFaceI < oldInternal.size())
{
patchFld[i] = oldInternal[oldFaceI];
if (map.flipFaceFlux().found(faceI))
{
patchFld[i] = flipOp()(patchFld[i]);
}
}
}
}
}
}
template<class GeoField, class PatchFieldType>
void Foam::fvMeshDistribute::initPatchFields
(
const typename GeoField::value_type& initVal
)
{
// Init patch fields of certain type
HashTable<GeoField*> flds
(
mesh_.objectRegistry::lookupClass<GeoField>()
@ -161,8 +253,7 @@ void Foam::fvMeshDistribute::initPatchFields
{
GeoField& fld = *iter();
typename GeoField::Boundary& bfld =
fld.boundaryFieldRef();
typename GeoField::Boundary& bfld = fld.boundaryFieldRef();
forAll(bfld, patchi)
{
@ -175,10 +266,11 @@ void Foam::fvMeshDistribute::initPatchFields
}
// correctBoundaryConditions patch fields of certain type
template<class GeoField>
void Foam::fvMeshDistribute::correctBoundaryConditions()
{
// correctBoundaryConditions patch fields of certain type
HashTable<GeoField*> flds
(
mesh_.objectRegistry::lookupClass<GeoField>()
@ -192,24 +284,6 @@ void Foam::fvMeshDistribute::correctBoundaryConditions()
}
// Send fields. Note order supplied so we can receive in exactly the same order.
// Note that field gets written as entry in dictionary so we
// can construct from subdictionary.
// (since otherwise the reading as-a-dictionary mixes up entries from
// consecutive fields)
// The dictionary constructed is:
// volScalarField
// {
// p {internalField ..; boundaryField ..;}
// k {internalField ..; boundaryField ..;}
// }
// volVectorField
// {
// U {internalField ... }
// }
// volVectorField {U {internalField ..; boundaryField ..;}}
//
template<class GeoField>
void Foam::fvMeshDistribute::sendFields
(
@ -219,6 +293,25 @@ void Foam::fvMeshDistribute::sendFields
Ostream& toNbr
)
{
// Send fields. Note order supplied so we can receive in exactly the same
// order.
// Note that field gets written as entry in dictionary so we
// can construct from subdictionary.
// (since otherwise the reading as-a-dictionary mixes up entries from
// consecutive fields)
// The dictionary constructed is:
// volScalarField
// {
// p {internalField ..; boundaryField ..;}
// k {internalField ..; boundaryField ..;}
// }
// volVectorField
// {
// U {internalField ... }
// }
// volVectorField {U {internalField ..; boundaryField ..;}}
toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
@ -244,7 +337,6 @@ void Foam::fvMeshDistribute::sendFields
}
// Opposite of sendFields
template<class GeoField>
void Foam::fvMeshDistribute::receiveFields
(

1562
src/dynamicMesh/fvMeshSubset/fvMeshSubset.C Normal file
View File

File diff suppressed because it is too large Load Diff

338
src/dynamicMesh/fvMeshSubset/fvMeshSubset.H Normal file
View File

@ -0,0 +1,338 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fvMeshSubset
Description
Post-processing mesh subset tool. Given the original mesh and the
list of selected cells, it creates the mesh consisting only of the
desired cells, with the mapping list for points, faces, and cells.
Puts all exposed internal faces into either
- a user supplied patch
- a newly created patch "oldInternalFaces"
- setCellSubset is for small subsets. Uses Maps to minimize memory.
- setLargeCellSubset is for largish subsets (>10% of mesh).
Uses labelLists instead.
- setLargeCellSubset does coupled patch subsetting as well. If it detects
a face on a coupled patch 'losing' its neighbour it will move the
face into the oldInternalFaces patch.
- if a user supplied patch is used it is up to the destination
patchField to handle exposed internal faces (mapping from face -1).
If not provided the default is to assign the internalField. All the
basic patch field types (e.g. fixedValue) will give a warning and
preferably derived patch field types should be used that know how to
handle exposed faces (e.g. use uniformFixedValue instead of fixedValue)
SourceFiles
fvMeshSubset.C
\*---------------------------------------------------------------------------*/
#ifndef fvMeshSubset_H
#define fvMeshSubset_H
#include "fvMesh.H"
#include "pointMesh.H"
#include "GeometricField.H"
#include "HashSet.H"
#include "surfaceMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class fvMeshSubset Declaration
\*---------------------------------------------------------------------------*/
class fvMeshSubset
{
private:
// Private data
//- Mesh to subset from
const fvMesh& baseMesh_;
//- Subset mesh pointer
autoPtr<fvMesh> fvMeshSubsetPtr_;
//- Point mapping array
labelList pointMap_;
//- Face mapping array
labelList faceMap_;
//- Cell mapping array
labelList cellMap_;
//- Patch mapping array
labelList patchMap_;
//- Optional face mapping array with flip encoded
mutable autoPtr<labelList> faceFlipMapPtr_;
// Private Member Functions
//- Check if subset has been performed
bool checkCellSubset() const;
//- Mark points in Map
static void markPoints(const labelList&, Map<label>&);
//- Mark points (with 0) in labelList
static void markPoints(const labelList&, labelList&);
//- Adapt nCellsUsingFace for coupled faces becoming 'uncoupled'.
void doCoupledPatches
(
const bool syncPar,
labelList& nCellsUsingFace
) const;
//- Subset of subset
static labelList subset
(
const label nElems,
const labelList& selectedElements,
const labelList& subsetMap
);
//- Create zones for submesh
void subsetZones();
//- Helper: extract cells-to-remove from cells-to-keep
labelList getCellsToRemove
(
const labelList& region,
const label currentRegion
) const;
//- Disallow default bitwise copy construct
fvMeshSubset(const fvMeshSubset&);
//- Disallow default bitwise assignment
void operator=(const fvMeshSubset&);
public:
// Constructors
//- Construct given a mesh to subset
explicit fvMeshSubset(const fvMesh&);
// Member Functions
// Edit
//- Set the subset. Create "oldInternalFaces" patch for exposed
// internal faces (patchID==-1) or use supplied patch.
// Does not handle coupled patches correctly if only one side
// gets deleted.
void setCellSubset
(
const labelHashSet& globalCellMap,
const label patchID = -1
);
//- Set the subset from all cells with region == currentRegion.
// Create "oldInternalFaces" patch for exposed
// internal faces (patchID==-1) or use supplied patch.
// Handles coupled patches by if necessary making coupled patch
// face part of patchID (so uncoupled)
void setLargeCellSubset
(
const labelList& region,
const label currentRegion,
const label patchID = -1,
const bool syncCouples = true
);
//- setLargeCellSubset but with labelHashSet.
void setLargeCellSubset
(
const labelHashSet& globalCellMap,
const label patchID = -1,
const bool syncPar = true
);
//- Two step subsetting
//- Get labels of exposed faces.
// These are
// - internal faces that become boundary faces
// - coupled faces that become uncoupled (since one of the
// sides gets deleted)
labelList getExposedFaces
(
const labelList& region,
const label currentRegion,
const bool syncCouples = true
) const;
//- For every exposed face (from above getExposedFaces)
// used supplied (existing!) patch
void setLargeCellSubset
(
const labelList& region,
const label currentRegion,
const labelList& exposedFaces,
const labelList& patchIDs,
const bool syncCouples = true
);
// Access
//- Original mesh
const fvMesh& baseMesh() const
{
return baseMesh_;
}
//- Have subMesh?
bool hasSubMesh() const;
//- Return reference to subset mesh
const fvMesh& subMesh() const;
fvMesh& subMesh();
//- Return point map
const labelList& pointMap() const;
//- Return face map
const labelList& faceMap() const;
//- Return face map with sign to encode flipped faces
const labelList& faceFlipMap() const;
//- Return cell map
const labelList& cellMap() const;
//- Return patch map
const labelList& patchMap() const;
// Field mapping
//- Map volume field
template<class Type>
static tmp<GeometricField<Type, fvPatchField, volMesh>>
interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
const fvMesh& sMesh,
const labelList& patchMap,
const labelList& cellMap,
const labelList& faceMap
);
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh>>
interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&
) const;
//- Map surface field. Optionally negates value if flipping
// a face (from exposing an internal face)
template<class Type>
static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
interpolate
(
const GeometricField<Type, fvsPatchField, surfaceMesh>&,
const fvMesh& sMesh,
const labelList& patchMap,
const labelList& cellMap,
const labelList& faceMap,
const bool negateIfFlipped = true
);
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
interpolate
(
const GeometricField<Type, fvsPatchField, surfaceMesh>&,
const bool negateIfFlipped = true
) const;
//- Map point field
template<class Type>
static tmp<GeometricField<Type, pointPatchField, pointMesh>>
interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>&,
const pointMesh& sMesh,
const labelList& patchMap,
const labelList& pointMap
);
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh>>
interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Map dimensioned field
template<class Type>
static tmp<DimensionedField<Type, volMesh>>
interpolate
(
const DimensionedField<Type, volMesh>&,
const fvMesh& sMesh,
const labelList& cellMap
);
template<class Type>
tmp<DimensionedField<Type, volMesh>>
interpolate(const DimensionedField<Type, volMesh>&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "fvMeshSubsetInterpolate.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

553
src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C Normal file
View File

@ -0,0 +1,553 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "fvMeshSubset.H"
#include "emptyFvsPatchField.H"
#include "emptyPointPatchField.H"
#include "emptyFvPatchFields.H"
#include "directFvPatchFieldMapper.H"
#include "directPointPatchFieldMapper.H"
#include "flipOp.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh>> fvMeshSubset::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
const fvMesh& sMesh,
const labelList& patchMap,
const labelList& cellMap,
const labelList& faceMap
)
{
// 1. Create the complete field with dummy patch fields
PtrList<fvPatchField<Type>> patchFields(patchMap.size());
forAll(patchFields, patchi)
{
// Set the first one by hand as it corresponds to the
// exposed internal faces. Additional interpolation can be put here
// as necessary.
if (patchMap[patchi] == -1)
{
patchFields.set
(
patchi,
new emptyFvPatchField<Type>
(
sMesh.boundary()[patchi],
DimensionedField<Type, volMesh>::null()
)
);
}
else
{
patchFields.set
(
patchi,
fvPatchField<Type>::New
(
calculatedFvPatchField<Type>::typeName,
sMesh.boundary()[patchi],
DimensionedField<Type, volMesh>::null()
)
);
}
}
tmp<GeometricField<Type, fvPatchField, volMesh>> tresF
(
new GeometricField<Type, fvPatchField, volMesh>
(
IOobject
(
"subset"+vf.name(),
sMesh.time().timeName(),
sMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
vf.dimensions(),
Field<Type>(vf.primitiveField(), cellMap),
patchFields
)
);
GeometricField<Type, fvPatchField, volMesh>& resF = tresF.ref();
// 2. Change the fvPatchFields to the correct type using a mapper
// constructor (with reference to the now correct internal field)
typename GeometricField<Type, fvPatchField, volMesh>::
Boundary& bf = resF.boundaryFieldRef();
forAll(bf, patchi)
{
if (patchMap[patchi] != -1)
{
// Construct addressing
const fvPatch& subPatch = sMesh.boundary()[patchi];
const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchi]];
const label baseStart = basePatch.start();
const label baseSize = basePatch.size();
labelList directAddressing(subPatch.size());
forAll(directAddressing, i)
{
label baseFacei = faceMap[subPatch.start()+i];
if (baseFacei >= baseStart && baseFacei < baseStart+baseSize)
{
directAddressing[i] = baseFacei-baseStart;
}
else
{
// Mapped from internal face. Do what? Leave up to
// fvPatchField
directAddressing[i] = -1;
}
}
bf.set
(
patchi,
fvPatchField<Type>::New
(
vf.boundaryField()[patchMap[patchi]],
subPatch,
resF(),
directFvPatchFieldMapper(directAddressing)
)
);
}
}
return tresF;
}
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh>> fvMeshSubset::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
return interpolate
(
vf,
subMesh(),
patchMap(),
cellMap(),
faceMap()
);
}
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvMeshSubset::interpolate
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& vf,
const fvMesh& sMesh,
const labelList& patchMap,
const labelList& cellMap,
const labelList& faceMap,
const bool negateIfFlipped
)
{
// 1. Create the complete field with dummy patch fields
PtrList<fvsPatchField<Type>> patchFields(patchMap.size());
forAll(patchFields, patchi)
{
// Set the first one by hand as it corresponds to the
// exposed internal faces. Additional interpolation can be put here
// as necessary.
if (patchMap[patchi] == -1)
{
patchFields.set
(
patchi,
new emptyFvsPatchField<Type>
(
sMesh.boundary()[patchi],
DimensionedField<Type, surfaceMesh>::null()
)
);
}
else
{
patchFields.set
(
patchi,
fvsPatchField<Type>::New
(
calculatedFvsPatchField<Type>::typeName,
sMesh.boundary()[patchi],
DimensionedField<Type, surfaceMesh>::null()
)
);
}
}
// Create the complete field from the pieces
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tresF
(
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
IOobject
(
"subset"+vf.name(),
sMesh.time().timeName(),
sMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
vf.dimensions(),
Field<Type>
(
vf.primitiveField(),
SubList<label>
(
faceMap,
sMesh.nInternalFaces()
)
),
patchFields
)
);
GeometricField<Type, fvsPatchField, surfaceMesh>& resF = tresF.ref();
// 2. Change the fvsPatchFields to the correct type using a mapper
// constructor (with reference to the now correct internal field)
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
Boundary& bf = resF.boundaryFieldRef();
forAll(bf, patchi)
{
if (patchMap[patchi] != -1)
{
// Construct addressing
const fvPatch& subPatch = sMesh.boundary()[patchi];
const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchi]];
const label baseStart = basePatch.start();
const label baseSize = basePatch.size();
labelList directAddressing(subPatch.size());
forAll(directAddressing, i)
{
label baseFacei = faceMap[subPatch.start()+i];
if (baseFacei >= baseStart && baseFacei < baseStart+baseSize)
{
directAddressing[i] = baseFacei-baseStart;
}
else
{
// Mapped from internal face. Do what? Leave up to
// patchField. This would require also to pass in
// original internal field so for now do as postprocessing
directAddressing[i] = -1;
}
}
bf.set
(
patchi,
fvsPatchField<Type>::New
(
vf.boundaryField()[patchMap[patchi]],
subPatch,
resF(),
directFvPatchFieldMapper(directAddressing)
)
);
// Postprocess patch field for exposed faces
fvsPatchField<Type>& pfld = bf[patchi];
const labelUList& fc = bf[patchi].patch().faceCells();
const labelList& own = vf.mesh().faceOwner();
forAll(pfld, i)
{
label baseFacei = faceMap[subPatch.start()+i];
if (baseFacei < vf.primitiveField().size())
{
Type val = vf.internalField()[baseFacei];
if (cellMap[fc[i]] == own[baseFacei] || !negateIfFlipped)
{
pfld[i] = val;
}
else
{
pfld[i] = flipOp()(val);
}
}
else
{
// Exposed face from other patch.
// Only possible in case of a coupled boundary
label patchi = vf.mesh().boundaryMesh().whichPatch
(
baseFacei
);
const fvPatch& otherPatch = vf.mesh().boundary()[patchi];
label patchFacei = otherPatch.patch().whichFace(baseFacei);
pfld[i] = vf.boundaryField()[patchi][patchFacei];
}
}
}
}
return tresF;
}
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvMeshSubset::interpolate
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& sf,
const bool negateIfFlipped
) const
{
return interpolate
(
sf,
subMesh(),
patchMap(),
cellMap(),
faceMap(),
negateIfFlipped
);
}
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh>>
fvMeshSubset::interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>& vf,
const pointMesh& sMesh,
const labelList& patchMap,
const labelList& pointMap
)
{
// 1. Create the complete field with dummy patch fields
PtrList<pointPatchField<Type>> patchFields(patchMap.size());
forAll(patchFields, patchi)
{
// Set the first one by hand as it corresponds to the
// exposed internal faces. Additional interpolation can be put here
// as necessary.
if (patchMap[patchi] == -1)
{
patchFields.set
(
patchi,
new emptyPointPatchField<Type>
(
sMesh.boundary()[patchi],
DimensionedField<Type, pointMesh>::null()
)
);
}
else
{
patchFields.set
(
patchi,
pointPatchField<Type>::New
(
calculatedPointPatchField<Type>::typeName,
sMesh.boundary()[patchi],
DimensionedField<Type, pointMesh>::null()
)
);
}
}
// Create the complete field from the pieces
tmp<GeometricField<Type, pointPatchField, pointMesh>> tresF
(
new GeometricField<Type, pointPatchField, pointMesh>
(
IOobject
(
"subset"+vf.name(),
sMesh.time().timeName(),
sMesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
vf.dimensions(),
Field<Type>(vf.primitiveField(), pointMap),
patchFields
)
);
GeometricField<Type, pointPatchField, pointMesh>& resF = tresF.ref();
// 2. Change the pointPatchFields to the correct type using a mapper
// constructor (with reference to the now correct internal field)
typename GeometricField<Type, pointPatchField, pointMesh>::
Boundary& bf = resF.boundaryFieldRef();
forAll(bf, patchi)
{
// Set the first one by hand as it corresponds to the
// exposed internal faces. Additional interpolation can be put here
// as necessary.
if (patchMap[patchi] != -1)
{
// Construct addressing
const pointPatch& basePatch =
vf.mesh().boundary()[patchMap[patchi]];
const labelList& meshPoints = basePatch.meshPoints();
// Make addressing from mesh to patch point
Map<label> meshPointMap(2*meshPoints.size());
forAll(meshPoints, localI)
{
meshPointMap.insert(meshPoints[localI], localI);
}
// Find which subpatch points originate from which patch point
const pointPatch& subPatch = sMesh.boundary()[patchi];
const labelList& subMeshPoints = subPatch.meshPoints();
// If mapped from outside patch leave handling up to patchField
labelList directAddressing(subPatch.size(), -1);
forAll(subMeshPoints, localI)
{
// Get mesh point on original mesh.
label meshPointI = pointMap[subMeshPoints[localI]];
Map<label>::const_iterator iter = meshPointMap.find(meshPointI);
if (iter != meshPointMap.end())
{
directAddressing[localI] = iter();
}
}
bf.set
(
patchi,
pointPatchField<Type>::New
(
vf.boundaryField()[patchMap[patchi]],
subPatch,
resF(),
directPointPatchFieldMapper(directAddressing)
)
);
}
}
return tresF;
}
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh>> fvMeshSubset::interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>& sf
) const
{
return interpolate
(
sf,
pointMesh::New(subMesh()), // subsetted point mesh
patchMap(),
pointMap()
);
}
template<class Type>
tmp<DimensionedField<Type, volMesh>> fvMeshSubset::interpolate
(
const DimensionedField<Type, volMesh>& df,
const fvMesh& sMesh,
const labelList& cellMap
)
{
// Create the complete field from the pieces
tmp<DimensionedField<Type, volMesh>> tresF
(
new DimensionedField<Type, volMesh>
(
IOobject
(
"subset"+df.name(),
sMesh.time().timeName(),
sMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
df.dimensions(),
Field<Type>(df, cellMap)
)
);
return tresF;
}
template<class Type>
tmp<DimensionedField<Type, volMesh>> fvMeshSubset::interpolate
(
const DimensionedField<Type, volMesh>& df
) const
{
return interpolate(df, subMesh(), cellMap());
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

3
src/dynamicMesh/meshCut/meshModifiers/multiDirRefinement/multiDirRefinement.C
View File

@ -256,7 +256,8 @@ void Foam::multiDirRefinement::refineHex8
false
),
List<refinementHistory::splitCell8>(0),
labelList(0)
labelList(0),
false
) // refinement history
);

7
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C → src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
View File

@ -793,7 +793,7 @@ Foam::label Foam::hexRef8::findLevel
// Gets cell level such that the face has four points <= level.
Foam::label Foam::hexRef8::getAnchorLevel(const label facei) const
Foam::label Foam::hexRef8::faceLevel(const label facei) const
{
const face& f = mesh_.faces()[facei];
@ -2218,7 +2218,8 @@ Foam::hexRef8::hexRef8
IOobject::AUTO_WRITE
),
List<refinementHistory::splitCell8>(0),
labelList(0)
labelList(0),
false
),
faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible
savedPointLevel_(0),
@ -3475,7 +3476,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
for (label facei = 0; facei < mesh_.nFaces(); facei++)
{
faceAnchorLevel[facei] = getAnchorLevel(facei);
faceAnchorLevel[facei] = faceLevel(facei);
}
// -1 : no need to split face

2
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H → src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H
View File

@ -411,7 +411,7 @@ public:
// Refinement
//- Gets level such that the face has four points <= level.
label getAnchorLevel(const label facei) const;
label faceLevel(const label facei) const;
//- Given valid mesh and current cell level and proposed
// cells to refine calculate any clashes (due to 2:1) and return

339
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8Data.C Normal file
View File

@ -0,0 +1,339 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "IOobject.H"
#include "UList.H"
#include "hexRef8Data.H"
#include "mapPolyMesh.H"
#include "mapDistributePolyMesh.H"
#include "polyMesh.H"
#include "syncTools.H"
#include "refinementHistory.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::hexRef8Data::hexRef8Data(const IOobject& io)
{
{
IOobject rio(io);
rio.rename("cellLevel");
bool haveFile = returnReduce(rio.headerOk(), orOp<bool>());
if (haveFile)
{
Info<< "Reading hexRef8 data : " << rio.name() << endl;
cellLevelPtr_.reset(new labelIOList(rio));
}
}
{
IOobject rio(io);
rio.rename("pointLevel");
bool haveFile = returnReduce(rio.headerOk(), orOp<bool>());
if (haveFile)
{
Info<< "Reading hexRef8 data : " << rio.name() << endl;
pointLevelPtr_.reset(new labelIOList(rio));
}
}
{
IOobject rio(io);
rio.rename("level0Edge");
bool haveFile = returnReduce(rio.headerOk(), orOp<bool>());
if (haveFile)
{
Info<< "Reading hexRef8 data : " << rio.name() << endl;
level0EdgePtr_.reset(new uniformDimensionedScalarField(rio));
}
}
{
IOobject rio(io);
rio.rename("refinementHistory");
bool haveFile = returnReduce(rio.headerOk(), orOp<bool>());
if (haveFile)
{
Info<< "Reading hexRef8 data : " << rio.name() << endl;
refHistoryPtr_.reset(new refinementHistory(rio));
}
}
}
Foam::hexRef8Data::hexRef8Data
(
const IOobject& io,
const hexRef8Data& data,
const labelList& cellMap,
const labelList& pointMap
)
{
if (data.cellLevelPtr_.valid())
{
IOobject rio(io);
rio.rename(data.cellLevelPtr_().name());
cellLevelPtr_.reset
(
new labelIOList
(
rio,
UIndirectList<label>(data.cellLevelPtr_(), cellMap)()
)
);
}
if (data.pointLevelPtr_.valid())
{
IOobject rio(io);
rio.rename(data.pointLevelPtr_().name());
pointLevelPtr_.reset
(
new labelIOList
(
rio,
UIndirectList<label>(data.pointLevelPtr_(), pointMap)()
)
);
}
if (data.level0EdgePtr_.valid())
{
IOobject rio(io);
rio.rename(data.level0EdgePtr_().name());
level0EdgePtr_.reset
(
new uniformDimensionedScalarField(rio, data.level0EdgePtr_())
);
}
if (data.refHistoryPtr_.valid())
{
IOobject rio(io);
rio.rename(data.refHistoryPtr_().name());
refHistoryPtr_ = data.refHistoryPtr_().clone(rio, cellMap);
}
}
Foam::hexRef8Data::hexRef8Data
(
const IOobject& io,
const UPtrList<const labelList>& cellMaps,
const UPtrList<const labelList>& pointMaps,
const UPtrList<const hexRef8Data>& procDatas
)
{
const polyMesh& mesh = dynamic_cast<const polyMesh&>(io.db());
// cellLevel
if (procDatas[0].cellLevelPtr_.valid())
{
IOobject rio(io);
rio.rename(procDatas[0].cellLevelPtr_().name());
cellLevelPtr_.reset(new labelIOList(rio, mesh.nCells()));
labelList& cellLevel = cellLevelPtr_();
forAll(procDatas, procI)
{
const labelList& procCellLevel = procDatas[procI].cellLevelPtr_();
UIndirectList<label>(cellLevel, cellMaps[procI]) = procCellLevel;
}
}
// pointLevel
if (procDatas[0].pointLevelPtr_.valid())
{
IOobject rio(io);
rio.rename(procDatas[0].pointLevelPtr_().name());
pointLevelPtr_.reset(new labelIOList(rio, mesh.nPoints()));
labelList& pointLevel = pointLevelPtr_();
forAll(procDatas, procI)
{
const labelList& procPointLevel = procDatas[procI].pointLevelPtr_();
UIndirectList<label>(pointLevel, pointMaps[procI]) = procPointLevel;
}
}
// level0Edge
if (procDatas[0].level0EdgePtr_.valid())
{
IOobject rio(io);
rio.rename(procDatas[0].level0EdgePtr_().name());
level0EdgePtr_.reset
(
new uniformDimensionedScalarField
(
rio,
procDatas[0].level0EdgePtr_()
)
);
}
// refinementHistory
if (procDatas[0].refHistoryPtr_.valid())
{
IOobject rio(io);
rio.rename(procDatas[0].refHistoryPtr_().name());
UPtrList<const refinementHistory> procRefs(procDatas.size());
forAll(procDatas, i)
{
procRefs.set(i, &procDatas[i].refHistoryPtr_());
}
refHistoryPtr_.reset
(
new refinementHistory
(
rio,
cellMaps,
procRefs
)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::hexRef8Data::~hexRef8Data()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::hexRef8Data::sync(const IOobject& io)
{
const polyMesh& mesh = dynamic_cast<const polyMesh&>(io.db());
bool hasCellLevel = returnReduce(cellLevelPtr_.valid(), orOp<bool>());
if (hasCellLevel && !cellLevelPtr_.valid())
{
IOobject rio(io);
rio.rename("cellLevel");
rio.readOpt() = IOobject::NO_READ;
cellLevelPtr_.reset(new labelIOList(rio, labelList(mesh.nCells(), 0)));
}
bool hasPointLevel = returnReduce(pointLevelPtr_.valid(), orOp<bool>());
if (hasPointLevel && !pointLevelPtr_.valid())
{
IOobject rio(io);
rio.rename("pointLevel");
rio.readOpt() = IOobject::NO_READ;
pointLevelPtr_.reset
(
new labelIOList(rio, labelList(mesh.nPoints(), 0))
);
}
bool hasLevel0Edge = returnReduce(level0EdgePtr_.valid(), orOp<bool>());
if (hasLevel0Edge)
{
// Get master length
scalar masterLen = level0EdgePtr_().value();
Pstream::scatter(masterLen);
if (!level0EdgePtr_.valid())
{
IOobject rio(io);
rio.rename("level0Edge");
rio.readOpt() = IOobject::NO_READ;
level0EdgePtr_.reset
(
new uniformDimensionedScalarField
(
rio,
dimensionedScalar("zero", dimLength, masterLen)
)
);
}
}
bool hasHistory = returnReduce(refHistoryPtr_.valid(), orOp<bool>());
if (hasHistory && !refHistoryPtr_.valid())
{
IOobject rio(io);
rio.rename("refinementHistory");
rio.readOpt() = IOobject::NO_READ;
refHistoryPtr_.reset(new refinementHistory(rio, mesh.nCells(), true));
}
}
void Foam::hexRef8Data::distribute(const mapDistributePolyMesh& map)
{
if (cellLevelPtr_.valid())
{
map.cellMap().distribute(cellLevelPtr_());
}
if (pointLevelPtr_.valid())
{
map.pointMap().distribute(pointLevelPtr_());
}
// No need to distribute the level0Edge
if (refHistoryPtr_.valid() && refHistoryPtr_().active())
{
refHistoryPtr_().distribute(map);
}
}
bool Foam::hexRef8Data::write() const
{
bool ok = true;
if (cellLevelPtr_.valid())
{
ok = ok && cellLevelPtr_().write();
}
if (pointLevelPtr_.valid())
{
ok = ok && pointLevelPtr_().write();
}
if (level0EdgePtr_.valid())
{
ok = ok && level0EdgePtr_().write();
}
if (refHistoryPtr_.valid())
{
ok = ok && refHistoryPtr_().write();
}
return ok;
}
// ************************************************************************* //

136
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8Data.H Normal file
View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::hexRef8Data
Description
Various for reading/decomposing/reconstructing/distributing refinement
data.
SourceFiles
hexRef8Data.C
\*---------------------------------------------------------------------------*/
#ifndef hexRef8Data_H
#define hexRef8Data_H
#include "labelIOList.H"
#include "uniformDimensionedFields.H"
#include "UPtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class mapPolyMesh;
class mapDistributePolyMesh;
class refinementHistory;
class fvMesh;
/*---------------------------------------------------------------------------*\
Class hexRef8Data Declaration
\*---------------------------------------------------------------------------*/
class hexRef8Data
{
private:
// Private data
autoPtr<labelIOList> cellLevelPtr_;
autoPtr<labelIOList> pointLevelPtr_;
autoPtr<uniformDimensionedScalarField> level0EdgePtr_;
autoPtr<refinementHistory> refHistoryPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
hexRef8Data(const hexRef8Data&);
//- Disallow default bitwise assignment
void operator=(const hexRef8Data&);
public:
// Constructors
//- Construct read. Has special provision for only some processors
// having the files so can be used in redistribution.
hexRef8Data(const IOobject& io);
//- Construct as subset
hexRef8Data
(
const IOobject& io,
const hexRef8Data&,
const labelList& cellMap,
const labelList& pointMap
);
//- Construct from multiple hexRef8Data
hexRef8Data
(
const IOobject& io,
const UPtrList<const labelList>& cellMaps,
const UPtrList<const labelList>& pointMaps,
const UPtrList<const hexRef8Data>&
);
//- Destructor
~hexRef8Data();
// Member Functions
//- Parallel synchronise. This enforces valid objects on all processors
// (even if they don't have a mesh). Used by redistributePar.
void sync(const IOobject& io);
//- In-place distribute
void distribute(const mapDistributePolyMesh&);
//- Write
bool write() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

660
src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C → src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/refinementHistory.C
View File

@ -23,12 +23,11 @@ License
\*---------------------------------------------------------------------------*/
#include "DynamicList.H"
#include "refinementHistory.H"
#include "ListOps.H"
#include "mapPolyMesh.H"
#include "mapDistributePolyMesh.H"
#include "polyMesh.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -138,7 +137,59 @@ Foam::refinementHistory::splitCell8::splitCell8(const splitCell8& sc)
{}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
void Foam::refinementHistory::splitCell8::operator=(const splitCell8& s)
{
//- Assignment operator since autoPtr otherwise 'steals' storage.
// Check for assignment to self
if (this == &s)
{
FatalErrorIn("splitCell8::operator=(const Foam::splitCell8&)")
<< "Attempted assignment to self"
<< abort(FatalError);
}
parent_ = s.parent_;
addedCellsPtr_.reset
(
s.addedCellsPtr_.valid()
? new FixedList<label, 8>(s.addedCellsPtr_())
: NULL
);
}
bool Foam::refinementHistory::splitCell8::operator==(const splitCell8& s) const
{
if (addedCellsPtr_.valid() != s.addedCellsPtr_.valid())
{
return false;
}
else if (parent_ != s.parent_)
{
return false;
}
else if (addedCellsPtr_.valid())
{
return addedCellsPtr_() == s.addedCellsPtr_();
}
else
{
return true;
}
}
bool Foam::refinementHistory::splitCell8::operator!=(const splitCell8& s) const
{
return !operator==(s);
}
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>(Istream& is, refinementHistory::splitCell8& sc)
{
@ -183,6 +234,8 @@ Foam::Ostream& Foam::operator<<
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::refinementHistory::checkIndices() const
{
// Check indices.
@ -319,11 +372,192 @@ void Foam::refinementHistory::markSplit
}
// Mark index and all its descendants
void Foam::refinementHistory::mark
(
const label val,
const label index,
labelList& splitToVal
) const
{
splitToVal[index] = val;
const splitCell8& split = splitCells_[index];
if (split.addedCellsPtr_.valid())
{
const FixedList<label, 8>& splits = split.addedCellsPtr_();
forAll(splits, i)
{
if (splits[i] >= 0)
{
mark(val, splits[i], splitToVal);
}
}
}
}
Foam::label Foam::refinementHistory::markCommonCells
(
labelList& cellToCluster
) const
{
label clusterI = 0;
labelList splitToCluster(splitCells_.size(), -1);
// Pass1: find top of all clusters
forAll(visibleCells_, cellI)
{
label index = visibleCells_[cellI];
if (index >= 0)
{
// Find highest ancestor
while (splitCells_[index].parent_ != -1)
{
index = splitCells_[index].parent_;
}
// Mark tree with clusterI
if (splitToCluster[index] == -1)
{
mark(clusterI, index, splitToCluster);
clusterI++;
}
}
}
// Pass2: mark all cells with cluster
cellToCluster.setSize(visibleCells_.size(), -1);
forAll(visibleCells_, cellI)
{
label index = visibleCells_[cellI];
if (index >= 0)
{
cellToCluster[cellI] = splitToCluster[index];
}
}
return clusterI;
}
void Foam::refinementHistory::add
(
boolList& blockedFace,
PtrList<labelList>& specifiedProcessorFaces,
labelList& specifiedProcessor,
List<labelPair>& explicitConnections
) const
{
const polyMesh& mesh = dynamic_cast<const polyMesh&>(db());
blockedFace.setSize(mesh.nFaces(), true);
// Find common parent for all cells
labelList cellToCluster;
markCommonCells(cellToCluster);
// Unblock all faces inbetween same cluster
label nUnblocked = 0;
forAll(mesh.faceNeighbour(), faceI)
{
label ownCluster = cellToCluster[mesh.faceOwner()[faceI]];
label neiCluster = cellToCluster[mesh.faceNeighbour()[faceI]];
if (ownCluster != -1 && ownCluster == neiCluster)
{
if (blockedFace[faceI])
{
blockedFace[faceI] = false;
nUnblocked++;
}
}
}
if (refinementHistory::debug)
{
reduce(nUnblocked, sumOp<label>());
Info<< type() << " : unblocked " << nUnblocked << " faces" << endl;
}
syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
}
void Foam::refinementHistory::apply
(
const boolList& blockedFace,
const PtrList<labelList>& specifiedProcessorFaces,
const labelList& specifiedProcessor,
const List<labelPair>& explicitConnections,
labelList& decomposition
) const
{
const polyMesh& mesh = dynamic_cast<const polyMesh&>(db());
// Find common parent for all cells
labelList cellToCluster;
label nClusters = markCommonCells(cellToCluster);
// Unblock all faces inbetween same cluster
labelList clusterToProc(nClusters, -1);
label nChanged = 0;
forAll(mesh.faceNeighbour(), faceI)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
label ownCluster = cellToCluster[own];
label neiCluster = cellToCluster[nei];
if (ownCluster != -1 && ownCluster == neiCluster)
{
if (clusterToProc[ownCluster] == -1)
{
clusterToProc[ownCluster] = decomposition[own];
}
if (decomposition[own] != clusterToProc[ownCluster])
{
decomposition[own] = clusterToProc[ownCluster];
nChanged++;
}
if (decomposition[nei] != clusterToProc[ownCluster])
{
decomposition[nei] = clusterToProc[ownCluster];
nChanged++;
}
}
}
if (refinementHistory::debug)
{
reduce(nChanged, sumOp<label>());
Info<< type() << " : changed decomposition on " << nChanged
<< " cells" << endl;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::refinementHistory::refinementHistory(const IOobject& io)
:
regIOobject(io)
regIOobject(io),
active_(false)
{
// Temporary warning
if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
@ -345,12 +579,18 @@ Foam::refinementHistory::refinementHistory(const IOobject& io)
close();
}
// When running in redistributPar + READ_IF_PRESENT it can happen
// that some processors do have refinementHistory and some don't so
// test for active has to be outside of above condition.
active_ = (returnReduce(visibleCells_.size(), sumOp<label>()) > 0);
if (debug)
{
Pout<< "refinementHistory::refinementHistory :"
<< " constructed history from IOobject :"
<< " splitCells:" << splitCells_.size()
<< " visibleCells:" << visibleCells_.size()
<< " active:" << active_
<< endl;
}
}
@ -360,10 +600,12 @@ Foam::refinementHistory::refinementHistory
(
const IOobject& io,
const List<splitCell8>& splitCells,
const labelList& visibleCells
const labelList& visibleCells,
const bool active
)
:
regIOobject(io),
active_(active),
splitCells_(splitCells),
freeSplitCells_(0),
visibleCells_(visibleCells)
@ -397,6 +639,66 @@ Foam::refinementHistory::refinementHistory
<< " constructed history from IOobject or components :"
<< " splitCells:" << splitCells_.size()
<< " visibleCells:" << visibleCells_.size()
<< " active:" << active_
<< endl;
}
}
Foam::refinementHistory::refinementHistory
(
const IOobject& io,
const label nCells
)
:
regIOobject(io),
active_(false),
freeSplitCells_(0)
{
// Temporary warning
if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
{
WarningInFunction
<< "Specified IOobject::MUST_READ_IF_MODIFIED but class"
<< " does not support automatic rereading."
<< endl;
}
if
(
io.readOpt() == IOobject::MUST_READ
|| io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
|| (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
readStream(typeName) >> *this;
close();
}
else
{
visibleCells_.setSize(nCells);
splitCells_.setCapacity(nCells);
for (label cellI = 0; cellI < nCells; cellI++)
{
visibleCells_[cellI] = cellI;
splitCells_.append(splitCell8());
}
}
active_ = (returnReduce(visibleCells_.size(), sumOp<label>()) > 0);
// Check indices.
checkIndices();
if (debug)
{
Pout<< "refinementHistory::refinementHistory :"
<< " constructed history from IOobject or initial size :"
<< " splitCells:" << splitCells_.size()
<< " visibleCells:" << visibleCells_.size()
<< " active:" << active_
<< endl;
}
}
@ -406,13 +708,15 @@ Foam::refinementHistory::refinementHistory
Foam::refinementHistory::refinementHistory
(
const IOobject& io,
const label nCells
const label nCells,
const bool active
)
:
regIOobject(io),
active_(active),
freeSplitCells_(0)
{
// Temporary warning
// Warn for MUST_READ_IF_MODIFIED
if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
{
WarningInFunction
@ -452,6 +756,7 @@ Foam::refinementHistory::refinementHistory
<< " constructed history from IOobject or initial size :"
<< " splitCells:" << splitCells_.size()
<< " visibleCells:" << visibleCells_.size()
<< " active:" << active_
<< endl;
}
}
@ -465,6 +770,7 @@ Foam::refinementHistory::refinementHistory
)
:
regIOobject(io),
active_(rh.active_),
splitCells_(rh.splitCells()),
freeSplitCells_(rh.freeSplitCells()),
visibleCells_(rh.visibleCells())
@ -477,6 +783,126 @@ Foam::refinementHistory::refinementHistory
}
// Construct from multiple
Foam::refinementHistory::refinementHistory
(
const IOobject& io,
const UPtrList<const labelList>& cellMaps,
const UPtrList<const refinementHistory>& refs
)
:
regIOobject(io),
active_(false)
{
if
(
io.readOpt() == IOobject::MUST_READ
|| io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
|| (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
WarningIn
(
"refinementHistory::refinementHistory(const IOobject&"
", const labelListList&, const PtrList<refinementHistory>&)"
) << "read option IOobject::MUST_READ, READ_IF_PRESENT or "
<< "MUST_READ_IF_MODIFIED"
<< " suggests that a read constructor would be more appropriate."
<< endl;
}
const polyMesh& mesh = dynamic_cast<const polyMesh&>(db());
// Determine offsets into splitCells
labelList offsets(refs.size()+1);
offsets[0] = 0;
forAll(refs, refI)
{
const DynamicList<splitCell8>& subSplits = refs[refI].splitCells();
offsets[refI+1] = offsets[refI]+subSplits.size();
}
// Construct merged splitCells
splitCells_.setSize(offsets.last());
forAll(refs, refI)
{
const DynamicList<splitCell8>& subSplits = refs[refI].splitCells();
forAll(subSplits, i)
{
splitCell8& newSplit = splitCells_[offsets[refI]+i];
// Copy
newSplit = subSplits[i];
// Offset indices
if (newSplit.parent_ >= 0)
{
newSplit.parent_ += offsets[refI];
}
if (newSplit.addedCellsPtr_.valid())
{
FixedList<label, 8>& splits = newSplit.addedCellsPtr_();
forAll(splits, i)
{
if (splits[i] >= 0)
{
splits[i] += offsets[refI];
}
}
}
}
}
// Construct merged visibleCells
visibleCells_.setSize(mesh.nCells(), -1);
forAll(refs, refI)
{
const labelList& cellMap = cellMaps[refI];
const labelList& subVis = refs[refI].visibleCells();
forAll(subVis, i)
{
label& newVis = visibleCells_[cellMap[i]];
newVis = subVis[i];
if (newVis >= 0)
{
newVis += offsets[refI];
}
}
}
// Is active if any of the refinementHistories is active (assumes active
// flag parallel synchronised)
active_ = false;
forAll(refs, refI)
{
if (refs[refI].active())
{
active_ = true;
break;
}
}
// Check indices.
checkIndices();
if (debug)
{
Pout<< "refinementHistory::refinementHistory :"
<< " constructed history from multiple refinementHistories :"
<< " splitCells:" << splitCells_.size()
<< " visibleCells:" << visibleCells_.size()
<< endl;
}
}
// Construct from Istream
Foam::refinementHistory::refinementHistory(const IOobject& io, Istream& is)
:
@ -485,6 +911,8 @@ Foam::refinementHistory::refinementHistory(const IOobject& io, Istream& is)
freeSplitCells_(0),
visibleCells_(is)
{
active_ = (returnReduce(visibleCells_.size(), sumOp<label>()) > 0);
// Check indices.
checkIndices();
@ -501,6 +929,192 @@ Foam::refinementHistory::refinementHistory(const IOobject& io, Istream& is)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::refinementHistory> Foam::refinementHistory::clone
(
const IOobject& io,
// Per visible cell the processor it is going to
const labelList& decomposition,
// Per splitCell entry the processor it moves to
const labelList& splitCellProc,
// Per splitCell entry the number of live cells that move to that processor
const labelList& splitCellNum,
const label procI,
// From old to new splitCells
labelList& oldToNewSplit
) const
{
oldToNewSplit.setSize(splitCells_.size());
oldToNewSplit = -1;
// Compacted splitCells
DynamicList<splitCell8> newSplitCells(splitCells_.size());
// Loop over all entries. Note: could recurse like countProc so only
// visit used entries but is probably not worth it.
forAll(splitCells_, index)
{
if (splitCellProc[index] == procI && splitCellNum[index] == 8)
{
// Entry moves in its whole to procI
oldToNewSplit[index] = newSplitCells.size();
newSplitCells.append(splitCells_[index]);
}
}
// Add live cells that are subsetted.
forAll(visibleCells_, cellI)
{
label index = visibleCells_[cellI];
if (index >= 0 && decomposition[cellI] == procI)
{
label parent = splitCells_[index].parent_;
// Create new splitCell with parent
oldToNewSplit[index] = newSplitCells.size();
newSplitCells.append(splitCell8(parent));
}
}
//forAll(oldToNewSplit, index)
//{
// Pout<< "old:" << index << " new:" << oldToNewSplit[index]
// << endl;
//}
newSplitCells.shrink();
// Renumber contents of newSplitCells
forAll(newSplitCells, index)
{
splitCell8& split = newSplitCells[index];
if (split.parent_ >= 0)
{
split.parent_ = oldToNewSplit[split.parent_];
}
if (split.addedCellsPtr_.valid())
{
FixedList<label, 8>& splits = split.addedCellsPtr_();
forAll(splits, i)
{
if (splits[i] >= 0)
{
splits[i] = oldToNewSplit[splits[i]];
}
}
}
}
// Count number of cells
label nSub = 0;
forAll(decomposition, cellI)
{
if (decomposition[cellI] == procI)
{
nSub++;
}
}
labelList newVisibleCells(nSub);
nSub = 0;
forAll(visibleCells_, cellI)
{
if (decomposition[cellI] == procI)
{
label index = visibleCells_[cellI];
if (index >= 0)
{
index = oldToNewSplit[index];
}
newVisibleCells[nSub++] = index;
}
}
return autoPtr<refinementHistory>
(
new refinementHistory
(
io,
newSplitCells,
newVisibleCells,
active_
)
);
}
Foam::autoPtr<Foam::refinementHistory> Foam::refinementHistory::clone
(
const IOobject& io,
const labelList& cellMap
) const
{
if (active_)
{
// Mark selected cells with '1'
labelList decomposition(visibleCells_.size(), 0);
forAll(cellMap, i)
{
decomposition[cellMap[i]] = 1;
}
// Per splitCell entry the processor it moves to
labelList splitCellProc(splitCells_.size(), -1);
// Per splitCell entry the number of live cells that move to that
// processor
labelList splitCellNum(splitCells_.size(), 0);
forAll(visibleCells_, cellI)
{
label index = visibleCells_[cellI];
if (index >= 0)
{
countProc
(
splitCells_[index].parent_,
decomposition[cellI],
splitCellProc,
splitCellNum
);
}
}
labelList oldToNewSplit;
return clone
(
io,
decomposition,
splitCellProc,
splitCellNum,
1, //procI,
oldToNewSplit
);
}
else
{
return autoPtr<refinementHistory>
(
new refinementHistory
(
io,
DynamicList<splitCell8>(0),
labelList(0),
false
)
);
}
}
void Foam::refinementHistory::resize(const label size)
{
label oldSize = visibleCells_.size();
@ -698,9 +1312,6 @@ void Foam::refinementHistory::distribute(const mapDistributePolyMesh& map)
}
}
//Pout<< "refinementHistory::distribute :"
// << " destination:" << destination << endl;
// Per splitCell entry the processor it moves to
labelList splitCellProc(splitCells_.size(), -1);
// Per splitCell entry the number of live cells that move to that processor
@ -746,21 +1357,11 @@ void Foam::refinementHistory::distribute(const mapDistributePolyMesh& map)
forAll(splitCells_, index)
{
// Pout<< "oldCell:" << index
// << " proc:" << splitCellProc[index]
// << " nCells:" << splitCellNum[index]
// << endl;
if (splitCellProc[index] == proci && splitCellNum[index] == 8)
{
// Entry moves in its whole to proci
oldToNew[index] = newSplitCells.size();
newSplitCells.append(splitCells_[index]);
//Pout<< "Added oldCell " << index
// << " info " << newSplitCells.last()
// << " at position " << newSplitCells.size()-1
// << endl;
}
}
@ -773,10 +1374,6 @@ void Foam::refinementHistory::distribute(const mapDistributePolyMesh& map)
{
label parent = splitCells_[index].parent_;
//Pout<< "Adding refined cell " << celli
// << " since moves to "
// << proci << " old parent:" << parent << endl;
// Create new splitCell with parent
oldToNew[index] = newSplitCells.size();
newSplitCells.append(splitCell8(parent));
@ -849,7 +1446,9 @@ void Foam::refinementHistory::distribute(const mapDistributePolyMesh& map)
// Remove all entries. Leave storage intact.
splitCells_.clear();
visibleCells_.setSize(map.mesh().nCells());
const polyMesh& mesh = dynamic_cast<const polyMesh&>(db());
visibleCells_.setSize(mesh.nCells());
visibleCells_ = -1;
for (label proci = 0; proci < Pstream::nProcs(); proci++)
@ -1138,6 +1737,17 @@ void Foam::refinementHistory::combineCells
}
bool Foam::refinementHistory::read()
{
bool ok = readData(readStream(typeName));
close();
active_ = (returnReduce(visibleCells_.size(), sumOp<label>()) > 0);
return ok;
}
bool Foam::refinementHistory::readData(Istream& is)
{
is >> *this;

153
src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.H → src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/refinementHistory.H
View File

@ -28,7 +28,7 @@ Description
All refinement history. Used in unrefinement.
- visibleCells: valid for the current mesh and contains per cell -1
(cell unrefined) or an index into splitCells_.
(cell unrefined) or an index into splitCells_.
- splitCells: for every split contains the parent (also index into
splitCells) and optionally a subsplit as 8 indices into splitCells.
Note that the numbers in splitCells are not cell labels, they are purely
@ -74,9 +74,10 @@ SourceFiles
#include "DynamicList.H"
#include "labelList.H"
#include "FixedList.H"
#include "SLList.H"
#include "autoPtr.H"
#include "regIOobject.H"
#include "boolList.H"
#include "labelPair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -122,50 +123,11 @@ public:
splitCell8(const splitCell8&);
//- Copy operator since autoPtr otherwise 'steals' storage.
void operator=(const splitCell8& s)
{
// Check for assignment to self
if (this == &s)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
void operator=(const splitCell8& s);
parent_ = s.parent_;
bool operator==(const splitCell8& s) const;
addedCellsPtr_.reset
(
s.addedCellsPtr_.valid()
? new FixedList<label, 8>(s.addedCellsPtr_())
: NULL
);
}
bool operator==(const splitCell8& s) const
{
if (addedCellsPtr_.valid() != s.addedCellsPtr_.valid())
{
return false;
}
else if (parent_ != s.parent_)
{
return false;
}
else if (addedCellsPtr_.valid())
{
return addedCellsPtr_() == s.addedCellsPtr_();
}
else
{
return true;
}
}
bool operator!=(const splitCell8& s) const
{
return !operator==(s);
}
bool operator!=(const splitCell8& s) const;
friend Istream& operator>>(Istream&, splitCell8&);
friend Ostream& operator<<(Ostream&, const splitCell8&);
@ -176,6 +138,9 @@ private:
// Private data
//- Is active?
bool active_;
//- Storage for splitCells
DynamicList<splitCell8> splitCells_;
@ -226,6 +191,15 @@ private:
labelList& splitCellNum
) const;
// For distribution:
//- Mark index and all its descendants
void mark(const label, const label, labelList&) const;
//- Mark cells according to top parent. Return number of clusters
// (set of cells originating from same parent)
label markCommonCells(labelList& cellToCluster) const;
public:
// Declare name of the class and its debug switch
@ -234,25 +208,47 @@ public:
// Constructors
//- Construct (read) given an IOobject
//- Construct (read) given an IOobject. If global number of visible
// cells > 0 becomes active
refinementHistory(const IOobject&);
//- Construct (read) or construct null
//- Construct (read) or construct from components
refinementHistory
(
const IOobject&,
const List<splitCell8>& splitCells,
const labelList& visibleCells
const labelList& visibleCells,
const bool active
);
//- Construct (read) or construct from initial number of cells
// (all visible)
// (all visible). If global number of visible
// cells > 0 becomes active
refinementHistory(const IOobject&, const label nCells);
//- Construct (read) or construct from initial number of cells
// (all visible) and active flag
refinementHistory
(
const IOobject&,
const label nCells,
const bool active
);
//- Construct as copy
refinementHistory(const IOobject&, const refinementHistory&);
//- Construct from Istream
//- Construct from multiple refinement histories. If global number of
// visible cells > 0 becomes active
refinementHistory
(
const IOobject&,
const UPtrList<const labelList>& cellMaps,
const UPtrList<const refinementHistory>&
);
//- Construct from Istream. If global number of
// visible cells > 0 becomes active
refinementHistory(const IOobject&, Istream&);
@ -278,12 +274,16 @@ public:
return freeSplitCells_;
}
//- Is there unrefinement history. Note that this will fall over if
// there are 0 cells in the mesh. But this gives problems with
// lots of other programs anyway.
//- Is there unrefinement history?
bool active() const
{
return visibleCells_.size() > 0;
return active_;
}
//- Is there unrefinement history?
bool& active()
{
return active_;
}
//- Get parent of cell
@ -314,6 +314,23 @@ public:
const labelList& combinedCells
);
//- Low level clone
autoPtr<refinementHistory> clone
(
const IOobject& io,
const labelList& decomposition,
const labelList& splitCellProc,
const labelList& splitCellNum,
const label procI,
labelList& oldToNewSplit
) const;
//- Create clone from subset
autoPtr<refinementHistory> clone
(
const IOobject& io,
const labelList& cellMap
) const;
//- Update numbering for mesh changes
void updateMesh(const mapPolyMesh&);
@ -343,20 +360,44 @@ public:
void writeDebug() const;
//- ReadData function required for regIOobject read operation
//- Read object. If global number of visible cells > 0 becomes active
virtual bool read();
//- ReadData function required for regIOobject read operation. Note:
// does not do a reduction - does not set active_ flag
virtual bool readData(Istream&);
//- WriteData function required for regIOobject write operation
virtual bool writeData(Ostream&) const;
// Helpers for decompositionConstraint
// Friend Functions
//- Add my decomposition constraints
void add
(
boolList& blockedFace,
PtrList<labelList>& specifiedProcessorFaces,
labelList& specifiedProcessor,
List<labelPair>& explicitConnections
) const;
//- Apply any additional post-decomposition constraints
void apply
(
const boolList& blockedFace,
const PtrList<labelList>& specifiedProcessorFaces,
const labelList& specifiedProcessor,
const List<labelPair>& explicitConnections,
labelList& decomposition
) const;
// Friend Operators
// IOstream Operators
//- Istream operator. Note: does not do a reduction - does not set
// active_ flag
friend Istream& operator>>(Istream&, refinementHistory&);
friend Ostream& operator<<(Ostream&, const refinementHistory&);
};