subsetting point fields

This commit is contained in:
mattijs
2008-05-23 11:06:43 +01:00
parent bc84848b6d
commit ffa939952c
23 changed files with 544 additions and 196 deletions

View File

@ -60,7 +60,7 @@ void subsetVolFields
Info<< "Subsetting field " << fieldName << endl;
GeometricField<Type, fvPatchField, volMesh> volField
GeometricField<Type, fvPatchField, volMesh> fld
(
IOobject
(
@ -73,7 +73,7 @@ void subsetVolFields
baseMesh
);
subFields.set(i, subsetter.interpolate(volField));
subFields.set(i, subsetter.interpolate(fld));
}
}
@ -94,7 +94,7 @@ void subsetSurfaceFields
Info<< "Subsetting field " << fieldName << endl;
GeometricField<Type, fvsPatchField, surfaceMesh> volField
GeometricField<Type, fvsPatchField, surfaceMesh> fld
(
IOobject
(
@ -107,7 +107,42 @@ void subsetSurfaceFields
baseMesh
);
subFields.set(i, subsetter.interpolate(volField));
subFields.set(i, subsetter.interpolate(fld));
}
}
template<class Type>
void subsetPointFields
(
const fvMeshSubset& subsetter,
const pointMesh& pMesh,
const wordList& fieldNames,
PtrList<GeometricField<Type, pointPatchField, pointMesh> >& subFields
)
{
const fvMesh& baseMesh = subsetter.baseMesh();
forAll(fieldNames, i)
{
const word& fieldName = fieldNames[i];
Info<< "Subsetting field " << fieldName << endl;
GeometricField<Type, pointPatchField, pointMesh> fld
(
IOobject
(
fieldName,
baseMesh.time().timeName(),
baseMesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
pMesh
);
subFields.set(i, subsetter.interpolate(fld));
}
}
@ -163,7 +198,9 @@ int main(int argc, char *argv[])
IOobjectList objects(mesh, runTime.timeName());
// Read vol fields and subset.
// Read vol fields and subset
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
wordList scalarNames(objects.names(volScalarField::typeName));
PtrList<volScalarField> scalarFlds(scalarNames.size());
@ -191,7 +228,9 @@ int main(int argc, char *argv[])
PtrList<volTensorField> tensorFlds(tensorNames.size());
subsetVolFields(subsetter, tensorNames, tensorFlds);
// Read surface fields and subset.
// Read surface fields and subset
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
wordList surfScalarNames(objects.names(surfaceScalarField::typeName));
PtrList<surfaceScalarField> surfScalarFlds(surfScalarNames.size());
@ -231,6 +270,59 @@ int main(int argc, char *argv[])
subsetSurfaceFields(subsetter, surfTensorNames, surfTensorFlds);
// Read point fields and subset
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pointMesh pMesh(mesh);
wordList pointScalarNames(objects.names(pointScalarField::typeName));
PtrList<pointScalarField> pointScalarFlds(pointScalarNames.size());
subsetPointFields(subsetter, pMesh, pointScalarNames, pointScalarFlds);
wordList pointVectorNames(objects.names(pointVectorField::typeName));
PtrList<pointVectorField> pointVectorFlds(pointVectorNames.size());
subsetPointFields(subsetter, pMesh, pointVectorNames, pointVectorFlds);
wordList pointSphericalTensorNames
(
objects.names(pointSphericalTensorField::typeName)
);
PtrList<pointSphericalTensorField> pointSphericalTensorFlds
(
pointSphericalTensorNames.size()
);
subsetPointFields
(
subsetter,
pMesh,
pointSphericalTensorNames,
pointSphericalTensorFlds
);
wordList pointSymmTensorNames
(
objects.names(pointSymmTensorField::typeName)
);
PtrList<pointSymmTensorField> pointSymmTensorFlds
(
pointSymmTensorNames.size()
);
subsetPointFields
(
subsetter,
pMesh,
pointSymmTensorNames,
pointSymmTensorFlds
);
wordList pointTensorNames(objects.names(pointTensorField::typeName));
PtrList<pointTensorField> pointTensorFlds(pointTensorNames.size());
subsetPointFields(subsetter, pMesh, pointTensorNames, pointTensorFlds);
// Write mesh and fields to new time
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
runTime++;
@ -252,6 +344,18 @@ int main(int argc, char *argv[])
vectorFlds[i].write();
}
forAll(sphericalTensorFlds, i)
{
sphericalTensorFlds[i].rename(sphericalTensorNames[i]);
sphericalTensorFlds[i].write();
}
forAll(symmTensorFlds, i)
{
symmTensorFlds[i].rename(symmTensorNames[i]);
symmTensorFlds[i].write();
}
forAll(tensorFlds, i)
{
tensorFlds[i].rename(tensorNames[i]);
@ -272,6 +376,18 @@ int main(int argc, char *argv[])
surfVectorFlds[i].write();
}
forAll(surfSphericalTensorFlds, i)
{
surfSphericalTensorFlds[i].rename(surfSphericalTensorNames[i]);
surfSphericalTensorFlds[i].write();
}
forAll(surfSymmTensorFlds, i)
{
surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
surfSymmTensorFlds[i].write();
}
forAll(surfTensorNames, i)
{
surfTensorFlds[i].rename(surfTensorNames[i]);
@ -279,6 +395,39 @@ int main(int argc, char *argv[])
surfTensorFlds[i].write();
}
// Point ones
forAll(pointScalarFlds, i)
{
pointScalarFlds[i].rename(pointScalarNames[i]);
pointScalarFlds[i].write();
}
forAll(pointVectorFlds, i)
{
pointVectorFlds[i].rename(pointVectorNames[i]);
pointVectorFlds[i].write();
}
forAll(pointSphericalTensorFlds, i)
{
pointSphericalTensorFlds[i].rename(pointSphericalTensorNames[i]);
pointSphericalTensorFlds[i].write();
}
forAll(pointSymmTensorFlds, i)
{
pointSymmTensorFlds[i].rename(pointSymmTensorNames[i]);
pointSymmTensorFlds[i].write();
}
forAll(pointTensorNames, i)
{
pointTensorFlds[i].rename(pointTensorNames[i]);
pointTensorFlds[i].write();
}
Info << nl << "End" << endl;
return 0;

