FEATURE: Update of dynamicTopoFvMesh/mesquiteMotionSolver. Author: Sandeep Menon. Merge: Dominik Christ.

This commit is contained in:
Dominik Christ 2015-04-27 10:04:06 +01:00
commit ae870dc409
51 changed files with 5764 additions and 4340 deletions

View file

@ -7,7 +7,9 @@ wmake libso dynamicMesh
# Make meshMotion solvers # Make meshMotion solvers
meshMotion/Allwmake meshMotion/Allwmake
wmake libso dynamicFvMesh # Make libraries in dynamicFvMesh
dynamicFvMesh/Allwmake
wmake libso topoChangerFvMesh wmake libso topoChangerFvMesh
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file

View file

@ -0,0 +1,12 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
# Clean dynamicFvMesh libraries
wclean
wclean dynamicTopoFvMesh
# Wipe out all lnInclude directories and re-link
wcleanLnIncludeAll
# ----------------------------------------------------------------- end-of-file

View file

@ -0,0 +1,9 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
# Make dynamicFvMesh libraries
wmake libso
wmake libso dynamicTopoFvMesh
# ----------------------------------------------------------------- end-of-file

View file

@ -17,42 +17,4 @@ turboFvMesh/turboFvMesh.C
fvMeshAdder/fvMeshAdder.C fvMeshAdder/fvMeshAdder.C
fvMeshDistribute/fvMeshDistribute.C fvMeshDistribute/fvMeshDistribute.C
tetMetrics = dynamicTopoFvMesh/tetMetrics
$(tetMetrics)/tetMetric.C
$(tetMetrics)/tetMetrics.C
lengthScaleEstimator = dynamicTopoFvMesh/lengthScaleEstimator
$(lengthScaleEstimator)/lengthScaleEstimator.C
eMesh = dynamicTopoFvMesh/eMesh
$(eMesh)/eMesh.C
$(eMesh)/eMeshDemandDrivenData.C
$(eMesh)/eBoundaryMesh/eBoundaryMesh.C
ePatches = $(eMesh)/ePatches
$(ePatches)/ePatch/ePatch.C
$(ePatches)/ePatch/newEPatch.C
dynamicTopoFvMesh/meshOps.C
dynamicTopoFvMesh/dynamicTopoFvMesh.C
dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C
dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
dynamicTopoFvMesh/dynamicTopoFvMeshCoupled.C
dynamicTopoFvMesh/edgeSwap.C
dynamicTopoFvMesh/edgeBisect.C
dynamicTopoFvMesh/edgeCollapse.C
dynamicTopoFvMesh/coupleMap.C
convexSetAlgorithm = dynamicTopoFvMesh/convexSetAlgorithm
$(convexSetAlgorithm)/convexSetAlgorithm.C
$(convexSetAlgorithm)/faceSetAlgorithm.C
$(convexSetAlgorithm)/cellSetAlgorithm.C
fieldMapping = dynamicTopoFvMesh/fieldMapping
$(fieldMapping)/topoMapper.C
$(fieldMapping)/fluxCorrector.C
$(fieldMapping)/topoCellMapper.C
$(fieldMapping)/topoPatchMapper.C
$(fieldMapping)/topoSurfaceMapper.C
LIB = $(FOAM_LIBBIN)/libdynamicFvMesh LIB = $(FOAM_LIBBIN)/libdynamicFvMesh

View file

@ -0,0 +1,39 @@
eMesh/eMesh.C
eMesh/eMeshDemandDrivenData.C
eMesh/eBoundaryMesh/eBoundaryMesh.C
ePatches = eMesh/ePatches
$(ePatches)/ePatch/ePatch.C
$(ePatches)/ePatch/newEPatch.C
dynamicTopoFvMesh.C
dynamicTopoFvMeshCheck.C
dynamicTopoFvMeshReOrder.C
dynamicTopoFvMeshMapping.C
edgeSwap.C
edgeBisect.C
edgeCollapse.C
coupledMesh/coupleMap.C
coupledMesh/dynamicTopoFvMeshCoupled.C
coupledMesh/subMeshProcessorPolyPatch.C
coupledMesh/subMeshProcessorFvPatch.C
convexSetAlgorithm/convexSetAlgorithm.C
convexSetAlgorithm/faceSetAlgorithm.C
convexSetAlgorithm/cellSetAlgorithm.C
fieldMapping/topoMapper.C
fieldMapping/fluxCorrector.C
fieldMapping/topoCellMapper.C
fieldMapping/topoPatchMapper.C
fieldMapping/topoSurfaceMapper.C
tetMetrics = tetMetrics
$(tetMetrics)/tetMetric.C
$(tetMetrics)/tetMetrics.C
lengthScaleEstimator = lengthScaleEstimator
$(lengthScaleEstimator)/lengthScaleEstimator.C
LIB = $(FOAM_LIBBIN)/libdynamicTopoFvMesh

View file

@ -0,0 +1,13 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/dynamicFvMesh \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/fvMeshDistribute \
-I$(LIB_SRC)/decompositionMethods/decompositionMethods/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
LIB_LIBS = \
-ldynamicMesh \
-ldynamicFvMesh \
-ldecompositionMethods \
-lfiniteVolume

View file

@ -34,7 +34,6 @@ Author
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "objectRegistry.H"
#include "Time.H" #include "Time.H"
#include "IOMap.H" #include "IOMap.H"
#include "meshOps.H" #include "meshOps.H"

View file

@ -1,647 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "coupledInfo.H"
#include "dynamicTopoFvMesh.H"
#include "emptyFvPatchFields.H"
#include "emptyFvsPatchFields.H"
#include "fixedValueFvPatchFields.H"
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Constructor for coupledInfo
coupledInfo::coupledInfo
(
const dynamicTopoFvMesh& mesh,
const coupleMap& cMap,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_(cMap),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
coupledInfo::coupledInfo
(
const dynamicTopoFvMesh& mesh,
const bool isTwoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label mPatch,
const label sPatch,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_
(
IOobject
(
"coupleMap_"
+ Foam::name(mPatch)
+ "_To_"
+ Foam::name(sPatch)
+ word(isLocal ? "_Local" : "_Proc")
+ word(isSend ? "_Send" : "_Recv"),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
true
),
isTwoDMesh,
isLocal,
isSend,
patchIndex,
mPatch,
sPatch
),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
//- Construct given addressing
coupledInfo::subMeshMapper::subMeshMapper
(
const coupledInfo& cInfo,
const label patchI
)
:
sizeBeforeMapping_(cInfo.baseMesh().boundary()[patchI].size()),
directAddressing_
(
SubList<label>
(
cInfo.map().faceMap(),
cInfo.subMesh().boundary()[patchI].size(),
cInfo.subMesh().boundary()[patchI].patch().start()
)
)
{
// Offset indices
label pStart = cInfo.baseMesh().boundary()[patchI].patch().start();
forAll(directAddressing_, faceI)
{
directAddressing_[faceI] -= pStart;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const dynamicTopoFvMesh& coupledInfo::baseMesh() const
{
return mesh_;
}
void coupledInfo::setMesh
(
label index,
dynamicTopoFvMesh* mesh
)
{
subMesh_.set(mesh);
}
dynamicTopoFvMesh& coupledInfo::subMesh()
{
if (!subMesh_.valid())
{
FatalErrorIn("dynamicTopoFvMesh& coupledInfo::subMesh()")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
const dynamicTopoFvMesh& coupledInfo::subMesh() const
{
if (!subMesh_.valid())
{
FatalErrorIn("const dynamicTopoFvMesh& coupledInfo::subMesh() const")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
bool coupledInfo::builtMaps() const
{
return builtMaps_;
}
void coupledInfo::setBuiltMaps()
{
builtMaps_ = true;
}
coupleMap& coupledInfo::map()
{
return map_;
}
const coupleMap& coupledInfo::map() const
{
return map_;
}
label coupledInfo::masterFaceZone() const
{
return masterFaceZone_;
}
label coupledInfo::slaveFaceZone() const
{
return slaveFaceZone_;
}
// Set subMesh centres
void coupledInfo::setCentres(PtrList<volVectorField>& centres) const
{
// Fetch reference to subMesh
const dynamicTopoFvMesh& mesh = subMesh();
// Set size
centres.setSize(1);
vectorField Cv(mesh.cellCentres());
vectorField Cf(mesh.faceCentres());
// Create and map the patch field values
label nPatches = mesh.boundary().size();
// Create field parts
PtrList<fvPatchField<vector> > volCentrePatches(nPatches);
// Over-ride and set all patches to fixedValue
for (label patchI = 0; patchI < nPatches; patchI++)
{
volCentrePatches.set
(
patchI,
new fixedValueFvPatchField<vector>
(
mesh.boundary()[patchI],
DimensionedField<vector, volMesh>::null()
)
);
// Slice field to patch (forced assignment)
volCentrePatches[patchI] ==
(
mesh.boundaryMesh()[patchI].patchSlice(Cf)
);
}
// Set the cell-centres pointer.
centres.set
(
0,
new volVectorField
(
IOobject
(
"cellCentres",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimLength,
SubField<vector>(Cv, mesh.nCells()),
volCentrePatches
)
);
}
// Subset volume field
template <class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
coupledInfo::subSetVolField
(
const GeometricField<Type, fvPatchField, volMesh>& fld
) const
{
// Create and map the internal-field values
Field<Type> internalField
(
fld.internalField(),
map().cellMap()
);
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
PtrList<fvPatchField<Type> > patchFields(nPatches);
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
new emptyFvPatchField<Type>
(
subMesh().boundary()[patchI],
DimensionedField<Type, volMesh>::null()
)
);
}
else
{
patchFields.set
(
patchI,
fvPatchField<Type>::New
(
fld.boundaryField()[patchI],
subMesh().boundary()[patchI],
DimensionedField<Type, volMesh>::null(),
subMeshMapper(*this, patchI)
)
);
}
}
// Create new field from pieces
tmp<GeometricField<Type, fvPatchField, volMesh> > subFld
(
new GeometricField<Type, fvPatchField, volMesh>
(
IOobject
(
"subField_" + fld.name(),
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fld.dimensions(),
internalField,
patchFields
)
);
return subFld;
}
// Subset surface field
template <class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
coupledInfo::subSetSurfaceField
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& fld
) const
{
// Create and map the internal-field values
Field<Type> internalField
(
fld.internalField(),
SubList<label>
(
map().faceMap(),
subMesh().nInternalFaces()
)
);
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
PtrList<fvsPatchField<Type> > patchFields(nPatches);
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
new emptyFvsPatchField<Type>
(
subMesh().boundary()[patchI],
DimensionedField<Type, surfaceMesh>::null()
)
);
}
else
{
patchFields.set
(
patchI,
fvsPatchField<Type>::New
(
fld.boundaryField()[patchI],
subMesh().boundary()[patchI],
DimensionedField<Type, surfaceMesh>::null(),
subMeshMapper(*this, patchI)
)
);
}
}
// Create new field from pieces
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > subFld
(
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
IOobject
(
"subField_" + fld.name(),
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fld.dimensions(),
internalField,
patchFields
)
);
return subFld;
}
template <class Type>
void coupledInfo::mapVolField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const
{
strStream
<< fieldType << token::NL
<< token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
const GeometricField<Type, fvPatchField, volMesh>& fld =
(
mesh_.lookupObject
<
GeometricField<Type, fvPatchField, volMesh>
>(fieldNames[i])
);
tmp<GeometricField<Type, fvPatchField, volMesh> > tsubFld =
(
subSetVolField(fld)
);
// Send field through stream
strStream
<< fieldNames[i]
<< token::NL << token::BEGIN_BLOCK
<< tsubFld
<< token::NL << token::END_BLOCK
<< token::NL;
}
strStream
<< token::END_BLOCK << token::NL;
}
template <class Type>
void coupledInfo::mapSurfaceField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const
{
strStream
<< fieldType << token::NL
<< token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
const GeometricField<Type, fvsPatchField, surfaceMesh>& fld =
(
mesh_.lookupObject
<
GeometricField<Type, fvsPatchField, surfaceMesh>
>(fieldNames[i])
);
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsubFld =
(
subSetSurfaceField(fld)
);
// Send field through stream
strStream
<< fieldNames[i]
<< token::NL << token::BEGIN_BLOCK
<< tsubFld
<< token::NL << token::END_BLOCK
<< token::NL;
}
strStream
<< token::END_BLOCK << token::NL;
}
// Set volume field pointer from input dictionary
template <class GeomField>
void coupledInfo::setField
(
const wordList& fieldNames,
const dictionary& fieldDicts,
PtrList<GeomField>& fields
) const
{
// Size up the pointer list
fields.setSize(fieldNames.size());
forAll(fieldNames, i)
{
fields.set
(
i,
new GeomField
(
IOobject
(
fieldNames[i],
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fieldDicts.subDict(fieldNames[i])
)
);
}
}
template <class GeomField>
void coupledInfo::resizeMap
(
const label srcIndex,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields,
GeomField& field
)
{
// autoMap the internal field
field.internalField().autoMap(internalMapper);
// Reverse map for additional cells
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.internalField().rmap
(
srcField.internalField(),
internalReverseMaps[pI]
);
}
// Map physical boundary-fields
forAll(boundaryMapper, patchI)
{
// autoMap the patchField
field.boundaryField()[patchI].autoMap(boundaryMapper[patchI]);
// Reverse map for additional patch faces
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.boundaryField()[patchI].rmap
(
srcField.boundaryField()[patchI],
boundaryReverseMaps[pI][patchI]
);
}
}
}
// Resize all fields in registry
template <class GeomField>
void coupledInfo::resizeMap
(
const wordList& names,
const objectRegistry& mesh,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields
)
{
forAll(names, indexI)
{
// Fetch field from registry
GeomField& field =
(
const_cast<GeomField&>
(
mesh.lookupObject<GeomField>(names[indexI])
)
);
// Map the field
coupledInfo::resizeMap
(
indexI,
internalMapper,
internalReverseMaps,
boundaryMapper,
boundaryReverseMaps,
srcFields,
field
);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void coupledInfo::operator=(const coupledInfo& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn
(
"void coupledInfo::operator=(const Foam::coupledInfo&)"
)
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -52,6 +52,7 @@ const char* coupleMap::names[coupleMap::INVALID + 1] =
"REMOVE_CELL", "REMOVE_CELL",
"MOVE_POINT", "MOVE_POINT",
"CONVERT_PATCH", "CONVERT_PATCH",
"CONVERT_PHYSICAL",
"INVALID" "INVALID"
}; };
@ -204,8 +205,8 @@ void coupleMap::makeFaces() const
face& f = faces[faceI]; face& f = faces[faceI];
labelList& fe = faceEdges[faceI]; labelList& fe = faceEdges[faceI];
// Fetch the buffer value for 2D meshes // Fetch the buffer value
label nfe = twoDMesh_ ? nfeBuffer[faceI] : 3; label nfe = nfeBuffer[faceI];
// Size up the lists // Size up the lists
f.setSize(nfe, -1); f.setSize(nfe, -1);
@ -279,6 +280,22 @@ void coupleMap::makeCellMap() const
} }
void coupleMap::makeInternalFaceMap() const
{
// It is an error to attempt to recalculate
// if the map is already calculated
if (internalFaceMap_.size())
{
FatalErrorIn("void coupleMap::makeInternalFaceMap() const")
<< "internal faceMap has already been calculated."
<< abort(FatalError);
}
// Slice for internal faces
internalFaceMap_ = SubList<label>(faceMap(), nEntities(INTERNAL_FACE));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
pointField& coupleMap::pointBuffer() const pointField& coupleMap::pointBuffer() const
@ -343,12 +360,8 @@ void coupleMap::allocateBuffers() const
entityBuffer(EDGE_SIZES).setSize(nEntities(NBDY)); entityBuffer(EDGE_SIZES).setSize(nEntities(NBDY));
entityBuffer(PATCH_ID).setSize(nEntities(NBDY)); entityBuffer(PATCH_ID).setSize(nEntities(NBDY));
// nFaceEdges buffer is required only for 2D, // Set face-sizes
// due to a mix of triangle / quad faces
if (twoDMesh_)
{
entityBuffer(NFE_BUFFER).setSize(nEntities(FACE)); entityBuffer(NFE_BUFFER).setSize(nEntities(FACE));
}
// Allocate for variable size face-lists // Allocate for variable size face-lists
entityBuffer(FACE).setSize(nEntities(NFE_SIZE)); entityBuffer(FACE).setSize(nEntities(NFE_SIZE));
@ -446,6 +459,43 @@ void coupleMap::mapMaster
} }
void coupleMap::pushOperation
(
const label index,
const opType oType
) const
{
entityIndices_.setSize(entityIndices_.size() + 1, index);
entityOperations_.setSize(entityOperations_.size() + 1, oType);
}
void coupleMap::pushOperation
(
const label index,
const opType oType,
const label pIndex
) const
{
if (oType == coupleMap::CONVERT_PHYSICAL)
{
entityIndices_.setSize(entityIndices_.size() + 1, index);
entityOperations_.setSize(entityOperations_.size() + 1, oType);
patchIndices_.setSize(patchIndices_.size() + 1, pIndex);
}
else
{
opType t = (oType < coupleMap::INVALID ? oType : coupleMap::INVALID);
FatalErrorIn("void coupleMap::pushOperation() const")
<< " Expected CONVERT_PHYSICAL" << nl
<< " Found: " << coupleMap::names[t]
<< abort(FatalError);
}
}
void coupleMap::pushOperation void coupleMap::pushOperation
( (
const label index, const label index,
@ -466,6 +516,15 @@ void coupleMap::pushOperation
moveNewPoints_.setSize(moveNewPoints_.size() + 1, newPoint); moveNewPoints_.setSize(moveNewPoints_.size() + 1, newPoint);
moveOldPoints_.setSize(moveOldPoints_.size() + 1, oldPoint); moveOldPoints_.setSize(moveOldPoints_.size() + 1, oldPoint);
} }
else
{
opType t = (oType < coupleMap::INVALID ? oType : coupleMap::INVALID);
FatalErrorIn("void coupleMap::pushOperation() const")
<< " Expected either MOVE_POINT or CONVERT_PATCH" << nl
<< " Found: " << coupleMap::names[t]
<< abort(FatalError);
}
} }
@ -485,6 +544,10 @@ void coupleMap::clearMaps() const
{ {
faceMap_.clear(); faceMap_.clear();
cellMap_.clear(); cellMap_.clear();
internalFaceMap_.clear();
subMeshPointMap_.clear();
subMeshEdgeMap_.clear();
forAll(entityMap_, mapI) forAll(entityMap_, mapI)
{ {
@ -510,6 +573,7 @@ void coupleMap::clearBuffers() const
entityIndices_.clear(); entityIndices_.clear();
entityOperations_.clear(); entityOperations_.clear();
patchIndices_.clear();
moveNewPoints_.clear(); moveNewPoints_.clear();
moveOldPoints_.clear(); moveOldPoints_.clear();
} }
@ -588,6 +652,17 @@ const labelList& coupleMap::cellMap() const
} }
const labelList& coupleMap::internalFaceMap() const
{
if (internalFaceMap_.empty())
{
makeInternalFaceMap();
}
return internalFaceMap_;
}
bool coupleMap::readData(Istream& is) bool coupleMap::readData(Istream& is)
{ {
Map<label> tmpMap(is); Map<label> tmpMap(is);

View file

@ -78,12 +78,14 @@ public:
FACE = 2, FACE = 2,
// Sizes only // Sizes only
CELL = 3, CELL = 3,
MAX_ENTITIES = 4,
INTERNAL_EDGE = 4, INTERNAL_EDGE = 4,
INTERNAL_FACE = 5, INTERNAL_FACE = 5,
SHARED_POINT = 6, SHARED_POINT = 6,
GLOBAL_POINT = 7, GLOBAL_POINT = 7,
NFE_SIZE = 8, NFE_SIZE = 8,
NBDY = 9, NBDY = 9,
MAX_SIZES = 10,
// Buffers only // Buffers only
OWNER = 3, OWNER = 3,
NEIGHBOUR = 4, NEIGHBOUR = 4,
@ -93,7 +95,8 @@ public:
FACE_SIZES = 8, FACE_SIZES = 8,
EDGE_STARTS = 9, EDGE_STARTS = 9,
EDGE_SIZES = 10, EDGE_SIZES = 10,
PATCH_ID = 11 PATCH_ID = 11,
MAX_BUFFERS = 12
}; };
//- Enumerants for operations //- Enumerants for operations
@ -106,9 +109,15 @@ public:
REMOVE_CELL = 4, REMOVE_CELL = 4,
MOVE_POINT = 5, MOVE_POINT = 5,
CONVERT_PATCH = 6, CONVERT_PATCH = 6,
CONVERT_PHYSICAL = 7,
INVALID INVALID
}; };
//- Public typedefs
typedef FixedList<label, MAX_SIZES> EntitySizeList;
typedef FixedList<Map<label>, MAX_ENTITIES> EntityMapList;
typedef FixedList<labelList, MAX_BUFFERS> EntityBufferList;
private: private:
// Private data // Private data
@ -136,14 +145,14 @@ private:
mutable List<labelPair> globalProcPoints_; mutable List<labelPair> globalProcPoints_;
// Entity sizes (as specified by the entityType enumerant) // Entity sizes (as specified by the entityType enumerant)
mutable FixedList<label,10> nEntities_; mutable EntitySizeList nEntities_;
// Maps for entities // Maps for entities
mutable FixedList<Map<label>,4> entityMap_; mutable EntityMapList entityMap_;
mutable FixedList<Map<label>,4> reverseEntityMap_; mutable EntityMapList reverseEntityMap_;
// Entity Buffers (as specified by the entityType enumerant) // Entity Buffers (as specified by the entityType enumerant)
mutable FixedList<labelList,12> entityBuffer_; mutable EntityBufferList entityBuffer_;
// List of entity indices with topological operations // List of entity indices with topological operations
mutable labelList entityIndices_; mutable labelList entityIndices_;
@ -151,6 +160,9 @@ private:
// List of operations performed on entities // List of operations performed on entities
mutable List<opType> entityOperations_; mutable List<opType> entityOperations_;
// Physical patch conversion indices
mutable labelList patchIndices_;
// List of point-locations to move points to // List of point-locations to move points to
mutable pointField moveNewPoints_; mutable pointField moveNewPoints_;
mutable pointField moveOldPoints_; mutable pointField moveOldPoints_;
@ -158,17 +170,25 @@ private:
// Addressing for field mapping // Addressing for field mapping
mutable labelList faceMap_; mutable labelList faceMap_;
mutable labelList cellMap_; mutable labelList cellMap_;
mutable labelList internalFaceMap_;
// Processor point and edge mapping
mutable Map<labelList> subMeshPointMap_;
mutable Map<labelList> subMeshEdgeMap_;
//- Demand-driven connectivity data. //- Demand-driven connectivity data.
mutable edgeList* edgesPtr_; mutable edgeList* edgesPtr_;
mutable faceList* facesPtr_; mutable faceList* facesPtr_;
mutable labelListList* faceEdgesPtr_; mutable labelListList* faceEdgesPtr_;
// Private Member Functions
void makeEdges() const; void makeEdges() const;
void makeFaces() const; void makeFaces() const;
void makeFaceMap() const; void makeFaceMap() const;
void makeCellMap() const; void makeCellMap() const;
void makeInternalFaceMap() const;
void clearAddressing() const; void clearAddressing() const;
@ -311,27 +331,44 @@ public:
const label master const label master
) const; ) const;
inline FixedList<label,10>& nEntities() const; inline EntitySizeList& nEntities() const;
inline label& nEntities(const label eType) const; inline label& nEntities(const label eType) const;
inline Map<label>& entityMap(const label eType) const; inline Map<label>& entityMap(const label eType) const;
inline Map<label>& reverseEntityMap(const label eType) const; inline Map<label>& reverseEntityMap(const label eType) const;
inline FixedList<labelList,12>& entityBuffer() const; inline EntityBufferList& entityBuffer() const;
inline labelList& entityBuffer(const label eType) const; inline labelList& entityBuffer(const label eType) const;
inline labelList& entityIndices() const; inline labelList& entityIndices() const;
inline List<opType>& entityOperations() const; inline List<opType>& entityOperations() const;
inline labelList& patchIndices() const;
inline pointField& moveNewPoints() const; inline pointField& moveNewPoints() const;
inline pointField& moveOldPoints() const; inline pointField& moveOldPoints() const;
inline Map<labelList>& subMeshPointMap() const;
inline Map<labelList>& subMeshEdgeMap() const;
void pushOperation
(
const label index,
const opType oType
) const;
void pushOperation void pushOperation
( (
const label index, const label index,
const opType oType, const opType oType,
const point& newPoint = vector::zero, const label pIndex
const point& oldPoint = vector::zero ) const;
void pushOperation
(
const label index,
const opType oType,
const point& newPoint,
const point& oldPoint
) const; ) const;
//- Demand-driven data //- Demand-driven data
@ -344,6 +381,7 @@ public:
const labelList& faceMap() const; const labelList& faceMap() const;
const labelList& cellMap() const; const labelList& cellMap() const;
const labelList& internalFaceMap() const;
void transferMaps void transferMaps
( (

View file

@ -222,7 +222,7 @@ inline bool coupleMap::isRecv() const
} }
inline FixedList<label,10>& coupleMap::nEntities() const inline coupleMap::EntitySizeList& coupleMap::nEntities() const
{ {
return nEntities_; return nEntities_;
} }
@ -246,7 +246,7 @@ inline Map<label>& coupleMap::reverseEntityMap(const label eType) const
} }
inline FixedList<labelList,12>& coupleMap::entityBuffer() const inline coupleMap::EntityBufferList& coupleMap::entityBuffer() const
{ {
return entityBuffer_; return entityBuffer_;
} }
@ -270,6 +270,12 @@ inline List<coupleMap::opType>& coupleMap::entityOperations() const
} }
inline labelList& coupleMap::patchIndices() const
{
return patchIndices_;
}
inline pointField& coupleMap::moveNewPoints() const inline pointField& coupleMap::moveNewPoints() const
{ {
return moveNewPoints_; return moveNewPoints_;
@ -282,6 +288,18 @@ inline pointField& coupleMap::moveOldPoints() const
} }
inline Map<labelList>& coupleMap::subMeshPointMap() const
{
return subMeshPointMap_;
}
inline Map<labelList>& coupleMap::subMeshEdgeMap() const
{
return subMeshEdgeMap_;
}
} // End namespace Foam } // End namespace Foam
// ************************************************************************* // // ************************************************************************* //

