zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen
2009-03-12 10:39:22 +01:00
parent 34da731896 471d45f39f
commit e011328d99
45 changed files with 1838 additions and 920 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
View File

@ -38,7 +38,8 @@ Foam::scalar Foam::compressibleCourantNo
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
if (mesh.nInternalFaces())
//- Can have fluid domains with 0 cells so do not test.
//if (mesh.nInternalFaces())
{
surfaceScalarField SfUfbyDelta =
mesh.surfaceInterpolation::deltaCoeffs()

21
applications/solvers/incompressible/pisoFoam/pisoFoam.C
View File

@ -23,7 +23,7 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
turbFoam
pisoFoam
Description
Transient solver for incompressible flow.
@ -40,15 +40,14 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
# include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -56,8 +55,8 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "CourantNo.H"
#include "readPISOControls.H"
#include "CourantNo.H"
// Pressure-velocity PISO corrector
{
@ -120,7 +119,7 @@ int main(int argc, char *argv[])
}
}
# include "continuityErrs.H"
#include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();

13
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
View File

@ -131,6 +131,12 @@ castellatedMeshControls
level (3 3);
}
}
// Optional angle to detect small-large cell situation perpendicular
// to the surface. Is the angle of face w.r.t the local surface
// normal. Use on flat(ish) surfaces only. Otherwise
// leave out or set to negative number.
//perpendicularAngle 10;
}
}
@ -317,9 +323,10 @@ meshQualityControls
minTriangleTwist -1;
//- if >0 : preserve single cells with all points on the surface if the
// resulting volume after snapping is larger than minVolFraction times old
// volume. If <0 : delete always.
minVolFraction 0.1;
// resulting volume after snapping (by approximation) is larger than
// minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
// If <0 : delete always.
//minVolCollapseRatio 0.5;
// Advanced

61
applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
View File