View File

@ -469,7 +469,6 @@ int main(int argc, char *argv[])
// Construct pointMesh only if nessecary since constructs edge
// addressing (expensive on polyhedral meshes)
autoPtr<pointMesh> pMeshPtr(NULL);
if (noPointValues)
{
Info<< " pointScalarFields : switched off"
@ -477,10 +476,6 @@ int main(int argc, char *argv[])
Info<< " pointVectorFields : switched off"
<< " (\"-noPointValues\" option)\n";
}
else
{
pMeshPtr.reset(new pointMesh(mesh));
}
PtrList<pointScalarField> psf;
PtrList<pointVectorField> pvf;
@ -488,21 +483,56 @@ int main(int argc, char *argv[])
PtrList<pointSymmTensorField> pSymmtf;
PtrList<pointTensorField> ptf;
if (pMeshPtr.valid())
if (!noPointValues)
{
readFields(pMeshPtr(), objects, selectedFields, psf);
readFields
(
vMesh,
vMesh.basePointMesh(),
objects,
selectedFields,
psf
);
print(" pointScalarFields :", Info, psf);
readFields(pMeshPtr(), objects, selectedFields, pvf);
readFields
(
vMesh,
vMesh.basePointMesh(),
objects,
selectedFields,
pvf
);
print(" pointVectorFields :", Info, pvf);
readFields(pMeshPtr(), objects, selectedFields, pSpheretf);
readFields
(
vMesh,
vMesh.basePointMesh(),
objects,
selectedFields,
pSpheretf
);
print(" pointSphericalTensorFields :", Info, pSpheretf);
readFields(pMeshPtr(), objects, selectedFields, pSymmtf);
readFields
(
vMesh,
vMesh.basePointMesh(),
objects,
selectedFields,
pSymmtf
);
print(" pointSymmTensorFields :", Info, pSymmtf);
readFields(pMeshPtr(), objects, selectedFields, ptf);
readFields
(
vMesh,
vMesh.basePointMesh(),
objects,
selectedFields,
ptf
);
print(" pointTensorFields :", Info, ptf);
}
Info<< endl;
@ -550,7 +580,7 @@ int main(int argc, char *argv[])
writer.write(vSymmtf);
writer.write(vtf);
if (pMeshPtr.valid())
if (!noPointValues)
{
writeFuns::writePointDataHeader
(
@ -567,7 +597,7 @@ int main(int argc, char *argv[])
writer.write(ptf);
// Interpolated volFields
volPointInterpolation pInterp(mesh, pMeshPtr());
volPointInterpolation pInterp(mesh, vMesh.pMesh());
writer.write(pInterp, vsf);
writer.write(pInterp, vvf);
writer.write(pInterp, vSpheretf);
@ -705,7 +735,7 @@ int main(int argc, char *argv[])
writer.write(vSymmtf);
writer.write(vtf);
if (pMeshPtr.valid())
if (!noPointValues)
{
writeFuns::writePointDataHeader
(
@ -785,7 +815,7 @@ int main(int argc, char *argv[])
writer.write(vSymmtf);
writer.write(vtf);
if (pMeshPtr.valid())
if (!noPointValues)
{
writeFuns::writePointDataHeader
(

View File

@ -34,50 +34,6 @@ namespace Foam
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
//template<class GeoField, class Mesh>
template<class GeoField>
void readFields
(
const typename GeoField::Mesh& mesh,
const IOobjectList& objects,
const HashSet<word>& selectedFields,
PtrList<GeoField>& fields
)
{
// Search list of objects for volScalarFields
IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
// Construct the vol scalar fields
fields.setSize(fieldObjects.size());
label nFields = 0;
for
(
IOobjectList::iterator iter = fieldObjects.begin();
iter != fieldObjects.end();
++iter
)
{
if (!selectedFields.size() || selectedFields.found(iter()->name()))
{
fields.set
(
nFields,
new GeoField
(
*iter(),
mesh
)
);
nFields++;
}
}
fields.setSize(nFields);
}
template<class GeoField>
void readFields
(

View File

@ -45,16 +45,6 @@ SourceFiles
namespace Foam
{
// Read the fields and put on the pointer list
template<class GeoField>
void readFields
(
const typename GeoField::Mesh& mesh,
const IOobjectList& objects,
const HashSet<word>& selectedFields,
PtrList<GeoField>& fields
);
// Read the fields and optionally subset and put on the pointer list
template<class GeoField>
void readFields

View File

@ -45,8 +45,7 @@ Foam::vtkMesh::vtkMesh
:
baseMesh_(baseMesh),
subsetter_(baseMesh),
setName_(setName),
topoPtr_(NULL)
setName_(setName)
{
if (setName.size() > 0)
{
@ -59,14 +58,6 @@ Foam::vtkMesh::vtkMesh
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::vtkMesh::~vtkMesh()
{
deleteDemandDrivenData(topoPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::polyMesh::readUpdateState Foam::vtkMesh::readUpdate()
@ -78,11 +69,11 @@ Foam::polyMesh::readUpdateState Foam::vtkMesh::readUpdate()
// Note: since fvMeshSubset has no movePoints() functionality reconstruct
// the subset even if only movement.
deleteDemandDrivenData(topoPtr_);
topoPtr_.clear();
if (setName_.size() > 0)
{
Pout<< "Subsetting mesh based on cellSet " << setName_ << endl;
Info<< "Subsetting mesh based on cellSet " << setName_ << endl;
// Read cellSet using whole mesh
cellSet currentSet(baseMesh_, setName_);

View File

@ -39,6 +39,7 @@ SourceFiles
#include "vtkTopo.H"
#include "fvMeshSubset.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -59,14 +60,17 @@ class vtkMesh
//- Reference to mesh
fvMesh& baseMesh_;
//- Subsetting engine
//- Demand driven pointMesh
mutable autoPtr<pointMesh> pointMeshPtr_;
//- Subsetting engine + sub-fvMesh
fvMeshSubset subsetter_;
//- Current cellSet (or empty)
const word setName_;
//- Current decomposition of topology
mutable vtkTopo* topoPtr_;
mutable autoPtr<vtkTopo> topoPtr_;
@ -87,11 +91,6 @@ public:
vtkMesh(fvMesh& baseMesh, const word& setName = "");
// Destructor
~vtkMesh();
// Member Functions
// Access
@ -102,6 +101,15 @@ public:
return baseMesh_;
}
const pointMesh& basePointMesh() const
{
if (!pointMeshPtr_.valid())
{
pointMeshPtr_.reset(new pointMesh(baseMesh_));
}
return pointMeshPtr_();
}
const fvMeshSubset& subsetter() const
{
return subsetter_;
@ -116,11 +124,11 @@ public:
//- topology
const vtkTopo& topo() const
{
if (!topoPtr_)
if (!topoPtr_.valid())
{
topoPtr_ = new vtkTopo(mesh());
topoPtr_.reset(new vtkTopo(mesh()));
}
return *topoPtr_;
return topoPtr_();
}
//- Access either mesh or submesh
@ -136,6 +144,19 @@ public:
}
}
//- Access either pointMesh of base or pointMesh of subset
const pointMesh& pMesh() const
{
if (useSubMesh())
{
return subsetter_.subPointMesh();
}
else
{
return basePointMesh();
}
}
//- Number of field cells
label nFieldCells() const
{
@ -145,7 +166,7 @@ public:
//- Number of field points
label nFieldPoints() const
{
return mesh().nPoints() + topo().addPointCellLabels().size();
return pMesh().size() + topo().addPointCellLabels().size();
}
@ -171,7 +192,6 @@ public:
return fld;
}
}
};

View File

@ -81,14 +81,17 @@ public:
// Public static data
// this must be consistent with the enumeration in "vtkCell.H"
static const label VTK_TRIANGLE = 5;
static const label VTK_POLYGON = 7;
static const label VTK_QUAD = 9;
enum vtkTypes
{
VTK_TRIANGLE = 5,
VTK_POLYGON = 7,
VTK_QUAD = 9,
static const label VTK_TETRA = 10;
static const label VTK_PYRAMID = 14;
static const label VTK_WEDGE = 13;
static const label VTK_HEXAHEDRON = 12;
VTK_TETRA = 10,
VTK_PYRAMID = 14,
VTK_WEDGE = 13,
VTK_HEXAHEDRON = 12,
};
// Constructors

View File

@ -66,7 +66,7 @@ basicSymmetryPointPatchField<Type>::basicSymmetryPointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{}

View File

@ -72,7 +72,7 @@ calculatedPointPatchField<Type>::calculatedPointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{}

View File

@ -65,7 +65,7 @@ coupledPointPatchField<Type>::coupledPointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{}

View File

@ -114,7 +114,7 @@ valuePointPatchField<Type>::valuePointPatchField
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(ptf, iF),
pointPatchField<Type>(p, iF),
Field<Type>(ptf, mapper)
{}

View File

@ -65,7 +65,7 @@ zeroGradientPointPatchField<Type>::zeroGradientPointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{}

View File

@ -81,10 +81,10 @@ cyclicPointPatchField<Type>::cyclicPointPatchField
const cyclicPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
coupledPointPatchField<Type>(ptf, iF),
coupledPointPatchField<Type>(ptf, p, iF, mapper),
cyclicPatch_(refCast<const cyclicPointPatch>(p))
{
if (!isType<cyclicPointPatch>(this->patch()))

View File

@ -81,7 +81,7 @@ emptyPointPatchField<Type>::emptyPointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{
if (!isType<emptyPointPatch>(this->patch()))
{

View File

@ -66,10 +66,10 @@ processorPointPatchField<Type>::processorPointPatchField
const processorPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
coupledPointPatchField<Type>(ptf, iF),
coupledPointPatchField<Type>(ptf, p, iF, mapper),
procPatch_(refCast<const processorPointPatch>(ptf.patch()))
{}

View File

@ -78,10 +78,10 @@ symmetryPointPatchField<Type>::symmetryPointPatchField
const symmetryPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
basicSymmetryPointPatchField<Type>(ptf, iF)
basicSymmetryPointPatchField<Type>(ptf, p, iF, mapper)
{
if (!isType<symmetryPointPatch>(this->patch()))
{

View File

@ -82,7 +82,7 @@ wedgePointPatchField<Type>::wedgePointPatchField
const pointPatchFieldMapper&
)
:
pointPatchField<Type>(ptf, iF)
pointPatchField<Type>(p, iF)
{
if (!isType<wedgePointPatch>(this->patch()))
{

View File

@ -82,10 +82,10 @@ globalPointPatchField<Type>::globalPointPatchField
const globalPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
coupledPointPatchField<Type>(ptf, iF),
coupledPointPatchField<Type>(ptf, p, iF, mapper),
globalPointPatch_(refCast<const globalPointPatch>(ptf.patch()))
{
if (!isType<globalPointPatch>(this->patch()))

View File

@ -62,10 +62,10 @@ slipPointPatchField<Type>::slipPointPatchField
const slipPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
basicSymmetryPointPatchField<Type>(ptf, iF)
basicSymmetryPointPatchField<Type>(ptf, p, iF, mapper)
{}

View File

@ -48,7 +48,7 @@ namespace Foam
bool Foam::fvMeshSubset::checkCellSubset() const
{
if (!fvMeshSubsetPtr_)
if (!fvMeshSubsetPtr_.valid())
{
FatalErrorIn("bool fvMeshSubset::checkCellSubset() const")
<< "Mesh subset not set. Please set the cell map using "
@ -267,7 +267,7 @@ void Foam::fvMeshSubset::subsetZones()
pz.name(),
subset(baseMesh().nPoints(), pz, pointMap()),
i,
fvMeshSubsetPtr_->pointZones()
fvMeshSubsetPtr_().pointZones()
);
}
@ -315,7 +315,7 @@ void Foam::fvMeshSubset::subsetZones()
subAddressing,
subFlipStatus,
i,
fvMeshSubsetPtr_->faceZones()
fvMeshSubsetPtr_().faceZones()
);
}
@ -333,13 +333,13 @@ void Foam::fvMeshSubset::subsetZones()
cz.name(),
subset(baseMesh().nCells(), cz, cellMap()),
i,
fvMeshSubsetPtr_->cellZones()
fvMeshSubsetPtr_().cellZones()
);
}
// Add the zones
fvMeshSubsetPtr_->addZones(pZonePtrs, fZonePtrs, cZonePtrs);
fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
}
@ -357,14 +357,6 @@ Foam::fvMeshSubset::fvMeshSubset(const fvMesh& baseMesh)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::fvMeshSubset::~fvMeshSubset()
{
deleteDemandDrivenData(fvMeshSubsetPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::fvMeshSubset::setCellSubset
@ -671,7 +663,9 @@ void Foam::fvMeshSubset::setCellSubset
// Make a new mesh
fvMeshSubsetPtr_ = new fvMesh
fvMeshSubsetPtr_.reset
(
new fvMesh
(
IOobject
(
@ -684,7 +678,9 @@ void Foam::fvMeshSubset::setCellSubset
newPoints,
newFaces,
newCells
)
);
pointMeshSubsetPtr_.clear();
// Add old patches
@ -700,7 +696,7 @@ void Foam::fvMeshSubset::setCellSubset
// Patch still exists. Add it
newBoundary[nNewPatches] = oldPatches[patchI].clone
(
fvMeshSubsetPtr_->boundaryMesh(),
fvMeshSubsetPtr_().boundaryMesh(),
nNewPatches,
boundaryPatchSizes[patchI],
patchStart
@ -723,7 +719,7 @@ void Foam::fvMeshSubset::setCellSubset
boundaryPatchSizes[oldInternalPatchID],
patchStart,
nNewPatches,
fvMeshSubsetPtr_->boundaryMesh()
fvMeshSubsetPtr_().boundaryMesh()
);
// The index for the first patch is -1 as it originates from
@ -738,7 +734,7 @@ void Foam::fvMeshSubset::setCellSubset
patchMap_.setSize(nNewPatches);
// Add the fvPatches
fvMeshSubsetPtr_->addFvPatches(newBoundary);
fvMeshSubsetPtr_().addFvPatches(newBoundary);
// Subset and add any zones
subsetZones();
@ -1166,7 +1162,9 @@ void Foam::fvMeshSubset::setLargeCellSubset
// not proper but cannot be avoided since otherwise surfaceInterpolation
// cannot find its fvSchemes (it will try to read e.g.
// system/region0SubSet/fvSchemes)
fvMeshSubsetPtr_ = new fvMesh
fvMeshSubsetPtr_.reset
(
new fvMesh
(
IOobject
(
@ -1180,8 +1178,9 @@ void Foam::fvMeshSubset::setLargeCellSubset
newFaces,
newCells,
syncPar // parallel synchronisation
)
);
pointMeshSubsetPtr_.clear();
// Add old patches
List<polyPatch*> newBoundary(nbSize);
@ -1251,7 +1250,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
// Clone (even if 0 size)
newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
(
fvMeshSubsetPtr_->boundaryMesh(),
fvMeshSubsetPtr_().boundaryMesh(),
nNewPatches,
newSize,
patchStart
@ -1282,7 +1281,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
boundaryPatchSizes[oldInternalPatchID],
patchStart,
nNewPatches,
fvMeshSubsetPtr_->boundaryMesh()
fvMeshSubsetPtr_().boundaryMesh()
);
//Pout<< " oldInternalFaces : "
@ -1310,7 +1309,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
// Patch still exists. Add it
newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
(
fvMeshSubsetPtr_->boundaryMesh(),
fvMeshSubsetPtr_().boundaryMesh(),
nNewPatches,
newSize,
patchStart
@ -1331,7 +1330,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
// Add the fvPatches
fvMeshSubsetPtr_->addFvPatches(newBoundary);
fvMeshSubsetPtr_().addFvPatches(newBoundary);
// Subset and add any zones
subsetZones();
@ -1359,7 +1358,7 @@ const fvMesh& Foam::fvMeshSubset::subMesh() const
{
checkCellSubset();
return *fvMeshSubsetPtr_;
return fvMeshSubsetPtr_();
}
@ -1367,7 +1366,27 @@ fvMesh& Foam::fvMeshSubset::subMesh()
{
checkCellSubset();
return *fvMeshSubsetPtr_;
return fvMeshSubsetPtr_();
}
const pointMesh& Foam::fvMeshSubset::subPointMesh() const
{
if (!pointMeshSubsetPtr_.valid())
{
pointMeshSubsetPtr_.reset(new pointMesh(subMesh()));
}
return pointMeshSubsetPtr_();
}
pointMesh& Foam::fvMeshSubset::subPointMesh()
{
if (!pointMeshSubsetPtr_.valid())
{
pointMeshSubsetPtr_.reset(new pointMesh(subMesh()));
}
return pointMeshSubsetPtr_();
}

View File

@ -56,11 +56,11 @@ SourceFiles
#define fvMeshSubset_H
#include "fvMesh.H"
#include "pointMesh.H"
#include "fvPatchFieldMapper.H"
#include "pointPatchFieldMapper.H"
#include "GeometricField.H"
#include "emptyFvPatchFields.H"
#include "labelHashSet.H"
#include "SubField.H"
#include "surfaceMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -119,6 +119,48 @@ public:
};
//- Patch-field subset interpolation class
class pointPatchFieldSubset
:
public pointPatchFieldMapper
{
const labelList& directAddressing_;
public:
// Constructors
//- Construct given addressing
pointPatchFieldSubset(const labelList& directAddressing)
:
directAddressing_(directAddressing)
{}
// Destructor
virtual ~pointPatchFieldSubset()
{}
// Member Functions
label size() const
{
return directAddressing_.size();
}
bool direct() const
{
return true;
}
const unallocLabelList& directAddressing() const
{
return directAddressing_;
}
};
private:
// Private data
@ -127,7 +169,9 @@ private:
const fvMesh& baseMesh_;
//- Subset mesh pointer
fvMesh* fvMeshSubsetPtr_;
autoPtr<fvMesh> fvMeshSubsetPtr_;
mutable autoPtr<pointMesh> pointMeshSubsetPtr_;
//- Point mapping array
labelList pointMap_;
@ -185,11 +229,6 @@ public:
explicit fvMeshSubset(const fvMesh&);
// Destructor
~fvMeshSubset();
// Member Functions
// Edit
@ -236,8 +275,14 @@ public:
//- Return reference to subset mesh
const fvMesh& subMesh() const;
fvMesh& subMesh();
//- Return reference to demand-driven subset pointMesh
const pointMesh& subPointMesh() const;
pointMesh& subPointMesh();
//- Return point map
const labelList& pointMap() const;
@ -289,6 +334,25 @@ public:
(
const GeometricField<Type, fvsPatchField, surfaceMesh>&
) const;
//- Map point field
template<class Type>
static tmp<GeometricField<Type, pointPatchField, pointMesh> >
interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>&,
const pointMesh& sMesh,
const objectRegistry& reg,
const labelList& patchMap,
const labelList& pointMap
);
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh> >
interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>&
) const;
};

View File

@ -26,6 +26,8 @@ License
#include "fvMeshSubset.H"
#include "emptyFvsPatchField.H"
#include "emptyPointPatchField.H"
#include "emptyFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -134,6 +136,23 @@ tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
}
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
return interpolate
(
vf,
subMesh(),
patchMap(),
cellMap(),
faceMap()
);
}
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
(
@ -258,23 +277,6 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
}
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
return interpolate
(
vf,
subMesh(),
patchMap(),
cellMap(),
faceMap()
);
}
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
(
@ -291,6 +293,130 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
}
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh> >
fvMeshSubset::interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>& vf,
const pointMesh& sMesh,
const objectRegistry& reg,
const labelList& patchMap,
const labelList& pointMap
)
{
// Create and map the internal-field values
Field<Type> internalField(vf.internalField(), pointMap);
// Create and map the patch field values
PtrList<pointPatchField<Type> > patchFields(patchMap.size());
forAll (patchFields, patchI)
{
// Set the first one by hand as it corresponds to the
// exposed internal faces. Additional interpolation can be put here
// as necessary.
if (patchMap[patchI] == -1)
{
patchFields.set
(
patchI,
new emptyPointPatchField<Type>
(
sMesh.boundary()[patchI],
DimensionedField<Type, pointMesh>::null()
)
);
}
else
{
// Construct addressing
const pointPatch& basePatch =
vf.mesh().boundary()[patchMap[patchI]];
const labelList& meshPoints = basePatch.meshPoints();
// Make addressing from mesh to patch point
Map<label> meshPointMap(2*meshPoints.size());
forAll(meshPoints, localI)
{
meshPointMap.insert(meshPoints[localI], localI);
}
// Find which subpatch points originate from which patch point
const pointPatch& subPatch = sMesh.boundary()[patchI];
const labelList& subMeshPoints = subPatch.meshPoints();
// If mapped from outside patch use point 0 for lack of better.
labelList directAddressing(subPatch.size(), 0);
forAll(subMeshPoints, localI)
{
// Get mesh point on original mesh.
label meshPointI = pointMap[subMeshPoints[localI]];
Map<label>::const_iterator iter = meshPointMap.find(meshPointI);
if (iter != meshPointMap.end())
{
directAddressing[localI] = iter();
}
}
patchFields.set
(
patchI,
pointPatchField<Type>::New
(
vf.boundaryField()[patchMap[patchI]],
subPatch,
DimensionedField<Type, pointMesh>::null(),
pointPatchFieldSubset(directAddressing)
)
);
}
}
// Create the complete field from the pieces
tmp<GeometricField<Type, pointPatchField, pointMesh> > tresF
(
new GeometricField<Type, pointPatchField, pointMesh>
(
IOobject
(
"subset"+vf.name(),
vf.time().timeName(),
reg,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
vf.dimensions(),
internalField,
patchFields
)
);
return tresF;
}
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh> > fvMeshSubset::interpolate
(
const GeometricField<Type, pointPatchField, pointMesh>& sf
) const
{
return interpolate
(
sf,
subPointMesh(), // subsetted point mesh
subMesh(), // registry (pointfields are stored on the polyMesh)
patchMap(),
pointMap()
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -77,7 +77,7 @@ surfaceSlipDisplacementPointPatchVectorField
const dictionary& dict
)
:
pointPatchVectorField(p, iF),
pointPatchVectorField(p, iF, dict),
surfaceNames_(dict.lookup("projectSurfaces")),
projectMode_(followModeNames_.read(dict.lookup("followMode"))),
projectDir_(dict.lookup("projectDirection")),