View file

@ -0,0 +1,747 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "coupledInfo.H"
#include "emptyPolyPatch.H"
#include "processorPolyPatch.H"
#include "fixedValueFvPatchFields.H"
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct given mesh, coupleMap and master / slave indices
template <class MeshType>
coupledInfo<MeshType>::coupledInfo
(
const MeshType& mesh,
const coupleMap& cMap,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_(cMap),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
// Construct from components
template <class MeshType>
coupledInfo<MeshType>::coupledInfo
(
const MeshType& mesh,
const bool isTwoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label mPatch,
const label sPatch,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_
(
IOobject
(
"coupleMap_"
+ Foam::name(mPatch)
+ "_To_"
+ Foam::name(sPatch)
+ word(isLocal ? "_Local" : "_Proc")
+ word(isSend ? "_Send" : "_Recv"),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
true
),
isTwoDMesh,
isLocal,
isSend,
patchIndex,
mPatch,
sPatch
),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
//- Construct given addressing
template <class MeshType>
coupledInfo<MeshType>::subMeshMapper::subMeshMapper
(
const coupledInfo& cInfo,
const label patchI
)
:
sizeBeforeMapping_(cInfo.baseMesh().boundary()[patchI].size()),
directAddressing_
(
SubList<label>
(
cInfo.map().faceMap(),
cInfo.subMesh().boundary()[patchI].size(),
cInfo.subMesh().boundary()[patchI].patch().start()
)
)
{
// Offset indices
label pStart = cInfo.baseMesh().boundary()[patchI].patch().start();
forAll(directAddressing_, faceI)
{
directAddressing_[faceI] -= pStart;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Return a const reference to the parent mesh
template <class MeshType>
const MeshType&
coupledInfo<MeshType>::baseMesh() const
{
return mesh_;
}
// Set a new subMesh
template <class MeshType>
void coupledInfo<MeshType>::setMesh
(
label index,
MeshType* mesh
)
{
subMesh_.set(mesh);
}
// Return a reference to the subMesh
template <class MeshType>
MeshType& coupledInfo<MeshType>::subMesh()
{
if (!subMesh_.valid())
{
FatalErrorIn("MeshType& coupledInfo::subMesh()")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
// Return a const reference to the subMesh
template <class MeshType>
const MeshType& coupledInfo<MeshType>::subMesh() const
{
if (!subMesh_.valid())
{
FatalErrorIn("const MeshType& coupledInfo::subMesh() const")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
// Return if maps have been built
template <class MeshType>
bool coupledInfo<MeshType>::builtMaps() const
{
return builtMaps_;
}
// Set internal state of maps as built
template <class MeshType>
void coupledInfo<MeshType>::setBuiltMaps()
{
builtMaps_ = true;
}
// Return a reference to the coupleMap
template <class MeshType>
coupleMap& coupledInfo<MeshType>::map()
{
return map_;
}
// Return a const reference to the coupleMap
template <class MeshType>
const coupleMap& coupledInfo<MeshType>::map() const
{
return map_;
}
// Return the master face zone ID
template <class MeshType>
label coupledInfo<MeshType>::masterFaceZone() const
{
return masterFaceZone_;
}
// Return the slave face zone ID
template <class MeshType>
label coupledInfo<MeshType>::slaveFaceZone() const
{
return slaveFaceZone_;
}
// Subset geometric field
template <class MeshType>
template <class GeomField, class ZeroType>
tmp<GeomField>
coupledInfo<MeshType>::subSetField
(
const GeomField& f,
const ZeroType& zeroValue,
const labelList& internalMapper
) const
{
typedef typename GeomField::InternalField InternalField;
typedef typename GeomField::PatchFieldType PatchFieldType;
typedef typename GeomField::GeometricBoundaryField GeomBdyFieldType;
typedef typename GeomField::DimensionedInternalField DimInternalField;
// Create and map the internal-field values
InternalField internalField(f.internalField(), internalMapper);
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
PtrList<PatchFieldType> patchFields(nPatches);
// Define patch type names, assumed to be
// common for volume and surface fields
word emptyType(emptyPolyPatch::typeName);
word processorType(processorPolyPatch::typeName);
// Create dummy types for initial field creation
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
PatchFieldType::New
(
emptyType,
subMesh().boundary()[patchI],
DimInternalField::null()
)
);
}
else
{
patchFields.set
(
patchI,
PatchFieldType::New
(
PatchFieldType::calculatedType(),
subMesh().boundary()[patchI],
DimInternalField::null()
)
);
}
}
// Create new field from pieces
tmp<GeomField> subFld
(
new GeomField
(
IOobject
(
"subField_" + f.name(),
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
f.dimensions(),
internalField,
patchFields
)
);
// Set correct references for patch internal fields,
// and map values from the supplied geometric field
GeomBdyFieldType& bf = subFld().boundaryField();
forAll(bf, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
bf.set
(
patchI,
PatchFieldType::New
(
emptyType,
subMesh().boundary()[patchI],
subFld().dimensionedInternalField()
)
);
}
else
if (isA<processorPolyPatch>(subMesh().boundary()[patchI].patch()))
{
bf.set
(
patchI,
PatchFieldType::New
(
processorType,
subMesh().boundary()[patchI],
subFld().dimensionedInternalField()
)
);
// Avoid dealing with uninitialised values
// by artificially assigning to zero
bf[patchI] == zeroValue;
}
else
{
bf.set
(
patchI,
PatchFieldType::New
(
f.boundaryField()[patchI],
subMesh().boundary()[patchI],
subFld().dimensionedInternalField(),
subMeshMapper(*this, patchI)
)
);
}
}
return subFld;
}
// Subset geometric fields from registry to output stream
template <class MeshType>
template <class GeomField, class ZeroType>
void coupledInfo<MeshType>::send
(
const wordList& fieldNames,
const word& fieldType,
const ZeroType& zeroValue,
const labelList& internalMapper,
OSstream& strStream
) const
{
strStream
<< fieldType << token::NL
<< token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
// Fetch object from registry
const objectRegistry& db = mesh_.thisDb();
const GeomField& fld = db.lookupObject<GeomField>(fieldNames[i]);
// Subset the field
tmp<GeomField> tsubFld = subSetField(fld, zeroValue, internalMapper);
// Send field subset through stream
strStream
<< fieldNames[i]
<< token::NL << token::BEGIN_BLOCK
<< tsubFld
<< token::NL << token::END_BLOCK
<< token::NL;
}
strStream
<< token::END_BLOCK << token::NL;
}
// Set geometric field pointers from input dictionary
template <class MeshType>
template <class GeomField>
void coupledInfo<MeshType>::setField
(
const wordList& fieldNames,
const dictionary& fieldDicts,
const label internalSize,
PtrList<GeomField>& fields
) const
{
typedef typename GeomField::InternalField InternalField;
typedef typename GeomField::PatchFieldType PatchFieldType;
typedef typename GeomField::GeometricBoundaryField GeomBdyFieldType;
typedef typename GeomField::DimensionedInternalField DimInternalField;
// Size up the pointer list
fields.setSize(fieldNames.size());
// Define patch type names, assumed to be
// common for volume and surface fields
word emptyType(emptyPolyPatch::typeName);
word processorType(processorPolyPatch::typeName);
forAll(fieldNames, i)
{
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
// Create field parts
PtrList<PatchFieldType> patchFields(nPatches);
// Read dimensions
dimensionSet dimSet
(
fieldDicts.subDict(fieldNames[i]).lookup("dimensions")
);
// Read the internal field
InternalField internalField
(
"internalField",
fieldDicts.subDict(fieldNames[i]),
internalSize
);
// Create dummy types for initial field creation
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
PatchFieldType::New
(
emptyType,
subMesh().boundary()[patchI],
DimInternalField::null()
)
);
}
else
{
patchFields.set
(
patchI,
PatchFieldType::New
(
PatchFieldType::calculatedType(),
subMesh().boundary()[patchI],
DimInternalField::null()
)
);
}
}
// Create field with dummy patches
fields.set
(
i,
new GeomField
(
IOobject
(
fieldNames[i],
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
dimSet,
internalField,
patchFields
)
);
// Set correct references for patch internal fields,
// and fetch values from the supplied geometric field dictionaries
GeomBdyFieldType& bf = fields[i].boundaryField();
forAll(bf, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
bf.set
(
patchI,
PatchFieldType::New
(
emptyType,
subMesh().boundary()[patchI],
fields[i].dimensionedInternalField()
)
);
}
else
if (isA<processorPolyPatch>(subMesh().boundary()[patchI].patch()))
{
bf.set
(
patchI,
PatchFieldType::New
(
processorType,
subMesh().boundary()[patchI],
fields[i].dimensionedInternalField()
)
);
}
else
{
bf.set
(
patchI,
PatchFieldType::New
(
subMesh().boundary()[patchI],
fields[i].dimensionedInternalField(),
fieldDicts.subDict
(
fieldNames[i]
).subDict("boundaryField").subDict
(
subMesh().boundary()[patchI].name()
)
)
);
}
}
}
}
// Resize map for individual field
template <class MeshType>
template <class GeomField>
void coupledInfo<MeshType>::resizeMap
(
const label srcIndex,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields,
GeomField& field
)
{
// autoMap the internal field
field.internalField().autoMap(internalMapper);
// Reverse map for additional cells
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.internalField().rmap
(
srcField.internalField(),
internalReverseMaps[pI]
);
}
// Map physical boundary-fields
forAll(boundaryMapper, patchI)
{
// autoMap the patchField
field.boundaryField()[patchI].autoMap(boundaryMapper[patchI]);
// Reverse map for additional patch faces
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.boundaryField()[patchI].rmap
(
srcField.boundaryField()[patchI],
boundaryReverseMaps[pI][patchI]
);
}
}
}
// Resize all fields in registry
template <class MeshType>
template <class GeomField>
void coupledInfo<MeshType>::resizeMap
(
const wordList& names,
const objectRegistry& mesh,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields
)
{
forAll(names, indexI)
{
// Fetch field from registry
GeomField& field =
(
const_cast<GeomField&>
(
mesh.lookupObject<GeomField>(names[indexI])
)
);
// Map the field
coupledInfo<MeshType>::resizeMap
(
indexI,
internalMapper,
internalReverseMaps,
boundaryMapper,
boundaryReverseMaps,
srcFields,
field
);
}
}
// Resize boundaryFields for all fields in the registry
template <class MeshType>
template <class GeomField>
void coupledInfo<MeshType>::resizeBoundaries
(
const objectRegistry& mesh,
const fvBoundaryMesh& boundary
)
{
typedef typename GeomField::PatchFieldType PatchFieldType;
typedef typename GeomField::GeometricBoundaryField GeomBoundaryType;
HashTable<const GeomField*> fields(mesh.lookupClass<GeomField>());
forAllIter(typename HashTable<const GeomField*>, fields, fIter)
{
// Fetch field from registry
GeomField& field = const_cast<GeomField&>(*fIter());
GeomBoundaryType& bf = field.boundaryField();
// Resize boundary
label nPatches = boundary.size();
label nOldPatches = field.boundaryField().size();
// Create a new list of boundaries
PtrList<PatchFieldType> newbf(nPatches);
// Existing fields are mapped with new fvBoundaryMesh references
for (label patchI = 0; patchI < nOldPatches; patchI++)
{
label oldPatchSize = bf[patchI].size();
newbf.set
(
patchI,
PatchFieldType::New
(
bf[patchI],
boundary[patchI],
field,
subMeshMapper(oldPatchSize, identity(oldPatchSize))
)
);
}
// Size up new patches
for (label patchI = nOldPatches; patchI < nPatches; patchI++)
{
newbf.set
(
patchI,
PatchFieldType::New
(
boundary[patchI].type(),
boundary[patchI],
field
)
);
}
// Transfer contents with new patches
bf.transfer(newbf);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// Disallow default bitwise assignment
template <class MeshType>
void coupledInfo<MeshType>::operator=(const coupledInfo& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn
(
"void coupledInfo::operator=(const Foam::coupledInfo&)"
)
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -43,26 +43,29 @@ SourceFiles
#include "autoPtr.H" #include "autoPtr.H"
#include "coupleMap.H" #include "coupleMap.H"
#include "volFieldsFwd.H" #include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
#include "fvPatchFieldMapper.H" #include "fvPatchFieldMapper.H"
namespace Foam namespace Foam
{ {
// Class forward declarations // Class forward declarations
class dynamicTopoFvMesh; class fvBoundaryMesh;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class coupledInfo Declaration Class coupledInfo Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
template <class MeshType>
class coupledInfo class coupledInfo
{ {
// Private data
// Reference to the parent mesh // Reference to the parent mesh
const dynamicTopoFvMesh& mesh_; const MeshType& mesh_;
// Auto pointer to a subMesh // Auto pointer to a subMesh
autoPtr<dynamicTopoFvMesh> subMesh_; autoPtr<MeshType> subMesh_;
// Flag to determine whether maps have been built. // Flag to determine whether maps have been built.
bool builtMaps_; bool builtMaps_;
@ -75,23 +78,28 @@ class coupledInfo
label masterFaceZone_; label masterFaceZone_;
label slaveFaceZone_; label slaveFaceZone_;
//- Disallow default bitwise assignment // Private member functions
// Disallow default bitwise assignment
inline void operator=(const coupledInfo&); inline void operator=(const coupledInfo&);
public: public:
// Constructor //- Constructors
// Construct given mesh, coupleMap and master / slave indices
inline coupledInfo inline coupledInfo
( (
const dynamicTopoFvMesh& mesh, const MeshType& mesh,
const coupleMap& cMap, const coupleMap& cMap,
const label mfzIndex = -1, const label mfzIndex = -1,
const label sfzIndex = -1 const label sfzIndex = -1
); );
// Construct from components
inline coupledInfo inline coupledInfo
( (
const dynamicTopoFvMesh& mesh, const MeshType& mesh,
const bool isTwoDMesh, const bool isTwoDMesh,
const bool isLocal, const bool isLocal,
const bool isSend, const bool isSend,
@ -105,21 +113,21 @@ public:
//- Access //- Access
// Return a const reference to the parent mesh // Return a const reference to the parent mesh
inline const dynamicTopoFvMesh& baseMesh() const; inline const MeshType& baseMesh() const;
// Set a new subMesh // Set a new subMesh
inline void setMesh(label index, dynamicTopoFvMesh* mesh); inline void setMesh(label index, MeshType* mesh);
// Return a reference to the subMesh // Return a reference to the subMesh
inline dynamicTopoFvMesh& subMesh(); inline MeshType& subMesh();
// Return a const reference to the subMesh // Return a const reference to the subMesh
inline const dynamicTopoFvMesh& subMesh() const; inline const MeshType& subMesh() const;
// Have maps been built? // Return if maps have been built
inline bool builtMaps() const; inline bool builtMaps() const;
// Change the flag // Set internal state of maps as built
inline void setBuiltMaps(); inline void setBuiltMaps();
// Return a reference to the coupleMap // Return a reference to the coupleMap
@ -128,8 +136,10 @@ public:
// Return a const reference to the coupleMap // Return a const reference to the coupleMap
inline const coupleMap& map() const; inline const coupleMap& map() const;
// Return the master / slave face zone IDs // Return the master face zone ID
inline label masterFaceZone() const; inline label masterFaceZone() const;
// Return the slave face zone ID
inline label slaveFaceZone() const; inline label slaveFaceZone() const;
//- Interpolation //- Interpolation
@ -179,7 +189,7 @@ public:
label sizeBeforeMapping() const label sizeBeforeMapping() const
{ {
return directAddressing_.size(); return sizeBeforeMapping_;
} }
bool direct() const bool direct() const
@ -193,49 +203,34 @@ public:
} }
}; };
// Set subMesh centres // Subset geometric field
inline void setCentres(PtrList<volVectorField>& centres) const; template <class GeomField, class ZeroType>
tmp<GeomField>
// Subset volume field subSetField
template <class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
subSetVolField
( (
const GeometricField<Type, fvPatchField, volMesh>& field const GeomField& f,
const ZeroType& zeroValue,
const labelList& internalMapper
) const; ) const;
// Subset surface field // Subset geometric fields from registry to output stream
template <class Type> template <class GeomField, class ZeroType>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > void send
subSetSurfaceField
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& field
) const;
// Subset the volume fields from registry to output stream
template <class Type>
void mapVolField
( (
const wordList& fieldNames, const wordList& fieldNames,
const word& fieldType, const word& fieldType,
const ZeroType& zeroValue,
const labelList& internalMapper,
OSstream& strStream OSstream& strStream
) const; ) const;
// Subset the surface field from registry to output stream // Set geometric field pointers from input dictionary
template <class Type>
void mapSurfaceField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const;
// Set volume field pointer from input dictionary
template <class GeomField> template <class GeomField>
void setField void setField
( (
const wordList& fieldNames, const wordList& fieldNames,
const dictionary& fieldDicts, const dictionary& fieldDicts,
const label internalSize,
PtrList<GeomField>& fields PtrList<GeomField>& fields
) const; ) const;
@ -264,6 +259,14 @@ public:
const List<labelListList>& boundaryReverseMaps, const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields const List<PtrList<GeomField> >& srcFields
); );
// Resize boundaryFields for all fields in the registry
template <class GeomField>
static void resizeBoundaries
(
const objectRegistry& mesh,
const fvBoundaryMesh& boundary
);
}; };
} // End namespace Foam } // End namespace Foam