@ -102,6 +102,67 @@ Foam::label Foam::checkTopology
}
}
if (allTopology)
{
labelList nInternalFaces(mesh.nCells(), 0);
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
nInternalFaces[mesh.faceOwner()[faceI]]++;
nInternalFaces[mesh.faceNeighbour()[faceI]]++;
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
if (patches[patchI].coupled())
{
const unallocLabelList& owners = patches[patchI].faceCells();
forAll(owners, i)
{
nInternalFaces[owners[i]]++;
}
}
}
faceSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
faceSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);
forAll(nInternalFaces, cellI)
{
if (nInternalFaces[cellI] <= 1)
{
oneCells.insert(cellI);
}
else if (nInternalFaces[cellI] == 2)
{
twoCells.insert(cellI);
}
}
label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
if (nOneCells > 0)
{
Info<< " <<Writing " << nOneCells
<< " cells with with single non-boundary face to set "
<< oneCells.name()
<< endl;
oneCells.write();
}
label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
if (nTwoCells > 0)
{
Info<< " <<Writing " << nTwoCells
<< " cells with with single non-boundary face to set "
<< twoCells.name()
<< endl;
twoCells.write();
}
}
{
regionSplit rs(mesh);

6
applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
View File

@ -65,10 +65,11 @@ SourceFiles
#include "className.H"
#include "fileName.H"
#include "volPointInterpolation.H"
#include "stringList.H"
#include "wordList.H"
#include "primitivePatch.H"
#include "PrimitivePatchInterpolation.H"
#include "volPointInterpolation.H"
// * * * * * * * * * * * * * Forward Declarations * * * * * * * * * * * * * //
@ -496,7 +497,6 @@ class vtkPV3Foam
void convertVolFields
(
const fvMesh&,
const volPointInterpolation&,
const PtrList<PrimitivePatchInterpolation<primitivePatch> >&,
const IOobjectList&,
vtkMultiBlockDataSet* output
@ -507,7 +507,7 @@ class vtkPV3Foam
void convertVolFieldBlock
(
const GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, pointPatchField, pointMesh>&,
autoPtr<GeometricField<Type, pointPatchField, pointMesh> >&,
vtkMultiBlockDataSet* output,
const partInfo& selector,
const List<polyDecomp>& decompLst

12
applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFields.C
View File

@ -106,8 +106,6 @@ void Foam::vtkPV3Foam::convertVolFields
printMemory();
}
// Construct interpolation on the raw mesh
volPointInterpolation pInterp(mesh);
PtrList<PrimitivePatchInterpolation<primitivePatch> >
ppInterpList(mesh.boundaryMesh().size());
@ -127,23 +125,23 @@ void Foam::vtkPV3Foam::convertVolFields
convertVolFields<scalar>
(
mesh, pInterp, ppInterpList, objects, output
mesh, ppInterpList, objects, output
);
convertVolFields<vector>
(
mesh, pInterp, ppInterpList, objects, output
mesh, ppInterpList, objects, output
);
convertVolFields<sphericalTensor>
(
mesh, pInterp, ppInterpList, objects, output
mesh, ppInterpList, objects, output
);
convertVolFields<symmTensor>
(
mesh, pInterp, ppInterpList, objects, output
mesh, ppInterpList, objects, output
);
convertVolFields<tensor>
(
mesh, pInterp, ppInterpList, objects, output
mesh, ppInterpList, objects, output
);
if (debug)

82
applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
View File

@ -34,6 +34,8 @@ InClass
#include "emptyFvPatchField.H"
#include "wallPolyPatch.H"
#include "faceSet.H"
#include "volPointInterpolation.H"
#include "vtkPV3FoamFaceField.H"
#include "vtkPV3FoamPatchField.H"
@ -43,7 +45,6 @@ template<class Type>
void Foam::vtkPV3Foam::convertVolFields
(
const fvMesh& mesh,
const volPointInterpolation& pInterp,
const PtrList<PrimitivePatchInterpolation<primitivePatch> >& ppInterpList,
const IOobjectList& objects,
vtkMultiBlockDataSet* output
@ -63,23 +64,22 @@ void Foam::vtkPV3Foam::convertVolFields
continue;
}
// Load field
GeometricField<Type, fvPatchField, volMesh> tf
(
*iter(),
mesh
);
tmp<GeometricField<Type, pointPatchField, pointMesh> > tptf
(
pInterp.interpolate(tf)
);
// Interpolated field (demand driven)
autoPtr<GeometricField<Type, pointPatchField, pointMesh> > ptfPtr;
// Convert activated internalMesh regions
convertVolFieldBlock
(
tf,
tptf(),
ptfPtr,
output,
partInfoVolume_,
regionPolyDecomp_
@ -89,7 +89,7 @@ void Foam::vtkPV3Foam::convertVolFields
convertVolFieldBlock
(
tf,
tptf(),
ptfPtr,
output,
partInfoCellZones_,
zonePolyDecomp_
@ -99,7 +99,7 @@ void Foam::vtkPV3Foam::convertVolFields
convertVolFieldBlock
(
tf,
tptf(),
ptfPtr,
output,
partInfoCellSets_,
csetPolyDecomp_
@ -258,46 +258,60 @@ void Foam::vtkPV3Foam::convertVolFields
}
}
template<class Type>
void Foam::vtkPV3Foam::convertVolFieldBlock
(
const GeometricField<Type, fvPatchField, volMesh>& tf,
const GeometricField<Type, pointPatchField, pointMesh>& ptf,
autoPtr<GeometricField<Type, pointPatchField, pointMesh> >& ptfPtr,
vtkMultiBlockDataSet* output,
const partInfo& selector,
const List<polyDecomp>& decompLst
)
{
for (int partId = selector.start(); partId < selector.end(); ++partId)
{
const label datasetNo = partDataset_[partId];
for (int partId = selector.start(); partId < selector.end(); ++partId)
{
const label datasetNo = partDataset_[partId];
if (datasetNo >= 0 && partStatus_[partId])
{
convertVolField
(
tf,
output,
selector,
datasetNo,
decompLst[datasetNo]
);
if (datasetNo >= 0 && partStatus_[partId])
{
convertVolField
(
tf,
output,
selector,
datasetNo,
decompLst[datasetNo]
);
convertPointField
(
ptf,
tf,
output,
selector,
datasetNo,
decompLst[datasetNo]
);
}
}
if (!ptfPtr.valid())
{
if (debug)
{
Info<< "convertVolFieldBlock interpolating:" << tf.name()
<< endl;
}
ptfPtr.reset
(
volPointInterpolation::New(tf.mesh()).interpolate(tf).ptr()
);
}
convertPointField
(
ptfPtr(),
tf,
output,
selector,
datasetNo,
decompLst[datasetNo]
);
}
}
}
template<class Type>
void Foam::vtkPV3Foam::convertVolField
(

2
applications/utilities/surface/surfaceConvert/surfaceConvert.C
View File

@ -115,7 +115,7 @@ int main(int argc, char *argv[])
}
Info<< "writing " << exportName;
if (scaleFactor > 0)
if (scaleFactor <= 0)
{
Info<< " without scaling" << endl;
}

2
bin/foamEbrowse
View File

@ -47,7 +47,7 @@ cd $WM_PROJECT_DIR
mkdir .tags 2>/dev/null
cd .tags
find -H .. \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) -print > $sourcesFile
find -H .. \( -name "*.[HC]" -o -name lnInclude -prune -o -name Doxygen -prune \) -print > $sourcesFile
ebrowse --files=$sourcesFile --output-file=ebrowse
#------------------------------------------------------------------------------

2
bin/foamPackSource
View File

@ -77,9 +77,9 @@ find -H $packDir \
-a ! -name "core.[1-9]*" \
-a ! -name "log[0-9]*" \
-a ! -name "libccmio*" \
-a ! -name "\.ebrowse" \
| sed \
-e "\@$packDir/.git/@d" \
-e "\@/.tags/@d" \
-e "\@$packDir/lib/@d" \
-e "\@libccmio.*/@d" \
-e '\@applications/bin/@d' \

21
bin/foamTags
View File

@ -42,15 +42,24 @@ fi
cd $WM_PROJECT_DIR
mkdir .tags 2>/dev/null
find -H . \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) | \
etags --declarations -l c++ -o .tags/etags -
find -H . \( -name "*.[HC]" -not -name "lnInclude" -not -name "Doxygen" \) | \
etags -l c++ -o .tags/etagsDef -
find -H . \( -name "*.H" -not -name "lnInclude" -not -name "Doxygen" \) | \
etags --declarations -l c++ -o .tags/etagsDec -
#etagsCmd="etags --declarations -l c++ -o .tags/etags -"
#etagsDefCmd="etags -l c++ -o .tags/etagsDef -"
#etagsDecCmd="etags --declarations -l c++ -o .tags/etagsDec -"
etagsCmd="ectags -e --extra=+fq --file-scope=no --c-kinds=+p -o .tags/etags -L -"
etagsDefCmd="ectags -e --extra=+fq --file-scope=no -o .tags/etagsDef -L -"
etagsDecCmd="ectags -e --extra=+fq --file-scope=no --c-kinds=+p -o .tags/etagsDec -L -"
ectagsDecCmd="ectags -o .tags/ectagsDec --file-scope=no --c-kinds=+p --excmd=n --extra=+fq --fields=+afiKmnsSzt -L -"
find -H . \( -name "*.[HC]" -o -name lnInclude -prune -o -name Doxygen -prune \) | $etagsCmd
find -H . \( -name "*.[HC]" -o -name lnInclude -prune -o -name Doxygen -prune \) | $etagsDefCmd
find -H . \( -name "*.H" -o -name lnInclude -prune -o -name Doxygen -prune \) | $etagsDecCmd
find -H . \( -name "*.H" -o -name lnInclude -prune -o -name Doxygen -prune \) | $ectagsDecCmd
gtags -i --gtagsconf bin/tools/gtagsrc .tags
foamEbrowse
#------------------------------------------------------------------------------

2
src/OpenFOAM/containers/Lists/FixedList/FixedListI.H
View File

@ -162,7 +162,7 @@ inline void Foam::FixedList<T, Size>::checkIndex(const label i) const
<< "attempt to access element from zero-sized list"
<< abort(FatalError);
}
else if (i < 0 || i >= Size)
else if (i < 0 || i >= label(Size))
{
FatalErrorIn("FixedList<T, Size>::checkIndex(const label)")
<< "index " << i << " out of range 0 ... " << (Size-1)

42
src/OpenFOAM/db/objectRegistry/objectRegistry.C
View File

@ -55,7 +55,8 @@ Foam::objectRegistry::objectRegistry
HashTable<regIOobject*>(nIoObjects),
time_(t),
parent_(t),
dbDir_(name())
dbDir_(name()),
event_(1)
{}
@ -69,7 +70,8 @@ Foam::objectRegistry::objectRegistry
HashTable<regIOobject*>(nIoObjects),
time_(io.time()),
parent_(io.db()),
dbDir_(parent_.dbDir()/local()/name())
dbDir_(parent_.dbDir()/local()/name()),
event_(1)
{
writeOpt() = IOobject::AUTO_WRITE;
}
@ -135,6 +137,42 @@ const Foam::objectRegistry& Foam::objectRegistry::subRegistry
}
Foam::label Foam::objectRegistry::getEvent() const
{
label curEvent = event_++;
if (event_ == labelMax)
{
WarningIn("objectRegistry::getEvent() const")
<< "Event counter has overflowed. Resetting counter on all"
<< " dependent objects." << endl
<< "This might cause extra evaluations." << endl;
// Reset event counter
curEvent = 1;
event_ = 2;
for (const_iterator iter = begin(); iter != end(); ++iter)
{
const regIOobject& io = *iter();
if (objectRegistry::debug)
{
Pout<< "objectRegistry::getEvent() : "
<< "resetting count on " << iter.key() << endl;
}
if (io.eventNo() != 0)
{
const_cast<regIOobject&>(io).eventNo() = curEvent;
}
}
}
return curEvent;
}
bool Foam::objectRegistry::checkIn(regIOobject& io) const
{
if (objectRegistry::debug)

6
src/OpenFOAM/db/objectRegistry/objectRegistry.H
View File

@ -64,6 +64,9 @@ class objectRegistry
//- Local directory path of this objectRegistry relative to time
fileName dbDir_;
//- Current event
mutable label event_;
// Private Member Functions
@ -151,6 +154,9 @@ public:
template<class Type>
const Type& lookupObject(const word& name) const;
//- Return new event number.
label getEvent() const;
// Edit

88
src/OpenFOAM/db/regIOobject/regIOobject.C
View File

@ -47,6 +47,7 @@ Foam::regIOobject::regIOobject(const IOobject& io)
registered_(false),
ownedByRegistry_(false),
lastModified_(0),
eventNo_(db().getEvent()),
isPtr_(NULL)
{
// Register with objectRegistry if requested
@ -64,6 +65,7 @@ Foam::regIOobject::regIOobject(const regIOobject& rio)
registered_(false),
ownedByRegistry_(false),
lastModified_(rio.lastModified_),
eventNo_(db().getEvent()),
isPtr_(NULL)
{
// Do not register copy with objectRegistry
@ -78,6 +80,7 @@ Foam::regIOobject::regIOobject(const regIOobject& rio, bool registerCopy)
registered_(false),
ownedByRegistry_(false),
lastModified_(rio.lastModified_),
eventNo_(db().getEvent()),
isPtr_(NULL)
{
if (registerCopy && rio.registered_)
@ -153,6 +156,91 @@ bool Foam::regIOobject::checkOut()
}
bool Foam::regIOobject::uptodate(const word& a) const
{
if (db().lookupObject<regIOobject>(a).eventNo() >= eventNo_)
{
return false;
}
else
{
return true;
}
}
bool Foam::regIOobject::uptodate(const word& a, const word& b) const
{
if
(
db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
)
{
return false;
}
else
{
return true;
}
}
bool Foam::regIOobject::uptodate
(
const word& a,
const word& b,
const word& c
) const
{
if
(
db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(c).eventNo() >= eventNo_
)
{
return false;
}
else
{
return true;
}
}
bool Foam::regIOobject::uptodate
(
const word& a,
const word& b,
const word& c,
const word& d
) const
{
if
(
db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(c).eventNo() >= eventNo_
|| db().lookupObject<regIOobject>(d).eventNo() >= eventNo_
)
{
return false;
}
else
{
return true;
}
}
//- Flag me as uptodate
void Foam::regIOobject::setUptodate()
{
eventNo_ = db().getEvent();
}
// Rename object and re-register with objectRegistry under new name
void Foam::regIOobject::rename(const word& newName)
{

22
src/OpenFOAM/db/regIOobject/regIOobject.H
View File

@ -45,6 +45,7 @@ SourceFiles
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
@ -70,6 +71,9 @@ private:
//- Time of last modification
mutable time_t lastModified_;
//- eventNo of last update
label eventNo_;
//- Istream for reading
Istream* isPtr_;
@ -141,6 +145,24 @@ public:
inline void release();
// Dependency checking
//- Event number at last update.
inline label eventNo() const;
//- Event number at last update.
inline label& eventNo();
//- Am I uptodate with respect to other regIOobjects
bool uptodate(const word&) const;
bool uptodate(const word&, const word&) const;
bool uptodate(const word&, const word&, const word&) const;
bool uptodate(const word&, const word&, const word&, const word&)
const;
//- Flag me as uptodate
void setUptodate();
// Edit
//- Rename

11
src/OpenFOAM/db/regIOobject/regIOobjectI.H
View File

@ -80,4 +80,15 @@ void Foam::regIOobject::release()
}
Foam::label Foam::regIOobject::eventNo() const
{
return eventNo_;
}
Foam::label& Foam::regIOobject::eventNo()
{
return eventNo_;
}
// ************************************************************************* //

120
src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
View File

@ -97,6 +97,68 @@ slicedBoundaryField
}
template
<
class Type,
template<class> class PatchField,
template<class> class SlicedPatchField,
class GeoMesh
>
Foam::tmp<Foam::FieldField<PatchField, Type> >
Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
slicedBoundaryField
(
const Mesh& mesh,
const FieldField<PatchField, Type>& bField,
const bool preserveCouples
)
{
tmp<FieldField<PatchField, Type> > tbf
(
new FieldField<PatchField, Type>(mesh.boundary().size())
);
FieldField<PatchField, Type>& bf = tbf();
forAll (mesh.boundary(), patchi)
{
if (preserveCouples && mesh.boundary()[patchi].coupled())
{
// For coupled patched construct the correct patch field type
bf.set
(
patchi,
PatchField<Type>::New
(
mesh.boundary()[patchi].type(),
mesh.boundary()[patchi],
*this
)
);
// Assign field
bf[patchi] == bField[patchi];
}
else
{
// Create unallocated copy of patch field
bf.set
(
patchi,
new SlicedPatchField<Type>
(
mesh.boundary()[patchi],
DimensionedField<Type, GeoMesh>::null()
)
);
bf[patchi].UList<Type>::operator=(bField[patchi]);
}
}
return tbf;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template
@ -204,6 +266,64 @@ SlicedGeometricField
}
template
<
class Type,
template<class> class PatchField,
template<class> class SlicedPatchField,
class GeoMesh
>
Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
SlicedGeometricField
(
const IOobject& io,
const GeometricField<Type, PatchField, GeoMesh>& gf,
const bool preserveCouples
)
:
GeometricField<Type, PatchField, GeoMesh>
(
io,
gf.mesh(),
gf.dimensions(),
Field<Type>(),
slicedBoundaryField(gf.mesh(), gf.boundaryField(), preserveCouples)
)
{
// Set the internalField to the supplied internal field
UList<Type>::operator=(gf.internalField());
correctBoundaryConditions();
}
template
<
class Type,
template<class> class PatchField,
template<class> class SlicedPatchField,
class GeoMesh
>
Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
SlicedGeometricField
(
const SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>& gf
)
:
GeometricField<Type, PatchField, GeoMesh>
(
gf,
gf.mesh(),
gf.dimensions(),
Field<Type>(),
slicedBoundaryField(gf.mesh(), gf.boundaryField(), true)
)
{
// Set the internalField to the supplied internal field
UList<Type>::operator=(gf.internalField());
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template

34
src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
View File

@ -87,8 +87,17 @@ private:
const bool preserveCouples
);
//- Disallow default bitwise copy construct
SlicedGeometricField(const SlicedGeometricField&);
//- Slice the given field and a create a PtrList of SlicedPatchField
// from which the boundary field is built
tmp<FieldField<PatchField, Type> > slicedBoundaryField
(
const Mesh& mesh,
const FieldField<PatchField, Type>& bField,
const bool preserveCouples
);
////- Disallow default bitwise copy construct
//SlicedGeometricField(const SlicedGeometricField&);
//- Disallow default bitwise assignment
void operator=(const SlicedGeometricField&);
@ -128,6 +137,27 @@ public:
const bool preserveCouples=true
);
//- Construct from GeometricField. Reuses full internal and
// patch fields except on couples (preserveCouples=true).
SlicedGeometricField
(
const IOobject&,
const GeometricField<Type, PatchField, GeoMesh>&,
const bool preserveCouples=true
);
//- Construct as copy
SlicedGeometricField
(
const SlicedGeometricField
<
Type,
PatchField,
SlicedPatchField,
GeoMesh
>&
);
// Destructor

12
src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
View File

@ -268,7 +268,12 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
// Derived info
bounds_ = boundBox(points_);
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Zones
pointZoneMesh newPointZones
(
IOobject
@ -418,6 +423,13 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
false
)
);
// Derived info
bounds_ = boundBox(points_);
// Rotation can cause direction vector to change
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
return polyMesh::POINTS_MOVED;
}

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