View file

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
subMeshLduAddressing
Description
Customized lduAddressing for subMeshes
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#ifndef subMeshLduAddressing_H
#define subMeshLduAddressing_H
#include "lduAddressing.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class subMeshLduAddressing Declaration
\*---------------------------------------------------------------------------*/
class subMeshLduAddressing
:
public lduAddressing
{
// Private data
//- Lower as a subList of allOwner
labelList::subList lowerAddr_;
//- Upper as a reference to neighbour
const labelList& upperAddr_;
//- Number of patches
label nPatches_;
// Private Member Functions
//- Disallow default bitwise copy construct
subMeshLduAddressing(const subMeshLduAddressing&);
//- Disallow default bitwise assignment
void operator=(const subMeshLduAddressing&);
public:
// Constructors
//- Construct from components
subMeshLduAddressing(const fvMesh& mesh)
:
lduAddressing(mesh.nCells()),
lowerAddr_
(
labelList::subList
(
mesh.faceOwner(),
mesh.nInternalFaces()
)
),
upperAddr_(mesh.faceNeighbour()),
nPatches_(mesh.boundary().size())
{}
// Destructor
virtual ~subMeshLduAddressing()
{}
// Member Functions
//- Return number of interfaces
virtual label nPatches() const
{
return nPatches_;
}
//- Return lower addressing (i.e. lower label = upper triangle)
virtual const unallocLabelList& lowerAddr() const
{
return lowerAddr_;
}
//- Return upper addressing (i.e. upper label)
virtual const unallocLabelList& upperAddr() const
{
return upperAddr_;
}
//- Return patch addressing
virtual const unallocLabelList& patchAddr(const label i) const
{
FatalErrorIn
(
"void subMeshLduAddressing::patchAddr(const label i)"
)
<< " Illegal request. "
<< abort(FatalError);
return unallocLabelList::null();
}
// Return patch field evaluation schedule
virtual const lduSchedule& patchSchedule() const
{
FatalErrorIn
(
"void subMeshLduAddressing::patchSchedule()"
)
<< " Illegal request. "
<< abort(FatalError);
return lduSchedule::null();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,72 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
subMeshProcessorFvPatch
Description
Functions specific to subMeshProcessorFvPatch
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "subMeshProcessorFvPatch.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(subMeshProcessorFvPatch, 0);
addToRunTimeSelectionTable(fvPatch, subMeshProcessorFvPatch, polyPatch);
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void subMeshProcessorFvPatch::makeWeights(scalarField& w) const
{
w = 1.0;
}
void subMeshProcessorFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
dc = 1.0/(nf() & fvPatch::delta());
}
tmp<vectorField> subMeshProcessorFvPatch::delta() const
{
return fvPatch::delta();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,110 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
subMeshProcessorFvPatch
Description
Customized processor type for subMeshes.
The intent of this class is to provide processor patch information
on subMeshes, but avoiding all forms of communication.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#ifndef subMeshProcessorFvPatch_H
#define subMeshProcessorFvPatch_H
#include "processorFvPatch.H"
#include "subMeshProcessorPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class subMeshProcessorFvPatch Declaration
\*---------------------------------------------------------------------------*/
class subMeshProcessorFvPatch
:
public processorFvPatch
{
// Private Data
const subMeshProcessorPolyPatch& procPolyPatch_;
protected:
// Protected Member Functions
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
public:
//- Runtime type information
TypeName(subMeshProcessorPolyPatch::typeName_());
// Constructors
//- Construct from components
subMeshProcessorFvPatch
(
const polyPatch& patch,
const fvBoundaryMesh& bm
)
:
processorFvPatch(patch, bm),
procPolyPatch_(refCast<const subMeshProcessorPolyPatch>(patch))
{}
// Destructor
virtual ~subMeshProcessorFvPatch()
{}
// Member functions
//- Return delta (P to N) vectors across coupled patch
virtual tmp<vectorField> delta() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
subMeshProcessorPolyPatch
Description
Functions specific to subMeshProcessorPolyPatch
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "subMeshProcessorPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(subMeshProcessorPolyPatch, 0);
addToRunTimeSelectionTable(polyPatch, subMeshProcessorPolyPatch, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
subMeshProcessorPolyPatch::subMeshProcessorPolyPatch
(
const word& name,
const label size,
const label start,
const label index,
const polyBoundaryMesh& bm,
const int myProcNo,
const int neighbProcNo
)
:
processorPolyPatch
(
name,
size,
start,
index,
bm,
myProcNo,
neighbProcNo
)
{}
subMeshProcessorPolyPatch::subMeshProcessorPolyPatch
(
const word& name,
const dictionary& dict,
const label index,
const polyBoundaryMesh& bm
)
:
processorPolyPatch(name, dict, index, bm)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
subMeshProcessorPolyPatch::~subMeshProcessorPolyPatch()
{}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,100 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
subMeshProcessorPolyPatch
Description
Customized processor type for subMeshes.
The intent of this class is to provide processor patch information
on subMeshes, but avoiding all forms of communication.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#ifndef subMeshProcessorPolyPatch_H
#define subMeshProcessorPolyPatch_H
#include "processorPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class subMeshProcessorPolyPatch Declaration
\*---------------------------------------------------------------------------*/
class subMeshProcessorPolyPatch
:
public processorPolyPatch
{
public:
//- Runtime type information
TypeName("subMeshProcessor");
// Constructors
//- Construct from components
subMeshProcessorPolyPatch
(
const word& name,
const label size,
const label start,
const label index,
const polyBoundaryMesh& bm,
const int myProcNo,
const int neighbProcNo
);
//- Construct from dictionary
subMeshProcessorPolyPatch
(
const word& name,
const dictionary& dict,
const label index,
const polyBoundaryMesh&
);
// Destructor
virtual ~subMeshProcessorPolyPatch();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -52,6 +52,7 @@ Author
#include "motionSolver.H" #include "motionSolver.H"
#include "fvPatchFields.H" #include "fvPatchFields.H"
#include "fvsPatchFields.H" #include "fvsPatchFields.H"
#include "subMeshLduAddressing.H"
#include "lengthScaleEstimator.H" #include "lengthScaleEstimator.H"
#include "conservativeMapFields.H" #include "conservativeMapFields.H"
@ -98,6 +99,7 @@ dynamicTopoFvMesh::dynamicTopoFvMesh(const IOobject& io)
loadMotionSolver_(true), loadMotionSolver_(true),
bandWidthReduction_(false), bandWidthReduction_(false),
coupledModification_(false), coupledModification_(false),
lduPtr_(NULL),
interval_(1), interval_(1),
eMeshPtr_(NULL), eMeshPtr_(NULL),
mapper_(NULL), mapper_(NULL),
@ -132,7 +134,7 @@ dynamicTopoFvMesh::dynamicTopoFvMesh(const IOobject& io)
nOldInternalEdges_(0), nOldInternalEdges_(0),
nInternalEdges_(0), nInternalEdges_(0),
maxModifications_(-1), maxModifications_(-1),
statistics_(0), statistics_(TOTAL_OP_TYPES, 0),
sliverThreshold_(0.1), sliverThreshold_(0.1),
slicePairs_(0), slicePairs_(0),
maxTetsPerEdge_(-1), maxTetsPerEdge_(-1),
@ -175,12 +177,6 @@ dynamicTopoFvMesh::dynamicTopoFvMesh(const IOobject& io)
// Load the field-mapper // Load the field-mapper
loadFieldMapper(); loadFieldMapper();
// Set sizes for the reverse maps
reversePointMap_.setSize(nPoints_, -7);
reverseEdgeMap_.setSize(nEdges_, -7);
reverseFaceMap_.setSize(nFaces_, -7);
reverseCellMap_.setSize(nCells_, -7);
// Load the length-scale estimator, // Load the length-scale estimator,
// and read refinement options // and read refinement options
loadLengthScaleEstimator(); loadLengthScaleEstimator();
@ -227,6 +223,7 @@ dynamicTopoFvMesh::dynamicTopoFvMesh
loadMotionSolver_(mesh.loadMotionSolver_), loadMotionSolver_(mesh.loadMotionSolver_),
bandWidthReduction_(mesh.bandWidthReduction_), bandWidthReduction_(mesh.bandWidthReduction_),
coupledModification_(false), coupledModification_(false),
lduPtr_(NULL),
interval_(1), interval_(1),
eMeshPtr_(NULL), eMeshPtr_(NULL),
mapper_(NULL), mapper_(NULL),
@ -261,7 +258,7 @@ dynamicTopoFvMesh::dynamicTopoFvMesh
nOldInternalEdges_(edgeStarts[0]), nOldInternalEdges_(edgeStarts[0]),
nInternalEdges_(edgeStarts[0]), nInternalEdges_(edgeStarts[0]),
maxModifications_(mesh.maxModifications_), maxModifications_(mesh.maxModifications_),
statistics_(0), statistics_(TOTAL_OP_TYPES, 0),
sliverThreshold_(mesh.sliverThreshold_), sliverThreshold_(mesh.sliverThreshold_),
slicePairs_(0), slicePairs_(0),
maxTetsPerEdge_(mesh.maxTetsPerEdge_), maxTetsPerEdge_(mesh.maxTetsPerEdge_),
@ -326,7 +323,7 @@ dynamicTopoFvMesh::dynamicTopoFvMesh
// Now build edgeFaces and pointEdges information. // Now build edgeFaces and pointEdges information.
edgeFaces_ = invertManyToMany<labelList, labelList>(nEdges_, faceEdges_); edgeFaces_ = invertManyToMany<labelList, labelList>(nEdges_, faceEdges_);
if (!twoDMesh_) if (is3D())
{ {
pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_); pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
} }
@ -336,7 +333,9 @@ dynamicTopoFvMesh::dynamicTopoFvMesh
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
dynamicTopoFvMesh::~dynamicTopoFvMesh() dynamicTopoFvMesh::~dynamicTopoFvMesh()
{} {
deleteDemandDrivenData(lduPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -469,6 +468,7 @@ label dynamicTopoFvMesh::insertFace
const face& newFace, const face& newFace,
const label newOwner, const label newOwner,
const label newNeighbour, const label newNeighbour,
const labelList& newFaceEdges,
const label zoneID const label zoneID
) )
{ {
@ -480,6 +480,7 @@ label dynamicTopoFvMesh::insertFace
faces_.append(newFace); faces_.append(newFace);
owner_.append(newOwner); owner_.append(newOwner);
neighbour_.append(newNeighbour); neighbour_.append(newNeighbour);
faceEdges_.append(newFaceEdges);
if (debug > 2) if (debug > 2)
{ {
@ -487,7 +488,8 @@ label dynamicTopoFvMesh::insertFace
<< newFaceIndex << ": " << newFaceIndex << ": "
<< newFace << newFace
<< " Owner: " << newOwner << " Owner: " << newOwner
<< " Neighbour: " << newNeighbour; << " Neighbour: " << newNeighbour
<< " faceEdges: " << newFaceEdges;
const polyBoundaryMesh& boundary = boundaryMesh(); const polyBoundaryMesh& boundary = boundaryMesh();
@ -510,7 +512,7 @@ label dynamicTopoFvMesh::insertFace
// Keep track of added boundary faces in a separate hash-table // Keep track of added boundary faces in a separate hash-table
// This information will be required at the reordering stage // This information will be required at the reordering stage
addedFacePatches_.insert(newFaceIndex,patch); addedFacePatches_.insert(newFaceIndex, patch);
if (newNeighbour == -1) if (newNeighbour == -1)
{ {
@ -814,7 +816,7 @@ label dynamicTopoFvMesh::insertEdge
} }
// Size-up the pointEdges list // Size-up the pointEdges list
if (!twoDMesh_) if (is3D())
{ {
meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[0]]); meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[0]]);
meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[1]]); meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[1]]);
@ -833,7 +835,7 @@ void dynamicTopoFvMesh::removeEdge
const label eIndex const label eIndex
) )
{ {
if (!twoDMesh_) if (is3D())
{ {
const edge& rEdge = edges_[eIndex]; const edge& rEdge = edges_[eIndex];
@ -978,7 +980,7 @@ label dynamicTopoFvMesh::insertPoint
// Add an empty entry to pointEdges as well. // Add an empty entry to pointEdges as well.
// This entry can be sized-up appropriately at a later stage. // This entry can be sized-up appropriately at a later stage.
if (!twoDMesh_) if (is3D())
{ {
pointEdges_.append(labelList(0)); pointEdges_.append(labelList(0));
} }
@ -1015,7 +1017,7 @@ void dynamicTopoFvMesh::removePoint
// oldPoints_[pIndex] = point(); // oldPoints_[pIndex] = point();
// Remove pointEdges as well // Remove pointEdges as well
if (!twoDMesh_) if (is3D())
{ {
pointEdges_[pIndex].clear(); pointEdges_[pIndex].clear();
} }
@ -1394,7 +1396,7 @@ void dynamicTopoFvMesh::handleMeshSlicing()
bool available = true; bool available = true;
if (twoDMesh_) if (is2D())
{ {
forAll(pairToCheck, indexI) forAll(pairToCheck, indexI)
{ {
@ -1544,7 +1546,7 @@ scalar dynamicTopoFvMesh::testProximity
scalar proxDistance = GREAT, testStep = 0.0; scalar proxDistance = GREAT, testStep = 0.0;
vector gCentre = vector::zero, gNormal = vector::zero; vector gCentre = vector::zero, gNormal = vector::zero;
if (twoDMesh_) if (is2D())
{ {
// Obtain the face-normal. // Obtain the face-normal.
gNormal = faces_[index].normal(points_); gNormal = faces_[index].normal(points_);
@ -1613,7 +1615,7 @@ scalar dynamicTopoFvMesh::testProximity
labelPair proxPoints(-1, -1); labelPair proxPoints(-1, -1);
bool foundPoint = false; bool foundPoint = false;
if (twoDMesh_) if (is2D())
{ {
proxPoints.first() = index; proxPoints.first() = index;
proxPoints.second() = proximityFace; proxPoints.second() = proximityFace;
@ -1744,16 +1746,58 @@ void dynamicTopoFvMesh::readOptionalParameters(bool reRead)
// Enable/disable run-time debug level // Enable/disable run-time debug level
if (meshSubDict.found("debug") || mandatory_) if (meshSubDict.found("debug") || mandatory_)
{ {
// Set an attachment point for debugging
if
(
!reRead &&
meshSubDict.found("parDebugAttach") &&
readBool(meshSubDict.lookup("parDebugAttach"))
)
{
label wait;
if (Pstream::master())
{
Info<< nl << "Waiting for debugger attachment..." << endl;
cin >> wait;
if (Pstream::parRun())
{
for (label i = 1; i < Pstream::nProcs(); i++)
{
meshOps::pWrite(i, wait);
}
}
}
else
{
meshOps::pRead(Pstream::masterNo(), wait);
}
}
// Look for a processor-specific debug flag // Look for a processor-specific debug flag
if (Pstream::parRun() && meshSubDict.found("parDebugProcs")) if
(
Pstream::parRun() &&
meshSubDict.found("parDebugProcs")
)
{
if (!reRead)
{ {
labelList procs(meshSubDict.lookup("parDebugProcs")); labelList procs(meshSubDict.lookup("parDebugProcs"));
label pdebug = readLabel(meshSubDict.lookup("debug"));
forAll(procs, procI) if (findIndex(procs, Pstream::myProcNo()) > -1)
{ {
if (Pstream::myProcNo() == procs[procI]) debug = pdebug;
}
else
if (pdebug)
{ {
debug = readLabel(meshSubDict.lookup("debug")); // Minimal debug information,
// so that synchronizations are successful
debug = 1;
} }
} }
} }
@ -1840,7 +1884,7 @@ void dynamicTopoFvMesh::readOptionalParameters(bool reRead)
} }
// For tetrahedral meshes... // For tetrahedral meshes...
if (!twoDMesh_) if (is3D())
{ {
// Check if swapping is to be avoided on any patches // Check if swapping is to be avoided on any patches
if (meshSubDict.found("noSwapPatches") || mandatory_) if (meshSubDict.found("noSwapPatches") || mandatory_)
@ -1892,8 +1936,14 @@ void dynamicTopoFvMesh::readOptionalParameters(bool reRead)
// Check for load-balancing in parallel // Check for load-balancing in parallel
if (reRead && (meshSubDict.found("loadBalancing") || mandatory_)) if (reRead && (meshSubDict.found("loadBalancing") || mandatory_))
{ {
// Fetch non-const reference
dictionary& balanceDict =
(
dict_.subDict("dynamicTopoFvMesh").subDict("loadBalancing")
);
// Execute balancing // Execute balancing
executeLoadBalancing(meshSubDict.subDict("loadBalancing")); executeLoadBalancing(balanceDict);
} }
} }
@ -1921,7 +1971,7 @@ void dynamicTopoFvMesh::initEdges()
edgeFaces_ = eMeshPtr_->edgeFaces(); edgeFaces_ = eMeshPtr_->edgeFaces();
faceEdges_ = eMeshPtr_->faceEdges(); faceEdges_ = eMeshPtr_->faceEdges();
if (!twoDMesh_) if (is3D())
{ {
// Invert edges to obtain pointEdges // Invert edges to obtain pointEdges
pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_); pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
@ -1938,7 +1988,7 @@ void dynamicTopoFvMesh::initEdges()
// Load the mesh-quality metric from the library // Load the mesh-quality metric from the library
void dynamicTopoFvMesh::loadMetric() void dynamicTopoFvMesh::loadMetric()
{ {
if (twoDMesh_) if (is2D())
{ {
return; return;
} }
@ -2169,7 +2219,7 @@ void dynamicTopoFvMesh::swap2DEdges(void *argument)
); );
Info<< " Swap Progress: " << percent << "% :" Info<< " Swap Progress: " << percent << "% :"
<< " Total: " << mesh.status(1) << " Total: " << mesh.status(TOTAL_SWAPS)
<< " \r" << " \r"
<< flush; << flush;
@ -2206,7 +2256,7 @@ void dynamicTopoFvMesh::swap2DEdges(void *argument)
if (reported) if (reported)
{ {
Info<< " Swap Progress: 100% :" Info<< " Swap Progress: 100% :"
<< " Total: " << mesh.status(1) << " Total: " << mesh.status(TOTAL_SWAPS)
<< " \r" << " \r"
<< endl; << endl;
} }
@ -2275,8 +2325,8 @@ void dynamicTopoFvMesh::swap3DEdges
); );
Info<< " Swap Progress: " << percent << "% :" Info<< " Swap Progress: " << percent << "% :"
<< " Surface: " << mesh.status(2) << " Surface: " << mesh.status(SURFACE_SWAPS)
<< ", Total: " << mesh.status(1) << ", Total: " << mesh.status(TOTAL_SWAPS)
<< " \r" << " \r"
<< flush; << flush;
@ -2333,8 +2383,8 @@ void dynamicTopoFvMesh::swap3DEdges
if (reported) if (reported)
{ {
Info<< " Swap Progress: 100% :" Info<< " Swap Progress: 100% :"
<< " Surface: " << mesh.status(2) << " Surface: " << mesh.status(SURFACE_SWAPS)
<< ", Total: " << mesh.status(1) << ", Total: " << mesh.status(TOTAL_SWAPS)
<< " \r" << " \r"
<< endl; << endl;
} }
@ -2395,9 +2445,9 @@ void dynamicTopoFvMesh::edgeRefinementEngine
); );
Info<< " Refinement Progress: " << percent << "% :" Info<< " Refinement Progress: " << percent << "% :"
<< " Bisections: " << mesh.status(3) << " Bisections: " << mesh.status(TOTAL_BISECTIONS)
<< ", Collapses: " << mesh.status(4) << ", Collapses: " << mesh.status(TOTAL_COLLAPSES)
<< ", Total: " << mesh.status(0) << ", Total: " << mesh.status(TOTAL_MODIFICATIONS)
<< " \r" << " \r"
<< flush; << flush;
@ -2445,9 +2495,9 @@ void dynamicTopoFvMesh::edgeRefinementEngine
if (reported) if (reported)
{ {
Info<< " Refinement Progress: 100% :" Info<< " Refinement Progress: 100% :"
<< " Bisections: " << mesh.status(3) << " Bisections: " << mesh.status(TOTAL_BISECTIONS)
<< ", Collapses: " << mesh.status(4) << ", Collapses: " << mesh.status(TOTAL_COLLAPSES)
<< ", Total: " << mesh.status(0) << ", Total: " << mesh.status(TOTAL_MODIFICATIONS)
<< " \r" << " \r"
<< endl; << endl;
} }
@ -2620,7 +2670,7 @@ void dynamicTopoFvMesh::remove2DSlivers()
{ {
if (collapseQuadFace(firstFace).type() > 0) if (collapseQuadFace(firstFace).type() > 0)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
continue; continue;
@ -2630,7 +2680,7 @@ void dynamicTopoFvMesh::remove2DSlivers()
{ {
if (collapseQuadFace(secondFace).type() > 0) if (collapseQuadFace(secondFace).type() > 0)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
continue; continue;
@ -2654,7 +2704,7 @@ void dynamicTopoFvMesh::remove2DSlivers()
{ {
if (collapseQuadFace(aF[faceI].index()).type() > 0) if (collapseQuadFace(aF[faceI].index()).type() > 0)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
break; break;
@ -2682,7 +2732,7 @@ void dynamicTopoFvMesh::remove2DSlivers()
{ {
if (collapseQuadFace(aF[faceI].index()).type() > 0) if (collapseQuadFace(aF[faceI].index()).type() > 0)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
break; break;
@ -2710,7 +2760,7 @@ void dynamicTopoFvMesh::remove2DSlivers()
{ {
if (collapseQuadFace(aF[faceI].index()).type() > 0) if (collapseQuadFace(aF[faceI].index()).type() > 0)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
break; break;
@ -3125,7 +3175,7 @@ void dynamicTopoFvMesh::removeSlivers()
} }
// Invoke the 2D sliver removal routine // Invoke the 2D sliver removal routine
if (twoDMesh_) if (is2D())
{ {
remove2DSlivers(); remove2DSlivers();
return; return;
@ -3185,7 +3235,7 @@ void dynamicTopoFvMesh::removeSlivers()
{ {
label proc = procIndices_[procI]; label proc = procIndices_[procI];
if (proc < Pstream::myProcNo()) if (priority(proc, lessOp<label>(), Pstream::myProcNo()))
{ {
Map<label>& rCellMap = Map<label>& rCellMap =
( (
@ -3334,13 +3384,13 @@ void dynamicTopoFvMesh::removeSlivers()
case 2: case 2:
{ {
// Cap cell. // Cap cell.
label opposingFace = readLabel(map.lookup("opposingFace")); //label opposingFace = readLabel(map.lookup("opposingFace"));
// Force trisection of the opposing face. // Force trisection of the opposing face.
changeMap faceMap = changeMap faceMap;
( //(
trisectFace(opposingFace, false, true) // trisectFace(opposingFace, false, true)
); //);
// Collapse the intermediate edge. // Collapse the intermediate edge.
// Since we don't know which edge it is, search // Since we don't know which edge it is, search
@ -3476,7 +3526,7 @@ void dynamicTopoFvMesh::removeSlivers()
// Increment the count for successful sliver removal // Increment the count for successful sliver removal
if (success) if (success)
{ {
statistics_[7]++; status(TOTAL_SLIVERS)++;
} }
} }
@ -3813,6 +3863,12 @@ scalar dynamicTopoFvMesh::edgeLengthScale
// MultiThreaded topology modifier // MultiThreaded topology modifier
void dynamicTopoFvMesh::threadedTopoModifier() void dynamicTopoFvMesh::threadedTopoModifier()
{ {
// Set sizes for the reverse maps
reversePointMap_.setSize(nPoints_, -7);
reverseEdgeMap_.setSize(nEdges_, -7);
reverseFaceMap_.setSize(nFaces_, -7);
reverseCellMap_.setSize(nCells_, -7);
// Remove sliver cells first. // Remove sliver cells first.
removeSlivers(); removeSlivers();
@ -3863,7 +3919,7 @@ void dynamicTopoFvMesh::threadedTopoModifier()
// Execute threads // Execute threads
if (threader_->multiThreaded()) if (threader_->multiThreaded())
{ {
if (twoDMesh_) if (is2D())
{ {
executeThreads(topoSequence, handlerPtr_, &swap2DEdges); executeThreads(topoSequence, handlerPtr_, &swap2DEdges);
} }
@ -3874,7 +3930,7 @@ void dynamicTopoFvMesh::threadedTopoModifier()
} }
// Set the master thread to implement modifications // Set the master thread to implement modifications
if (twoDMesh_) if (is2D())
{ {
swap2DEdges(&(handlerPtr_[0])); swap2DEdges(&(handlerPtr_[0]));
} }
@ -3901,31 +3957,36 @@ bool dynamicTopoFvMesh::resetMesh()
// Reduce across processors. // Reduce across processors.
reduce(topoChangeFlag_, orOp<bool>()); reduce(topoChangeFlag_, orOp<bool>());
const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
if (topoChangeFlag_) if (topoChangeFlag_)
{ {
// Write out statistics // Write out statistics
if (Pstream::parRun() && debug) if (Pstream::parRun())
{ {
Pout<< " Bisections :: Total: " << status(3) if (debug)
<< ", Surface: " << status(5) << nl
<< " Collapses :: Total: " << status(4)
<< ", Surface: " << status(6) << nl
<< " Swaps :: Total: " << status(1)
<< ", Surface: " << status(2) << endl;
}
else
{ {
Info<< " Bisections :: Total: " << status(3) Pout<< " Bisections :: Total: " << status(TOTAL_BISECTIONS)
<< ", Surface: " << status(5) << nl << ", Surface: " << status(SURFACE_BISECTIONS) << nl
<< " Collapses :: Total: " << status(4) << " Collapses :: Total: " << status(TOTAL_COLLAPSES)
<< ", Surface: " << status(6) << nl << ", Surface: " << status(SURFACE_COLLAPSES) << nl
<< " Swaps :: Total: " << status(1) << " Swaps :: Total: " << status(TOTAL_SWAPS)
<< ", Surface: " << status(2) << endl; << ", Surface: " << status(SURFACE_SWAPS) << endl;
} }
if (status(7)) reduce(statistics_, sumOp<labelList>());
}
Info<< " Bisections :: Total: " << status(TOTAL_BISECTIONS)
<< ", Surface: " << status(SURFACE_BISECTIONS) << nl
<< " Collapses :: Total: " << status(TOTAL_COLLAPSES)
<< ", Surface: " << status(SURFACE_COLLAPSES) << nl
<< " Swaps :: Total: " << status(TOTAL_SWAPS)
<< ", Surface: " << status(SURFACE_SWAPS) << endl;
if (status(TOTAL_SLIVERS))
{ {
Pout<< " Slivers :: " << status(7) << endl; Pout<< " Slivers :: " << status(TOTAL_SLIVERS) << endl;
} }
// Fetch reference to mapper // Fetch reference to mapper
@ -3935,6 +3996,17 @@ bool dynamicTopoFvMesh::resetMesh()
// - Must be done prior to field-transfers and mesh reset // - Must be done prior to field-transfers and mesh reset
fieldMapper.storeMeshInformation(); fieldMapper.storeMeshInformation();
// Obtain the number of patches before
// any possible boundary reset
label nOldPatches = boundaryMesh().size();
// If the number of patches have changed
// at run-time, reset boundaries first
if (nPatches_ != nOldPatches)
{
resetBoundaries();
}
// Set up field-transfers before dealing with mapping // Set up field-transfers before dealing with mapping
wordList fieldTypes; wordList fieldTypes;
List<wordList> fieldNames; List<wordList> fieldNames;
@ -3957,8 +4029,6 @@ bool dynamicTopoFvMesh::resetMesh()
// Determine if mapping is to be skipped // Determine if mapping is to be skipped
// Optionally skip mapping for remeshing-only / pre-processing // Optionally skip mapping for remeshing-only / pre-processing
const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
bool skipMapping = false; bool skipMapping = false;
if (meshSubDict.found("skipMapping") || mandatory_) if (meshSubDict.found("skipMapping") || mandatory_)
@ -4052,9 +4122,6 @@ bool dynamicTopoFvMesh::resetMesh()
<< reOrderingTimer.elapsedTime() << " s" << reOrderingTimer.elapsedTime() << " s"
<< endl; << endl;
// Obtain the number of patches before mesh reset
label nOldPatches = boundaryMesh().size();
// Obtain the patch-point maps before resetting the mesh // Obtain the patch-point maps before resetting the mesh
List<Map<label> > oldPatchPointMaps(nOldPatches); List<Map<label> > oldPatchPointMaps(nOldPatches);
@ -4077,13 +4144,6 @@ bool dynamicTopoFvMesh::resetMesh()
xferMove(cellCentres_) xferMove(cellCentres_)
); );
// If the number of patches have changed
// at run-time, reset boundaries first
if (nPatches_ != nOldPatches)
{
resetBoundaries();
}
// Reset the mesh, and specify a non-valid // Reset the mesh, and specify a non-valid
// boundary to avoid globalData construction // boundary to avoid globalData construction
polyMesh::resetPrimitives polyMesh::resetPrimitives
@ -4276,11 +4336,11 @@ bool dynamicTopoFvMesh::resetMesh()
// Clear flipFaces // Clear flipFaces
flipFaces_.clear(); flipFaces_.clear();
// Set new sizes for the reverse maps // Clear reverse maps
reversePointMap_.setSize(nPoints_, -7); reversePointMap_.clear();
reverseEdgeMap_.setSize(nEdges_, -7); reverseEdgeMap_.clear();
reverseFaceMap_.setSize(nFaces_, -7); reverseFaceMap_.clear();
reverseCellMap_.setSize(nCells_, -7); reverseCellMap_.clear();
// Update "old" information // Update "old" information
nOldPoints_ = nPoints_; nOldPoints_ = nPoints_;
@ -4303,6 +4363,7 @@ bool dynamicTopoFvMesh::resetMesh()
if (Pstream::parRun()) if (Pstream::parRun())
{ {
procIndices_.clear(); procIndices_.clear();
procPriority_.clear();
sendMeshes_.clear(); sendMeshes_.clear();
recvMeshes_.clear(); recvMeshes_.clear();
} }
@ -4349,6 +4410,31 @@ bool dynamicTopoFvMesh::resetMesh()
// Dump length-scale to disk, if requested. // Dump length-scale to disk, if requested.
calculateLengthScale(true); calculateLengthScale(true);
// Optionally check if decomposition is to be written out
if
(
Pstream::parRun() &&
time().outputTime() &&
meshSubDict.found("writeDecomposition") &&
readBool(meshSubDict.lookup("writeDecomposition"))
)
{
volScalarField
(
IOobject
(
"decomposition",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
dimensionedScalar("rank", dimless, Pstream::myProcNo())
).write();
}
// Reset and return flag // Reset and return flag
if (topoChangeFlag_) if (topoChangeFlag_)
{ {
@ -4361,6 +4447,26 @@ bool dynamicTopoFvMesh::resetMesh()
} }
// Return custom lduAddressing
const lduAddressing& dynamicTopoFvMesh::lduAddr() const
{
// For subMeshes, avoid globalData construction
// due to lduAddressing::patchSchedule, using a
// customized approach.
if (isSubMesh_)
{
if (!lduPtr_)
{
lduPtr_ = new subMeshLduAddressing(*this);
}
return *lduPtr_;
}
return fvMesh::lduAddr();
}
// Update mesh corresponding to the given map // Update mesh corresponding to the given map
void dynamicTopoFvMesh::updateMesh(const mapPolyMesh& mpm) void dynamicTopoFvMesh::updateMesh(const mapPolyMesh& mpm)
{ {
@ -4375,14 +4481,8 @@ void dynamicTopoFvMesh::updateMesh(const mapPolyMesh& mpm)
// Delete oldPoints in polyMesh // Delete oldPoints in polyMesh
polyMesh::resetMotion(); polyMesh::resetMotion();
// Clear-out fvMesh geometry and addressing // Update fvMesh and polyMesh, and map all fields
fvMesh::clearOut(); fvMesh::updateMesh(mpm);
// Update polyMesh.
polyMesh::updateMesh(mpm);
// Map all fields
mapFields(mpm);
} }
@ -4409,6 +4509,10 @@ void dynamicTopoFvMesh::mapFields(const mapPolyMesh& mpm) const
// Set the mapPolyMesh object in the mapper // Set the mapPolyMesh object in the mapper
fieldMapper.setMapper(mpm); fieldMapper.setMapper(mpm);
// Deregister gradient fields and centres,
// but retain for mapping
fieldMapper.deregisterMeshInformation();
// Conservatively map scalar/vector volFields // Conservatively map scalar/vector volFields
conservativeMapVolFields<scalar>(fieldMapper); conservativeMapVolFields<scalar>(fieldMapper);
conservativeMapVolFields<vector>(fieldMapper); conservativeMapVolFields<vector>(fieldMapper);

View file

@ -65,10 +65,13 @@ class eMesh;
class Stack; class Stack;
class changeMap; class changeMap;
class objectMap; class objectMap;
class coupledInfo;
class motionSolver; class motionSolver;
class convexSetAlgorithm; class convexSetAlgorithm;
class lengthScaleEstimator; class lengthScaleEstimator;
class subMeshLduAddressing;
template <class MeshType>
class coupledInfo;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class dynamicTopoFvMesh Declaration Class dynamicTopoFvMesh Declaration
@ -78,6 +81,26 @@ class dynamicTopoFvMesh
: :
public dynamicFvMesh public dynamicFvMesh
{ {
// Private typedefs
typedef coupledInfo<dynamicTopoFvMesh> coupledMesh;
typedef threadHandler<dynamicTopoFvMesh> meshHandler;
// Private enumerants
enum opType
{
TOTAL_MODIFICATIONS = 0,
TOTAL_SWAPS = 1,
TOTAL_BISECTIONS = 2,
TOTAL_COLLAPSES = 3,
SURFACE_SWAPS = 4,
SURFACE_BISECTIONS = 5,
SURFACE_COLLAPSES = 6,
TOTAL_SLIVERS = 7,
TOTAL_OP_TYPES
};
// Private data // Private data
//- Reference to base mesh //- Reference to base mesh
@ -90,7 +113,7 @@ class dynamicTopoFvMesh
bool topoChangeFlag_; bool topoChangeFlag_;
//- Is this a subMesh? //- Is this a subMesh?
bool isSubMesh_; Switch isSubMesh_;
//- Dynamic mesh dictionary //- Dynamic mesh dictionary
IOdictionary dict_; IOdictionary dict_;
@ -100,7 +123,7 @@ class dynamicTopoFvMesh
Switch mandatory_; Switch mandatory_;
//- Mesh characteristics [2D/3D] //- Mesh characteristics [2D/3D]
Switch twoDMesh_; bool twoDMesh_;
//- Edge refinement switch //- Edge refinement switch
Switch edgeRefinement_; Switch edgeRefinement_;
@ -114,6 +137,9 @@ class dynamicTopoFvMesh
//- Coupled modification switch //- Coupled modification switch
mutable Switch coupledModification_; mutable Switch coupledModification_;
//- Sub-Mesh lduAddressing
mutable subMeshLduAddressing* lduPtr_;
//- Specify the re-meshing interval //- Specify the re-meshing interval
label interval_; label interval_;
@ -218,7 +244,7 @@ class dynamicTopoFvMesh
//- Run-time statistics //- Run-time statistics
label maxModifications_; label maxModifications_;
FixedList<label, 8> statistics_; labelList statistics_;
//- Sliver exudation //- Sliver exudation
scalar sliverThreshold_; scalar sliverThreshold_;
@ -240,8 +266,6 @@ class dynamicTopoFvMesh
autoPtr<IOmultiThreader> threader_; autoPtr<IOmultiThreader> threader_;
// Pointer-list of thread-handlers // Pointer-list of thread-handlers
typedef threadHandler<dynamicTopoFvMesh> meshHandler;
PtrList<meshHandler> handlerPtr_; PtrList<meshHandler> handlerPtr_;
// Entity mutexes used for synchronization // Entity mutexes used for synchronization
@ -249,14 +273,17 @@ class dynamicTopoFvMesh
FixedList<Mutex, 4> entityMutex_; FixedList<Mutex, 4> entityMutex_;
// Local coupled patch information // Local coupled patch information
PtrList<coupledInfo> patchCoupling_; PtrList<coupledMesh> patchCoupling_;
// List of processors that this sub-domain talks to // List of processors that this sub-domain talks to
labelList procIndices_; labelList procIndices_;
// Processor priority list
labelList procPriority_;
// Sub-Mesh pointers // Sub-Mesh pointers
PtrList<coupledInfo> sendMeshes_; PtrList<coupledMesh> sendMeshes_;
PtrList<coupledInfo> recvMeshes_; PtrList<coupledMesh> recvMeshes_;
// Private Member Functions // Private Member Functions
@ -334,6 +361,12 @@ class dynamicTopoFvMesh
// Initialize the threading environment // Initialize the threading environment
void initializeThreadingEnvironment(const label specThreads = -1); void initializeThreadingEnvironment(const label specThreads = -1);
// Return if mesh is 2D
inline bool is2D() const;
// Return if mesh is 3D
inline bool is3D() const;
// Return a non-const reference to the lengthScaleEstimator // Return a non-const reference to the lengthScaleEstimator
inline lengthScaleEstimator& lengthEstimator(); inline lengthScaleEstimator& lengthEstimator();
@ -409,8 +442,8 @@ class dynamicTopoFvMesh
// Method to determine the boundary patch index for a given edge // Method to determine the boundary patch index for a given edge
inline label whichEdgePatch(const label index) const; inline label whichEdgePatch(const label index) const;
// Report topo-modification status // Report or alter topo-modification status
inline label status(const label type) const; inline label& status(const label type);
// Method for the swapping of a quad-face in 2D // Method for the swapping of a quad-face in 2D
const changeMap const changeMap
@ -448,15 +481,6 @@ class dynamicTopoFvMesh
bool forceOp = false bool forceOp = false
); );
// Method for the trisection of a face in 3D
const changeMap
trisectFace
(
const label fIndex,
bool checkOnly = false,
bool forceOp = false
);
// Method for the collapse of an edge in 3D // Method for the collapse of an edge in 3D
const changeMap const changeMap
collapseEdge collapseEdge
@ -560,6 +584,7 @@ class dynamicTopoFvMesh
const face& newFace, const face& newFace,
const label newOwner, const label newOwner,
const label newNeighbour, const label newNeighbour,
const labelList& newFaceEdges,
const label zoneID = -1 const label zoneID = -1
); );
@ -648,10 +673,6 @@ class dynamicTopoFvMesh
// around an edge after bisection. // around an edge after bisection.
scalar computeBisectionQuality(const label eIndex) const; scalar computeBisectionQuality(const label eIndex) const;
// Utility method to compute the quality of cells
// around a face after trisection.
scalar computeTrisectionQuality(const label fIndex) const;
// Check whether the given edge is on a bounding curve // Check whether the given edge is on a bounding curve
bool checkBoundingCurve bool checkBoundingCurve
( (
@ -945,6 +966,9 @@ class dynamicTopoFvMesh
PtrList<labelListList>& triangulations PtrList<labelListList>& triangulations
) const; ) const;
// Initialize coupled weights calculation
void initCoupledWeights();
// Additional mapping contributions for coupled entities // Additional mapping contributions for coupled entities
void computeCoupledWeights void computeCoupledWeights
( (
@ -986,7 +1010,7 @@ class dynamicTopoFvMesh
void readOptionalParameters(bool reRead = false); void readOptionalParameters(bool reRead = false);
// Execute load balancing operations on the mesh // Execute load balancing operations on the mesh
void executeLoadBalancing(const dictionary& balanceDict); void executeLoadBalancing(dictionary& balanceDict);
// Identify coupled patches. // Identify coupled patches.
bool identifyCoupledPatches(); bool identifyCoupledPatches();
@ -1035,7 +1059,7 @@ class dynamicTopoFvMesh
// Build patch sub-mesh for a specified processor // Build patch sub-mesh for a specified processor
void buildProcessorPatchMesh void buildProcessorPatchMesh
( (
coupledInfo& subMesh, coupledMesh& subMesh,
Map<labelList>& commonCells Map<labelList>& commonCells
); );
@ -1089,12 +1113,24 @@ class dynamicTopoFvMesh
labelListList& rotations labelListList& rotations
) const; ) const;
// Set processor rank priority
void initProcessorPriority();
// Check for processor priority
template <class BinaryOp>
inline bool priority
(
const label first,
const BinaryOp& bop,
const label second
) const;
// Dump cell-quality statistics // Dump cell-quality statistics
bool meshQuality(bool outputOption); bool meshQuality(bool outputOption);
public: public:
// Declare the name of the class and its debug switch //- Runtime type information
TypeName("dynamicTopoFvMesh"); TypeName("dynamicTopoFvMesh");
// Constructors // Constructors
@ -1108,14 +1144,17 @@ public:
// Member Functions // Member Functions
//- Update mesh corresponding to the given map // Return custom lduAddressing
virtual const lduAddressing& lduAddr() const;
// Update mesh corresponding to the given map
virtual void updateMesh(const mapPolyMesh& mpm); virtual void updateMesh(const mapPolyMesh& mpm);
// Map all fields in time using a customized mapper // Map all fields in time using a customized mapper
virtual void mapFields(const mapPolyMesh& mpm) const; virtual void mapFields(const mapPolyMesh& mpm) const;
// Update the mesh for motion / topology changes // Update the mesh for motion / topology changes
// Return true if topology changes have occurred // - Return true if topology changes have occurred
virtual bool update(); virtual bool update();
}; };

View file

@ -57,8 +57,7 @@ bool dynamicTopoFvMesh::meshQuality
label nCells = 0, minCell = -1; label nCells = 0, minCell = -1;
scalar maxQuality = -GREAT; scalar maxQuality = -GREAT;
scalar minQuality = GREAT; scalar minQuality = GREAT;
scalar cQuality = 0.0; scalar cQuality, meanQuality = 0.0;
scalar meanQuality = 0.0;
// Track slivers // Track slivers
bool sliversAbsent = true; bool sliversAbsent = true;
@ -89,7 +88,7 @@ bool dynamicTopoFvMesh::meshQuality
continue; continue;
} }
if (twoDMesh_) if (is2D())
{ {
// Assume XY plane here // Assume XY plane here
vector n = vector(0,0,1); vector n = vector(0,0,1);
@ -198,25 +197,20 @@ bool dynamicTopoFvMesh::meshQuality
// Output statistics: // Output statistics:
if (outputOption || (debug > 0)) if (outputOption || (debug > 0))
{ {
if (minQuality < 0.0)
{
Pout<< nl
<< " Warning: Minimum cell quality is: " << minQuality
<< " at cell: " << minCell
<< endl;
}
// Reduce statistics across processors. // Reduce statistics across processors.
reduce(minQuality, minOp<scalar>()); reduce(minQuality, minOp<scalar>());
reduce(maxQuality, maxOp<scalar>()); reduce(maxQuality, maxOp<scalar>());
reduce(meanQuality, sumOp<scalar>()); reduce(meanQuality, sumOp<scalar>());
reduce(nCells, sumOp<label>()); reduce(nCells, sumOp<label>());
if (minQuality < 0.0)
{
WarningIn
(
"bool dynamicTopoFvMesh::meshQuality"
"(bool outputOption)"
)
<< nl
<< "Minimum cell quality is: " << minQuality
<< " at cell: " << minCell
<< endl;
}
Info<< nl Info<< nl
<< " ~~~ Mesh Quality Statistics ~~~ " << nl << " ~~~ Mesh Quality Statistics ~~~ " << nl
<< " Min: " << minQuality << nl << " Min: " << minQuality << nl
@ -600,8 +594,7 @@ bool dynamicTopoFvMesh::checkTriangulationVolumes
) const ) const
{ {
label m = hullVertices.size(); label m = hullVertices.size();
scalar oldTetVol = 0.0; scalar oldTetVol = 0.0, newTetVol = 0.0;
scalar newTetVol = 0.0;
const edge& edgeToCheck = edges_[eIndex]; const edge& edgeToCheck = edges_[eIndex];
@ -760,7 +753,7 @@ void dynamicTopoFvMesh::writeEdgeConnectivity
3, false, true 3, false, true
); );
if (twoDMesh_) if (is2D())
{ {
return; return;
} }
@ -1132,6 +1125,25 @@ void dynamicTopoFvMesh::checkConnectivity(const label maxErrors) const
if (curFace.empty()) if (curFace.empty())
{ {
// This might be a deleted face
if (faceI < nOldFaces_)
{
if (reverseFaceMap_[faceI] == -1)
{
continue;
}
}
else
if (deletedFaces_.found(faceI))
{
continue;
}
Pout<< "Face " << faceI
<< " has no vertex labels."
<< endl;
nFailedChecks++;
continue; continue;
} }
@ -1196,7 +1208,7 @@ void dynamicTopoFvMesh::checkConnectivity(const label maxErrors) const
if (nCommon != 1) if (nCommon != 1)
{ {
Pout<< "Cells: " << nl Pout<< "Cells for face: " << faceI << "::" << curFace << nl
<< '\t' << owner_[faceI] << ":: " << ownCell << nl << '\t' << owner_[faceI] << ":: " << ownCell << nl
<< '\t' << neighbour_[faceI] << " :: " << neiCell << nl << '\t' << neighbour_[faceI] << " :: " << neiCell << nl
<< " share multiple faces. " << " share multiple faces. "
@ -1691,7 +1703,7 @@ void dynamicTopoFvMesh::checkConnectivity(const label maxErrors) const
} }
} }
if (!twoDMesh_) if (is3D())
{ {
Pout<< "Checking point-edge connectivity..."; Pout<< "Checking point-edge connectivity...";
@ -1818,8 +1830,8 @@ void dynamicTopoFvMesh::checkConnectivity(const label maxErrors) const
if if
( (
(cellToNode[cellI].size() != 6 && twoDMesh_) || (cellToNode[cellI].size() != 6 && is2D()) ||
(cellToNode[cellI].size() != 4 && !twoDMesh_) (cellToNode[cellI].size() != 4 && is3D())
) )
{ {
Pout<< nl << "Warning: Cell: " Pout<< nl << "Warning: Cell: "
@ -2183,9 +2195,7 @@ bool dynamicTopoFvMesh::checkCollapse
) const ) const
{ {
label faceIndex = -1; label faceIndex = -1;
scalar cQuality = 0; scalar cQuality = 0.0, oldVolume = 0.0, newVolume = 0.0;
scalar oldVolume = 0;
scalar newVolume = 0;
const cell& cellToCheck = cells_[cellIndex]; const cell& cellToCheck = cells_[cellIndex];
// Look for a face that doesn't contain 'pointIndex' // Look for a face that doesn't contain 'pointIndex'

View file

@ -47,6 +47,20 @@ namespace Foam
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Return if mesh is 2D
inline bool dynamicTopoFvMesh::is2D() const
{
return twoDMesh_;
}
// Return if mesh is 3D
inline bool dynamicTopoFvMesh::is3D() const
{
return !twoDMesh_;
}
// Return a non-const reference to the lengthScaleEstimator // Return a non-const reference to the lengthScaleEstimator
inline lengthScaleEstimator& inline lengthScaleEstimator&
dynamicTopoFvMesh::lengthEstimator() dynamicTopoFvMesh::lengthEstimator()
@ -121,7 +135,7 @@ inline label dynamicTopoFvMesh::getEdgeIndex
const edge& edgeToCheck const edge& edgeToCheck
) const ) const
{ {
if (twoDMesh_) if (is2D())
{ {
// No efficient search method. Loop through all edges. // No efficient search method. Loop through all edges.
forAll(edges_, edgeI) forAll(edges_, edgeI)
@ -206,7 +220,7 @@ inline bool dynamicTopoFvMesh::checkBisection
{ {
scalar length = 0.0, scale = 0.0; scalar length = 0.0, scale = 0.0;
if (twoDMesh_) if (is2D())
{ {
// If this entity was deleted, skip it. // If this entity was deleted, skip it.
if (faces_[index].empty()) if (faces_[index].empty())
@ -274,7 +288,7 @@ inline bool dynamicTopoFvMesh::checkCollapse
{ {
scalar length = 0.0, scale = 0.0; scalar length = 0.0, scale = 0.0;
if (twoDMesh_) if (is2D())
{ {
// If this entity was deleted, skip it. // If this entity was deleted, skip it.
if (faces_[index].empty()) if (faces_[index].empty())
@ -390,7 +404,7 @@ inline void dynamicTopoFvMesh::initStacks
tID = 0; tID = 0;
} }
if (twoDMesh_) if (is2D())
{ {
forAll(faces_, faceI) forAll(faces_, faceI)
{ {
@ -447,7 +461,7 @@ inline void dynamicTopoFvMesh::initStacks
stackElements[elemI] = stack(0)[elemI]; stackElements[elemI] = stack(0)[elemI];
} }
label elemType = twoDMesh_ ? 2 : 1; label elemType = is2D() ? 2 : 1;
writeVTK writeVTK
( (
@ -568,34 +582,20 @@ inline label dynamicTopoFvMesh::whichEdgePatch
} }
// Report topo-modification status // Report or alter topo-modification status
// - Enumerants are as follows inline label& dynamicTopoFvMesh::status(const label type)
// [0] nModifications
// [1] Total swaps (internal and surface)
// [2] Surface swaps
// [3] Total bisections
// [4] Total collapses
// [5] Surface bisections
// [6] Surface collapses
// [7] Slivers removed
inline label dynamicTopoFvMesh::status(const label type) const
{ {
if (type < statistics_.size()) if (type < 0 || type >= TOTAL_OP_TYPES)
{
return statistics_[type];
}
else
{ {
FatalErrorIn FatalErrorIn
( (
"inline label dynamicTopoFvMesh::status" "inline label& dynamicTopoFvMesh::status(const label type)"
"(const label type) const"
) )
<< " Unknown type: " << type << nl << " Unknown type: " << type << nl
<< abort(FatalError); << abort(FatalError);
} }
return 0; return statistics_[type];
} }
@ -618,6 +618,19 @@ inline void dynamicTopoFvMesh::setFlip(const label fIndex)
} }
// Check for processor priority
template <class BinaryOp>
inline bool dynamicTopoFvMesh::priority
(
const label first,
const BinaryOp& bop,
const label second
) const
{
return bop(procPriority_[first], procPriority_[second]);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam } // End namespace Foam

View file

@ -131,7 +131,7 @@ void dynamicTopoFvMesh::computeMapping
// for mild convexity. // for mild convexity.
const cell& cellToCheck = cells_[cIndex]; const cell& cellToCheck = cells_[cIndex];
if (twoDMesh_) if (is2D())
{ {
const labelList& parents = cellParents_[cIndex]; const labelList& parents = cellParents_[cIndex];
@ -603,6 +603,9 @@ void dynamicTopoFvMesh::threadedMapping
boundary[patchI].faceFaces(); boundary[patchI].faceFaces();
} }
// Force calculation of demand-driven data on subMeshes
initCoupledWeights();
// Execute threads in linear sequence // Execute threads in linear sequence
executeThreads(identity(nThreads), hdl, &computeMappingThread); executeThreads(identity(nThreads), hdl, &computeMappingThread);
} }
@ -668,7 +671,7 @@ void dynamicTopoFvMesh::setCellMapping
masterCells.append(mapCells[cellI]); masterCells.append(mapCells[cellI]);
} }
} }
else
if (cellParents_.found(mapCells[cellI])) if (cellParents_.found(mapCells[cellI]))
{ {
const labelList& nParents = cellParents_[mapCells[cellI]]; const labelList& nParents = cellParents_[mapCells[cellI]];

View file

@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "objectMap.H" #include "objectMap.H"
#include "objectRegistry.H"
#include "coupledInfo.H" #include "coupledInfo.H"
#include "dynamicTopoFvMesh.H" #include "dynamicTopoFvMesh.H"
@ -119,10 +118,6 @@ void dynamicTopoFvMesh::reOrderPoints
<< abort(FatalError); << abort(FatalError);
} }
// Wipe out pointsFromPoints, since this is not
// really used in this context for mapping.
pointsFromPoints_.clear();
// Renumber all maps. // Renumber all maps.
forAll(pointsFromPoints_, indexI) forAll(pointsFromPoints_, indexI)
{ {
@ -534,7 +529,7 @@ void dynamicTopoFvMesh::reOrderEdges
entityMutex_[2].unlock(); entityMutex_[2].unlock();
} }
if (!twoDMesh_) if (is3D())
{ {
// Invert edges to obtain pointEdges // Invert edges to obtain pointEdges
pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_); pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
@ -1725,7 +1720,30 @@ void dynamicTopoFvMesh::reOrderMesh
{ {
if (debug) if (debug)
{ {
Info<< nl // Buffered print-out of mesh-stats
if (Pstream::parRun())
{
if (Pstream::master())
{
for (label i = Pstream::nProcs() - 1; i >= 1; i--)
{
// Notify slave
meshOps::pWrite(i, i);
// Wait for response from slave
label j = 0;
meshOps::pRead(i, j);
}
}
else
{
// Wait for master
label i = 0;
meshOps::pRead(Pstream::masterNo(), i);
}
}
Pout<< nl
<< "=================" << nl << "=================" << nl
<< " Mesh reOrdering " << nl << " Mesh reOrdering " << nl
<< "=================" << nl << "=================" << nl
@ -1740,13 +1758,13 @@ void dynamicTopoFvMesh::reOrderMesh
if (debug > 1) if (debug > 1)
{ {
Info<< " Patch Starts [Face]: " << oldPatchStarts_ << nl Pout<< " Patch Starts [Face]: " << oldPatchStarts_ << nl
<< " Patch Sizes [Face]: " << oldPatchSizes_ << nl << " Patch Sizes [Face]: " << oldPatchSizes_ << nl
<< " Patch Starts [Edge]: " << oldEdgePatchStarts_ << nl << " Patch Starts [Edge]: " << oldEdgePatchStarts_ << nl
<< " Patch Sizes [Edge]: " << oldEdgePatchSizes_ << nl; << " Patch Sizes [Edge]: " << oldEdgePatchSizes_ << nl;
} }
Info<< "=================" << nl Pout<< "=================" << nl
<< " Mesh Info [n+1]:" << nl << " Mesh Info [n+1]:" << nl
<< " Points: " << nPoints_ << nl << " Points: " << nPoints_ << nl
<< " Edges: " << nEdges_ << nl << " Edges: " << nEdges_ << nl
@ -1758,13 +1776,22 @@ void dynamicTopoFvMesh::reOrderMesh
if (debug > 1) if (debug > 1)
{ {
Info<< " Patch Starts [Face]: " << patchStarts_ << nl Pout<< " Patch Starts [Face]: " << patchStarts_ << nl
<< " Patch Sizes: [Face]: " << patchSizes_ << nl << " Patch Sizes: [Face]: " << patchSizes_ << nl
<< " Patch Starts [Edge]: " << edgePatchStarts_ << nl << " Patch Starts [Edge]: " << edgePatchStarts_ << nl
<< " Patch Sizes: [Edge]: " << edgePatchSizes_ << nl; << " Patch Sizes: [Edge]: " << edgePatchSizes_ << nl;
} }
Info<< "=================" << endl; Pout<< "=================" << endl;
// Notify master
if (Pstream::parRun() && !Pstream::master())
{
meshOps::pWrite(Pstream::masterNo(), Pstream::myProcNo());
}
bool checkPoint = false;
reduce(checkPoint, andOp<bool>());
// Check connectivity structures for consistency // Check connectivity structures for consistency
// before entering the reOrdering phase. // before entering the reOrdering phase.

View file

@ -63,8 +63,6 @@ void eMesh::clearAddressing() const
edges_.clear(); edges_.clear();
deleteDemandDrivenData(pePtr_);
deleteDemandDrivenData(epPtr_);
deleteDemandDrivenData(fePtr_); deleteDemandDrivenData(fePtr_);
deleteDemandDrivenData(efPtr_); deleteDemandDrivenData(efPtr_);
} }
@ -147,8 +145,6 @@ eMesh::eMesh(const polyMesh& pMesh, const word& subDir)
), ),
*this *this
), ),
pePtr_(NULL),
epPtr_(NULL),
fePtr_(NULL), fePtr_(NULL),
efPtr_(NULL) efPtr_(NULL)
{ {
@ -391,28 +387,6 @@ void eMesh::setInstance(const fileName& inst)
} }
const labelListList& eMesh::pointEdges() const
{
if (!pePtr_)
{
calcPointEdges();
}
return *pePtr_;
}
const labelListList& eMesh::edgePoints() const
{
if (!epPtr_)
{
calcEdgePoints();
}
return *epPtr_;
}
const labelListList& eMesh::faceEdges() const const labelListList& eMesh::faceEdges() const
{ {
if (!fePtr_) if (!fePtr_)

View file

@ -35,7 +35,6 @@ SourceFiles
#ifndef eMesh_H #ifndef eMesh_H
#define eMesh_H #define eMesh_H
#include "objectRegistry.H"
#include "Time.H" #include "Time.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "edgeIOList.H" #include "edgeIOList.H"
@ -80,12 +79,6 @@ class eMesh
//- Mapping between ordered edge-data and regular edges. //- Mapping between ordered edge-data and regular edges.
labelList reverseEdgeMap_; labelList reverseEdgeMap_;
//- Point-edges
mutable labelListList* pePtr_;
//- Edge-points
mutable labelListList* epPtr_;
//- Face-edges //- Face-edges
mutable labelListIOList* fePtr_; mutable labelListIOList* fePtr_;
@ -111,12 +104,6 @@ class eMesh
//- Calculate ordered edges (and edgeFaces) //- Calculate ordered edges (and edgeFaces)
void calcOrderedEdgeList(); void calcOrderedEdgeList();
//- Calculate ordered edgePoints
void calcEdgePoints() const;
//- Calculate point-edges
void calcPointEdges() const;
//- Calculate face-edges //- Calculate face-edges
void calcFaceEdges() const; void calcFaceEdges() const;
@ -197,13 +184,6 @@ public:
// Demand-driven data // Demand-driven data
//- Return constant reference to the pointEdges list
const labelListList& pointEdges() const;
//- Return constant reference to the edgePoints list
//- Valid for 3D tetrahedral meshes only.
const labelListList& edgePoints() const;
//- Return constant reference to the faceEdges list //- Return constant reference to the faceEdges list
const labelListList& faceEdges() const; const labelListList& faceEdges() const;

View file

@ -145,201 +145,6 @@ void eMesh::calcOrderedEdgeList()
calcFaceEdges(); calcFaceEdges();
} }
void eMesh::calcEdgePoints() const
{
if (debug)
{
Info<< "void eMesh::calcEdgePoints() const : "
<< "Calculating ordered edgePoints" << endl;
}
if (epPtr_)
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgePoints already allocated."
<< abort(FatalError);
}
if (!efPtr_->size())
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgeFaces doesn't exist."
<< abort(FatalError);
}
// Size the list
epPtr_ = new labelListList(nEdges_);
labelListList& edgePoints = *epPtr_;
// EdgePoints are ordered such that points appear in
// counter-clockwise direction around point[0] of the edge.
// If a boundary point is present, that must begin the list.
// NOTE: Will work only on tetrahedral meshes!
bool found;
label faceIndex = -1, cellIndex = -1;
const labelList& owner = mesh_.faceOwner();
const labelList& neighbour = mesh_.faceNeighbour();
const cellList& cells = mesh_.cells();
const faceList& faces = mesh_.faces();
const labelListList& eFaces = this->edgeFaces();
for (label eIndex=0; eIndex < nEdges_; eIndex++)
{
const edge& e = edges_[eIndex];
const labelList& eFace = eFaces[eIndex];
// Size the list
edgePoints[eIndex].setSize(eFace.size(), -1);
if (eIndex < nInternalEdges_)
{
// Pick the first face and start with that
faceIndex = eFace[0];
}
else
{
// Need to find a properly oriented start-face
forAll(eFace, faceI)
{
if
(
(!mesh_.isInternalFace(eFace[faceI]))
&& (edgeDirection(faces[eFace[faceI]], e) == 1)
)
{
faceIndex = eFace[faceI];
break;
}
}
}
// Shuffle vertices to appear in CCW order
for (label j=0; j < eFace.size(); j++)
{
const face& f = faces[faceIndex];
// Add the isolated point
edgePoints[eIndex][j] = findIsolatedPoint(f, e);
// Figure out how this edge is oriented.
if (edgeDirection(f, e) == 1)
{
// Counter-clockwise. Pick the owner.
cellIndex = owner[faceIndex];
}
else if (mesh_.isInternalFace(faceIndex))
{
// Clockwise. Pick the neighbour.
cellIndex = neighbour[faceIndex];
}
else
{
// Looks like we've hit a boundary face. Break out.
break;
}
const cell& cellToCheck = cells[cellIndex];
found = false;
// Assuming tet-cells,
// Loop through edgeFaces and get the next face
forAll(eFace, faceI)
{
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[0]
)
{
faceIndex = cellToCheck[0];
found = true;
break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[1]
)
{
faceIndex = cellToCheck[1];
found = true;
break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[2]
)
{
faceIndex = cellToCheck[2];
found = true;
break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[3]
)
{
faceIndex = cellToCheck[3];
found = true;
break;
}
}
# ifdef FULLDEBUG
if (!found)
{
// Something's terribly wrong
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
)
<< " Failed to determine a vertex ring. " << nl
<< " edgeFaces connectivity is inconsistent. " << nl
<< "Edge: " << eIndex << ":: " << e << nl
<< "edgeFaces: " << eFace
<< abort(FatalError);
}
# else
// dummy statement to quech compiler warning
found = found;
# endif
}
}
}
void eMesh::calcPointEdges() const
{
if (debug)
{
Info<< "void eMesh::calcPointEdges() const : "
<< "Calculating PointEdges" << endl;
}
if (pePtr_)
{
FatalErrorIn
(
"void eMesh::calcPointEdges() const"
) << "pePtr_ already allocated"
<< abort(FatalError);
}
pePtr_ = new labelListList(mesh_.nPoints());
invertManyToMany(mesh_.nPoints(), edges(), *pePtr_);
}
void eMesh::calcFaceEdges() const void eMesh::calcFaceEdges() const
{ {

View file

@ -79,6 +79,8 @@ ePatch::ePatch
size_(readLabel(dict.lookup("size"))) size_(readLabel(dict.lookup("size")))
{} {}
//- Construct as copy, resetting the boundary mesh
ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm) ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm)
: :
patchIdentifier(p, p.index()), patchIdentifier(p, p.index()),
@ -88,6 +90,16 @@ ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm)
{} {}
//- Construct as copy
ePatch::ePatch(const ePatch& p)
:
patchIdentifier(p),
boundaryMesh_(p.boundaryMesh_),
start_(p.start_),
size_(p.size_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
ePatch::~ePatch() ePatch::~ePatch()

View file

@ -68,14 +68,6 @@ private:
//- Size of the patch //- Size of the patch
label size_; label size_;
// Demand-driven private data
// Private Member Functions
//- Disallow construct as copy
ePatch(const ePatch&);
protected: protected:
// The ePatch geometry initialisation is called by eBoundaryMesh // The ePatch geometry initialisation is called by eBoundaryMesh
@ -159,7 +151,6 @@ public:
const eBoundaryMesh& bm const eBoundaryMesh& bm
); );
//- Construct from dictionary //- Construct from dictionary
ePatch ePatch
( (
@ -172,6 +163,9 @@ public:
//- Construct as copy, resetting the boundary mesh //- Construct as copy, resetting the boundary mesh
ePatch(const ePatch&, const eBoundaryMesh&); ePatch(const ePatch&, const eBoundaryMesh&);
//- Construct as copy
ePatch(const ePatch&);
// Selectors // Selectors

View file

@ -65,7 +65,7 @@ autoPtr<ePatch> ePatch::New
) << "Unknown ePatch type " << patchType << " for patch " << name ) << "Unknown ePatch type " << patchType << " for patch " << name
<< endl << endl << endl << endl
<< "Valid ePatch types are :" << endl << "Valid ePatch types are :" << endl
<< wordConstructorTablePtr_->sortedToc() << wordConstructorTablePtr_->toc()
<< exit(FatalError); << exit(FatalError);
} }
@ -102,7 +102,7 @@ autoPtr<ePatch> ePatch::New
dict dict
) << "Unknown ePatch type " << patchType << endl << endl ) << "Unknown ePatch type " << patchType << endl << endl
<< "Valid ePatch types are :" << endl << "Valid ePatch types are :" << endl
<< dictionaryConstructorTablePtr_->sortedToc() << dictionaryConstructorTablePtr_->toc()
<< exit(FatalIOError); << exit(FatalIOError);
} }

File diff suppressed because it is too large Load diff

View file

@ -160,7 +160,7 @@ bool dynamicTopoFvMesh::testDelaunay
forAll(procIndices_, pI) forAll(procIndices_, pI)
{ {
// Fetch non-const reference to subMeshes // Fetch non-const reference to subMeshes
const coupledInfo& recvMesh = recvMeshes_[pI]; const coupledMesh& recvMesh = recvMeshes_[pI];
const coupleMap& cMap = recvMesh.map(); const coupleMap& cMap = recvMesh.map();
const dynamicTopoFvMesh& sMesh = recvMesh.subMesh(); const dynamicTopoFvMesh& sMesh = recvMesh.subMesh();
@ -1073,7 +1073,7 @@ const changeMap dynamicTopoFvMesh::swapQuadFace
topoChangeFlag_ = true; topoChangeFlag_ = true;
// Increment the counter // Increment the counter
statistics_[1]++; status(TOTAL_SWAPS)++;
// Return a successful operation. // Return a successful operation.
map.type() = 1; map.type() = 1;
@ -1853,7 +1853,7 @@ const changeMap dynamicTopoFvMesh::removeEdgeFlips
map.removeEdge(eIndex); map.removeEdge(eIndex);
// Increment the counter // Increment the counter
statistics_[1]++; status(TOTAL_SWAPS)++;
// Set the flag // Set the flag
topoChangeFlag_ = true; topoChangeFlag_ = true;
@ -2337,7 +2337,7 @@ const changeMap dynamicTopoFvMesh::swap23
if (debug > 2) if (debug > 2)
{ {
Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl; Pout<< " On SubMesh: " << isSubMesh_ << nl;
Pout<< " coupledModification: " << coupledModification_ << nl; Pout<< " coupledModification: " << coupledModification_ << nl;
label bPatch = whichEdgePatch(eIndex); label bPatch = whichEdgePatch(eIndex);
@ -2379,7 +2379,7 @@ const changeMap dynamicTopoFvMesh::swap23
+ "_beforeSwap_" + "_beforeSwap_"
+ Foam::name(numTriangulations) + "_" + Foam::name(numTriangulations) + "_"
+ Foam::name(triangulationIndex), + Foam::name(triangulationIndex),
cellsForRemoval, labelList(cellsForRemoval),
3, false, true 3, false, true
); );
} }
@ -2487,7 +2487,8 @@ const changeMap dynamicTopoFvMesh::swap23
-1, -1,
tmpTriFace, tmpTriFace,
newCellIndex[0], newCellIndex[0],
newCellIndex[1] newCellIndex[1],
labelList(0)
) )
); );
@ -2509,7 +2510,8 @@ const changeMap dynamicTopoFvMesh::swap23
-1, -1,
tmpTriFace, tmpTriFace,
newCellIndex[1], newCellIndex[1],
newCellIndex[2] newCellIndex[2],
labelList(0)
) )
); );
@ -2528,19 +2530,14 @@ const changeMap dynamicTopoFvMesh::swap23
-1, -1,
tmpTriFace, tmpTriFace,
newCellIndex[0], newCellIndex[0],
newCellIndex[2] newCellIndex[2],
labelList(0)
) )
); );
// Add this face to the map. // Add this face to the map.
map.addFace(newFaceIndex[2]); map.addFace(newFaceIndex[2]);
// Append three dummy faceEdges entries.
for (label i = 0; i < 3; i++)
{
faceEdges_.append(labelList(0));
}
// Add an entry to edgeFaces // Add an entry to edgeFaces
labelList newEdgeFaces(3); labelList newEdgeFaces(3);
newEdgeFaces[0] = newFaceIndex[0]; newEdgeFaces[0] = newFaceIndex[0];
@ -2901,7 +2898,7 @@ const changeMap dynamicTopoFvMesh::swap23
+ "_afterSwap_" + "_afterSwap_"
+ Foam::name(numTriangulations) + "_" + Foam::name(numTriangulations) + "_"
+ Foam::name(triangulationIndex), + Foam::name(triangulationIndex),
newCellIndex, labelList(newCellIndex),
3, false, true 3, false, true
); );
} }
@ -2996,7 +2993,7 @@ const changeMap dynamicTopoFvMesh::swap32
if (debug > 2) if (debug > 2)
{ {
Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl; Pout<< " On SubMesh: " << isSubMesh_ << nl;
Pout<< " coupledModification: " << coupledModification_ << nl; Pout<< " coupledModification: " << coupledModification_ << nl;
label bPatch = whichEdgePatch(eIndex); label bPatch = whichEdgePatch(eIndex);
@ -3082,7 +3079,8 @@ const changeMap dynamicTopoFvMesh::swap32
-1, -1,
newTriFace, newTriFace,
newCellIndex[0], newCellIndex[0],
newCellIndex[1] newCellIndex[1],
labelList(3, -1)
) )
); );
@ -3092,9 +3090,6 @@ const changeMap dynamicTopoFvMesh::swap32
// Note the triangulation face in index() // Note the triangulation face in index()
map.index() = newFaceIndex; map.index() = newFaceIndex;
// Add faceEdges for the new face as well.
faceEdges_.append(labelList(3));
// Define the three edges to check while building faceEdges: // Define the three edges to check while building faceEdges:
FixedList<edge,3> check; FixedList<edge,3> check;
@ -3208,7 +3203,8 @@ const changeMap dynamicTopoFvMesh::swap32
facePatch, facePatch,
newBdyTriFace[0], newBdyTriFace[0],
newCellIndex[1], newCellIndex[1],
-1 -1,
labelList(3, -1)
) )
); );
@ -3223,7 +3219,8 @@ const changeMap dynamicTopoFvMesh::swap32
facePatch, facePatch,
newBdyTriFace[1], newBdyTriFace[1],
newCellIndex[0], newCellIndex[0],
-1 -1,
labelList(3, -1)
) )
); );
@ -3286,11 +3283,11 @@ const changeMap dynamicTopoFvMesh::swap32
} }
// Add faceEdges for the two new boundary faces // Add faceEdges for the two new boundary faces
faceEdges_.append(bdyFaceEdges[0]); faceEdges_[newBdyFaceIndex[0]] = bdyFaceEdges[0];
faceEdges_.append(bdyFaceEdges[1]); faceEdges_[newBdyFaceIndex[1]] = bdyFaceEdges[1];
// Update the number of surface swaps. // Update the number of surface swaps.
statistics_[2]++; status(SURFACE_SWAPS)++;
} }
newTetCell[0][nF0++] = newFaceIndex; newTetCell[0][nF0++] = newFaceIndex;
@ -3520,7 +3517,7 @@ const changeMap dynamicTopoFvMesh::swap32
+ "_afterSwap_" + "_afterSwap_"
+ Foam::name(numTriangulations) + "_" + Foam::name(numTriangulations) + "_"
+ Foam::name(triangulationIndex), + Foam::name(triangulationIndex),
newCellIndex, labelList(newCellIndex),
3, false, true 3, false, true
); );
} }

View file

@ -76,7 +76,7 @@ void conservativeMapVolFields(const topoMapper& mapper)
// Map patch fields // Map patch fields
forAll(bMap, patchI) forAll(bMap, patchI)
{ {
bMap[patchI].mapPatchField bMap[patchI].mapFvPatchField
( (
field.name(), field.name(),
field.boundaryField()[patchI] field.boundaryField()[patchI]
@ -131,7 +131,7 @@ void conservativeMapSurfaceFields(const topoMapper& mapper)
// Map patch fields // Map patch fields
forAll(bMap, patchI) forAll(bMap, patchI)
{ {
bMap[patchI].mapPatchField bMap[patchI].mapFvsPatchField
( (
field.name(), field.name(),
field.boundaryField()[patchI] field.boundaryField()[patchI]

View file

@ -90,7 +90,7 @@ autoPtr<fluxCorrector> fluxCorrector::New
) << "Unknown fluxCorrector type " << correctorTypeName ) << "Unknown fluxCorrector type " << correctorTypeName
<< endl << endl << endl << endl
<< "Valid fluxCorrector types are: " << endl << "Valid fluxCorrector types are: " << endl
<< meshConstructorTablePtr_->sortedToc() << meshConstructorTablePtr_->toc()
<< exit(FatalError); << exit(FatalError);
} }

View file

@ -39,7 +39,6 @@ Author
#include "topoCellMapper.H" #include "topoCellMapper.H"
#include "topoSurfaceMapper.H" #include "topoSurfaceMapper.H"
#include "topoBoundaryMeshMapper.H" #include "topoBoundaryMeshMapper.H"
#include "fixedValueFvPatchFields.H"
namespace Foam namespace Foam
{ {
@ -49,13 +48,14 @@ namespace Foam
//- Store gradients prior to mesh reset //- Store gradients prior to mesh reset
void topoMapper::storeGradients() const void topoMapper::storeGradients() const
{ {
storeGradients<scalar>(sGrads_); // Go in order of highest to lowest rank,
storeGradients<vector>(vGrads_); // to avoid double storage
storeGradients<vector>(gradTable_, vGradPtrs_);
storeGradients<scalar>(gradTable_, sGradPtrs_);
if (fvMesh::debug) if (fvMesh::debug)
{ {
Info<< "Registered volScalarFields: " << scalarGrads() << endl; Info<< "Registered gradients: " << gradientTable() << endl;
Info<< "Registered volVectorFields: " << vectorGrads() << endl;
} }
} }
@ -63,6 +63,10 @@ void topoMapper::storeGradients() const
//- Store geometric information //- Store geometric information
void topoMapper::storeGeometry() const void topoMapper::storeGeometry() const
{ {
typedef volVectorField::PatchFieldType PatchFieldType;
typedef volVectorField::GeometricBoundaryField GeomBdyFieldType;
typedef volVectorField::DimensionedInternalField DimInternalField;
// Wipe out existing information // Wipe out existing information
deleteDemandDrivenData(cellCentresPtr_); deleteDemandDrivenData(cellCentresPtr_);
@ -73,27 +77,42 @@ void topoMapper::storeGeometry() const
label nPatches = mesh_.boundary().size(); label nPatches = mesh_.boundary().size();
// Create field parts // Create field parts
PtrList<fvPatchField<vector> > volCentrePatches(nPatches); PtrList<PatchFieldType> volCentrePatches(nPatches);
// Over-ride and set all patches to fixedValue // Define patch type names
for (label patchI = 0; patchI < nPatches; patchI++) word emptyType("empty");
word fixedValueType("fixedValue");
// Create dummy types for initial field creation
forAll(volCentrePatches, patchI)
{
if (mesh_.boundary()[patchI].type() == emptyType)
{ {
volCentrePatches.set volCentrePatches.set
( (
patchI, patchI,
new fixedValueFvPatchField<vector> PatchFieldType::New
( (
emptyType,
mesh_.boundary()[patchI], mesh_.boundary()[patchI],
DimensionedField<vector, volMesh>::null() DimInternalField::null()
) )
); );
}
// Slice field to patch (forced assignment) else
volCentrePatches[patchI] == {
volCentrePatches.set
( (
mesh_.boundaryMesh()[patchI].patchSlice(Cf) patchI,
PatchFieldType::New
(
fixedValueType,
mesh_.boundary()[patchI],
DimInternalField::null()
)
); );
} }
}
// Set the cell-centres pointer. // Set the cell-centres pointer.
cellCentresPtr_ = cellCentresPtr_ =
@ -107,7 +126,7 @@ void topoMapper::storeGeometry() const
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE, IOobject::NO_WRITE,
false true
), ),
mesh_, mesh_,
dimLength, dimLength,
@ -115,6 +134,48 @@ void topoMapper::storeGeometry() const
volCentrePatches volCentrePatches
) )
); );
// Alias for convenience
volVectorField& centres = *cellCentresPtr_;
// Set correct references for patch internal fields
GeomBdyFieldType& bf = centres.boundaryField();
forAll(bf, patchI)
{
if (mesh_.boundary()[patchI].type() == emptyType)
{
bf.set
(
patchI,
PatchFieldType::New
(
emptyType,
mesh_.boundary()[patchI],
centres.dimensionedInternalField()
)
);
}
else
{
bf.set
(
patchI,
PatchFieldType::New
(
fixedValueType,
mesh_.boundary()[patchI],
centres.dimensionedInternalField()
)
);
// Slice field to patch (forced assignment)
bf[patchI] == mesh_.boundaryMesh()[patchI].patchSlice(Cf);
}
}
// Set the cell-volumes pointer
cellVolumesPtr_ = new scalarField(mesh_.cellVolumes());
} }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -132,6 +193,7 @@ topoMapper::topoMapper
surfaceMap_(NULL), surfaceMap_(NULL),
boundaryMap_(NULL), boundaryMap_(NULL),
fluxCorrector_(fluxCorrector::New(mesh, dict)), fluxCorrector_(fluxCorrector::New(mesh, dict)),
cellVolumesPtr_(NULL),
cellCentresPtr_(NULL) cellCentresPtr_(NULL)
{} {}
@ -339,6 +401,22 @@ const vectorField& topoMapper::internalCentres() const
} }
//- Return non-const access to cell volumes
scalarField& topoMapper::internalVolumes() const
{
if (!cellVolumesPtr_)
{
FatalErrorIn
(
"scalarField& topoMapper::internalVolumes() const"
) << nl << " Pointer has not been set. "
<< abort(FatalError);
}
return *cellVolumesPtr_;
}
//- Return stored patch centre information //- Return stored patch centre information
const vectorField& topoMapper::patchCentres(const label i) const const vectorField& topoMapper::patchCentres(const label i) const
{ {
@ -356,17 +434,10 @@ const vectorField& topoMapper::patchCentres(const label i) const
} }
//- Return names of stored scalar gradients //- Return names of stored gradients
const wordList topoMapper::scalarGrads() const const wordList topoMapper::gradientTable() const
{ {
return sGrads_.toc(); return gradTable_.toc();
}
//- Return names of stored vector gradients
const wordList topoMapper::vectorGrads() const
{
return vGrads_.toc();
} }
@ -374,7 +445,7 @@ const wordList topoMapper::vectorGrads() const
template <> template <>
volVectorField& topoMapper::gradient(const word& name) const volVectorField& topoMapper::gradient(const word& name) const
{ {
if (!sGrads_.found(name)) if (!gradTable_.found(name))
{ {
FatalErrorIn FatalErrorIn
( (
@ -384,7 +455,7 @@ volVectorField& topoMapper::gradient(const word& name) const
<< abort(FatalError); << abort(FatalError);
} }
return sGrads_[name](); return sGradPtrs_[gradTable_[name].second()];
} }
@ -392,7 +463,7 @@ volVectorField& topoMapper::gradient(const word& name) const
template <> template <>
volTensorField& topoMapper::gradient(const word& name) const volTensorField& topoMapper::gradient(const word& name) const
{ {
if (!vGrads_.found(name)) if (!gradTable_.found(name))
{ {
FatalErrorIn FatalErrorIn
( (
@ -402,7 +473,28 @@ volTensorField& topoMapper::gradient(const word& name) const
<< abort(FatalError); << abort(FatalError);
} }
return vGrads_[name](); return vGradPtrs_[gradTable_[name].second()];
}
//- Deregister gradient fields and centres,
// but retain for mapping
void topoMapper::deregisterMeshInformation() const
{
// Check out scalar gradients
forAll(sGradPtrs_, fieldI)
{
mesh_.objectRegistry::checkOut(sGradPtrs_[fieldI]);
}
// Check out vector gradients
forAll(vGradPtrs_, fieldI)
{
mesh_.objectRegistry::checkOut(vGradPtrs_[fieldI]);
}
// Check out cell centres
mesh_.objectRegistry::checkOut(*cellCentresPtr_);
} }
@ -495,11 +587,15 @@ void topoMapper::clear() const
surfaceMap_.clear(); surfaceMap_.clear();
boundaryMap_.clear(); boundaryMap_.clear();
// Clear index maps
gradTable_.clear();
// Clear stored gradients // Clear stored gradients
sGrads_.clear(); sGradPtrs_.clear();
vGrads_.clear(); vGradPtrs_.clear();
// Wipe out geomtry information // Wipe out geomtry information
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_); deleteDemandDrivenData(cellCentresPtr_);
// Clear maps // Clear maps

View file

@ -41,7 +41,9 @@ SourceFiles
#ifndef topoMapper_H #ifndef topoMapper_H
#define topoMapper_H #define topoMapper_H
#include "Tuple2.H"
#include "autoPtr.H" #include "autoPtr.H"
#include "PtrList.H"
#include "IOmanip.H" #include "IOmanip.H"
#include "volFields.H" #include "volFields.H"
@ -62,6 +64,11 @@ class topoBoundaryMeshMapper;
class topoMapper class topoMapper
{ {
// Private typedefs
typedef Tuple2<word, label> GradientMap;
typedef HashTable<GradientMap> GradientTable;
// Private data // Private data
//- Reference to fvMesh //- Reference to fvMesh
@ -84,11 +91,15 @@ class topoMapper
//- Flux corrector //- Flux corrector
mutable autoPtr<fluxCorrector> fluxCorrector_; mutable autoPtr<fluxCorrector> fluxCorrector_;
// Index map for gradients
mutable GradientTable gradTable_;
// Stored gradients for mapping // Stored gradients for mapping
mutable HashTable<autoPtr<volVectorField> > sGrads_; mutable PtrList<volVectorField> sGradPtrs_;
mutable HashTable<autoPtr<volTensorField> > vGrads_; mutable PtrList<volTensorField> vGradPtrs_;
//- Geometric information on the old mesh //- Geometric information on the old mesh
mutable scalarField* cellVolumesPtr_;
mutable volVectorField* cellCentresPtr_; mutable volVectorField* cellCentresPtr_;
//- Intersection weights //- Intersection weights
@ -120,7 +131,8 @@ class topoMapper
template <class Type, class gradType> template <class Type, class gradType>
void storeGradients void storeGradients
( (
HashTable<autoPtr<gradType> >& gradTable GradientTable& gradTable,
PtrList<gradType>& gradList
) const; ) const;
//- Store gradients prior to mesh reset //- Store gradients prior to mesh reset
@ -209,20 +221,24 @@ public:
//- Store mesh information for the mapping stage //- Store mesh information for the mapping stage
void storeMeshInformation() const; void storeMeshInformation() const;
//- Deregister gradient fields and centres,
// but retain for mapping
void deregisterMeshInformation() const;
//- Return non-const access to cell centres //- Return non-const access to cell centres
volVectorField& volCentres() const; volVectorField& volCentres() const;
//- Return non-const access to cell volumes
scalarField& internalVolumes() const;
//- Return stored cell centre information //- Return stored cell centre information
const vectorField& internalCentres() const; const vectorField& internalCentres() const;
//- Return stored patch centre information //- Return stored patch centre information
const vectorField& patchCentres(const label i) const; const vectorField& patchCentres(const label i) const;
//- Return names of stored scalar gradients //- Return names of stored gradients
const wordList scalarGrads() const; const wordList gradientTable() const;
//- Return names of stored vector gradients
const wordList vectorGrads() const;
//- Fetch the gradient field //- Fetch the gradient field
template <class Type> template <class Type>

View file

@ -35,68 +35,106 @@ namespace Foam
template <class Type, class gradType> template <class Type, class gradType>
void topoMapper::storeGradients void topoMapper::storeGradients
( (
HashTable<autoPtr<gradType> >& gradTable GradientTable& gradTable,
PtrList<gradType>& gradList
) const ) const
{ {
// Define a few typedefs for convenience // Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType; typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gVolType; typedef const GeometricField<Type, fvPatchField, volMesh> constVolType;
typedef HashTable<constVolType*> volTypeTable;
// Fetch all fields from registry // Fetch all fields from registry
HashTable<const volType*> fields volTypeTable fields(mesh_.objectRegistry::lookupClass<volType>());
(
mesh_.objectRegistry::lookupClass<volType>() // Track field count
); label nFields = 0;
// Store old-times before gradient computation // Store old-times before gradient computation
forAllIter(typename HashTable<const volType*>, fields, fIter) for
(
typename volTypeTable::iterator fIter = fields.begin();
fIter != fields.end();
++fIter
)
{ {
fIter()->storeOldTimes(); fIter()->storeOldTimes();
nFields++;
} }
forAllConstIter(typename HashTable<const volType*>, fields, fIter) // Size up the list
gradList.setSize(nFields);
label fieldIndex = 0;
for
(
typename volTypeTable::const_iterator fIter = fields.begin();
fIter != fields.end();
++fIter
)
{ {
const volType& field = *fIter(); const volType& field = *fIter();
// Compute the gradient. // Compute the gradient.
tmp<gVolType> tGrad;
// If the fvSolution dictionary contains an entry, // If the fvSolution dictionary contains an entry,
// use that, otherwise, default to leastSquares // use that, otherwise, default to leastSquares
word gradName("grad(" + field.name() + ')'); word gradName("grad(" + field.name() + ')');
// Register field under a name that's unique
word registerName("remapGradient(" + field.name() + ')');
// Make a new entry
if (mesh_.schemesDict().subDict("gradSchemes").found(gradName)) if (mesh_.schemesDict().subDict("gradSchemes").found(gradName))
{ {
tGrad = fvc::grad(field); gradList.set
}
else
{
tGrad = fv::leastSquaresGrad<Type>(mesh_).grad(field);
}
// Make a new entry, but don't register the field.
gradTable.insert
(
field.name(),
autoPtr<gradType>
( (
fieldIndex,
new gradType new gradType
( (
IOobject IOobject
( (
tGrad().name(), registerName,
mesh_.time().timeName(), mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE, IOobject::NO_WRITE,
false true
), ),
tGrad() fvc::grad(field, gradName)()
) )
);
}
else
{
gradList.set
(
fieldIndex,
new gradType
(
IOobject
(
registerName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
true
),
fv::leastSquaresGrad<Type>(mesh_).grad(field)()
) )
); );
} }
// Add a map entry
gradTable.insert
(
field.name(),
GradientMap(registerName, fieldIndex++)
);
}
} }

View file

@ -313,7 +313,7 @@ void topoPatchMapper::calcAddressing() const
// Do mapped faces. Note that this can already be set by insertedFaces // Do mapped faces. Note that this can already be set by insertedFaces
// so check if addressing size still zero. // so check if addressing size still zero.
const labelList& fm = patch_.patch().patchSlice(mpm_.faceMap()); const labelList::subList fm = patch_.patch().patchSlice(mpm_.faceMap());
forAll(fm, faceI) forAll(fm, faceI)
{ {
@ -367,6 +367,8 @@ void topoPatchMapper::calcAddressing() const
continue; continue;
} }
label oldFace = (faceI < fm.size() ? fm[faceI] : -1);
FatalErrorIn FatalErrorIn
( (
"void topoPatchMapper::calcAddressing() const" "void topoPatchMapper::calcAddressing() const"
@ -374,8 +376,11 @@ void topoPatchMapper::calcAddressing() const
<< "Addressing is missing." << nl << "Addressing is missing." << nl
<< " Patch face index: " << faceI << nl << " Patch face index: " << faceI << nl
<< " nInsertedFaces: " << insertedFaces.size() << nl << " nInsertedFaces: " << insertedFaces.size() << nl
<< " faceMap: " << fm[faceI] << nl << " faceMap[faceI]: " << oldFace << nl
<< " patch faceMap size: " << fm.size() << nl
<< " Patch: " << patch_.name() << nl << " Patch: " << patch_.name() << nl
<< " polyPatch: " << nl << patch_.patch() << nl
<< " faceMap size: " << mpm_.faceMap().size() << nl
<< abort(FatalError); << abort(FatalError);
} }
} }

View file

@ -187,10 +187,18 @@ public:
//- Map the patch field //- Map the patch field
template <class Type> template <class Type>
void mapPatchField void mapFvPatchField
( (
const word& fieldName, const word& fieldName,
Field<Type>& pF fvPatchField<Type>& pF
) const;
//- Map the patch field
template <class Type>
void mapFvsPatchField
(
const word& fieldName,
fvsPatchField<Type>& pF
) const; ) const;
}; };

View file

@ -32,15 +32,49 @@ namespace Foam
//- Map the patch field //- Map the patch field
template <class Type> template <class Type>
void topoPatchMapper::mapPatchField void topoPatchMapper::mapFvPatchField
( (
const word& fieldName, const word& fieldName,
Field<Type>& pF fvPatchField<Type>& pF
) const ) const
{ {
// Specify that mapping is conservative // Specify that mapping is conservative
conservative_ = true; conservative_ = true;
if (fvMesh::debug)
{
Pout<< " Field: " << fieldName
<< " Mapping patch: " << pF.patch().name()
<< " size: " << this->size()
<< " sizeBeforeMapping: " << this->sizeBeforeMapping()
<< endl;
}
// Map patchField onto itself
pF.autoMap(*this);
}
//- Map the patch field
template <class Type>
void topoPatchMapper::mapFvsPatchField
(
const word& fieldName,
fvsPatchField<Type>& pF
) const
{
// Specify that mapping is conservative
conservative_ = true;
if (fvMesh::debug)
{
Pout<< " Field: " << fieldName
<< " Mapping patch: " << pF.patch().name()
<< " size: " << this->size()
<< " sizeBeforeMapping: " << this->sizeBeforeMapping()
<< endl;
}
// Map patchField onto itself // Map patchField onto itself
pF.autoMap(*this); pF.autoMap(*this);
} }

View file

@ -34,6 +34,7 @@ Author
\*----------------------------------------------------------------------------*/ \*----------------------------------------------------------------------------*/
#include "fvc.H"
#include "volFields.H" #include "volFields.H"
#include "lengthScaleEstimator.H" #include "lengthScaleEstimator.H"
#include "processorPolyPatch.H" #include "processorPolyPatch.H"
@ -43,7 +44,14 @@ namespace Foam
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(lengthScaleEstimator,0); defineTypeNameAndDebug(lengthScaleEstimator, 0);
lengthScaleEstimator::ScaleFnPair lengthScaleEstimator::methods_[] =
{
ScaleFnPair("constant", &lengthScaleEstimator::constantScale),
ScaleFnPair("direct", &lengthScaleEstimator::directScale),
ScaleFnPair("inverse", &lengthScaleEstimator::inverseScale)
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -66,6 +74,8 @@ lengthScaleEstimator::lengthScaleEstimator
sliceHoldOff_(0), sliceHoldOff_(0),
sliceBoxes_(0), sliceBoxes_(0),
field_("none"), field_("none"),
scale_(NULL),
gradient_(false),
fieldLength_(0.0), fieldLength_(0.0),
lowerRefineLevel_(0.001), lowerRefineLevel_(0.001),
upperRefineLevel_(0.999), upperRefineLevel_(0.999),
@ -303,7 +313,7 @@ void lengthScaleEstimator::writeLengthScaleInfo
// First perform a blocking send/receive of the number of faces. // First perform a blocking send/receive of the number of faces.
OPstream::write OPstream::write
( (
Pstream::blocking, Pstream::nonBlocking,
neiProcNo, neiProcNo,
reinterpret_cast<const char*>(&(nSendFaces[pI])), reinterpret_cast<const char*>(&(nSendFaces[pI])),
sizeof(label) sizeof(label)
@ -311,7 +321,7 @@ void lengthScaleEstimator::writeLengthScaleInfo
IPstream::read IPstream::read
( (
Pstream::blocking, Pstream::nonBlocking,
neiProcNo, neiProcNo,
reinterpret_cast<char*>(&(nRecvFaces[pI])), reinterpret_cast<char*>(&(nRecvFaces[pI])),
sizeof(label) sizeof(label)
@ -319,6 +329,10 @@ void lengthScaleEstimator::writeLengthScaleInfo
} }
} }
// Wait for all transfers to complete.
OPstream::waitRequests();
IPstream::waitRequests();
// Send info to neighbouring processors. // Send info to neighbouring processors.
forAll(boundary, patchI) forAll(boundary, patchI)
{ {
@ -548,6 +562,69 @@ void lengthScaleEstimator::readLengthScaleInfo
} }
// Use a constant length scale for field-based refinement
scalar lengthScaleEstimator::constantScale(const scalar fieldValue) const
{
return fieldLength_;
}
// Use a direct-proportion length scale for field-based refinement
scalar lengthScaleEstimator::directScale(const scalar fieldValue) const
{
if (fieldValue > upperRefineLevel_)
{
return meanScale_;
}
scalar value = fieldValue;
// Bracket the value within limits
value = min(value, upperRefineLevel_);
value = max(value, lowerRefineLevel_);
scalar ratio
(
(value - lowerRefineLevel_) / (upperRefineLevel_ - lowerRefineLevel_)
);
scalar scale
(
minLengthScale_ + ((maxLengthScale_ - minLengthScale_) * ratio)
);
return scale;
}
// Use an inverse-proportion length scale for field-based refinement
scalar lengthScaleEstimator::inverseScale(const scalar fieldValue) const
{
if (fieldValue < lowerRefineLevel_)
{
return meanScale_;
}
scalar value = fieldValue;
// Bracket the value within limits
value = min(value, upperRefineLevel_);
value = max(value, lowerRefineLevel_);
scalar ratio
(
(value - lowerRefineLevel_) / (upperRefineLevel_ - lowerRefineLevel_)
);
scalar scale
(
minLengthScale_ + ((maxLengthScale_ - minLengthScale_) * (1.0 - ratio))
);
return scale;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Read edge refinement options from the dictionary // Read edge refinement options from the dictionary
@ -764,6 +841,43 @@ void lengthScaleEstimator::readRefinementOptions
{ {
field_ = word(refineDict.lookup("fieldRefinement")); field_ = word(refineDict.lookup("fieldRefinement"));
word scaleMethod(refineDict.lookup("fieldScaleMethod"));
const label nMethods(sizeof(methods_) / sizeof(methods_[0]));
bool invalidMethod = true;
for (label i = 0; i < nMethods; i++)
{
if (scaleMethod == methods_[i].first())
{
scale_ = methods_[i].second();
invalidMethod = false;
break;
}
}
if (invalidMethod)
{
Info << nl << " Available methods: " << endl;
for (label i = 0; i < nMethods; i++)
{
Info << " " << methods_[i].first() << nl;
}
FatalErrorIn
(
"void lengthScaleEstimator::readRefinementOptions"
"(const dictionary&, bool, bool)"
)
<< " Invalid method: " << scaleMethod << " specified." << nl
<< abort(FatalError);
}
// Define whether the gradient must be evaluated
gradient_ = readBool(refineDict.lookup("fieldGradient"));
// Lookup a specified length-scale // Lookup a specified length-scale
fieldLength_ = readScalar(refineDict.lookup("fieldLengthScale")); fieldLength_ = readScalar(refineDict.lookup("fieldLengthScale"));
@ -799,7 +913,7 @@ void lengthScaleEstimator::readRefinementOptions
} }
//- Set explicitly coupled patch information // Set explicitly coupled patch information
void lengthScaleEstimator::setCoupledPatches void lengthScaleEstimator::setCoupledPatches
( (
const dictionary& coupledPatches const dictionary& coupledPatches
@ -844,7 +958,7 @@ void lengthScaleEstimator::setCoupledPatches
} }
//- Calculate the length scale field // Calculate the length scale field
void lengthScaleEstimator::calculateLengthScale void lengthScaleEstimator::calculateLengthScale
( (
UList<scalar>& lengthScale UList<scalar>& lengthScale
@ -915,11 +1029,103 @@ void lengthScaleEstimator::calculateLengthScale
// If a field has been specified, use that. // If a field has been specified, use that.
if (field_ != "none") if (field_ != "none")
{ {
const volScalarField& vFld = const volScalarField* vFldPtr = NULL;
// Temporary pointer for gradient
volScalarField* magGradPtr = NULL;
if (gradient_)
{
// Lookup various field types, and evaluate the gradient
bool invalidObject = true;
// Evaluate using gradient scheme
word gradName("grad(" + field_ + ')');
// Register field under a name that's unique
word registerName("lengthScaleGradient(" + field_ + ')');
if (mesh_.objectRegistry::foundObject<volScalarField>(field_))
{
const volScalarField& field =
( (
mesh_.objectRegistry::lookupObject<volScalarField>(field_) mesh_.objectRegistry::lookupObject<volScalarField>(field_)
); );
magGradPtr =
(
new volScalarField
(
IOobject
(
registerName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mag(fvc::grad(field, gradName))
)
);
invalidObject = false;
}
if (mesh_.objectRegistry::foundObject<volVectorField>(field_))
{
const volVectorField& field =
(
mesh_.objectRegistry::lookupObject<volVectorField>(field_)
);
magGradPtr =
(
new volScalarField
(
IOobject
(
registerName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mag(fvc::grad(field, gradName))
)
);
invalidObject = false;
}
if (invalidObject)
{
FatalErrorIn
(
"void lengthScaleEstimator::calculateLengthScale"
"(UList<scalar>& lengthScale)"
)
<< " Invalid field: " << field_
<< " specified for gradient." << nl
<< abort(FatalError);
}
vFldPtr = magGradPtr;
}
else
{
// Lookup a scalar field by default
const volScalarField& field =
(
mesh_.objectRegistry::lookupObject<volScalarField>(field_)
);
vFldPtr = &field;
}
const volScalarField& vFld = *vFldPtr;
const labelList& own = mesh_.faceOwner(); const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour(); const labelList& nei = mesh_.faceNeighbour();
@ -935,7 +1141,7 @@ void lengthScaleEstimator::calculateLengthScale
if (!cellLevels[ownCell]) if (!cellLevels[ownCell])
{ {
cellLevels[ownCell] = level; cellLevels[ownCell] = level;
lengthScale[ownCell] = fieldLength_; lengthScale[ownCell] = (this->*scale_)(fAvg);
levelCells.insert(ownCell); levelCells.insert(ownCell);
@ -945,7 +1151,7 @@ void lengthScaleEstimator::calculateLengthScale
if (!cellLevels[neiCell]) if (!cellLevels[neiCell])
{ {
cellLevels[neiCell] = level; cellLevels[neiCell] = level;
lengthScale[neiCell] = fieldLength_; lengthScale[neiCell] = (this->*scale_)(fAvg);
levelCells.insert(neiCell); levelCells.insert(neiCell);
@ -953,6 +1159,12 @@ void lengthScaleEstimator::calculateLengthScale
} }
} }
} }
// Clean up gradient, if necessary
if (gradient_)
{
deleteDemandDrivenData(magGradPtr);
}
} }
bool doneWithSweeps = false; bool doneWithSweeps = false;

View file

@ -43,6 +43,7 @@ SourceFiles
#ifndef lengthScaleEstimator_H #ifndef lengthScaleEstimator_H
#define lengthScaleEstimator_H #define lengthScaleEstimator_H
#include "Tuple2.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "dictionary.H" #include "dictionary.H"
@ -108,8 +109,14 @@ class lengthScaleEstimator
// are to be avoided // are to be avoided
labelList noModPatchIDs_; labelList noModPatchIDs_;
//- Typedefs for field refinement function pointers
typedef scalar (lengthScaleEstimator::*ScaleFn)(const scalar) const;
typedef Tuple2<const char*, ScaleFn> ScaleFnPair;
//- Field-based refinement //- Field-based refinement
word field_; word field_;
ScaleFn scale_;
bool gradient_;
scalar fieldLength_; scalar fieldLength_;
scalar lowerRefineLevel_; scalar lowerRefineLevel_;
scalar upperRefineLevel_; scalar upperRefineLevel_;
@ -118,6 +125,9 @@ class lengthScaleEstimator
scalar meanScale_; scalar meanScale_;
label maxRefineLevel_; label maxRefineLevel_;
//- Function pointers for scale methods
static ScaleFnPair methods_[];
// Private Member Functions // Private Member Functions
// Check for legitimacy of patches // Check for legitimacy of patches
@ -153,9 +163,19 @@ class lengthScaleEstimator
labelHashSet& levelCells labelHashSet& levelCells
) const; ) const;
// Use a constant length scale for field-based refinement
scalar constantScale(const scalar fieldValue) const;
// Use a direct-proportion length scale for field-based refinement
scalar directScale(const scalar fieldValue) const;
// Use an inverse-proportion length scale for field-based refinement
scalar inverseScale(const scalar fieldValue) const;
public: public:
// Declare the name of the class and its debug switch // Declare the name of the class and its debug switch
TypeName("lengthScaleEstimator"); TypeName("lengthScaleEstimator");
// Constructors // Constructors

View file

@ -34,17 +34,15 @@ Author
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "objectRegistry.H"
#include "Time.H" #include "Time.H"
#include "meshOps.H" #include "meshOps.H"
#include "ListOps.H" #include "ListOps.H"
#include "Pstream.H" #include "Pstream.H"
#include "triFace.H" #include "triFace.H"
#include "IOmanip.H" #include "IOmanip.H"
#include "HashSet.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "OFstream.H"
#include "triPointRef.H" #include "triPointRef.H"
#include "tetPointRef.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -56,7 +54,7 @@ namespace meshOps
// Utility method to build a hull of cells // Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes) // connected to the edge (for 2D simplical meshes)
void constructPrismHull inline void constructPrismHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -126,7 +124,7 @@ void constructPrismHull
// Utility method to build a hull of cells (and faces) // Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes) // around an edge (for 3D simplical meshes)
void constructHull inline void constructHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -333,7 +331,7 @@ void constructHull
// Renaud Waldura // Renaud Waldura
// Dijkstra's Shortest Path Algorithm in Java // Dijkstra's Shortest Path Algorithm in Java
// http://renaud.waldura.com/ // http://renaud.waldura.com/
bool Dijkstra inline bool Dijkstra
( (
const Map<point>& points, const Map<point>& points,
const Map<edge>& edges, const Map<edge>& edges,
@ -463,13 +461,47 @@ bool Dijkstra
} }
} }
/*
// Write out the path
if (debug > 3)
{
if (foundEndPoint)
{
DynamicList<label> pathNodes(50);
label currentPoint = endPoint;
while (currentPoint != startPoint)
{
pathNodes.append(currentPoint);
currentPoint = pi[currentPoint];
}
pathNodes.append(startPoint);
pathNodes.shrink();
writeVTK
(
"DijkstraPath_"
+ Foam::name(startPoint)
+ '_'
+ Foam::name(endPoint),
pathNodes,
0
);
}
}
*/
return foundEndPoint; return foundEndPoint;
} }
// Select a list of elements from connectivity, // Select a list of elements from connectivity,
// and output to a VTK format // and output to a VTK format
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,
@ -847,7 +879,7 @@ void writeVTK
// Actual routine to write out the VTK file // Actual routine to write out the VTK file
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,