@ -58,35 +58,6 @@ defineTypeNameAndDebug(autoSnapDriver, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::autoSnapDriver::getZonedSurfaces
(
labelList& zonedSurfaces,
labelList& unzonedSurfaces
) const
{ // Surfaces with zone information
const wordList& faceZoneNames = meshRefiner_.surfaces().faceZoneNames();
zonedSurfaces.setSize(faceZoneNames.size());
label zonedI = 0;
unzonedSurfaces.setSize(faceZoneNames.size());
label unzonedI = 0;
forAll(faceZoneNames, surfI)
{
if (faceZoneNames[surfI].size())
{
zonedSurfaces[zonedI++] = surfI;
}
else
{
unzonedSurfaces[unzonedI++] = surfI;
}
}
zonedSurfaces.setSize(zonedI);
unzonedSurfaces.setSize(unzonedI);
}
// Get faces to repatch. Returns map from face to patch.
Foam::Map<Foam::label> Foam::autoSnapDriver::getZoneBafflePatches
(
@ -711,9 +682,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::createZoneBaffles
List<labelPair>& baffles
)
{
labelList zonedSurfaces;
labelList unzonedSurfaces;
getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
autoPtr<mapPolyMesh> map;
@ -798,9 +767,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
const List<labelPair>& baffles
)
{
labelList zonedSurfaces;
labelList unzonedSurfaces;
getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
autoPtr<mapPolyMesh> map;
@ -1048,9 +1015,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
labelList snapSurf(localPoints.size(), -1);
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces;
labelList unzonedSurfaces;
getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
labelList zonedSurfaces =
meshRefiner_.surfaces().getNamedSurfaces();
labelList unzonedSurfaces =
meshRefiner_.surfaces().getUnnamedSurfaces();
// 1. All points to non-interface surfaces
@ -1368,9 +1336,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
indirectPrimitivePatch& pp = ppPtr();
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces;
labelList unzonedSurfaces;
getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
// Faces that do not move

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

@ -84,9 +84,6 @@ class autoSnapDriver
// Snapping
//- Split surfaces into non-zoned and zones ones
void getZonedSurfaces(labelList&, labelList&) const;
//- Get faces to repatch. Returns map from face to patch.
Map<label> getZoneBafflePatches(const bool allowBoundary) const;

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

@ -374,11 +374,28 @@ private:
const labelList&
) const;
bool isCollapsedFace
(
const pointField&,
const pointField& neiCc,
const scalar minFaceArea,
const scalar maxNonOrtho,
const label faceI
) const;
bool isCollapsedCell
(
const pointField&,
const scalar volFraction,
const label cellI
) const;
//- Returns list with for every internal face -1 or the patch
// they should be baffled into. If removeEdgeConnectedCells is set
// removes cells based on perpendicularAngle.
labelList markFacesOnProblemCells
(
const dictionary& motionDict,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const labelList& globalToPatch

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

@ -1503,6 +1503,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
(
markFacesOnProblemCells
(
motionDict,
removeEdgeConnectedCells,
perpendicularAngle,
globalToPatch
@ -1548,6 +1549,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
(
markFacesOnProblemCells
(
motionDict,
removeEdgeConnectedCells,
perpendicularAngle,
globalToPatch

727
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -278,13 +278,105 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
}
// Check if moving face to new points causes a 'collapsed' face.
// Uses new point position only for the face, not the neighbouring
// cell centres
bool Foam::meshRefinement::isCollapsedFace
(
const pointField& points,
const pointField& neiCc,
const scalar minFaceArea,
const scalar maxNonOrtho,
const label faceI
) const
{
vector s = mesh_.faces()[faceI].normal(points);
scalar magS = mag(s);
// Check face area
if (magS < minFaceArea)
{
return true;
}
// Check orthogonality
const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
vector d = ownCc - mesh_.cellCentres()[nei];
scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
if (dDotS < maxNonOrtho)
{
return true;
}
else
{
return false;
}
}
else
{
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
if (mesh_.boundaryMesh()[patchI].coupled())
{
vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
if (dDotS < maxNonOrtho)
{
return true;
}
else
{
return false;
}
}
else
{
// Collapsing normal boundary face does not cause problems with
// non-orthogonality
return true;
}
}
}
// Check if moving cell to new points causes it to collapse.
bool Foam::meshRefinement::isCollapsedCell
(
const pointField& points,
const scalar volFraction,
const label cellI
) const
{
scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
if (vol/mesh_.cellVolumes()[cellI] < volFraction)
{
return true;
}
else
{
return false;
}
}
// Returns list with for every internal face -1 or the patch they should
// be baffled into. Gets run after createBaffles so all the surface
// intersections have already been turned into baffles. Used to remove cells
// by baffling all their faces and have the splitMeshRegions chuck away non
// used regions.
// be baffled into. Gets run after createBaffles so all the unzoned surface
// intersections have already been turned into baffles. (Note: zoned surfaces
// are not baffled at this stage)
// Used to remove cells by baffling all their faces and have the
// splitMeshRegions chuck away non used regions.
Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
(
const dictionary& motionDict,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const labelList& globalToPatch
@ -294,6 +386,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
const labelList& pointLevel = meshCutter_.pointLevel();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Swap neighbouring cell centres and cell level
labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
// Per internal face (boundary faces not used) the patch that the
// baffle should get (or -1)
@ -334,6 +430,9 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
// Count of faces marked for baffling
label nBaffleFaces = 0;
// Count of faces not baffled since would not cause a collapse
label nPrevented = 0;
if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
{
Info<< "markFacesOnProblemCells :"
@ -420,6 +519,74 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
);
// See if checking for collapse
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Collapse checking parameters
scalar volFraction = -1;
if (motionDict.found("minVolCollapseRatio"))
{
volFraction = readScalar(motionDict.lookup("minVolCollapseRatio"));
}
const bool checkCollapse = (volFraction > 0);
scalar minArea = -1;
scalar maxNonOrtho = -1;
// Find nearest (non-baffle) surface
pointField newPoints;
if (checkCollapse)
{
minArea = readScalar(motionDict.lookup("minArea"));
maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
Info<< "markFacesOnProblemCells :"
<< " Deleting all-anchor surface cells only if"
<< "snapping them violates mesh quality constraints:" << nl
<< " snapped/original cell volume < " << volFraction << nl
<< " face area < " << minArea << nl
<< " non-orthogonality > " << maxNonOrtho << nl
<< endl;
// Construct addressing engine.
autoPtr<indirectPrimitivePatch> ppPtr
(
meshRefinement::makePatch
(
mesh_,
adaptPatchIDs
)
);
const indirectPrimitivePatch& pp = ppPtr();
const pointField& localPoints = pp.localPoints();
const labelList& meshPoints = pp.meshPoints();
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces_.findNearest
(
surfaces_.getUnnamedSurfaces(),
localPoints,
scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
hitSurface,
hitInfo
);
// Start of from current points
newPoints = mesh_.points();
forAll(hitInfo, i)
{
if (hitInfo[i].hit())
{
newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
}
}
}
// For each cell count the number of anchor points that are on
// the boundary:
// 8 : check the number of (baffle) boundary faces. If 3 or more block
@ -430,6 +597,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
// 3 or more of these cells. (note that on a flat surface a boundary
// point will only have 4 cells connected to it)
// Does cell have exactly 7 of its 8 anchor points on the boundary?
PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
// If so what is the remaining non-boundary anchor point?
@ -494,23 +662,45 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
// Eugene: delete cell no matter what.
//if (nBfaces > 1)
{
forAll(cFaces, cf)
if
(
checkCollapse
&& !isCollapsedCell(newPoints, volFraction, cellI)
)
{
label faceI = cFaces[cf];
if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
nPrevented++;
//Pout<< "Preventing collapse of 8 anchor point cell "
// << cellI << " at " << mesh_.cellCentres()[cellI]
// << " since new volume "
// << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
// << " old volume " << mesh_.cellVolumes()[cellI]
// << endl;
}
else
{
// Block all faces of cell
forAll(cFaces, cf)
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
label faceI = cFaces[cf];
// Mark face as a 'boundary'
markBoundaryFace
if
(
faceI,
isBoundaryFace,
isBoundaryEdge,
isBoundaryPoint
);
facePatch[faceI] == -1
&& mesh_.isInternalFace(faceI)
)
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
// Mark face as a 'boundary'
markBoundaryFace
(
faceI,
isBoundaryFace,
isBoundaryEdge,
isBoundaryPoint
);
}
}
}
}
@ -556,29 +746,52 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
{
const cell& cFaces = mesh_.cells()[cellI];
forAll(cFaces, cf)
if
(
checkCollapse
&& !isCollapsedCell(newPoints, volFraction, cellI)
)
{
label faceI = cFaces[cf];
nPrevented++;
//Pout<< "Preventing collapse of 7 anchor cell "
// << cellI
// << " at " << mesh_.cellCentres()[cellI]
// << " since new volume "
// << mesh_.cells()[cellI].mag
// (newPoints, mesh_.faces())
// << " old volume " << mesh_.cellVolumes()[cellI]
// << endl;
}
else
{
const cell& cFaces = mesh_.cells()[cellI];
if
(
facePatch[faceI] == -1
&& mesh_.isInternalFace(faceI)
)
forAll(cFaces, cf)
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
label faceI = cFaces[cf];
// Mark face as a 'boundary'
markBoundaryFace
if
(
faceI,
isBoundaryFace,
isBoundaryEdge,
isBoundaryPoint
);
facePatch[faceI] == -1
&& mesh_.isInternalFace(faceI)
)
{
facePatch[faceI] = getBafflePatch
(
facePatch,
faceI
);
nBaffleFaces++;
// Mark face as a 'boundary'
markBoundaryFace
(
faceI,
isBoundaryFace,
isBoundaryEdge,
isBoundaryPoint
);
}
}
}
}
@ -635,11 +848,33 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (nFaceBoundaryEdges == fEdges.size())
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
if
(
checkCollapse
&& !isCollapsedFace
(
newPoints,
neiCc,
minArea,
maxNonOrtho,
faceI
)
)
{
nPrevented++;
//Pout<< "Preventing collapse of face " << faceI
// << " with all boundary edges "
// << " at " << mesh_.faceCentres()[faceI]
// << endl;
}
else
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
// Do NOT update boundary data since this would grow blocked
// faces across gaps.
// Do NOT update boundary data since this would grow blocked
// faces across gaps.
}
}
}
}
@ -669,11 +904,34 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (nFaceBoundaryEdges == fEdges.size())
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
if
(
checkCollapse
&& !isCollapsedFace
(
newPoints,
neiCc,
minArea,
maxNonOrtho,
faceI
)
)
{
nPrevented++;
//Pout<< "Preventing collapse of coupled face "
// << faceI
// << " with all boundary edges "
// << " at " << mesh_.faceCentres()[faceI]
// << endl;
}
else
{
facePatch[faceI] = getBafflePatch(facePatch, faceI);
nBaffleFaces++;
// Do NOT update boundary data since this would grow
// blocked faces across gaps.
// Do NOT update boundary data since this would grow
// blocked faces across gaps.
}
}
}
@ -687,266 +945,181 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
<< " additional internal faces to be converted into baffles."
<< endl;
return facePatch;
}
//XXXXXXXXXXXXXX
// Mark faces to be baffled to prevent snapping problems. Does
// test to find nearest surface and checks which faces would get squashed.
Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
(
const dictionary& motionDict,
const labelList& globalToPatch
) const
{
// Get the labels of added patches.
labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
// Construct addressing engine.
autoPtr<indirectPrimitivePatch> ppPtr
(
meshRefinement::makePatch
(
mesh_,
adaptPatchIDs
)
);
const indirectPrimitivePatch& pp = ppPtr();
const pointField& localPoints = pp.localPoints();
const labelList& meshPoints = pp.meshPoints();
// Find nearest (non-baffle) surface
pointField newPoints(mesh_.points());
if (checkCollapse)
{
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces_.findNearest
(
surfaces_.getUnnamedSurfaces(),
localPoints,
scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
hitSurface,
hitInfo
);
forAll(hitInfo, i)
{
if (hitInfo[i].hit())
{
//label pointI = meshPoints[i];
//Pout<< " " << pointI << " moved from "
// << mesh_.points()[pointI] << " by "
// << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
// << endl;
newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
}
}
Info<< "markFacesOnProblemCells : prevented "
<< returnReduce(nPrevented, sumOp<label>())
<< " internal faces fom getting converted into baffles."
<< endl;
}
// Per face (internal or coupled!) the patch that the
// baffle should get (or -1).
labelList facePatch(mesh_.nFaces(), -1);
// Count of baffled faces
label nBaffleFaces = 0;
return facePatch;
}
// // Sync position? Or not since same face on both side so just sync
// // result of baffle.
//// Mark faces to be baffled to prevent snapping problems. Does
//// test to find nearest surface and checks which faces would get squashed.
//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
//(
// const dictionary& motionDict,
// const labelList& globalToPatch
//) const
//{
// // Get the labels of added patches.
// labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
//
// const scalar minArea(readScalar(motionDict.lookup("minArea")));
// // Construct addressing engine.
// autoPtr<indirectPrimitivePatch> ppPtr
// (
// meshRefinement::makePatch
// (
// mesh_,
// adaptPatchIDs
// )
// );
// const indirectPrimitivePatch& pp = ppPtr();
// const pointField& localPoints = pp.localPoints();
// const labelList& meshPoints = pp.meshPoints();
//
// Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
// << minArea << endl;
//
// pointField facePoints;
// for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
// // Find nearest (non-baffle) surface
// pointField newPoints(mesh_.points());
// {
// const face& f = mesh_.faces()[faceI];
// List<pointIndexHit> hitInfo;
// labelList hitSurface;
// surfaces_.findNearest
// (
// surfaces_.getUnnamedSurfaces(),
// localPoints,
// scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
// hitSurface,
// hitInfo
// );
//
// bool usesPatchPoint = false;
//
// facePoints.setSize(f.size());
// forAll(f, fp)
// forAll(hitInfo, i)
// {
// Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
//
// if (iter != pp.meshPointMap().end())
// if (hitInfo[i].hit())
// {
// facePoints[fp] = newPosition[iter()];
// usesPatchPoint = true;
// //label pointI = meshPoints[i];
// //Pout<< " " << pointI << " moved from "
// // << mesh_.points()[pointI] << " by "
// // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
// // << endl;
// newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
// }
// else
// }
// }
//
// // Per face (internal or coupled!) the patch that the
// // baffle should get (or -1).
// labelList facePatch(mesh_.nFaces(), -1);
// // Count of baffled faces
// label nBaffleFaces = 0;
//
// {
// pointField oldPoints(mesh_.points());
// mesh_.movePoints(newPoints);
// faceSet wrongFaces(mesh_, "wrongFaces", 100);
// {
// //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
//
// // Just check the errors from squashing
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// const labelList allFaces(identity(mesh_.nFaces()));
// label nWrongFaces = 0;
//
// scalar minArea(readScalar(motionDict.lookup("minArea")));
// if (minArea > -SMALL)
// {
// facePoints[fp] = mesh_.points()[f[fp]];
// polyMeshGeometry::checkFaceArea
// (
// false,
// minArea,
// mesh_,
// mesh_.faceAreas(),
// allFaces,
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce
// (
// wrongFaces.size(),
// sumOp<label>()
// );
//
// Info<< " faces with area < "
// << setw(5) << minArea
// << " m^2 : "
// << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
// }
//
//// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
// scalar minDet = 0.01;
// if (minDet > -1)
// {
// polyMeshGeometry::checkCellDeterminant
// (
// false,
// minDet,
// mesh_,
// mesh_.faceAreas(),
// allFaces,
// polyMeshGeometry::affectedCells(mesh_, allFaces),
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce
// (
// wrongFaces.size(),
// sumOp<label>()
// );
//
// Info<< " faces on cells with determinant < "
// << setw(5) << minDet << " : "
// << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
// }
// }
//
// if (usesPatchPoint)
// {
// // Check area of face wrt original area
// face identFace(identity(f.size()));
//
// if (identFace.mag(facePoints) < minArea)
// forAllConstIter(faceSet, wrongFaces, iter)
// {
// label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
//
// if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
// {
// facePatch[faceI] = getBafflePatch(facePatch, faceI);
// facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
// nBaffleFaces++;
//
// //Pout<< " " << iter.key()
// // //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
// // << " is destined for patch " << facePatch[iter.key()]
// // << endl;
// }
// }
// // Restore points.
// mesh_.movePoints(oldPoints);
// }
//
//
// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
// Info<< "markFacesOnProblemCellsGeometric : marked "
// << returnReduce(nBaffleFaces, sumOp<label>())
// << " additional internal and coupled faces"
// << " to be converted into baffles." << endl;
//
// if (pp.coupled())
// {
// forAll(pp, i)
// {
// label faceI = pp.start()+i;
// syncTools::syncFaceList
// (
// mesh_,
// facePatch,
// maxEqOp<label>(),
// false // no separation
// );
//
// const face& f = mesh_.faces()[faceI];
//
// bool usesPatchPoint = false;
//
// facePoints.setSize(f.size());
// forAll(f, fp)
// {
// Map<label>::const_iterator iter =
// pp.meshPointMap().find(f[fp]);
//
// if (iter != pp.meshPointMap().end())
// {
// facePoints[fp] = newPosition[iter()];
// usesPatchPoint = true;
// }
// else
// {
// facePoints[fp] = mesh_.points()[f[fp]];
// }
// }
//
// if (usesPatchPoint)
// {
// // Check area of face wrt original area
// face identFace(identity(f.size()));
//
// if (identFace.mag(facePoints) < minArea)
// {
// facePatch[faceI] = getBafflePatch(facePatch, faceI);
// nBaffleFaces++;
// }
// }
// }
// }
// }
{
pointField oldPoints(mesh_.points());
mesh_.movePoints(newPoints);
faceSet wrongFaces(mesh_, "wrongFaces", 100);
{
//motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
// Just check the errors from squashing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const labelList allFaces(identity(mesh_.nFaces()));
label nWrongFaces = 0;
scalar minArea(readScalar(motionDict.lookup("minArea")));
if (minArea > -SMALL)
{
polyMeshGeometry::checkFaceArea
(
false,
minArea,
mesh_,
mesh_.faceAreas(),
allFaces,
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces with area < "
<< setw(5) << minArea
<< " m^2 : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
scalar minDet = 0.01;
if (minDet > -1)
{
polyMeshGeometry::checkCellDeterminant
(
false,
minDet,
mesh_,
mesh_.faceAreas(),
allFaces,
polyMeshGeometry::affectedCells(mesh_, allFaces),
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces on cells with determinant < "
<< setw(5) << minDet << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
}
forAllConstIter(faceSet, wrongFaces, iter)
{
label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
{
facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
nBaffleFaces++;
//Pout<< " " << iter.key()
// //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
// << " is destined for patch " << facePatch[iter.key()]
// << endl;
}
}
// Restore points.
mesh_.movePoints(oldPoints);
}
Info<< "markFacesOnProblemCellsGeometric : marked "
<< returnReduce(nBaffleFaces, sumOp<label>())
<< " additional internal and coupled faces"
<< " to be converted into baffles." << endl;
syncTools::syncFaceList
(
mesh_,
facePatch,
maxEqOp<label>(),
false // no separation
);
return facePatch;
}
//XXXXXXXX
// return facePatch;
//}
// ************************************************************************* //

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

@ -154,7 +154,7 @@ public:
return cellZoneNames_;
}
//- Get indices of named surfaces (surfaces with faceZoneName)
//- Get indices of unnamed surfaces (surfaces without faceZoneName)
labelList getUnnamedSurfaces() const;
//- Get indices of named surfaces (surfaces with faceZoneName)

19
src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
View File

@ -165,6 +165,11 @@ void Foam::hierarchGeomDecomp::findBinary
// (one beyond) index of highValue
label high = values.size();
// Safeguards to avoid infinite loop.
label lowPrev = -1;
label midPrev = -1;
label highPrev = -1;
//while (low <= high)
while (true)
{
@ -197,6 +202,20 @@ void Foam::hierarchGeomDecomp::findBinary
// Update mid, midValue
midValue = 0.5*(lowValue+highValue);
mid = findLower(values, midValue, low, high);
// Safeguard if same as previous.
bool hasNotChanged = (mid == midPrev) && (low == lowPrev) && (high == highPrev);
if (returnReduce(hasNotChanged, andOp<bool>()))
{
WarningIn("hierarchGeomDecomp::findBinary(..)")
<< "unable to find desired deomposition split, making do!"
<< endl;
break;
}
midPrev = mid;
lowPrev = low;
highPrev = high;
}
}

374
src/finiteVolume/cfdTools/general/MRF/MRFZone.C
View File

@ -32,53 +32,96 @@ License
#include "syncTools.H"
#include "faceSet.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::MRFZone, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::MRFZone::setMRFFaces
(
labelList& faceType,
const labelList& excludedPatchIDs
)
void Foam::MRFZone::setMRFFaces()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Knock out coupled patches
forAll(patches, patchi)
// Type per face:
// 0:not in zone
// 1:moving with frame
// 2:other
labelList faceType(mesh_.nFaces(), 0);
// Determine faces in cell zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// (without constructing cells)
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
// Cells in zone
boolList zoneCell(mesh_.nCells(), false);
if (cellZoneID_ != -1)
{
if (patches[patchi].coupled())
const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
forAll(cellLabels, i)
{
const polyPatch& pp = patches[patchi];
zoneCell[cellLabels[i]] = true;
}
}
forAll(pp, j)
label nZoneFaces = 0;
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
{
faceType[faceI] = 1;
nZoneFaces++;
}
}
labelHashSet excludedPatches(excludedPatchLabels_);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled() || excludedPatches.found(patchI))
{
forAll(pp, i)
{
label faceI = pp.start()+j;
label faceI = pp.start()+i;
if (faceType[faceI] == 1)
if (zoneCell[own[faceI]])
{
faceType[faceI] = 2;
nZoneFaces++;
}
}
}
else if (!isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label faceI = pp.start()+i;
if (zoneCell[own[faceI]])
{
faceType[faceI] = 1;
nZoneFaces++;
}
}
}
}
// All explicitly provided exclusions
forAll(excludedPatchIDs, i)
{
const polyPatch& pp = patches[excludedPatchIDs[i]];
// Now we have for faceType:
// 0 : face not in cellZone
// 1 : internal face or normal patch face
// 2 : coupled patch face or excluded patch face
forAll(pp, j)
{
label faceI = pp.start()+j;
// Sort into lists per patch.
if (faceType[faceI] == 1)
{
faceType[faceI] = 2;
}
}
}
// Collect into lists per patch.
internalFaces_.setSize(mesh_.nFaces());
label nInternal = 0;
@ -142,149 +185,43 @@ void Foam::MRFZone::setMRFFaces
}
}
//if (debug)
//{
// faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
// Pout<< "Writing internal faces in MRF zone to faceSet "
// << internalFaces.name() << endl;
// internalFaces.write();
//}
//{
// faceSet MRFFaces(mesh_, "includedFaces", 100);
// forAll(includedFaces_, patchi)
// {
// forAll(includedFaces_[patchi], i)
// {
// label patchFacei = includedFaces_[patchi][i];
// MRFFaces.insert(patches[patchi].start()+patchFacei);
// }
// }
// Pout<< "Writing patch faces in MRF zone to faceSet "
// << MRFFaces.name() << endl;
// MRFFaces.write();
//}
//{
// faceSet excludedFaces(mesh_, "excludedFaces", 100);
// forAll(excludedFaces_, patchi)
// {
// forAll(excludedFaces_[patchi], i)
// {
// label patchFacei = excludedFaces_[patchi][i];
// excludedFaces.insert(patches[patchi].start()+patchFacei);
// }
// }
// Pout<< "Writing faces in MRF zone with special handling to faceSet "
// << excludedFaces.name() << endl;
// excludedFaces.write();
//}
}
void Foam::MRFZone::setMRFFaces()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Type per face:
// 0:not in zone
// 1:moving with frame
// 2:other
labelList faceType(mesh_.nFaces(), 0);
bool faceZoneFound = (faceZoneID_ != -1);
reduce(faceZoneFound, orOp<bool>());
if (faceZoneFound)
if (debug)
{
// Explicitly provided faces.
if (faceZoneID_ != -1)
faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
Pout<< "Writing " << internalFaces.size()
<< " internal faces in MRF zone to faceSet "
<< internalFaces.name() << endl;
internalFaces.write();
faceSet MRFFaces(mesh_, "includedFaces", 100);
forAll(includedFaces_, patchi)
{
const labelList& zoneFaces = mesh_.faceZones()[faceZoneID_];
forAll(zoneFaces, i)
forAll(includedFaces_[patchi], i)
{
faceType[zoneFaces[i]] = 1;
}
if (allPatchesMove_)
{
// Explicitly provided patches that do not move
setMRFFaces(faceType, patchLabels_);
}
else
{
setMRFFaces(faceType, labelList(0));
label patchFacei = includedFaces_[patchi][i];
MRFFaces.insert(patches[patchi].start()+patchFacei);
}
}
Pout<< "Writing " << MRFFaces.size()
<< " patch faces in MRF zone to faceSet "
<< MRFFaces.name() << endl;
MRFFaces.write();
faceSet excludedFaces(mesh_, "excludedFaces", 100);
forAll(excludedFaces_, patchi)
{
forAll(excludedFaces_[patchi], i)
{
label patchFacei = excludedFaces_[patchi][i];
excludedFaces.insert(patches[patchi].start()+patchFacei);
}
}
Pout<< "Writing " << excludedFaces.size()
<< " faces in MRF zone with special handling to faceSet "
<< excludedFaces.name() << endl;
excludedFaces.write();
}
else
{
// Determine faces in cell zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// (without constructing cells)
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
// Cells in zone
boolList zoneCell(mesh_.nCells(), false);
if (cellZoneID_ != -1)
{
const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
forAll(cellLabels, i)
{
zoneCell[cellLabels[i]] = true;
}
}
label nZoneFaces = 0;
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
{
faceType[faceI] = 1;
nZoneFaces++;
}
}
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (!isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label faceI = pp.start()+i;
if (zoneCell[own[faceI]])
{
faceType[faceI] = 1;
nZoneFaces++;
}
}
}
}
syncTools::syncFaceList(mesh_, faceType, maxEqOp<label>(), false);
Info<< nl
<< "MRFZone " << name_ << " : did not find a faceZone; using "
<< returnReduce(nZoneFaces, sumOp<label>())
<< " faces from the cellZone instead." << endl;
if (allPatchesMove_)
{
// Explicitly provided excluded patches
setMRFFaces(faceType, patchLabels_);
}
else
{
setMRFFaces(faceType, labelList(0));
}
}
}
@ -296,36 +233,41 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
name_(is),
dict_(is),
cellZoneID_(mesh_.cellZones().findZoneID(name_)),
faceZoneID_(mesh_.faceZones().findZoneID(name_)),
allPatchesMove_(dict_.found("nonRotatingPatches")),
patchNames_
excludedPatchNames_
(
allPatchesMove_
? dict_.lookup("nonRotatingPatches")
: dict_.lookup("patches")
dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
),
origin_(dict_.lookup("origin")),
axis_(dict_.lookup("axis")),
omega_(dict_.lookup("omega")),
Omega_("Omega", omega_*axis_)
{
if (dict_.found("patches"))
{
WarningIn("MRFZone(const fvMesh&, Istream&)")
<< "Ignoring entry 'patches'\n"
<< " By default all patches within the rotating region rotate.\n"
<< " Optionally supply excluded patches using 'nonRotatingPatches'."
<< endl;
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
axis_ = axis_/mag(axis_);
Omega_ = omega_*axis_;
patchLabels_.setSize(patchNames_.size());
excludedPatchLabels_.setSize(excludedPatchNames_.size());
forAll(patchNames_, i)
forAll(excludedPatchNames_, i)
{
patchLabels_[i] = patches.findPatchID(patchNames_[i]);
excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
if (patchLabels_[i] == -1)
if (excludedPatchLabels_[i] == -1)
{
FatalErrorIn
(
"Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
) << "cannot find MRF patch " << patchNames_[i]
) << "cannot find MRF patch " << excludedPatchNames_[i]
<< exit(FatalError);
}
}
@ -364,7 +306,46 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
forAll(cells, i)
{
Usource[cells[i]] -= V[cells[i]]*(Omega ^ U[cells[i]]);
label celli = cells[i];
Usource[celli] -= V[celli]*(Omega ^ U[celli]);
}
}
void Foam::MRFZone::relativeVelocity(volVectorField& U) const
{
const volVectorField& C = mesh_.C();
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const labelList& cells = mesh_.cellZones()[cellZoneID_];
forAll(cells, i)
{
label celli = cells[i];
U[celli] -= (Omega ^ (C[celli] - origin));
}
// Included patches
forAll(includedFaces_, patchi)
{
forAll(includedFaces_[patchi], i)
{
label patchFacei = includedFaces_[patchi][i];
U.boundaryField()[patchi][patchFacei] = vector::zero;
}
}
// Excluded patches
forAll(excludedFaces_, patchi)
{
forAll(excludedFaces_[patchi], i)
{
label patchFacei = excludedFaces_[patchi][i];
U.boundaryField()[patchi][patchFacei] -=
(Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
}
}
}
@ -410,6 +391,49 @@ void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
}
void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
{
const surfaceVectorField& Cf = mesh_.Cf();
const surfaceVectorField& Sf = mesh_.Sf();
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
// Internal faces
forAll(internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
}
// Included patches
forAll(includedFaces_, patchi)
{
forAll(includedFaces_[patchi], i)
{
label patchFacei = includedFaces_[patchi][i];
phi.boundaryField()[patchi][patchFacei] +=
(Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
& Sf.boundaryField()[patchi][patchFacei];
}
}
// Excluded patches
forAll(excludedFaces_, patchi)
{
forAll(excludedFaces_[patchi], i)
{
label patchFacei = excludedFaces_[patchi][i];
phi.boundaryField()[patchi][patchFacei] +=
(Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
& Sf.boundaryField()[patchi][patchFacei];
}
}
}
void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
{
const vector& origin = origin_.value();

31
src/finiteVolume/cfdTools/general/MRF/MRFZone.H
View File

@ -26,7 +26,7 @@ Class
Foam::MRFZone
Description
MRF zone definition based on cell zone and optional face zone and parameters
MRF zone definition based on cell zone and parameters
obtained from a control dictionary constructed from the given stream.
The rotation of the MRF region is defined by an origin and axis of
@ -74,15 +74,8 @@ class MRFZone
label cellZoneID_;
//- label of face zone with faces on outside of cell zone.
label faceZoneID_;
//- Do patches move with frame (true) or are explicitly provided (false,
// old behaviour)
Switch allPatchesMove_;
const wordList patchNames_;
labelList patchLabels_;
const wordList excludedPatchNames_;
labelList excludedPatchLabels_;
//- Internal faces that are part of MRF
labelList internalFaces_;
@ -101,17 +94,9 @@ class MRFZone
// Private Member Functions
//- Divide faces in frame according to patch
void setMRFFaces
(
labelList& faceType,
const labelList& excludedPatchIDs
);
//- Divide faces in frame according to patch
void setMRFFaces();
//- Disallow default bitwise copy construct
MRFZone(const MRFZone&);
@ -121,6 +106,10 @@ class MRFZone
public:
// Declare name of the class and its debug switch
ClassName("MRFZone");
// Constructors
//- Construct from fvMesh and Istream
@ -165,9 +154,15 @@ public:
//- Add the Coriolis force contribution to the momentum equation
void addCoriolis(fvVectorMatrix& UEqn) const;
//- Make the given absolute velocity relative within the MRF region
void relativeVelocity(volVectorField& U) const;
//- Make the given absolute flux relative within the MRF region
void relativeFlux(surfaceScalarField& phi) const;
//- Make the given relative flux absolute within the MRF region
void absoluteFlux(surfaceScalarField& phi) const;
//- Correct the boundary velocity for the roation of the MRF region
void correctBoundaryVelocity(volVectorField& U) const;

18
src/finiteVolume/cfdTools/general/MRF/MRFZones.C
View File

@ -65,6 +65,15 @@ void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const
}
void Foam::MRFZones::relativeVelocity(volVectorField& U) const
{
forAll(*this, i)
{
operator[](i).relativeVelocity(U);
}
}
void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
{
forAll(*this, i)
@ -74,6 +83,15 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
}
void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
{
forAll(*this, i)
{
operator[](i).absoluteFlux(phi);
}
}
void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
{
forAll(*this, i)

6
src/finiteVolume/cfdTools/general/MRF/MRFZones.H
View File

@ -76,9 +76,15 @@ public:
//- Add the Coriolis force contribution to the momentum equation
void addCoriolis(fvVectorMatrix& UEqn) const;
//- Make the given absolute velocity relative within the MRF region
void relativeVelocity(volVectorField& U) const;
//- Make the given absolute flux relative within the MRF region
void relativeFlux(surfaceScalarField& phi) const;
//- Make the given relative flux absolute within the MRF region
void absoluteFlux(surfaceScalarField& phi) const;
//- Correct the boundary velocity for the roation of the MRF region
void correctBoundaryVelocity(volVectorField& U) const;
};

4
src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
View File

@ -22,10 +22,6 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Formulation based on relative velocities
\*---------------------------------------------------------------------------*/
#include "SRFModel.H"

6
src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
View File

@ -27,7 +27,7 @@ License
#include "adjustPhi.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "processorFvPatchFields.H"
#include "processorFvsPatchFields.H"
#include "inletOutletFvPatchFields.H"
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
@ -52,7 +52,7 @@ bool Foam::adjustPhi
const fvPatchVectorField& Up = U.boundaryField()[patchi];
const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
if (!isType<processorFvPatchScalarField>(phip))
if (!isType<processorFvsPatchScalarField>(phip))
{
if
(
@ -128,7 +128,7 @@ bool Foam::adjustPhi
const fvPatchVectorField& Up = U.boundaryField()[patchi];
fvsPatchScalarField& phip = phi.boundaryField()[patchi];
if (!isType<processorFvPatchScalarField>(phip))
if (!isType<processorFvsPatchScalarField>(phip))
{
if
(

8
src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
View File

@ -56,13 +56,7 @@ slicedFvPatchField<Type>::slicedFvPatchField
)
:
fvPatchField<Type>(p, iF)
{
notImplemented
(
"slicedFvPatchField<Type>::"
"slicedFvPatchField(const fvPatch&, const Field<Type>&)"
);
}
{}
template<class Type>

2
src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H
View File

@ -75,7 +75,7 @@ public:
const Field<Type>&
);
//- Construct from patch and internal field
//- Construct from patch and internal field. Assign value later.
slicedFvPatchField
(
const fvPatch&,

8
src/finiteVolume/fields/fvsPatchFields/basic/sliced/slicedFvsPatchField.C
View File

@ -56,13 +56,7 @@ slicedFvsPatchField<Type>::slicedFvsPatchField
)
:
fvsPatchField<Type>(p, iF)
{
notImplemented
(
"slicedFvsPatchField<Type>::"
"slicedFvsPatchField(const fvPatch&, const Field<Type>&)"
);
}
{}
template<class Type>

8
src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
View File

@ -62,13 +62,7 @@ Foam::distanceSurface::interpolateField
);
// Sample.
return surface().interpolate
(
cellDistancePtr_(),
pointDistance_,
volFld,
pointFld()
);
return surface().interpolate(volFld, pointFld());
}

560
src/sampling/sampledSurface/isoSurface/isoSurface.C
View File

@ -30,6 +30,7 @@ License
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
#include "slicedVolFields.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,6 +42,243 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::isoSurface::noTransform(const tensor& tt) const
{
return
(mag(tt.xx()-1) < mergeDistance_)
&& (mag(tt.yy()-1) < mergeDistance_)
&& (mag(tt.zz()-1) < mergeDistance_)
&& (mag(tt.xy()) < mergeDistance_)
&& (mag(tt.xz()) < mergeDistance_)
&& (mag(tt.yx()) < mergeDistance_)
&& (mag(tt.yz()) < mergeDistance_)
&& (mag(tt.zx()) < mergeDistance_)
&& (mag(tt.zy()) < mergeDistance_);
}
bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
return cpp.parallel() && !cpp.separated();
}
// Calculates per face whether couple is collocated.
Foam::PackedBoolList Foam::isoSurface::collocatedFaces
(
const coupledPolyPatch& pp
) const
{
// Initialise to false
PackedBoolList collocated(pp.size());
const vectorField& separation = pp.separation();
const tensorField& forwardT = pp.forwardT();
if (forwardT.size() == 0)
{
// Parallel.
if (separation.size() == 0)
{
collocated = 1u;
}
else if (separation.size() == 1)
{
// Fully separate. Do not synchronise.
}
else
{
// Per face separation.
forAll(pp, faceI)
{
if (mag(separation[faceI]) < mergeDistance_)
{
collocated[faceI] = 1u;
}
}
}
}
else if (forwardT.size() == 1)
{
// Fully transformed.
}
else
{
// Per face transformation.
forAll(pp, faceI)
{
if (noTransform(forwardT[faceI]))
{
collocated[faceI] = 1u;
}
}
}
return collocated;
}
// Insert the data for local point patchPointI into patch local values
// and/or into the shared values field.
void Foam::isoSurface::insertPointData
(
const processorPolyPatch& pp,
const Map<label>& meshToShared,
const pointField& pointValues,
const label patchPointI,
pointField& patchValues,
pointField& sharedValues
) const
{
label meshPointI = pp.meshPoints()[patchPointI];
// Store in local field
label nbrPointI = pp.neighbPoints()[patchPointI];
if (nbrPointI >= 0 && nbrPointI < patchValues.size())
{
minEqOp<point>()(patchValues[nbrPointI], pointValues[meshPointI]);
}
// Store in shared field
Map<label>::const_iterator iter = meshToShared.find(meshPointI);
if (iter != meshToShared.end())
{
minEqOp<point>()(sharedValues[iter()], pointValues[meshPointI]);
}
}
void Foam::isoSurface::syncUnseparatedPoints
(
pointField& pointValues,
const point& nullValue
) const
{
// Until syncPointList handles separated coupled patches with multiple
// transforms do our own synchronisation of non-separated patches only
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const globalMeshData& pd = mesh_.globalData();
const labelList& sharedPtAddr = pd.sharedPointAddr();
const labelList& sharedPtLabels = pd.sharedPointLabels();
// Create map from meshPoint to globalShared index.
Map<label> meshToShared(2*sharedPtLabels.size());
forAll(sharedPtLabels, i)
{
meshToShared.insert(sharedPtLabels[i], sharedPtAddr[i]);
}
// Values on shared points.
pointField sharedInfo(pd.nGlobalPoints(), nullValue);
if (Pstream::parRun())
{
// Send
forAll(patches, patchI)
{
if
(
isA<processorPolyPatch>(patches[patchI])
&& patches[patchI].nPoints() > 0
)
{
const processorPolyPatch& pp =
refCast<const processorPolyPatch>(patches[patchI]);
const labelList& meshPts = pp.meshPoints();
pointField patchInfo(meshPts.size(), nullValue);
PackedBoolList isCollocated(collocatedFaces(pp));
forAll(isCollocated, faceI)
{
if (isCollocated[faceI])
{
const face& f = pp.localFaces()[faceI];
forAll(f, fp)
{
label pointI = f[fp];
insertPointData
(
pp,
meshToShared,
pointValues,
pointI,
patchInfo,
sharedInfo
);
}
}
}
OPstream toNbr(Pstream::blocking, pp.neighbProcNo());
toNbr << patchInfo;
}
}
// Receive and combine.
forAll(patches, patchI)
{
if
(
isA<processorPolyPatch>(patches[patchI])
&& patches[patchI].nPoints() > 0
)
{
const processorPolyPatch& pp =
refCast<const processorPolyPatch>(patches[patchI]);
pointField nbrPatchInfo(pp.nPoints());
{
// We do not know the number of points on the other side
// so cannot use Pstream::read.
IPstream fromNbr(Pstream::blocking, pp.neighbProcNo());
fromNbr >> nbrPatchInfo;
}
// Null any value which is not on neighbouring processor
nbrPatchInfo.setSize(pp.nPoints(), nullValue);
const labelList& meshPts = pp.meshPoints();
forAll(meshPts, pointI)
{
label meshPointI = meshPts[pointI];
minEqOp<point>()
(
pointValues[meshPointI],
nbrPatchInfo[pointI]
);
}
}
}
}
// Don't do cyclics for now. Are almost always separated anyway.
// Shared points
// Combine on master.
Pstream::listCombineGather(sharedInfo, minEqOp<point>());
Pstream::listCombineScatter(sharedInfo);
// Now we will all have the same information. Merge it back with
// my local information. (Note assignment and not combine operator)
forAll(sharedPtLabels, i)
{
label meshPointI = sharedPtLabels[i];
pointValues[meshPointI] = sharedInfo[sharedPtAddr[i]];
}
}
Foam::scalar Foam::isoSurface::isoFraction
(
const scalar s0,
@ -89,6 +327,7 @@ bool Foam::isoSurface::isEdgeOfFaceCut
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const label cellI,
const label faceI,
@ -98,13 +337,12 @@ void Foam::isoSurface::getNeighbour
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const surfaceScalarField& weights = mesh_.weights();
if (mesh_.isInternalFace(faceI))
{
label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
nbrValue = cVals[nbr];
nbrPoint = mesh_.cellCentres()[nbr];
nbrPoint = meshC[nbr];
}
else
{
@ -113,38 +351,8 @@ void Foam::isoSurface::getNeighbour
const polyPatch& pp = mesh_.boundaryMesh()[patchI];
label patchFaceI = faceI-pp.start();
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient
nbrValue = cVals[own[faceI]];
nbrPoint = mesh_.faceCentres()[faceI];
}
else if (pp.coupled())
{
if (!refCast<const coupledPolyPatch>(pp).separated())
{
// Behave as internal face:
// other side value
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
// other side cell centre
nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
}
else
{
// Do some interpolation for now
const scalarField& w = weights.boundaryField()[patchI];
nbrPoint = mesh_.faceCentres()[faceI];
nbrValue =
(1-w[patchFaceI])*cVals[own[faceI]]
+ w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
}
}
else
{
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = mesh_.faceCentres()[faceI];
}
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
}
}
@ -153,6 +361,7 @@ void Foam::isoSurface::getNeighbour
void Foam::isoSurface::calcCutTypes
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals
)
@ -174,6 +383,7 @@ void Foam::isoSurface::calcCutTypes
getNeighbour
(
boundaryRegion,
meshC,
cVals,
own[faceI],
faceI,
@ -215,6 +425,7 @@ void Foam::isoSurface::calcCutTypes
getNeighbour
(
boundaryRegion,
meshC,
cVals,
own[faceI],
faceI,
@ -408,6 +619,7 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
void Foam::isoSurface::calcSnappedCc
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals,
@ -448,6 +660,7 @@ void Foam::isoSurface::calcSnappedCc
getNeighbour
(
boundaryRegion,
meshC,
cVals,
cellI,
faceI,
@ -574,6 +787,7 @@ void Foam::isoSurface::calcSnappedPoint
(
const PackedBoolList& isBoundaryPoint,
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals,
@ -639,6 +853,7 @@ void Foam::isoSurface::calcSnappedPoint
getNeighbour
(
boundaryRegion,
meshC,
cVals,
own,
faceI,
@ -747,14 +962,22 @@ void Foam::isoSurface::calcSnappedPoint
}
}
syncTools::syncPointList
(
mesh_,
collapsedPoint,
minEqOp<point>(),
greatPoint,
true // are coordinates so separate
);
//Pout<< "**hack" << endl;
//pointField oldCollaped(collapsedPoint);
syncUnseparatedPoints(collapsedPoint, greatPoint);
//forAll(collapsedPoint, pointI)
//{
// if (collapsedPoint[pointI] != oldCollaped[pointI])
// {
// Pout<< "**Synced point " << pointI
// << " coord:" << mesh_.points()[pointI]
// << " from " << oldCollaped[pointI]
// << " to " << collapsedPoint[pointI]
// << endl;
// }
//}
snappedPoint.setSize(mesh_.nPoints());
snappedPoint = -1;
@ -1534,6 +1757,7 @@ Foam::isoSurface::isoSurface
)
:
mesh_(cVals.mesh()),
pVals_(pVals),
iso_(iso),
mergeDistance_(mergeTol*mesh_.bounds().mag())
{
@ -1545,8 +1769,8 @@ Foam::isoSurface::isoSurface
<< min(cVals.internalField()) << " / "
<< max(cVals.internalField()) << nl
<< " point min/max : "
<< min(pVals) << " / "
<< max(pVals) << nl
<< min(pVals_) << " / "
<< max(pVals_) << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
@ -1555,6 +1779,89 @@ Foam::isoSurface::isoSurface
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Rewrite input field
// ~~~~~~~~~~~~~~~~~~~
// Rewrite input volScalarField to have interpolated values
// on separated patches.
cValsPtr_.reset(adaptPatchFields(cVals).ptr());
// Construct cell centres field consistent with cVals
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Generate field to interpolate. This is identical to the mesh.C()
// except on separated coupled patches and on empty patches.
slicedVolVectorField meshC
(
IOobject
(
"C",
mesh_.pointsInstance(),
mesh_.meshSubDir,
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimLength,
mesh_.cellCentres(),
mesh_.faceCentres()
);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
// Adapt separated coupled (proc and cyclic) patches
if (isA<coupledPolyPatch>(pp) && !collocatedPatch(pp))
{
fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
(
meshC.boundaryField()[patchI]
);
PackedBoolList isCollocated
(
collocatedFaces(refCast<const coupledPolyPatch>(pp))
);
forAll(isCollocated, i)
{
if (!isCollocated[i])
{
pfld[i] = mesh_.faceCentres()[pp.start()+i];
}
}
}
else if (isA<emptyPolyPatch>(pp))
{
typedef slicedVolVectorField::GeometricBoundaryField bType;
bType& bfld = const_cast<bType&>(meshC.boundaryField());
// Clear old value. Cannot resize it since is a slice.
bfld.set(patchI, NULL);
// Set new value we can change
bfld.set
(
patchI,
new calculatedFvPatchField<vector>
(
mesh_.boundary()[patchI],
meshC
)
);
// Change to face centres
bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
}
}
// Pre-calculate patch-per-face to avoid whichPatch call.
labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
@ -1572,8 +1879,9 @@ Foam::isoSurface::isoSurface
}
// Determine if any cut through face/cell
calcCutTypes(boundaryRegion, cVals, pVals);
calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
DynamicList<point> snappedPoints(nCutCells_);
@ -1585,8 +1893,9 @@ Foam::isoSurface::isoSurface
calcSnappedCc
(
boundaryRegion,
cVals,
pVals,
meshC,
cValsPtr_(),
pVals_,
snappedPoints,
snappedCc
@ -1609,48 +1918,67 @@ Foam::isoSurface::isoSurface
label nCellSnaps = snappedPoints.size();
// Determine if point is on boundary. Points on boundaries are never
// snapped. Coupled boundaries are handled explicitly so not marked here.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
forAll(patches, patchI)
{
// Mark all boundary points that are not physically coupled (so anything
// but collocated coupled patches)
const polyPatch& pp = patches[patchI];
if
(
!pp.coupled()
|| refCast<const coupledPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
{
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
faceI++;
}
}
}
// Per point -1 or a point inside snappedPoints.
labelList snappedPoint;
if (regularise)
{
// Determine if point is on boundary.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
forAll(patches, patchI)
{
// Mark all boundary points that are not physically coupled
// (so anything but collocated coupled patches)
if (patches[patchI].coupled())
{
if (!collocatedPatch(patches[patchI]))
{
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>
(
patches[patchI]
);
PackedBoolList isCollocated(collocatedFaces(cpp));
forAll(isCollocated, i)
{
if (!isCollocated[i])
{
const face& f = mesh_.faces()[cpp.start()+i];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
}
}
}
}
else
{
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
const face& f = mesh_.faces()[pp.start()+i];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
}
}
}
calcSnappedPoint
(
isBoundaryPoint,
boundaryRegion,
cVals,
pVals,
meshC,
cValsPtr_(),
pVals_,
snappedPoints,
snappedPoint
@ -1669,80 +1997,14 @@ Foam::isoSurface::isoSurface
}
// Generate field to interpolate. This is identical to the mesh.C()
// except on separated coupled patches and on empty patches.
slicedVolVectorField meshC
(
IOobject
(
"C",
mesh_.pointsInstance(),
mesh_.meshSubDir,
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimLength,
mesh_.cellCentres(),
mesh_.faceCentres()
);
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
pp.coupled()
&& refCast<const coupledPolyPatch>(pp).separated()
)
{
fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
(
meshC.boundaryField()[patchI]
);
pfld.operator==
(
pp.patchSlice(mesh_.faceCentres())
);
}
else if (isA<emptyPolyPatch>(pp))
{
typedef slicedVolVectorField::GeometricBoundaryField bType;
bType& bfld = const_cast<bType&>(meshC.boundaryField());
// Clear old value. Cannot resize it since slice.
bfld.set(patchI, NULL);
// Set new value we can change
bfld.set
(
patchI,
new calculatedFvPatchField<vector>
(
mesh_.boundary()[patchI],
meshC
)
);
// Change to face centres
bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
}
}
}
DynamicList<point> triPoints(nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
generateTriPoints
(
cVals,
pVals,
cValsPtr_(),
pVals_,
meshC,
mesh_.points(),

58
src/sampling/sampledSurface/isoSurface/isoSurface.H
View File

@ -56,6 +56,7 @@ SourceFiles
#include "pointIndexHit.H"
#include "PackedBoolList.H"
#include "volFields.H"
#include "slicedVolFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -92,6 +93,11 @@ class isoSurface
//- Reference to mesh
const fvMesh& mesh_;
const scalarField& pVals_;
//- Input volScalarField with separated coupled patches rewritten
autoPtr<slicedVolScalarField> cValsPtr_;
//- Isosurface value
const scalar iso_;
@ -117,6 +123,38 @@ class isoSurface
// Private Member Functions
// Point synchronisation
//- Does tensor differ (to within mergeTolerance) from identity
bool noTransform(const tensor& tt) const;
//- Is patch a collocated (i.e. non-separated) coupled patch?
static bool collocatedPatch(const polyPatch&);
//- Per face whether is collocated
PackedBoolList collocatedFaces(const coupledPolyPatch&) const;
//- Take value at local point patchPointI and assign it to its
// correct place in patchValues (for transferral) and sharedValues
// (for reduction)
void insertPointData
(
const processorPolyPatch& pp,
const Map<label>& meshToShared,
const pointField& pointValues,
const label patchPointI,
pointField& patchValues,
pointField& sharedValues
) const;
//- Synchonise points on all non-separated coupled patches
void syncUnseparatedPoints
(
pointField& collapsedPoint,
const point& nullValue
) const;
//- Get location of iso value as fraction inbetween s0,s1
scalar isoFraction
(
@ -136,6 +174,7 @@ class isoSurface
void getNeighbour
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const label cellI,
const label faceI,
@ -147,6 +186,7 @@ class isoSurface
void calcCutTypes
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals
);
@ -170,6 +210,7 @@ class isoSurface
void calcSnappedCc
(
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals,
DynamicList<point>& snappedPoints,
@ -182,12 +223,23 @@ class isoSurface
(
const PackedBoolList& isBoundaryPoint,
const labelList& boundaryRegion,
const volVectorField& meshC,
const volScalarField& cVals,
const scalarField& pVals,
DynamicList<point>& snappedPoints,
labelList& snappedPoint
) const;
//- Return input field with coupled (and empty) patch values rewritten
template<class Type>
tmp<SlicedGeometricField
<Type, fvPatchField, slicedFvPatchField, volMesh> >
adaptPatchFields
(
const GeometricField<Type, fvPatchField, volMesh>& fld
) const;
//- Generate single point by interpolation or snapping
template<class Type>
Type generatePoint
@ -345,8 +397,8 @@ public:
// Constructors
//- Construct from cell values and point values. Uses boundaryField
// for boundary values. Requires on coupled patchfields to be set
// to the opposite cell value.
// for boundary values. Holds reference to cellIsoVals and
// pointIsoVals.
isoSurface
(
const volScalarField& cellIsoVals,
@ -376,8 +428,6 @@ public:
template <class Type>
tmp<Field<Type> > interpolate
(
const volScalarField& cellIsoVals,
const scalarField& pointIsoVals,
const GeometricField<Type, fvPatchField, volMesh>& cCoords,
const Field<Type>& pCoords
) const;

224
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View File

@ -27,9 +27,114 @@ License
#include "isoSurface.H"
#include "polyMesh.H"
#include "syncTools.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::SlicedGeometricField
<
Type,
Foam::fvPatchField,
Foam::slicedFvPatchField,
Foam::volMesh
> >
Foam::isoSurface::adaptPatchFields
(
const GeometricField<Type, fvPatchField, volMesh>& fld
) const
{
typedef SlicedGeometricField
<
Type,
fvPatchField,
slicedFvPatchField,
volMesh
> FieldType;
tmp<FieldType> tsliceFld
(
new FieldType
(
IOobject
(
fld.name(),
fld.instance(),
fld.db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
fld, // internal field
true // preserveCouples
)
);
FieldType& sliceFld = tsliceFld();
const fvMesh& mesh = fld.mesh();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<emptyPolyPatch>(pp))
{
// Clear old value. Cannot resize it since is a slice.
sliceFld.boundaryField().set(patchI, NULL);
// Set new value we can change
sliceFld.boundaryField().set
(
patchI,
new calculatedFvPatchField<Type>
(
mesh.boundary()[patchI],
sliceFld
)
);
sliceFld.boundaryField()[patchI] ==
mesh.boundary()[patchI].patchInternalField
(
sliceFld
);
}
else if (isA<cyclicPolyPatch>(pp))
{
// Already has interpolate as value
}
else if (isA<processorPolyPatch>(pp) && !collocatedPatch(pp))
{
fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
(
fld.boundaryField()[patchI]
);
const scalarField& w = mesh.weights().boundaryField()[patchI];
tmp<Field<Type> > f =
w*pfld.patchInternalField()
+ (1.0-w)*pfld.patchNeighbourField();
PackedBoolList isCollocated
(
collocatedFaces(refCast<const processorPolyPatch>(pp))
);
forAll(isCollocated, i)
{
if (!isCollocated[i])
{
pfld[i] = f()[i];
}
}
}
}
return tsliceFld;
}
template<class Type>
Type Foam::isoSurface::generatePoint
(
@ -389,7 +494,6 @@ void Foam::isoSurface::generateTriPoints
}
// Generate triangle points
triPoints.clear();
@ -456,16 +560,16 @@ void Foam::isoSurface::generateTriPoints
syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& !refCast<const processorPolyPatch>(pp).separated()
)
if (isA<processorPolyPatch>(pp) && collocatedPatch(pp))
{
// Coincident processor patches. Boundary field of
// cVals and cCoords is opposite cell.
//if (refCast<const processorPolyPatch>(pp).owner())
{
label faceI = pp.start();
@ -474,18 +578,6 @@ void Foam::isoSurface::generateTriPoints
{
if (faceCutType_[faceI] != NOTCUT)
{
label bFaceI = faceI-mesh_.nInternalFaces();
if
(
neiSnapped[bFaceI]
&& (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
)
{
FatalErrorIn("isoSurface::generateTriPoints(..)")
<< "problem:" << abort(FatalError);
}
generateFaceTriPoints
(
cVals,
@ -512,84 +604,6 @@ void Foam::isoSurface::generateTriPoints
}
}
}
else if (isA<emptyPolyPatch>(pp))
{
// Check if field is there (when generating geometry the
// empty patches have been rewritten to be the face centres),
// otherwise use zero-gradient.
label faceI = pp.start();
const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
// Owner value of cVals
scalarField internalVals;
if (fvp.size() == 0)
{
internalVals.setSize(pp.size());
forAll(pp, i)
{
internalVals[i] = cVals[own[pp.start()+i]];
}
}
const scalarField& bVals =
(
fvp.size() > 0
? static_cast<const scalarField&>(fvp)
: internalVals
);
const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
// Owner value of cCoords
Field<Type> internalCoords;
if (pc.size() == 0)
{
internalCoords.setSize(pp.size());
forAll(pp, i)
{
internalCoords[i] = cCoords[own[pp.start()+i]];
}
}
const Field<Type>& bCoords =
(
pc.size() > 0
? static_cast<const Field<Type>&>(pc)
: internalCoords
);
forAll(pp, i)
{
if (faceCutType_[faceI] != NOTCUT)
{
generateFaceTriPoints
(
cVals,
pVals,
cCoords,
pCoords,
snappedPoints,
snappedCc,
snappedPoint,
faceI,
bVals[i],
bCoords[i],
false, // fc not snapped
pTraits<Type>::zero,
triPoints,
triMeshCells
);
}
faceI++;
}
}
else
{
label faceI = pp.start();
@ -642,12 +656,20 @@ template <class Type>
Foam::tmp<Foam::Field<Type> >
Foam::isoSurface::interpolate
(
const volScalarField& cVals,
const scalarField& pVals,
const GeometricField<Type, fvPatchField, volMesh>& cCoords,
const Field<Type>& pCoords
) const
{
// Recalculate boundary values
tmp<SlicedGeometricField
<
Type,
fvPatchField,
slicedFvPatchField,
volMesh
> > c2(adaptPatchFields(cCoords));
DynamicList<Type> triPoints(nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
@ -658,10 +680,10 @@ Foam::isoSurface::interpolate
generateTriPoints
(
cVals,
pVals,
cValsPtr_(),
pVals_,
cCoords,
c2(),
pCoords,
snappedPoints,

16
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
View File

@ -71,13 +71,7 @@ Foam::sampledIsoSurface::interpolateField
volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
// Sample.
return surface().interpolate
(
*volSubFieldPtr_,
*pointSubFieldPtr_,
volSubFld,
tpointSubFld()
);
return surface().interpolate(volSubFld, tpointSubFld());
}
else
{
@ -85,13 +79,7 @@ Foam::sampledIsoSurface::interpolateField
volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
// Sample.
return surface().interpolate
(
*volFieldPtr_,
*pointFieldPtr_,
volFld,
tpointFld()
);
return surface().interpolate(volFld, tpointFld());
}
}

16
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
View File

@ -67,13 +67,7 @@ Foam::sampledCuttingPlane::interpolateField
volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
// Sample.
return surface().interpolate
(
cellDistancePtr_(),
pointDistance_,
volSubFld,
tpointSubFld()
);
return surface().interpolate(volSubFld, tpointSubFld());
}
else
{
@ -83,13 +77,7 @@ Foam::sampledCuttingPlane::interpolateField
);
// Sample.
return surface().interpolate
(
cellDistancePtr_(),
pointDistance_,
volFld,
tpointFld()
);
return surface().interpolate(volFld, tpointFld());
}
}

4
tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
View File

@ -26,7 +26,7 @@ Application
MRFSimpleFoam
Description
Steady-state solver for incompressible, turbulent flow of non-Newtonian
Steady-state solver for incompressible, turbulent flow of non-Newtonian
fluids with MRF regions.
\*---------------------------------------------------------------------------*/
@ -63,10 +63,10 @@ int main(int argc, char *argv[])
// Pressure-velocity SIMPLE corrector
{
// Momentum predictor
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);
mrfZones.addCoriolis(UEqn());

6
tutorials/incompressible/simpleSRFFoam/simpleSRFFoam/UEqn.H → tutorials/incompressible/simpleSRFFoam/simpleSRFFoam/UrelEqn.H
View File

@ -1,9 +1,9 @@
// Momentum predictor
// Relative momentum predictor
tmp<fvVectorMatrix> UrelEqn
(
fvm::div(phi, Urel)
+ turbulence->divDevReff(Urel)
+ SRF->Su()
+ turbulence->divDevReff(Urel)
+ SRF->Su()
);
UrelEqn().relax();

28
tutorials/incompressible/simpleSRFFoam/simpleSRFFoam/simpleSRFFoam.C
View File

@ -27,7 +27,7 @@ Application
Description
Steady-state solver for incompressible, turbulent flow of non-Newtonian
fluids with single rotating frame.
fluids in a single rotating frame.
\*---------------------------------------------------------------------------*/
@ -40,16 +40,13 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
//mesh.clearPrimitives();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -57,20 +54,19 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSIMPLEControls.H"
# include "initConvergenceCheck.H"
#include "readSIMPLEControls.H"
#include "initConvergenceCheck.H"
p.storePrevIter();
// Pressure-velocity SIMPLE corrector
{
# include "UEqn.H"
# include "pEqn.H"
#include "UrelEqn.H"
#include "pEqn.H"
}
turbulence->correct();
if (runTime.outputTime())
{
volVectorField Uabs
@ -93,7 +89,7 @@ int main(int argc, char *argv[])
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
# include "convergenceCheck.H"
#include "convergenceCheck.H"
}
Info<< "End\n" << endl;