View file

@ -33,8 +33,8 @@ Author
All rights reserved All rights reserved
SourceFiles SourceFiles
meshOpsI.H
meshOps.C meshOps.C
meshOpsTemplates.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -93,7 +93,7 @@ namespace meshOps
// Utility method to build a hull of cells // Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes) // connected to the edge (for 2D simplical meshes)
void constructPrismHull inline void constructPrismHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -107,7 +107,7 @@ namespace meshOps
// Utility method to build a hull of cells (and faces) // Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes) // around an edge (for 3D simplical meshes)
void constructHull inline void constructHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -181,7 +181,7 @@ namespace meshOps
); );
// Dijkstra's algorithm for the shortest path problem // Dijkstra's algorithm for the shortest path problem
bool Dijkstra inline bool Dijkstra
( (
const Map<point>& points, const Map<point>& points,
const Map<edge>& edges, const Map<edge>& edges,
@ -235,7 +235,8 @@ namespace meshOps
inline void pWrite inline void pWrite
( (
const label toID, const label toID,
const label& data const label& data,
bool blocking = true
); );
// Parallel send (for fixed lists) // Parallel send (for fixed lists)
@ -258,7 +259,8 @@ namespace meshOps
inline void pRead inline void pRead
( (
const label fromID, const label fromID,
label& data label& data,
bool blocking = true
); );
// Parallel receive (for fixed lists) // Parallel receive (for fixed lists)
@ -282,7 +284,7 @@ namespace meshOps
// Select a list of elements from connectivity, // Select a list of elements from connectivity,
// and output to a VTK format // and output to a VTK format
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,
@ -299,7 +301,7 @@ namespace meshOps
); );
// Actual routine to write out the VTK file // Actual routine to write out the VTK file
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,
@ -320,6 +322,23 @@ namespace meshOps
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Add additional operations missing in ops.H
template<class T>
class lessOp
{ public: T operator()(const T& x, const T& y) const { return x < y; } };
template<class T>
class lessEqOp
{ public: T operator()(const T& x, const T& y) const { return x <= y; } };
template<class T>
class greaterOp
{ public: T operator()(const T& x, const T& y) const { return x > y; } };
template<class T>
class greaterEqOp
{ public: T operator()(const T& x, const T& y) const { return x >= y; } };
} // End namespace Foam } // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -327,7 +346,7 @@ namespace meshOps
#include "meshOpsI.H" #include "meshOpsI.H"
#ifdef NoRepository #ifdef NoRepository
# include "meshOpsTemplates.C" # include "meshOps.C"
#endif #endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -462,12 +462,13 @@ inline bool pointInTriFace
inline void pWrite inline void pWrite
( (
const label toID, const label toID,
const label& data const label& data,
bool blocking
) )
{ {
OPstream::write OPstream::write
( (
Pstream::blocking, (blocking ? Pstream::blocking : Pstream::nonBlocking),
toID, toID,
reinterpret_cast<const char*>(&data), reinterpret_cast<const char*>(&data),
sizeof(label) sizeof(label)
@ -479,12 +480,13 @@ inline void pWrite
inline void pRead inline void pRead
( (
const label fromID, const label fromID,
label& data label& data,
bool blocking
) )
{ {
IPstream::read IPstream::read
( (
Pstream::blocking, (blocking ? Pstream::blocking : Pstream::nonBlocking),
fromID, fromID,
reinterpret_cast<char*>(&data), reinterpret_cast<char*>(&data),
sizeof(label) sizeof(label)
@ -492,6 +494,79 @@ inline void pRead
} }
// Parallel non-blocking send for fixed lists
template <class Type, unsigned Size>
inline void pWrite
(
const label toID,
const FixedList<Type, Size>& data
)
{
OPstream::write
(
Pstream::blocking,
toID,
reinterpret_cast<const char*>(&data[0]),
Size*sizeof(Type)
);
}
// Parallel non-blocking receive for fixed lists
template <class Type, unsigned Size>
inline void pRead
(
const label fromID,
FixedList<Type, Size>& data
)
{
IPstream::read
(
Pstream::blocking,
fromID,
reinterpret_cast<char*>(data.begin()),
Size*sizeof(Type)
);
}
// Parallel non-blocking send for lists
template <class Type>
inline void pWrite
(
const label toID,
const UList<Type>& data
)
{
OPstream::write
(
Pstream::nonBlocking,
toID,
reinterpret_cast<const char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Parallel non-blocking receive for lists
template <class Type>
inline void pRead
(
const label fromID,
UList<Type>& data
)
{
IPstream::read
(
Pstream::nonBlocking,
fromID,
reinterpret_cast<char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Wait for buffer transfer completion.
inline void waitForBuffers() inline void waitForBuffers()
{ {
if (Pstream::parRun()) if (Pstream::parRun())
@ -588,6 +663,76 @@ inline void replaceLabel
} }
// Utility method to size-up the list to include an item
template <class Type>
inline void sizeUpList
(
const Type& item,
List<Type>& list
)
{
list.setSize(list.size() + 1, item);
}
// Utility method to size-down the list to remove an item
template <class Type>
inline void sizeDownList
(
const Type& item,
List<Type>& list
)
{
label index = -1;
if ((index = findIndex(list, item)) > -1)
{
meshOps::removeIndex(index, list);
}
else
{
FatalErrorIn
(
"inline void meshOps::sizeDownList"
"(const Type& item, List<Type>& list)"
)
<< nl << "Item: " << item
<< " was not found in list. " << nl
<< " List: " << nl << list
<< abort(FatalError);
}
}
// Remove an item at a particular index in the list
template <class Type>
inline void removeIndex
(
const label index,
List<Type>& list
)
{
// Create a new list
List<Type> newList(list.size() - 1);
// Copy individual items
label n = 0;
forAll(list, itemI)
{
if (itemI == index)
{
continue;
}
newList[n++] = list[itemI];
}
// Overwrite
list.transfer(newList);
}
} // End namespace meshOps } // End namespace meshOps

View file

@ -1,195 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
meshOps
Description
Various utility functions that perform mesh-related operations.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "Pair.H"
#include "meshOps.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace meshOps
{
// Parallel non-blocking send for fixed lists
template <class Type, unsigned Size>
inline void pWrite
(
const label toID,
const FixedList<Type, Size>& data
)
{
OPstream::write
(
Pstream::blocking,
toID,
reinterpret_cast<const char*>(&data[0]),
Size*sizeof(Type)
);
}
// Parallel non-blocking receive for fixed lists
template <class Type, unsigned Size>
inline void pRead
(
const label fromID,
FixedList<Type, Size>& data
)
{
IPstream::read
(
Pstream::blocking,
fromID,
reinterpret_cast<char*>(data.begin()),
Size*sizeof(Type)
);
}
// Parallel non-blocking send for lists
template <class Type>
inline void pWrite
(
const label toID,
const UList<Type>& data
)
{
OPstream::write
(
Pstream::nonBlocking,
toID,
reinterpret_cast<const char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Parallel non-blocking receive for lists
template <class Type>
inline void pRead
(
const label fromID,
UList<Type>& data
)
{
IPstream::read
(
Pstream::nonBlocking,
fromID,
reinterpret_cast<char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Utility method to size-up the list to include an item
template <class Type>
inline void sizeUpList
(
const Type& item,
List<Type>& list
)
{
list.setSize(list.size() + 1, item);
}
// Utility method to size-down the list to remove an item
template <class Type>
inline void sizeDownList
(
const Type& item,
List<Type>& list
)
{
label index = -1;
if ((index = findIndex(list, item)) > -1)
{
meshOps::removeIndex(index, list);
}
else
{
FatalErrorIn
(
"inline void meshOps::sizeDownList"
"(const Type& item, List<Type>& list)"
)
<< nl << "Item: " << item
<< " was not found in list. " << nl
<< " List: " << nl << list
<< abort(FatalError);
}
}
// Remove an item at a particular index in the list
template <class Type>
inline void removeIndex
(
const label index,
List<Type>& list
)
{
// Create a new list
List<Type> newList(list.size() - 1);
// Copy individual items
label n = 0;
forAll(list, itemI)
{
if (itemI == index)
{
continue;
}
newList[n++] = list[itemI];
}
// Overwrite
list.transfer(newList);
}
} // End namespace meshOps
} // End namespace Foam
// ************************************************************************* //

View file

@ -49,10 +49,18 @@ SourceFiles
// Have gcc ignore certain warnings while including mesquite headers // Have gcc ignore certain warnings while including mesquite headers
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) #if defined(__GNUC__) && !defined(__INTEL_COMPILER)
# pragma GCC diagnostic ignored "-Wold-style-cast" # pragma GCC diagnostic ignored "-Wold-style-cast"
# pragma GCC diagnostic ignored "-Wnon-virtual-dtor"
# pragma GCC diagnostic ignored "-Wunused-but-set-variable"
#endif #endif
#include "Mesquite_all_headers.hpp" #include "Mesquite_all_headers.hpp"
#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
# pragma GCC diagnostic warning "-Wold-style-cast"
# pragma GCC diagnostic warning "-Wnon-virtual-dtor"
# pragma GCC diagnostic warning "-Wunused-but-set-variable"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam namespace Foam
@ -88,6 +96,17 @@ class mesquiteMotionSolver
//- Number of elements //- Number of elements
unsigned long nCells_; unsigned long nCells_;
//- Polyhedra decomposition type
// 1: for cell, 2: for face
label decompType_;
//- Number of points from poly decomposition
unsigned long nDecompPoints_;
//- List of decomposition centres
List<labelPair> decompCellCentres_;
List<labelPair> decompFaceCentres_;
//- Number of auxiliary points //- Number of auxiliary points
labelList nAuxPoints_; labelList nAuxPoints_;
@ -106,7 +125,7 @@ class mesquiteMotionSolver
//- Specify max volume correction iterations //- Specify max volume correction iterations
label volCorrMaxIter_; label volCorrMaxIter_;
// Specify tolerance for the CG solver //- Specify tolerance for the CG solver
scalar tolerance_; scalar tolerance_;
//- Specify multiple mesh-motion sweeps //- Specify multiple mesh-motion sweeps
@ -159,6 +178,7 @@ class mesquiteMotionSolver
List<vectorField> gradEdge_; List<vectorField> gradEdge_;
List<vectorField> localPts_; List<vectorField> localPts_;
List<scalarField> edgeMarker_; List<scalarField> edgeMarker_;
List<scalarField> edgeConstant_;
//- Data for the auxiliary entities //- Data for the auxiliary entities
labelList procIndices_; labelList procIndices_;
@ -219,7 +239,7 @@ class mesquiteMotionSolver
scalar cmptSumMag(const vectorField& field); scalar cmptSumMag(const vectorField& field);
// Transfer buffers for surface point fields // Transfer buffers for surface point fields
void transferBuffers(vectorField &field); void transferBuffers(vectorField &field, bool fix = false);
// Apply boundary conditions // Apply boundary conditions
void applyBCs(vectorField &field); void applyBCs(vectorField &field);
@ -253,6 +273,14 @@ class mesquiteMotionSolver
// Prepare mesquite connectivity for parallel runs // Prepare mesquite connectivity for parallel runs
void initMesquiteParallelArrays(); void initMesquiteParallelArrays();
// Compute centroids based on up-to-date point positions
void computeCentroids
(
const pointField& points,
pointField& faceCentres,
pointField& cellCentres
);
// Copy auxiliary points to/from buffers // Copy auxiliary points to/from buffers
void copyAuxiliaryPoints(bool firstCopy); void copyAuxiliaryPoints(bool firstCopy);
@ -286,6 +314,9 @@ class mesquiteMotionSolver
// Prepare point-normals with updated point positions // Prepare point-normals with updated point positions
void preparePointNormals(); void preparePointNormals();
// Prepare non-uniform edge constants with updated point positions
void prepareEdgeConstants(const vectorField& p);
// Utility method to check validity of cells connected to a point. // Utility method to check validity of cells connected to a point.
bool checkValidity bool checkValidity
( (

View file

@ -15,6 +15,9 @@ FoamFile
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Specify dynamicFvMesh library
dynamicFvMeshLibs ("libdynamicTopoFvMesh.so");
//- Select the type of dynamicFvMesh //- Select the type of dynamicFvMesh
dynamicFvMesh dynamicTopoFvMesh; dynamicFvMesh dynamicTopoFvMesh;