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 branch 'master' into dsmc

Browse Source
This commit is contained in:
graham
2009-11-20 11:31:30 +00:00
parent 35f0c0ab8d 8ce765e3f5
commit fabbd8bb47
116 changed files with 4333 additions and 5805 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H
View File

@ -61,8 +61,16 @@ namespace Foam
template<class T, class Container> class CompactListList;
template<class T, class Container> Istream& operator>>(Istream&, CompactListList<T, Container>&);
template<class T, class Container> Ostream& operator<<(Ostream&, const CompactListList<T, Container>&);
template<class T, class Container> Istream& operator>>
(
Istream&,
CompactListList<T, Container>&
);
template<class T, class Container> Ostream& operator<<
(
Ostream&,
const CompactListList<T, Container>&
);
/*---------------------------------------------------------------------------*\

77
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
View File

@ -52,20 +52,11 @@ Foam::tmp
typename Foam::GeometricField<Type, PatchField, GeoMesh>::
GeometricBoundaryField
>
Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
Foam::GeometricField<Type, PatchField, GeoMesh>::readField
(
const dictionary& fieldDict
)
{
if (is.version() < 2.0)
{
FatalIOErrorIn
(
"GeometricField<Type, PatchField, GeoMesh>::readField(Istream&)",
is
) << "IO versions < 2.0 are not supported."
<< exit(FatalIOError);
}
dictionary fieldDict(is);
DimensionedField<Type, GeoMesh>::readField(fieldDict, "internalField");
tmp<GeometricBoundaryField> tboundaryField
@ -96,6 +87,28 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
}
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::tmp
<
typename Foam::GeometricField<Type, PatchField, GeoMesh>::
GeometricBoundaryField
>
Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
{
if (is.version() < 2.0)
{
FatalIOErrorIn
(
"GeometricField<Type, PatchField, GeoMesh>::readField(Istream&)",
is
) << "IO versions < 2.0 are not supported."
<< exit(FatalIOError);
}
return readField(dictionary(is));
}
template<class Type, template<class> class PatchField, class GeoMesh>
bool Foam::GeometricField<Type, PatchField, GeoMesh>::readIfPresent()
{
@ -404,6 +417,44 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
}
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
(
const IOobject& io,
const Mesh& mesh,
const dictionary& dict
)
:
DimensionedField<Type, GeoMesh>(io, mesh, dimless),
timeIndex_(this->time().timeIndex()),
field0Ptr_(NULL),
fieldPrevIterPtr_(NULL),
boundaryField_(*this, readField(dict))
{
// Check compatibility between field and mesh
if (this->size() != GeoMesh::size(this->mesh()))
{
FatalErrorIn
(
"GeometricField<Type, PatchField, GeoMesh>::GeometricField"
"(const IOobject&, const Mesh&, const dictionary&)"
) << " number of field elements = " << this->size()
<< " number of mesh elements = " << GeoMesh::size(this->mesh())
<< exit(FatalIOError);
}
readOldTimeIfPresent();
if (debug)
{
Info<< "Finishing dictionary-construct of "
"GeometricField<Type, PatchField, GeoMesh>"
<< endl << this->info() << endl;
}
}
// construct as copy
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField

11
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
View File

@ -240,6 +240,9 @@ private:
// Private member functions
//- Read the field from the dictionary
tmp<GeometricBoundaryField> readField(const dictionary&);
//- Read the field from the given stream
tmp<GeometricBoundaryField> readField(Istream& is);
@ -327,6 +330,14 @@ public:
Istream&
);
//- Construct from dictionary
GeometricField
(
const IOobject&,
const Mesh&,
const dictionary&
);
//- Construct as copy
GeometricField
(

57
src/Pstream/mpi/UIPread.C
View File

@ -66,6 +66,14 @@ Foam::UIPstream::UIPstream
label wantedSize = externalBuf_.capacity();
if (debug)
{
Pout<< "UIPstream::UIPstream : read from:" << fromProcNo
<< " tag:" << tag << " wanted size:" << wantedSize
<< Foam::endl;
}
// If the buffer size is not specified, probe the incomming message
// and set it
if (!wantedSize)
@ -75,6 +83,12 @@ Foam::UIPstream::UIPstream
externalBuf_.setCapacity(messageSize_);
wantedSize = messageSize_;
if (debug)
{
Pout<< "UIPstream::UIPstream : probed size:" << wantedSize
<< Foam::endl;
}
}
messageSize_ = UIPstream::read
@ -127,6 +141,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
if (commsType() == UPstream::nonBlocking)
{
// Message is already received into externalBuf
messageSize_ = buffers.recvBuf_[fromProcNo].size();
}
else
{
@ -134,6 +149,14 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
label wantedSize = externalBuf_.capacity();
if (debug)
{
Pout<< "UIPstream::UIPstream PstreamBuffers :"
<< " read from:" << fromProcNo
<< " tag:" << tag_ << " wanted size:" << wantedSize
<< Foam::endl;
}
// If the buffer size is not specified, probe the incomming message
// and set it
if (!wantedSize)
@ -143,6 +166,12 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
externalBuf_.setCapacity(messageSize_);
wantedSize = messageSize_;
if (debug)
{
Pout<< "UIPstream::UIPstream PstreamBuffers : probed size:"
<< wantedSize << Foam::endl;
}
}
messageSize_ = UIPstream::read
@ -180,6 +209,14 @@ Foam::label Foam::UIPstream::read
const int tag
)
{
if (debug)
{
Pout<< "UIPstream::read : starting read from:" << fromProcNo
<< " tag:" << tag << " wanted size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< Foam::endl;
}
if (commsType == blocking || commsType == scheduled)
{
MPI_Status status;
@ -214,6 +251,14 @@ Foam::label Foam::UIPstream::read
label messageSize;
MPI_Get_count(&status, MPI_BYTE, &messageSize);
if (debug)
{
Pout<< "UIPstream::read : finished read from:" << fromProcNo
<< " tag:" << tag << " read size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< Foam::endl;
}
if (messageSize > bufSize)
{
FatalErrorIn
@ -256,6 +301,15 @@ Foam::label Foam::UIPstream::read
return 0;
}
if (debug)
{
Pout<< "UIPstream::read : started read from:" << fromProcNo
<< " tag:" << tag << " read size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< " request:" << PstreamGlobals::outstandingRequests_.size()
<< Foam::endl;
}
PstreamGlobals::outstandingRequests_.append(request);
// Assume the message is completely received.
@ -267,7 +321,8 @@ Foam::label Foam::UIPstream::read
(
"UIPstream::read"
"(const int fromProcNo, char* buf, std::streamsize bufSize)"
) << "Unsupported communications type " << commsType
) << "Unsupported communications type "
<< commsType
<< Foam::abort(FatalError);
return 0;

36
src/Pstream/mpi/UOPwrite.C
View File

@ -43,6 +43,14 @@ bool Foam::UOPstream::write
const int tag
)
{
if (debug)
{
Pout<< "UIPstream::write : starting write to:" << toProcNo
<< " tag:" << tag << " size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< Foam::endl;
}
bool transferFailed = true;
if (commsType == blocking)
@ -56,6 +64,14 @@ bool Foam::UOPstream::write
tag,
MPI_COMM_WORLD
);
if (debug)
{
Pout<< "UIPstream::write : finished write to:" << toProcNo
<< " tag:" << tag << " size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< Foam::endl;
}
}
else if (commsType == scheduled)
{
@ -68,6 +84,14 @@ bool Foam::UOPstream::write
tag,
MPI_COMM_WORLD
);
if (debug)
{
Pout<< "UIPstream::write : finished write to:" << toProcNo
<< " tag:" << tag << " size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< Foam::endl;
}
}
else if (commsType == nonBlocking)
{
@ -84,6 +108,15 @@ bool Foam::UOPstream::write
&request
);
if (debug)
{
Pout<< "UIPstream::write : started write to:" << toProcNo
<< " tag:" << tag << " size:" << label(bufSize)
<< " commsType:" << UPstream::commsTypeNames[commsType]
<< " request:" << PstreamGlobals::outstandingRequests_.size()
<< Foam::endl;
}
PstreamGlobals::outstandingRequests_.append(request);
}
else
@ -93,7 +126,8 @@ bool Foam::UOPstream::write
"UOPstream::write"
"(const int fromProcNo, char* buf, std::streamsize bufSize"
", const int)"
) << "Unsupported communications type " << commsType
) << "Unsupported communications type "
<< UPstream::commsTypeNames[commsType]
<< Foam::abort(FatalError);
}

45
src/Pstream/mpi/UPstream.C
View File

@ -70,6 +70,12 @@ bool Foam::UPstream::init(int& argc, char**& argv)
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myProcNo_);
if (debug)
{
Pout<< "UPstream::init : initialised with numProcs:" << numprocs
<< " myProcNo:" << myProcNo_ << endl;
}
if (numprocs <= 1)
{
FatalErrorIn("UPstream::init(int& argc, char**& argv)")
@ -124,6 +130,11 @@ bool Foam::UPstream::init(int& argc, char**& argv)
void Foam::UPstream::exit(int errnum)
{
if (debug)
{
Pout<< "UPstream::exit." << endl;
}
# ifndef SGIMPI
int size;
char* buff;
@ -164,6 +175,11 @@ void Foam::UPstream::abort()
void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
{
if (Pstream::debug)
{
Pout<< "Foam::reduce : value:" << Value << endl;
}
if (!UPstream::parRun())
{
return;
@ -433,11 +449,23 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
}
*/
}
if (Pstream::debug)
{
Pout<< "Foam::reduce : reduced value:" << Value << endl;
}
}
void Foam::UPstream::waitRequests()
{
if (debug)
{
Pout<< "UPstream::waitRequests : starting wait for "
<< PstreamGlobals::outstandingRequests_.size()
<< " outstanding requests." << endl;
}
if (PstreamGlobals::outstandingRequests_.size())
{
if
@ -458,11 +486,22 @@ void Foam::UPstream::waitRequests()
PstreamGlobals::outstandingRequests_.clear();
}
if (debug)
{
Pout<< "UPstream::waitRequests : finished wait." << endl;
}
}
bool Foam::UPstream::finishedRequest(const label i)
{
if (debug)
{
Pout<< "UPstream::waitRequests : starting wait for request:" << i
<< endl;
}
if (i >= PstreamGlobals::outstandingRequests_.size())
{
FatalErrorIn
@ -483,6 +522,12 @@ bool Foam::UPstream::finishedRequest(const label i)
MPI_STATUS_IGNORE
);
if (debug)
{
Pout<< "UPstream::waitRequests : finished wait for request:" << i
<< endl;
}
return flag != 0;
}

1
src/dynamicMesh/Make/files
View File

@ -49,7 +49,6 @@ perfectInterface/perfectInterface.C
boundaryMesh/boundaryMesh.C
boundaryPatch/boundaryPatch.C
boundaryMesh/octreeDataFaceList.C
setUpdater/setUpdater.C
meshModifiers = meshCut/meshModifiers

138
src/dynamicMesh/boundaryMesh/boundaryMesh.C
View File

@ -29,11 +29,12 @@ License
#include "polyMesh.H"
#include "repatchPolyTopoChanger.H"
#include "faceList.H"
#include "octree.H"
#include "octreeDataFaceList.H"
#include "indexedOctree.H"
#include "treeDataPrimitivePatch.H"
#include "triSurface.H"
#include "SortableList.H"
#include "OFstream.H"
#include "uindirectPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -892,6 +893,17 @@ Foam::labelList Foam::boundaryMesh::getNearest
<< endl;
}
uindirectPrimitivePatch leftPatch
(
UIndirectList<face>(mesh(), leftFaces),
mesh().points()
);
uindirectPrimitivePatch rightPatch
(
UIndirectList<face>(mesh(), rightFaces),
mesh().points()
);
// Overall bb
treeBoundBox overallBb(mesh().localPoints());
@ -911,29 +923,35 @@ Foam::labelList Foam::boundaryMesh::getNearest
bbMax.z() += 2*tol;
// Create the octrees
octree<octreeDataFaceList> leftTree
indexedOctree
<
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
> leftTree
(
overallBb,
octreeDataFaceList
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
(
mesh(),
leftFaces
false, // cacheBb
leftPatch
),
1,
20,
10
overallBb,
10, // maxLevel
10, // leafSize
3.0 // duplicity
);
octree<octreeDataFaceList> rightTree
indexedOctree
<
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
> rightTree
(
overallBb,
octreeDataFaceList
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
(
mesh(),
rightFaces
false, // cacheBb
rightPatch
),
1,
20,
10
overallBb,
10, // maxLevel
10, // leafSize
3.0 // duplicity
);
if (debug)
@ -953,7 +971,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
treeBoundBox tightest;
const scalar searchDim = mag(searchSpan);
const scalar searchDimSqr = magSqr(searchSpan);
forAll(nearestBFaceI, patchFaceI)
{
@ -982,50 +1000,25 @@ Foam::labelList Foam::boundaryMesh::getNearest
}
// Search right tree
tightest.min() = ctr - searchSpan;
tightest.max() = ctr + searchSpan;
scalar rightDist = searchDim;
label rightI = rightTree.findNearest(ctr, tightest, rightDist);
pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
// Search left tree. Note: could start from rightDist bounding box
// instead of starting from top.
tightest.min() = ctr - searchSpan;
tightest.max() = ctr + searchSpan;
scalar leftDist = searchDim;
label leftI = leftTree.findNearest(ctr, tightest, leftDist);
pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
if (rightI == -1)
if (rightInfo.hit())
{
// No face found in right tree.
if (leftI == -1)
{
// No face found in left tree.
nearestBFaceI[patchFaceI] = -1;
}
else
{
// Found in left but not in right. Choose left regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = leftFaces[leftI];
}
}
else
{
if (leftI == -1)
{
// Found in right but not in left. Choose right regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = rightFaces[rightI];
}
else
if (leftInfo.hit())
{
// Found in both trees. Compare normals.
label rightFaceI = rightFaces[rightInfo.index()];
label leftFaceI = leftFaces[leftInfo.index()];
scalar rightSign = n & ns[rightFaces[rightI]];
scalar leftSign = n & ns[leftFaces[leftI]];
label rightDist = mag(rightInfo.hitPoint()-ctr);
label leftDist = mag(leftInfo.hitPoint()-ctr);
scalar rightSign = n & ns[rightFaceI];
scalar leftSign = n & ns[leftFaceI];
if
(
@ -1036,11 +1029,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Both same sign. Choose nearest.
if (rightDist < leftDist)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
else
@ -1059,11 +1052,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Different sign and nearby. Choosing matching normal
if (rightSign > 0)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
else
@ -1071,15 +1064,38 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Different sign but faraway. Choosing nearest.
if (rightDist < leftDist)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
}
}
else
{
// Found in right but not in left. Choose right regardless if
// correct sign. Note: do we want this?
label rightFaceI = rightFaces[rightInfo.index()];
nearestBFaceI[patchFaceI] = rightFaceI;
}
}
else
{
// No face found in right tree.
if (leftInfo.hit())
{
// Found in left but not in right. Choose left regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()];
}
else
{
// No face found in left tree.
nearestBFaceI[patchFaceI] = -1;
}
}
}

573
src/dynamicMesh/boundaryMesh/octreeDataFaceList.C
View File

@ -1,573 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/
#include "octreeDataFaceList.H"
#include "octree.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeDataFaceList, 0);
Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::octreeDataFaceList::calcBb()
{
allBb_.setSize(faceLabels_.size());
allBb_ = treeBoundBox
(
vector(GREAT, GREAT, GREAT),
vector(-GREAT, -GREAT, -GREAT)
);
forAll (faceLabels_, faceLabelI)
{
label faceI = faceLabels_[faceLabelI];
// Update bb of face
treeBoundBox& myBb = allBb_[faceLabelI];
const face& f = mesh_.localFaces()[faceI];
forAll(f, fp)
{
const point& coord = mesh_.localPoints()[f[fp]];
myBb.min() = min(myBb.min(), coord);
myBb.max() = max(myBb.max(), coord);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from all faces in bMesh
Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh)
:
mesh_(mesh),
faceLabels_(mesh_.size()),
allBb_(mesh_.size())
{
forAll(faceLabels_, faceI)
{
faceLabels_[faceI] = faceI;
}
// Generate tight fitting bounding box
calcBb();
}
// Construct from selected faces in bMesh
Foam::octreeDataFaceList::octreeDataFaceList
(
const bMesh& mesh,
const labelList& faceLabels
)
:
mesh_(mesh),
faceLabels_(faceLabels),
allBb_(faceLabels.size())
{
// Generate tight fitting bounding box
calcBb();
}
// Construct as copy
Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes)
:
mesh_(shapes.mesh()),
faceLabels_(shapes.faceLabels()),
allBb_(shapes.allBb())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::octreeDataFaceList::~octreeDataFaceList()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataFaceList::getSampleType
(
const octree<octreeDataFaceList>& oc,
const point& sample
) const
{
// Need to determine whether sample is 'inside' or 'outside'
// Done by finding nearest face. This gives back a face which is
// guaranteed to contain nearest point. This point can be
// - in interior of face: compare to face normal
// - on edge of face: compare to edge normal
// - on point of face: compare to point normal
// Unfortunately the octree does not give us back the intersection point
// or where on the face it has hit so we have to recreate all that
// information.
// Find nearest face to sample
treeBoundBox tightest(treeBoundBox::greatBox);
scalar tightestDist = GREAT;
label index = oc.findNearest(sample, tightest, tightestDist);
if (index == -1)
{
FatalErrorIn
(
"octreeDataFaceList::getSampleType"
"(octree<octreeDataFaceList>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
label faceI = faceLabels_[index];
// Get actual intersection point on face
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
const face& f = mesh_.localFaces()[faceI];
const pointField& points = mesh_.localPoints();
pointHit curHit = f.nearestPoint(sample, points);
//
// 1] Check whether sample is above face
//
if (curHit.hit())
{
// Simple case. Compare to face normal.
if (debug & 2)
{
Pout<< " -> face hit:" << curHit.hitPoint()
<< " comparing to face normal " << mesh_.faceNormals()[faceI]
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.faceNormals()[faceI],
sample - curHit.hitPoint()
);
}
if (debug & 2)
{
Pout<< " -> face miss:" << curHit.missPoint();
}
//
// 2] Check whether intersection is on one of the face vertices or
// face centre
//
// typical dimension as sqrt of face area.
scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
forAll(f, fp)
{
if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face vertex fp
if (debug & 2)
{
Pout<< " -> face point hit :" << points[f[fp]]
<< " point normal:" << mesh_.pointNormals()[f[fp]]
<< " distance:"
<< mag(points[f[fp]] - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.pointNormals()[f[fp]],
sample - curHit.missPoint()
);
}
}
// Get face centre
point ctr(f.centre(points));
if ((mag(ctr - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
if (debug & 2)
{
Pout<< " -> centre hit:" << ctr
<< " distance:"
<< mag(ctr - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.faceNormals()[faceI],
sample - curHit.missPoint()
);
}
//
// 3] Get the 'real' edge the face intersection is on
//
const labelList& myEdges = mesh_.faceEdges()[faceI];
forAll(myEdges, myEdgeI)
{
const edge& e = mesh_.edges()[myEdges[myEdgeI]];
pointHit edgeHit =
line<point, const point&>
(
points[e.start()],
points[e.end()]
).nearestDist(sample);
point edgePoint;
if (edgeHit.hit())
{
edgePoint = edgeHit.hitPoint();
}
else
{
edgePoint = edgeHit.missPoint();
}
if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge e
// Calculate edge normal (wrong: uses face normals instead of
// triangle normals)
const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
vector edgeNormal(vector::zero);
forAll(myFaces, myFaceI)
{
edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
}
if (debug & 2)
{
Pout<< " -> real edge hit point:" << edgePoint
<< " comparing to edge normal:" << edgeNormal
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFaceList>::getVolType
(
edgeNormal,
sample - curHit.missPoint()
);
}
}
//
// 4] Get the internal edge (vertex - faceCentre) the face intersection
// is on
//
forAll(f, fp)
{
pointHit edgeHit =
line<point, const point&>
(
points[f[fp]],
ctr
).nearestDist(sample);
point edgePoint;
if (edgeHit.hit())
{
edgePoint = edgeHit.hitPoint();
}
else
{
edgePoint = edgeHit.missPoint();
}
if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge between two face triangles
// Calculate edge normal as average of the two triangle normals
label fpPrev = f.rcIndex(fp);
label fpNext = f.fcIndex(fp);
vector e = points[f[fp]] - ctr;
vector ePrev = points[f[fpPrev]] - ctr;
vector eNext = points[f[fpNext]] - ctr;
vector nLeft = ePrev ^ e;
nLeft /= mag(nLeft) + VSMALL;
vector nRight = e ^ eNext;
nRight /= mag(nRight) + VSMALL;
if (debug & 2)
{
Pout<< " -> internal edge hit point:" << edgePoint
<< " comparing to edge normal "
<< 0.5*(nLeft + nRight)
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFaceList>::getVolType
(
0.5*(nLeft + nRight),
sample - curHit.missPoint()
);
}
}
if (debug & 2)
{
Pout<< "Did not find sample " << sample
<< " anywhere related to nearest face " << faceI << endl
<< "Face:";
forAll(f, fp)
{
Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
<< endl;
}
}
// Can't determine status of sample with respect to nearest face.
// Either
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return octree<octreeDataFaceList>::UNKNOWN;
}
bool Foam::octreeDataFaceList::overlaps
(
const label index,
const treeBoundBox& sampleBb
) const
{
return sampleBb.overlaps(allBb_[index]);
}
bool Foam::octreeDataFaceList::contains
(
const label,
const point&
) const
{
notImplemented
(
"octreeDataFaceList::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataFaceList::intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
const vector dir(end - start);
// Disable picking up intersections behind us.
scalar oldTol = intersection::setPlanarTol(0.0);
pointHit inter =
f.ray
(
start,
dir,
mesh_.localPoints(),
intersection::HALF_RAY,
intersection::VECTOR
);
intersection::setPlanarTol(oldTol);
if (inter.hit() && inter.distance() <= mag(dir))
{
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
bool Foam::octreeDataFaceList::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// Get nearest and furthest away vertex
point myNear, myFar;
allBb_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataFaceList::calcSign
(
const label index,
const point& sample,
vector&
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
point ctr = f.centre(mesh_.localPoints());
vector vec = sample - ctr;
vec /= mag(vec) + VSMALL;
return mesh_.faceNormals()[faceI] & vec;
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataFaceList::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
if (nearHit.hit())
{
nearest = nearHit.hitPoint();
}
else
{
nearest = nearHit.missPoint();
}
if (debug & 1)
{
point ctr = f.centre(mesh_.localPoints());
scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
Pout<< "octreeDataFaceList::calcNearest : "
<< "sample:" << sample
<< " faceI:" << faceI
<< " ctr:" << ctr
<< " sign:" << sign
<< " nearest point:" << nearest
<< " distance to face:" << nearHit.distance()
<< endl;
}
return nearHit.distance();
}
void Foam::octreeDataFaceList::write
(
Ostream& os,
const label index
) const
{
os << faceLabels_[index] << " " << allBb_[index];
}
// ************************************************************************* //

220
src/dynamicMesh/boundaryMesh/octreeDataFaceList.H
View File

@ -1,220 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::octreeDataFaceList
Description
Holds data for octree to work on list of faces on a bMesh
(= PrimitivePatch which holds faces, not references them)
Same as octreeDataFace except for that.
SourceFiles
octreeDataFaceList.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataFaceList_H
#define octreeDataFaceList_H
#include "treeBoundBoxList.H"
#include "faceList.H"
#include "point.H"
#include "className.H"
#include "bMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataFaceList Declaration
\*---------------------------------------------------------------------------*/
class octreeDataFaceList
{
// Static data
//- tolerance on linear dimensions
static scalar tol;
// Static function
static inline label nexti(label max, label i)
{
return (i + 1) % max;
}
// Private data
//- the mesh
const bMesh& mesh_;
//- labels (in mesh indexing) of faces
labelList faceLabels_;
//- bbs for all above faces
treeBoundBoxList allBb_;
// Private Member Functions
//- Set allBb to tight fitting bounding box
void calcBb();
public:
// Declare name of the class and its debug switch
ClassName("octreeDataFaceList");
// Constructors
//- Construct from all faces in bMesh.
octreeDataFaceList(const bMesh& mesh);
//- Construct from selected faces in bMesh.
octreeDataFaceList(const bMesh& mesh, const labelList& faceLabels);
//- Construct as copy
octreeDataFaceList(const octreeDataFaceList&);
// Destructor
~octreeDataFaceList();
// Member Functions
// Access
const bMesh& mesh() const
{
return mesh_;
}
const labelList& faceLabels() const
{
return faceLabels_;
}
const treeBoundBoxList& allBb() const
{
return allBb_;
}
label size() const
{
return allBb_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataFaceList>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample). Returns GREAT if
// sample does not project onto (triangle decomposition) of face.
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
// Edit
// Write
//- Write shape at index
void write(Ostream& os, const label index) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

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

@ -330,8 +330,6 @@ void Foam::fvMeshAdder::MapVolFields
if (fieldsToAdd.found(fld.name()))
{
Pout<< "Mapping field " << fld.name() << endl;
const GeometricField<Type, fvPatchField, volMesh>& fldToAdd =
*fieldsToAdd[fld.name()];
@ -339,7 +337,7 @@ void Foam::fvMeshAdder::MapVolFields
}
else
{
WarningIn("fvMeshAdder::MapVolFields")
WarningIn("fvMeshAdder::MapVolFields(..)")
<< "Not mapping field " << fld.name()
<< " since not present on mesh to add"
<< endl;
@ -642,15 +640,13 @@ void Foam::fvMeshAdder::MapSurfaceFields
if (fieldsToAdd.found(fld.name()))
{
Pout<< "Mapping field " << fld.name() << endl;
const fldType& fldToAdd = *fieldsToAdd[fld.name()];
MapSurfaceField<Type>(meshMap, fld, fldToAdd);
}
else
{
WarningIn("fvMeshAdder::MapSurfaceFields")
WarningIn("fvMeshAdder::MapSurfaceFields(..)")
<< "Not mapping field " << fld.name()
<< " since not present on mesh to add"
<< endl;

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

@ -25,8 +25,6 @@ License
\*----------------------------------------------------------------------------*/
#include "fvMeshDistribute.H"
#include "ProcessorTopology.H"
#include "commSchedule.H"
#include "PstreamCombineReduceOps.H"
#include "fvMeshAdder.H"
#include "faceCoupleInfo.H"
@ -39,6 +37,8 @@ License
#include "mergePoints.H"
#include "mapDistributePolyMesh.H"
#include "surfaceFields.H"
#include "syncTools.H"
#include "CompactListList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,71 +47,6 @@ defineTypeNameAndDebug(Foam::fvMeshDistribute, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//Foam::List<Foam::labelPair> Foam::fvMeshDistribute::getSchedule
//(
// const labelList& distribution
//)
//{
// labelList nCellsPerProc(countCells(distribution));
//
// if (debug)
// {
// Pout<< "getSchedule : Wanted distribution:" << nCellsPerProc << endl;
// }
//
// // Processors I need to send data to
// labelListList mySendProcs(Pstream::nProcs());
//
// // Count
// label nSendProcs = 0;
// forAll(nCellsPerProc, sendProcI)
// {
// if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0)
// {
// nSendProcs++;
// }
// }
//
// // Fill
// mySendProcs[Pstream::myProcNo()].setSize(nSendProcs);
// nSendProcs = 0;
// forAll(nCellsPerProc, sendProcI)
// {
// if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0)
// {
// mySendProcs[Pstream::myProcNo()][nSendProcs++] = sendProcI;
// }
// }
//
// // Synchronise
// Pstream::gatherList(mySendProcs);
// Pstream::scatterList(mySendProcs);
//
// // Combine into list (same on all procs) giving sending and receiving
// // processor
// label nComms = 0;
// forAll(mySendProcs, procI)
// {
// nComms += mySendProcs[procI].size();
// }
//
// List<labelPair> schedule(nComms);
// nComms = 0;
//
// forAll(mySendProcs, procI)
// {
// const labelList& sendProcs = mySendProcs[procI];
//
// forAll(sendProcs, i)
// {
// schedule[nComms++] = labelPair(procI, sendProcs[i]);
// }
// }
//
// return schedule;
//}
Foam::labelList Foam::fvMeshDistribute::select
(
const bool selectEqual,
@ -144,14 +79,29 @@ Foam::labelList Foam::fvMeshDistribute::select
// Check all procs have same names and in exactly same order.
void Foam::fvMeshDistribute::checkEqualWordList(const wordList& lst)
void Foam::fvMeshDistribute::checkEqualWordList
(
const string& msg,
const wordList& lst
)
{
wordList myObjects(lst);
List<wordList> allNames(Pstream::nProcs());
allNames[Pstream::myProcNo()] = lst;
Pstream::gatherList(allNames);
Pstream::scatterList(allNames);
// Check that all meshes have same objects
Pstream::listCombineGather(myObjects, checkEqualType());
// Below scatter only needed to balance sends and receives.
Pstream::listCombineScatter(myObjects);
for (label procI = 1; procI < Pstream::nProcs(); procI++)
{
if (allNames[procI] != allNames[0])
{
FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)")
<< "When checking for equal " << msg.c_str() << " :" << endl
<< "processor0 has:" << allNames[0] << endl
<< "processor" << procI << " has:" << allNames[procI] << endl
<< msg.c_str() << " need to be synchronised on all processors."
<< exit(FatalError);
}
}
}
@ -664,21 +614,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
<< " newPointI:" << newPointI << abort(FatalError);
}
}
//- old: use pointToMaster map.
//forAll(constructMap, i)
//{
// label oldPointI = constructMap[i];
//
// // See if merged into other point
// Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
//
// if (iter != pointToMaster.end())
// {
// oldPointI = iter();
// }
//
// constructMap[i] = map().reversePointMap()[oldPointI];
//}
}
return map;
}
@ -693,46 +628,49 @@ void Foam::fvMeshDistribute::getNeighbourData
labelList& sourceNewProc
) const
{
sourceFace.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
sourceProc.setSize(sourceFace.size());
sourceNewProc.setSize(sourceFace.size());
label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
sourceFace.setSize(nBnd);
sourceProc.setSize(nBnd);
sourceNewProc.setSize(nBnd);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Send meshFace labels of processor patches and destination processor.
// Get neighbouring meshFace labels and new processor of coupled boundaries.
labelList nbrFaces(nBnd, -1);
labelList nbrNewProc(nBnd, -1);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
label offset = pp.start() - mesh_.nInternalFaces();
// Labels of faces on this side
labelList meshFaceLabels(pp.size());
forAll(meshFaceLabels, i)
// Mesh labels of faces on this side
forAll(pp, i)
{
meshFaceLabels[i] = pp.start()+i;
label bndI = offset + i;
nbrFaces[bndI] = pp.start()+i;
}
// Which processor they will end up on
const labelList newProc
SubList<label>(nbrNewProc, pp.size(), offset).assign
(
UIndirectList<label>(distribution, pp.faceCells())
UIndirectList<label>(distribution, pp.faceCells())()
);
OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
toNeighbour << meshFaceLabels << newProc;
}
}
// Receive meshFace labels and destination processors of processor faces.
// Exchange the boundary data
syncTools::swapBoundaryFaceList(mesh_, nbrFaces, false);
syncTools::swapBoundaryFaceList(mesh_, nbrNewProc, false);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
label offset = pp.start() - mesh_.nInternalFaces();
if (isA<processorPolyPatch>(pp))
@ -740,11 +678,6 @@ void Foam::fvMeshDistribute::getNeighbourData
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
// Receive the data
IPstream fromNeighbour(Pstream::blocking, procPatch.neighbProcNo());
labelList nbrFaces(fromNeighbour);
labelList nbrNewProc(fromNeighbour);
// Check which of the two faces we store.
if (Pstream::myProcNo() < procPatch.neighbProcNo())
@ -752,9 +685,10 @@ void Foam::fvMeshDistribute::getNeighbourData
// Use my local face labels
forAll(pp, i)
{
sourceFace[offset + i] = pp.start()+i;
sourceProc[offset + i] = Pstream::myProcNo();
sourceNewProc[offset + i] = nbrNewProc[i];
label bndI = offset + i;
sourceFace[bndI] = pp.start()+i;
sourceProc[bndI] = Pstream::myProcNo();
sourceNewProc[bndI] = nbrNewProc[bndI];
}
}
else
@ -762,9 +696,10 @@ void Foam::fvMeshDistribute::getNeighbourData
// Use my neighbours face labels
forAll(pp, i)
{
sourceFace[offset + i] = nbrFaces[i];
sourceProc[offset + i] = procPatch.neighbProcNo();
sourceNewProc[offset + i] = nbrNewProc[i];
label bndI = offset + i;
sourceFace[bndI] = nbrFaces[bndI];
sourceProc[bndI] = procPatch.neighbProcNo();
sourceNewProc[bndI] = nbrNewProc[bndI];
}
}
}
@ -773,9 +708,10 @@ void Foam::fvMeshDistribute::getNeighbourData
// Normal (physical) boundary
forAll(pp, i)
{
sourceFace[offset + i] = patchI;
sourceProc[offset + i] = -1;
sourceNewProc[offset + i] = -1;
label bndI = offset + i;
sourceFace[bndI] = patchI;
sourceProc[bndI] = -1;
sourceNewProc[bndI] = -1;
}
}
}
@ -1120,7 +1056,8 @@ void Foam::fvMeshDistribute::sendMesh
const labelList& sourceFace,
const labelList& sourceProc,
const labelList& sourceNewProc
const labelList& sourceNewProc,
UOPstream& toDomain
)
{
if (debug)
@ -1133,52 +1070,95 @@ void Foam::fvMeshDistribute::sendMesh
<< endl;
}
// Assume sparse point zones. Get contents in merged-zone indices.
labelListList zonePoints(pointZoneNames.size());
// Assume sparse, possibly overlapping point zones. Get contents
// in merged-zone indices.
CompactListList<label> zonePoints;
{
const pointZoneMesh& pointZones = mesh.pointZones();
labelList rowSizes(pointZoneNames.size(), 0);
forAll(pointZoneNames, nameI)
{
label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
if (myZoneID != -1)
{
zonePoints[nameI] = pointZones[myZoneID];
rowSizes[nameI] = pointZones[myZoneID].size();
}
}
zonePoints.setSize(rowSizes);
forAll(pointZoneNames, nameI)
{
label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
if (myZoneID != -1)
{
zonePoints[nameI].assign(pointZones[myZoneID]);
}
}
}
// Assume sparse face zones
labelListList zoneFaces(faceZoneNames.size());
boolListList zoneFaceFlip(faceZoneNames.size());
// Assume sparse, possibly overlapping face zones
CompactListList<label> zoneFaces;
CompactListList<bool> zoneFaceFlip;
{
const faceZoneMesh& faceZones = mesh.faceZones();
labelList rowSizes(faceZoneNames.size(), 0);
forAll(faceZoneNames, nameI)
{
label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
if (myZoneID != -1)
{
zoneFaces[nameI] = faceZones[myZoneID];
zoneFaceFlip[nameI] = faceZones[myZoneID].flipMap();
rowSizes[nameI] = faceZones[myZoneID].size();
}
}
zoneFaces.setSize(rowSizes);
zoneFaceFlip.setSize(rowSizes);
forAll(faceZoneNames, nameI)
{
label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
if (myZoneID != -1)
{
zoneFaces[nameI].assign(faceZones[myZoneID]);
zoneFaceFlip[nameI].assign(faceZones[myZoneID].flipMap());
}
}
}
// Assume sparse cell zones
labelListList zoneCells(cellZoneNames.size());
// Assume sparse, possibly overlapping cell zones
CompactListList<label> zoneCells;
{
const cellZoneMesh& cellZones = mesh.cellZones();
labelList rowSizes(pointZoneNames.size(), 0);
forAll(cellZoneNames, nameI)
{
label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
if (myZoneID != -1)
{
zoneCells[nameI] = cellZones[myZoneID];
rowSizes[nameI] = cellZones[myZoneID].size();
}
}
zoneCells.setSize(rowSizes);
forAll(cellZoneNames, nameI)
{
label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
if (myZoneID != -1)
{
zoneCells[nameI].assign(cellZones[myZoneID]);
}
}
}
@ -1198,10 +1178,9 @@ void Foam::fvMeshDistribute::sendMesh
//}
// Send
OPstream toDomain(Pstream::blocking, domain);
toDomain
<< mesh.points()
<< mesh.faces()
<< CompactListList<label, face>(mesh.faces())
<< mesh.faceOwner()
<< mesh.faceNeighbour()
<< mesh.boundaryMesh()
@ -1214,6 +1193,13 @@ void Foam::fvMeshDistribute::sendMesh
<< sourceFace
<< sourceProc
<< sourceNewProc;
if (debug)
{
Pout<< "Started sending mesh to domain " << domain
<< endl;
}
}
@ -1227,21 +1213,20 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
const Time& runTime,
labelList& domainSourceFace,
labelList& domainSourceProc,
labelList& domainSourceNewProc
labelList& domainSourceNewProc,
UIPstream& fromNbr
)
{
IPstream fromNbr(Pstream::blocking, domain);
pointField domainPoints(fromNbr);
faceList domainFaces(fromNbr);
faceList domainFaces = CompactListList<label, face>(fromNbr)();
labelList domainAllOwner(fromNbr);
labelList domainAllNeighbour(fromNbr);
PtrList<entry> patchEntries(fromNbr);
labelListList zonePoints(fromNbr);
labelListList zoneFaces(fromNbr);
boolListList zoneFaceFlip(fromNbr);
labelListList zoneCells(fromNbr);
CompactListList<label> zonePoints(fromNbr);
CompactListList<label> zoneFaces(fromNbr);
CompactListList<bool> zoneFaceFlip(fromNbr);
CompactListList<label> zoneCells(fromNbr);
fromNbr
>> domainSourceFace
@ -1442,6 +1427,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
// Local environment of all boundary faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -1481,36 +1467,39 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
mesh_.clearOut();
mesh_.resetMotion();
// Get data to send. Make sure is synchronised
const wordList volScalars(mesh_.names(volScalarField::typeName));
checkEqualWordList(volScalars);
checkEqualWordList("volScalarFields", volScalars);
const wordList volVectors(mesh_.names(volVectorField::typeName));
checkEqualWordList(volVectors);
checkEqualWordList("volVectorFields", volVectors);
const wordList volSphereTensors
(
mesh_.names(volSphericalTensorField::typeName)
);
checkEqualWordList(volSphereTensors);
checkEqualWordList("volSphericalTensorFields", volSphereTensors);
const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
checkEqualWordList(volSymmTensors);
checkEqualWordList("volSymmTensorFields", volSymmTensors);
const wordList volTensors(mesh_.names(volTensorField::typeName));
checkEqualWordList(volTensors);
checkEqualWordList("volTensorField", volTensors);
const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
checkEqualWordList(surfScalars);
checkEqualWordList("surfaceScalarFields", surfScalars);
const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
checkEqualWordList(surfVectors);
checkEqualWordList("surfaceVectorFields", surfVectors);
const wordList surfSphereTensors
(
mesh_.names(surfaceSphericalTensorField::typeName)
);
checkEqualWordList(surfSphereTensors);
checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
const wordList surfSymmTensors
(
mesh_.names(surfaceSymmTensorField::typeName)
);
checkEqualWordList(surfSymmTensors);
checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
checkEqualWordList(surfTensors);
checkEqualWordList("surfaceTensorFields", surfTensors);
// Find patch to temporarily put exposed and processor faces into.
@ -1589,6 +1578,9 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
Pstream::scatterList(nSendCells);
// Allocate buffers
PstreamBuffers pBufs(Pstream::nonBlocking);
// What to send to neighbouring domains
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -1612,6 +1604,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
<< nl << endl;
}
// Pstream for sending mesh and fields
//OPstream str(Pstream::blocking, recvProc);
UOPstream str(recvProc, pBufs);
// Mesh subsetting engine
fvMeshSubset subsetter(mesh_);
@ -1659,6 +1655,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
procSourceNewProc
);
// Send to neighbour
sendMesh
(
@ -1671,43 +1669,69 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
procSourceFace,
procSourceProc,
procSourceNewProc
procSourceNewProc,
str
);
sendFields<volScalarField>(recvProc, volScalars, subsetter);
sendFields<volVectorField>(recvProc, volVectors, subsetter);
sendFields<volScalarField>(recvProc, volScalars, subsetter, str);
sendFields<volVectorField>(recvProc, volVectors, subsetter, str);
sendFields<volSphericalTensorField>
(
recvProc,
volSphereTensors,
subsetter
subsetter,
str
);
sendFields<volSymmTensorField>
(
recvProc,
volSymmTensors,
subsetter
subsetter,
str
);
sendFields<volTensorField>(recvProc, volTensors, subsetter);
sendFields<volTensorField>(recvProc, volTensors, subsetter, str);
sendFields<surfaceScalarField>(recvProc, surfScalars, subsetter);
sendFields<surfaceVectorField>(recvProc, surfVectors, subsetter);
sendFields<surfaceScalarField>
(
recvProc,
surfScalars,
subsetter,
str
);
sendFields<surfaceVectorField>
(
recvProc,
surfVectors,
subsetter,
str
);
sendFields<surfaceSphericalTensorField>
(
recvProc,
surfSphereTensors,
subsetter
subsetter,
str
);
sendFields<surfaceSymmTensorField>
(
recvProc,
surfSymmTensors,
subsetter
subsetter,
str
);
sendFields<surfaceTensorField>
(
recvProc,
surfTensors,
subsetter,
str
);
sendFields<surfaceTensorField>(recvProc, surfTensors, subsetter);
}
}
// Start sending&receiving from buffers
pBufs.finishedSends();
// Subset the part that stays
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -1817,109 +1841,132 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
<< nl << endl;
}
// Pstream for receiving mesh and fields
UIPstream str(sendProc, pBufs);
// Receive from sendProc
labelList domainSourceFace;
labelList domainSourceProc;
labelList domainSourceNewProc;
// Opposite of sendMesh
autoPtr<fvMesh> domainMeshPtr = receiveMesh
(
sendProc,
pointZoneNames,
faceZoneNames,
cellZoneNames,
const_cast<Time&>(mesh_.time()),
domainSourceFace,
domainSourceProc,
domainSourceNewProc
);
fvMesh& domainMesh = domainMeshPtr();
// Receive fields
autoPtr<fvMesh> domainMeshPtr;
PtrList<volScalarField> vsf;
receiveFields<volScalarField>
(
sendProc,
volScalars,
domainMesh,
vsf
);
PtrList<volVectorField> vvf;
receiveFields<volVectorField>
(
sendProc,
volVectors,
domainMesh,
vvf
);
PtrList<volSphericalTensorField> vsptf;
receiveFields<volSphericalTensorField>
(
sendProc,
volSphereTensors,
domainMesh,
vsptf
);
PtrList<volSymmTensorField> vsytf;
receiveFields<volSymmTensorField>
(
sendProc,
volSymmTensors,
domainMesh,
vsytf
);
PtrList<volTensorField> vtf;
receiveFields<volTensorField>
(
sendProc,
volTensors,
domainMesh,
vtf
);
PtrList<surfaceScalarField> ssf;
receiveFields<surfaceScalarField>
(
sendProc,
surfScalars,
domainMesh,
ssf
);
PtrList<surfaceVectorField> svf;
receiveFields<surfaceVectorField>
(
sendProc,
surfVectors,
domainMesh,
svf
);
PtrList<surfaceSphericalTensorField> ssptf;
receiveFields<surfaceSphericalTensorField>
(
sendProc,
surfSphereTensors,
domainMesh,
ssptf
);
PtrList<surfaceSymmTensorField> ssytf;
receiveFields<surfaceSymmTensorField>
(
sendProc,
surfSymmTensors,
domainMesh,
ssytf
);
PtrList<surfaceTensorField> stf;
receiveFields<surfaceTensorField>
(
sendProc,
surfTensors,
domainMesh,
stf
);
// Opposite of sendMesh
{
domainMeshPtr = receiveMesh
(
sendProc,
pointZoneNames,
faceZoneNames,
cellZoneNames,
const_cast<Time&>(mesh_.time()),
domainSourceFace,
domainSourceProc,
domainSourceNewProc,
str
);
fvMesh& domainMesh = domainMeshPtr();
// Receive fields. Read as single dictionary because
// of problems reading consecutive fields from single stream.
dictionary fieldDicts(str);
receiveFields<volScalarField>
(
sendProc,
volScalars,
domainMesh,
vsf,
fieldDicts.subDict(volScalarField::typeName)
);
receiveFields<volVectorField>
(
sendProc,
volVectors,
domainMesh,
vvf,
fieldDicts.subDict(volVectorField::typeName)
);
receiveFields<volSphericalTensorField>
(
sendProc,
volSphereTensors,
domainMesh,
vsptf,
fieldDicts.subDict(volSphericalTensorField::typeName)
);
receiveFields<volSymmTensorField>
(
sendProc,
volSymmTensors,
domainMesh,
vsytf,
fieldDicts.subDict(volSymmTensorField::typeName)
);
receiveFields<volTensorField>
(
sendProc,
volTensors,
domainMesh,
vtf,
fieldDicts.subDict(volTensorField::typeName)
);
receiveFields<surfaceScalarField>
(
sendProc,
surfScalars,
domainMesh,
ssf,
fieldDicts.subDict(surfaceScalarField::typeName)
);
receiveFields<surfaceVectorField>
(
sendProc,
surfVectors,
domainMesh,
svf,
fieldDicts.subDict(surfaceVectorField::typeName)
);
receiveFields<surfaceSphericalTensorField>
(
sendProc,
surfSphereTensors,
domainMesh,
ssptf,
fieldDicts.subDict(surfaceSphericalTensorField::typeName)
);
receiveFields<surfaceSymmTensorField>
(
sendProc,
surfSymmTensors,
domainMesh,
ssytf,
fieldDicts.subDict(surfaceSymmTensorField::typeName)
);
receiveFields<surfaceTensorField>
(
sendProc,
surfTensors,
domainMesh,
stf,
fieldDicts.subDict(surfaceTensorField::typeName)
);
}
const fvMesh& domainMesh = domainMeshPtr();
constructCellMap[sendProc] = identity(domainMesh.nCells());

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

@ -82,38 +82,8 @@ class fvMeshDistribute
const scalar mergeTol_;
// Private classes
//- Check words are the same. Used in patch type checking of similarly
// named patches.
class checkEqualType
{
public:
void operator()(word& x, const word& y) const
{
if (x != y)
{
FatalErrorIn("checkEqualType()(word&, const word&) const")
<< "Patch type " << x << " possibly on processor "
<< Pstream::myProcNo()
<< " does not equal patch type " << y
<< " on some other processor." << nl
<< "Please check similarly named patches for"
<< " having exactly the same type."
<< abort(FatalError);
}
}
};
// Private Member Functions
//- Given distribution work out a communication schedule. Is list
// of pairs. First element of pair is send processor, second is
// receive processor.
//static List<labelPair> getSchedule(const labelList&);
//- Find indices with value
static labelList select
(
@ -123,7 +93,7 @@ class fvMeshDistribute
);
//- Check all procs have same names and in exactly same order.
static void checkEqualWordList(const wordList&);
static void checkEqualWordList(const string&, const wordList&);
//- Merge wordlists over all processors
static wordList mergeWordList(const wordList&);
@ -288,7 +258,8 @@ class fvMeshDistribute
const wordList& cellZoneNames,
const labelList& sourceFace,
const labelList& sourceProc,
const labelList& sourceNewProc
const labelList& sourceNewProc,
UOPstream& toDomain
);
//- Send subset of fields
template<class GeoField>
@ -296,7 +267,8 @@ class fvMeshDistribute
(
const label domain,
const wordList& fieldNames,
const fvMeshSubset&
const fvMeshSubset&,
UOPstream& toNbr
);
//- Receive mesh. Opposite of sendMesh
@ -309,7 +281,8 @@ class fvMeshDistribute
const Time& runTime,
labelList& domainSourceFace,
labelList& domainSourceProc,
labelList& domainSourceNewProc
labelList& domainSourceNewProc,
UIPstream& fromNbr
);
//- Receive fields. Opposite of sendFields
@ -319,7 +292,8 @@ class fvMeshDistribute
const label domain,
const wordList& fieldNames,
fvMesh&,
PtrList<GeoField>&
PtrList<GeoField>&,
const dictionary& fieldDicts
);
//- Disallow default bitwise copy construct

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

@ -184,7 +184,7 @@ void Foam::fvMeshDistribute::mapBoundaryFields
if (flds.size() != oldBflds.size())
{
FatalErrorIn("fvMeshDistribute::mapBoundaryFields") << "problem"
FatalErrorIn("fvMeshDistribute::mapBoundaryFields(..)") << "problem"
<< abort(FatalError);
}
@ -273,19 +273,40 @@ void Foam::fvMeshDistribute::initPatchFields
// Send fields. Note order supplied so we can receive in exactly the same order.
// Note that field gets written as entry in dictionary so we
// can construct from subdictionary.
// (since otherwise the reading as-a-dictionary mixes up entries from
// consecutive fields)
// The dictionary constructed is:
// volScalarField
// {
// p {internalField ..; boundaryField ..;}
// k {internalField ..; boundaryField ..;}
// }
// volVectorField
// {
// U {internalField ... }
// }
// volVectorField {U {internalField ..; boundaryField ..;}}
//
template<class GeoField>
void Foam::fvMeshDistribute::sendFields
(
const label domain,
const wordList& fieldNames,
const fvMeshSubset& subsetter
const fvMeshSubset& subsetter,
UOPstream& toNbr
)
{
toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
//Pout<< "Subsetting field " << fieldNames[i]
// << " for domain:" << domain
// << endl;
if (debug)
{
Pout<< "Subsetting field " << fieldNames[i]
<< " for domain:" << domain << endl;
}
// Send all fieldNames. This has to be exactly the same set as is
// being received!
@ -294,10 +315,12 @@ void Foam::fvMeshDistribute::sendFields
tmp<GeoField> tsubfld = subsetter.interpolate(fld);
// Send
OPstream toNbr(Pstream::blocking, domain);
toNbr << tsubfld();
toNbr
<< fieldNames[i] << token::NL << token::BEGIN_BLOCK
<< tsubfld
<< token::NL << token::END_BLOCK << token::NL;
}
toNbr << token::END_BLOCK << token::NL;
}
@ -308,18 +331,25 @@ void Foam::fvMeshDistribute::receiveFields
const label domain,
const wordList& fieldNames,
fvMesh& mesh,
PtrList<GeoField>& fields
PtrList<GeoField>& fields,
const dictionary& fieldDicts
)
{
if (debug)
{
Pout<< "Receiving fields " << fieldNames
<< " from domain:" << domain << endl;
}
fields.setSize(fieldNames.size());
forAll(fieldNames, i)
{
//Pout<< "Receiving field " << fieldNames[i]
// << " from domain:" << domain
// << endl;
IPstream fromNbr(Pstream::blocking, domain);
if (debug)
{
Pout<< "Constructing field " << fieldNames[i]
<< " from domain:" << domain << endl;
}
fields.set
(
@ -335,7 +365,7 @@ void Foam::fvMeshDistribute::receiveFields
IOobject::AUTO_WRITE
),
mesh,
fromNbr
fieldDicts.subDict(fieldNames[i])
)
);
}

1
src/finiteVolume/Make/files
View File

@ -122,6 +122,7 @@ $(derivedFvPatchFields)/inletOutletTotalTemperature/inletOutletTotalTemperatureF
$(derivedFvPatchFields)/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/movingWallVelocity/movingWallVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/oscillatingFixedValue/oscillatingFixedValueFvPatchFields.C
$(derivedFvPatchFields)/outletInlet/outletInletFvPatchFields.C
$(derivedFvPatchFields)/partialSlip/partialSlipFvPatchFields.C

2
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
View File

@ -154,7 +154,7 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
label sizeby2 = this->size()/2;
const unallocLabelList& faceCells = this->cyclicPatch().faceCells();
if (long(&psiInternal) == long(&this->internalField()))
if (&psiInternal == &this->internalField())
{
tmp<Field<scalar> > tjf = jump();
const Field<scalar>& jf = tjf();

142
src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C Normal file
View File

@ -0,0 +1,142 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "translatingWallVelocityFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
translatingWallVelocityFvPatchVectorField::
translatingWallVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchField<vector>(p, iF),
U_(vector::zero)
{}
translatingWallVelocityFvPatchVectorField::
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField& ptf,
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
U_(ptf.U_)
{}
translatingWallVelocityFvPatchVectorField::
translatingWallVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchField<vector>(p, iF),
U_(dict.lookup("U"))
{
// Evaluate the wall velocity
updateCoeffs();
}
translatingWallVelocityFvPatchVectorField::
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField& twvpvf
)
:
fixedValueFvPatchField<vector>(twvpvf),
U_(twvpvf.U_)
{}
translatingWallVelocityFvPatchVectorField::
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField& twvpvf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchField<vector>(twvpvf, iF),
U_(twvpvf.U_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void translatingWallVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
// Remove the component of U normal to the wall in case the wall is not flat
vectorField n = patch().nf();
vectorField::operator=(U_ - n*(n & U_));
fixedValueFvPatchVectorField::updateCoeffs();
}
void translatingWallVelocityFvPatchVectorField::write(Ostream& os) const
{
fvPatchVectorField::write(os);
os.writeKeyword("U") << U_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchVectorField,
translatingWallVelocityFvPatchVectorField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

155
src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H Normal file
View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::translatingWallVelocityFvPatchVectorField
Description
Foam::translatingWallVelocityFvPatchVectorField
SourceFiles
translatingWallVelocityFvPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef translatingWallVelocityFvPatchVectorField_H
#define translatingWallVelocityFvPatchVectorField_H
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class translatingWallVelocityFvPatchField Declaration
\*---------------------------------------------------------------------------*/
class translatingWallVelocityFvPatchVectorField
:
public fixedValueFvPatchVectorField
{
// Private data
//- Origin of the rotation
vector U_;
public:
//- Runtime type information
TypeName("translatingWallVelocity");
// Constructors
//- Construct from patch and internal field
translatingWallVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&
);
//- Construct from patch, internal field and dictionary
translatingWallVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const dictionary&
);
//- Construct by mapping given a
// translatingWallVelocityFvPatchVectorField onto a new patch
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField&,
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField&
);
//- Construct and return a clone
virtual tmp<fvPatchVectorField> clone() const
{
return tmp<fvPatchVectorField>
(
new translatingWallVelocityFvPatchVectorField(*this)
);
}
//- Construct as copy setting internal field reference
translatingWallVelocityFvPatchVectorField
(
const translatingWallVelocityFvPatchVectorField&,
const DimensionedField<vector, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchVectorField> clone
(
const DimensionedField<vector, volMesh>& iF
) const
{
return tmp<fvPatchVectorField>
(
new translatingWallVelocityFvPatchVectorField(*this, iF)
);
}
// Member functions
// Access functions
//- Return the velocity
const vector& U() const
{
return U_;
}
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

4
src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
View File

@ -146,6 +146,10 @@ void Foam::singleCellFvMesh::agglomerateMesh
// From new patch face back to agglomeration
patchFaceMap_.setSize(oldPatches.size());
// From fine face to coarse face (or -1)
reverseFaceMap_.setSize(mesh.nFaces());
reverseFaceMap_.labelList::operator=(-1);
// Face counter
coarseI = 0;

2
src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
View File

@ -329,7 +329,7 @@ surfaceDisplacementPointPatchVectorField
surfacesDict_(dict.subDict("geometry")),
projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
projectDir_(dict.lookup("projectDirection")),
wedgePlane_(dict.lookupOrDefault(dict.lookup("wedgePlane"), -1)),
wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
{
if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)

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

@ -1858,7 +1858,7 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
// Get local mesh bounding box. Single box for now.
List<treeBoundBox> meshBb(1);
treeBoundBox& bb = meshBb[0];
bb = boundBox(mesh_.points(), false);
bb = treeBoundBox(mesh_.points());
bb = bb.extend(rndGen, 1E-4);
// Distribute all geometry (so refinementSurfaces and shellSurfaces)

1
src/meshTools/Make/files
View File

@ -50,6 +50,7 @@ indexedOctree/treeDataCell.C
indexedOctree/treeDataEdge.C
indexedOctree/treeDataFace.C
indexedOctree/treeDataPoint.C
indexedOctree/treeDataPrimitivePatchName.C
indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface

567
src/meshTools/indexedOctree/treeDataPrimitivePatch.C Normal file
View File

@ -0,0 +1,567 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "treeDataPrimitivePatch.H"
#include "indexedOctree.H"
#include "triangleFuncs.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::scalar
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
tolSqr = sqr(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::treeBoundBox
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
calcBb
(
const pointField& points,
const face& f
)
{
treeBoundBox bb(points[f[0]], points[f[0]]);
for (label fp = 1; fp < f.size(); fp++)
{
const point& p = points[f[fp]];
bb.min() = min(bb.min(), p);
bb.max() = max(bb.max(), p);
}
return bb;
}
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
update()
{
if (cacheBb_)
{
bbs_.setSize(patch_.size());
forAll(patch_, i)
{
bbs_[i] = calcBb(patch_.points(), patch_[i]);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
treeDataPrimitivePatch
(
const bool cacheBb,
const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
)
:
patch_(patch),
cacheBb_(cacheBb)
{
update();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::pointField
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
points() const
{
pointField cc(patch_.size());
forAll(patch_, i)
{
cc[i] = patch_[i].centre(patch_.points());
}
return cc;
}
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::label
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
getVolumeType
(
const indexedOctree
<
treeDataPrimitivePatch
<
Face,
FaceList,
PointField,
PointType
>
>& oc,
const point& sample
) const
{
// Need to determine whether sample is 'inside' or 'outside'
// Done by finding nearest face. This gives back a face which is
// guaranteed to contain nearest point. This point can be
// - in interior of face: compare to face normal
// - on edge of face: compare to edge normal
// - on point of face: compare to point normal
// Unfortunately the octree does not give us back the intersection point
// or where on the face it has hit so we have to recreate all that
// information.
// Find nearest face to sample
pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
if (info.index() == -1)
{
FatalErrorIn
(
"treeDataPrimitivePatch::getSampleType"
"(indexedOctree<treeDataPrimitivePatch>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
// Get actual intersection point on face
label faceI = info.index();
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
const pointField& points = patch_.localPoints();
const face& f = patch_.localFaces()[faceI];
// Retest to classify where on face info is. Note: could be improved. We
// already have point.
pointHit curHit = f.nearestPoint(sample, points);
const vector area = f.normal(points);
const point& curPt = curHit.rawPoint();
//
// 1] Check whether sample is above face
//
if (curHit.hit())
{
// Nearest point inside face. Compare to face normal.
if (debug & 2)
{
Pout<< " -> face hit:" << curPt
<< " comparing to face normal " << area << endl;
}
return indexedOctree<treeDataPrimitivePatch>::getSide
(
area,
sample - curPt
);
}
if (debug & 2)
{
Pout<< " -> face miss:" << curPt;
}
//
// 2] Check whether intersection is on one of the face vertices or
// face centre
//
const scalar typDimSqr = mag(area) + VSMALL;
forAll(f, fp)
{
if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point equals face vertex fp
// Calculate point normal (wrong: uses face normals instead of
// triangle normals)
return indexedOctree<treeDataPrimitivePatch>::getSide
(
patch_.pointNormals()[f[fp]],
sample - curPt
);
}
}
const point fc(f.centre(points));
if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
if (debug & 2)
{
Pout<< " -> centre hit:" << fc
<< " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
}
return indexedOctree<treeDataPrimitivePatch>::getSide
(
area,
sample - curPt
);
}
//
// 3] Get the 'real' edge the face intersection is on
//
const labelList& fEdges = patch_.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
label edgeI = fEdges[fEdgeI];
const edge& e = patch_.edges()[edgeI];
pointHit edgeHit = e.line(points).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point lies on edge e
// Calculate edge normal (wrong: uses face normals instead of
// triangle normals)
const labelList& eFaces = patch_.edgeFaces()[edgeI];
vector edgeNormal(vector::zero);
forAll(eFaces, i)
{
edgeNormal += patch_.faceNormal()[eFaces[i]];
}
if (debug & 2)
{
Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal:" << edgeNormal
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return indexedOctree<treeDataPrimitivePatch>::getSide
(
edgeNormal,
sample - curPt
);
}
}
//
// 4] Get the internal edge the face intersection is on
//
forAll(f, fp)
{
pointHit edgeHit = linePointRef
(
points[f[fp]],
fc
).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point lies on edge between two face triangles
// Calculate edge normal as average of the two triangle normals
vector e = points[f[fp]] - fc;
vector ePrev = points[f[f.rcIndex(fp)]] - fc;
vector eNext = points[f[f.fcIndex(fp)]] - fc;
vector nLeft = ePrev ^ e;
nLeft /= mag(nLeft) + VSMALL;
vector nRight = e ^ eNext;
nRight /= mag(nRight) + VSMALL;
if (debug & 2)
{
Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal "
<< 0.5*(nLeft + nRight)
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return indexedOctree<treeDataPrimitivePatch>::getSide
(
0.5*(nLeft + nRight),
sample - curPt
);
}
}
if (debug & 2)
{
Pout<< "Did not find sample " << sample
<< " anywhere related to nearest face " << faceI << endl
<< "Face:";
forAll(f, fp)
{
Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
<< endl;
}
}
// Can't determine status of sample with respect to nearest face.
// Either
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
}
// Check if any point on shape is inside cubeBb.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
bool
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
overlaps
(
const label index,
const treeBoundBox& cubeBb
) const
{
// 1. Quick rejection: bb does not intersect face bb at all
if (cacheBb_)
{
if (!cubeBb.overlaps(bbs_[index]))
{
return false;
}
}
else
{
if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
{
return false;
}
}
// 2. Check if one or more face points inside
const pointField& points = patch_.points();
const face& f = patch_[index];
forAll(f, fp)
{
if (cubeBb.contains(points[f[fp]]))
{
return true;
}
}
// 3. Difficult case: all points are outside but connecting edges might
// go through cube. Use triangle-bounding box intersection.
const point fc = f.centre(points);
forAll(f, fp)
{
bool triIntersects = triangleFuncs::intersectBb
(
points[f[fp]],
points[f[f.fcIndex(fp)]],
fc,
cubeBb
);
if (triIntersects)
{
return true;
}
}
return false;
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
void
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
findNearest
(
const labelList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const
{
const pointField& points = patch_.points();
forAll(indices, i)
{
label index = indices[i];
const face& f = patch_[index];
pointHit nearHit = f.nearestPoint(sample, points);
scalar distSqr = sqr(nearHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = nearHit.rawPoint();
}
}
}
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
bool
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
// Do quick rejection test
if (cacheBb_)
{
const treeBoundBox& faceBb = bbs_[index];
if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
{
// start and end in same block outside of faceBb.
return false;
}
}
const pointField& points = patch_.points();
const face& f = patch_[index];
const point fc = f.centre(points);
const vector dir(end - start);
pointHit inter = patch_[index].intersection
(
start,
dir,
fc,
points,
intersection::HALF_RAY
);
if (inter.hit() && inter.distance() <= 1)
{
// Note: no extra test on whether intersection is in front of us
// since using half_ray
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
// ************************************************************************* //

210
src/meshTools/indexedOctree/treeDataPrimitivePatch.H Normal file
View File

@ -0,0 +1,210 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::treeDataPrimitivePatch
Description
Encapsulation of data needed to search on PrimitivePatches
SourceFiles
treeDataPrimitivePatch.C
\*---------------------------------------------------------------------------*/
#ifndef treeDataPrimitivePatch_H
#define treeDataPrimitivePatch_H
#include "PrimitivePatch.H"
//#include "indexedOctree.H"
#include "treeBoundBoxList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
Class treeDataPrimitivePatchName Declaration
\*---------------------------------------------------------------------------*/
TemplateName(treeDataPrimitivePatch);
/*---------------------------------------------------------------------------*\
Class treeDataPrimitivePatch Declaration
\*---------------------------------------------------------------------------*/
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType=point
>
class treeDataPrimitivePatch
:
public treeDataPrimitivePatchName
{
// Static data
//- tolerance on linear dimensions
static scalar tolSqr;
// Private data
//- Underlying geometry
const PrimitivePatch<Face, FaceList, PointField, PointType>& patch_;
//- Whether to precalculate and store face bounding box
const bool cacheBb_;
//- face bounding boxes (valid only if cacheBb_)
treeBoundBoxList bbs_;
// Private Member Functions
//- Calculate face bounding box
static treeBoundBox calcBb(const pointField&, const face&);
//- Initialise all member data
void update();
public:
// Constructors
//- Construct from patch.
treeDataPrimitivePatch
(
const bool cacheBb,
const PrimitivePatch<Face, FaceList, PointField, PointType>&
);
// Member Functions
// Access
label size() const
{
return patch_.size();
}
//- Get representative point cloud for all shapes inside
// (one point per shape)
pointField points() const;
// Search
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
(
const indexedOctree
<
treeDataPrimitivePatch
<
Face,
FaceList,
PointField,
PointType
>
>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataPrimitivePatch::findNearest"
"(const labelList&, const linePointRef&, ..)"
);
}
//- Calculate intersection of shape with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "treeDataPrimitivePatch.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

33
src/meshTools/indexedOctree/treeDataPrimitivePatchName.C Normal file
View File

@ -0,0 +1,33 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "treeDataPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::treeDataPrimitivePatchName, 0);
// ************************************************************************* //

2
src/meshTools/indexedOctree/treeDataTriSurface.C
View File

@ -276,7 +276,7 @@ bool Foam::treeDataTriSurface::overlaps
const point& p1 = points[f[1]];
const point& p2 = points[f[2]];
boundBox triBb(p0, p0);
treeBoundBox triBb(p0, p0);
triBb.min() = min(triBb.min(), p1);
triBb.min() = min(triBb.min(), p2);

2
src/meshTools/octree/octreeDataPoint.H
View File

@ -69,7 +69,7 @@ public:
// Constructors
//- Construct from components. Holds reference to points!
octreeDataPoint(const pointField&);
explicit octreeDataPoint(const pointField&);
// Member Functions

4
src/meshTools/octree/treeBoundBox.H
View File

@ -170,11 +170,11 @@ public:
inline treeBoundBox(const point& min, const point& max);
//- Construct from components
inline treeBoundBox(const boundBox& bb);
explicit inline treeBoundBox(const boundBox& bb);
//- Construct as the bounding box of the given pointField.
// Local processor domain only (no reduce as in boundBox)
treeBoundBox(const UList<point>&);
explicit treeBoundBox(const UList<point>&);
//- Construct as subset of points
// Local processor domain only (no reduce as in boundBox)

2
src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
View File

@ -957,7 +957,7 @@ bool Foam::distributedTriSurfaceMesh::overlaps
{
const treeBoundBox& bb = bbs[bbI];
boundBox triBb(p0, p0);
treeBoundBox triBb(p0, p0);
triBb.min() = min(triBb.min(), p1);
triBb.min() = min(triBb.min(), p2);

12
src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
View File

@ -757,9 +757,13 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
) const
{
// Build tree out of all samples.
//Note: cannot be done one the fly - gcc4.4 compiler bug.
treeBoundBox bb(samples);
octree<octreeDataPoint> ppTree
(
treeBoundBox(samples), // overall search domain
bb, // overall search domain
octreeDataPoint(samples), // all information needed to do checks
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
@ -857,10 +861,12 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
scalar maxSearch = max(maxDist);
vector span(maxSearch, maxSearch, maxSearch);
// octree.shapes holds reference!
//Note: cannot be done one the fly - gcc4.4 compiler bug.
treeBoundBox bb(samples);
octree<octreeDataPoint> ppTree
(
treeBoundBox(samples), // overall search domain
bb, // overall search domain
octreeDataPoint(samples), // all information needed to do checks
1, // min levels
20.0, // maximum ratio of cubes v.s. cells

2
src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
View File

@ -32,6 +32,8 @@ License
#include "error.H"
#include "IStringStream.H"
// For EOF only
#include <cstdio>
// flex input buffer size
int Foam::chemkinReader::yyBufSize = YY_BUF_SIZE;

8
src/triSurface/triSurface/stitchTriangles.C
View File

@ -46,10 +46,7 @@ bool triSurface::stitchTriangles
pointField newPoints;
bool hasMerged = mergePoints(rawPoints, tol, verbose, pointMap, newPoints);
pointField& ps = storedPoints();
// Set the coordinates to the merged ones
ps = newPoints;
if (hasMerged)
{
@ -59,6 +56,11 @@ bool triSurface::stitchTriangles
<< " points down to " << newPoints.size() << endl;
}
pointField& ps = storedPoints();
// Set the coordinates to the merged ones
ps = newPoints;
// Reset the triangle point labels to the unique points array
label newTriangleI = 0;
forAll(*this, i)

2
src/triSurface/triSurface/triSurface.C
View File

@ -814,7 +814,7 @@ void Foam::triSurface::scalePoints(const scalar scaleFactor)
void Foam::triSurface::cleanup(const bool verbose)
{
// Merge points (already done for STL, TRI)
stitchTriangles(pointField(points()), SMALL, verbose);
stitchTriangles(points(), SMALL, verbose);
// Merging points might have changed geometric factors
clearOut();

2
src/turbulenceModels/compressible/RAS/Make/files
View File

@ -38,8 +38,8 @@ derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFv
derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.C
derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C
derivedFvPatchFields/turbulentTemperatureCoupledBaffle/regionProperties.C
backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
LIB = $(FOAM_LIBBIN)/libcompressibleRASModels

132
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.C
View File

@ -33,16 +33,9 @@ License
#include "basicThermo.H"
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool turbulentTemperatureCoupledBaffleFvPatchScalarField::interfaceOwner
bool Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::interfaceOwner
(
const polyMesh& nbrRegion,
const polyPatch& nbrPatch
@ -103,6 +96,7 @@ bool turbulentTemperatureCoupledBaffleFvPatchScalarField::interfaceOwner
}
nbrIndex = props.fluidRegionNames().size() + i;
}
return myIndex < nbrIndex;
}
}
@ -110,25 +104,20 @@ bool turbulentTemperatureCoupledBaffleFvPatchScalarField::interfaceOwner
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
turbulentTemperatureCoupledBaffleFvPatchScalarField::
Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::
turbulentTemperatureCoupledBaffleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(p, iF),
fixedValueFvPatchScalarField(p, iF),
neighbourFieldName_("undefined-neighbourFieldName"),
KName_("undefined-K")
{
this->refValue() = 0.0;
this->refGrad() = 0.0;
this->valueFraction() = 1.0;
this->fixesValue_ = true;
}
{}
turbulentTemperatureCoupledBaffleFvPatchScalarField::
Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::
turbulentTemperatureCoupledBaffleFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleFvPatchScalarField& ptf,
@ -137,14 +126,13 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
mixedFvPatchScalarField(ptf, p, iF, mapper),
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
neighbourFieldName_(ptf.neighbourFieldName_),
KName_(ptf.KName_),
fixesValue_(ptf.fixesValue_)
KName_(ptf.KName_)
{}
turbulentTemperatureCoupledBaffleFvPatchScalarField::
Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::
turbulentTemperatureCoupledBaffleFvPatchScalarField
(
const fvPatch& p,
@ -152,7 +140,7 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField
const dictionary& dict
)
:
mixedFvPatchScalarField(p, iF),
fixedValueFvPatchScalarField(p, iF, dict),
neighbourFieldName_(dict.lookup("neighbourFieldName")),
KName_(dict.lookup("K"))
{
@ -174,46 +162,26 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField
<< " in file " << dimensionedInternalField().objectPath()
<< exit(FatalError);
}
fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
if (dict.found("refValue"))
{
// Full restart
refValue() = scalarField("refValue", dict, p.size());
refGrad() = scalarField("refGradient", dict, p.size());
valueFraction() = scalarField("valueFraction", dict, p.size());
fixesValue_ = readBool(dict.lookup("fixesValue"));
}
else
{
// Start from user entered data. Assume fixedValue.
refValue() = *this;
refGrad() = 0.0;
valueFraction() = 1.0;
fixesValue_ = true;
}
}
turbulentTemperatureCoupledBaffleFvPatchScalarField::
Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::
turbulentTemperatureCoupledBaffleFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleFvPatchScalarField& wtcsf,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(wtcsf, iF),
fixedValueFvPatchScalarField(wtcsf, iF),
neighbourFieldName_(wtcsf.neighbourFieldName_),
KName_(wtcsf.KName_),
fixesValue_(wtcsf.fixesValue_)
KName_(wtcsf.KName_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<scalarField>
turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const
Foam::tmp<Foam::scalarField>
Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const
{
const fvMesh& mesh = patch().boundaryMesh().mesh();
@ -260,7 +228,7 @@ turbulentTemperatureCoupledBaffleFvPatchScalarField::K() const
}
void turbulentTemperatureCoupledBaffleFvPatchScalarField::updateCoeffs()
void Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::updateCoeffs()
{
if (updated())
{
@ -353,18 +321,13 @@ void turbulentTemperatureCoupledBaffleFvPatchScalarField::updateCoeffs()
).fvPatchScalarField::operator=(Twall);
}
// Switch between fixed value (of harmonic avg) or gradient
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFixed = 0;
// Like snGrad but bypass switching on refValue/refGrad.
tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
if (debug)
{
scalar Q = gSum(K()*patch().magSf()*normalGradient());
//tmp<scalarField> normalGradient =
// (*this-intFld())
// * patch().deltaCoeffs();
scalar Q = gSum(K()*patch().magSf()*snGrad());
Info<< patch().boundaryMesh().mesh().name() << ':'
<< patch().name() << ':'
@ -380,74 +343,33 @@ void turbulentTemperatureCoupledBaffleFvPatchScalarField::updateCoeffs()
<< endl;
}
forAll(*this, i)
{
// if outgoing flux use fixed value.
if (normalGradient()[i] < 0.0)
{
this->refValue()[i] = operator[](i);
this->refGrad()[i] = 0.0; // not used
this->valueFraction()[i] = 1.0;
nFixed++;
}
else
{
this->refValue()[i] = 0.0; // not used
this->refGrad()[i] = normalGradient()[i];
this->valueFraction()[i] = 0.0;
}
}
reduce(nFixed, sumOp<label>());
fixesValue_ = (nFixed > 0);
if (debug)
{
label nTotSize = returnReduce(this->size(), sumOp<label>());
Info<< patch().boundaryMesh().mesh().name() << ':'
<< patch().name() << ':'
<< this->dimensionedInternalField().name() << " -> "
<< nbrMesh.name() << ':'
<< nbrPatch.name() << ':'
<< this->dimensionedInternalField().name() << " :"
<< " patch:" << patch().name()
<< " out of:" << nTotSize
<< " fixedBC:" << nFixed
<< " gradient:" << nTotSize-nFixed << endl;
}
mixedFvPatchScalarField::updateCoeffs();
fixedValueFvPatchScalarField::updateCoeffs();
}
void turbulentTemperatureCoupledBaffleFvPatchScalarField::write
void Foam::turbulentTemperatureCoupledBaffleFvPatchScalarField::write
(
Ostream& os
) const
{
mixedFvPatchScalarField::write(os);
fixedValueFvPatchScalarField::write(os);
os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
<< token::END_STATEMENT << nl;
os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
turbulentTemperatureCoupledBaffleFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace compressible
} // End namespace Foam
// ************************************************************************* //

32
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
View File

@ -23,16 +23,15 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::compressible::turbulentTemperatureCoupledBaffleFvPatchScalarField
turbulentTemperatureCoupledBaffleFvPatchScalarField
Description
Mixed boundary condition for temperature, to be used for heat-transfer
on back-to-back baffles.
Harmonic fixed value boundary condition for temperature, to be used
for heat-transfer on back-to-back baffles.
If my temperature is T1, neighbour is T2:
If my temperature is T1, heat conductivity K1 and neighbour is T2,K2
T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
both sides get fixedValue (K1/dx1*T1 + K2/dx2*T2)/(K1/dx1+K2/dx2)
Example usage:
myInterfacePatchName
@ -65,16 +64,12 @@ SourceFiles
#ifndef turbulentTemperatureCoupledBaffleFvPatchScalarField_H
#define turbulentTemperatureCoupledBaffleFvPatchScalarField_H
//#include "fvPatchFields.H"
#include "mixedFvPatchFields.H"
//#include "fvPatch.H"
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
/*---------------------------------------------------------------------------*\
Class turbulentTemperatureCoupledBaffleFvPatchScalarField Declaration
@ -82,7 +77,7 @@ namespace compressible
class turbulentTemperatureCoupledBaffleFvPatchScalarField
:
public mixedFvPatchScalarField
public fixedValueFvPatchScalarField
{
// Private data
@ -92,8 +87,6 @@ class turbulentTemperatureCoupledBaffleFvPatchScalarField
//- Name of thermal conductivity field
const word KName_;
bool fixesValue_;
// Private Member Functions
@ -104,7 +97,7 @@ class turbulentTemperatureCoupledBaffleFvPatchScalarField
public:
//- Runtime type information
TypeName("compressible::turbulentTemperatureCoupledBaffle");
TypeName("turbulentTemperatureCoupledBaffle");
// Constructors
@ -172,14 +165,6 @@ public:
//- Get corresponding K field
tmp<scalarField> K() const;
//- Return true if this patch field fixes a value.
// Needed to check if a level has to be specified while solving
// Poissons equations.
virtual bool fixesValue() const
{
return fixesValue_;
}
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
@ -190,7 +175,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace compressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

458
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C Normal file
View File

@ -0,0 +1,458 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "directMappedPatchBase.H"
#include "regionProperties.H"
#include "basicThermo.H"
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::interfaceOwner
(
const polyMesh& nbrRegion,
const polyPatch& nbrPatch
) const
{
const fvMesh& myRegion = patch().boundaryMesh().mesh();
if (nbrRegion.name() == myRegion.name())
{
return patch().index() < nbrPatch.index();
}
else
{
const regionProperties& props =
myRegion.objectRegistry::parent().lookupObject<regionProperties>
(
"regionProperties"
);
label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
if (myIndex == -1)
{
label i = findIndex(props.solidRegionNames(), myRegion.name());
if (i == -1)
{
FatalErrorIn
(
"turbulentTemperatureCoupledBaffleMixedFvPatchScalarField"
"::interfaceOwner(const polyMesh&"
", const polyPatch&)const"
) << "Cannot find region " << myRegion.name()
<< " neither in fluids " << props.fluidRegionNames()
<< " nor in solids " << props.solidRegionNames()
<< exit(FatalError);
}
myIndex = props.fluidRegionNames().size() + i;
}
label nbrIndex = findIndex
(
props.fluidRegionNames(),
nbrRegion.name()
);
if (nbrIndex == -1)
{
label i = findIndex(props.solidRegionNames(), nbrRegion.name());
if (i == -1)
{
FatalErrorIn
(
"coupleManager::interfaceOwner"
"(const polyMesh&, const polyPatch&) const"
) << "Cannot find region " << nbrRegion.name()
<< " neither in fluids " << props.fluidRegionNames()
<< " nor in solids " << props.solidRegionNames()
<< exit(FatalError);
}
nbrIndex = props.fluidRegionNames().size() + i;
}
return myIndex < nbrIndex;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(p, iF),
neighbourFieldName_("undefined-neighbourFieldName"),
KName_("undefined-K")
{
this->refValue() = 0.0;
this->refGrad() = 0.0;
this->valueFraction() = 1.0;
this->fixesValue_ = true;
}
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
mixedFvPatchScalarField(ptf, p, iF, mapper),
neighbourFieldName_(ptf.neighbourFieldName_),
KName_(ptf.KName_),
fixesValue_(ptf.fixesValue_)
{}
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
mixedFvPatchScalarField(p, iF),
neighbourFieldName_(dict.lookup("neighbourFieldName")),
KName_(dict.lookup("K"))
{
if (!isA<directMappedPatchBase>(this->patch().patch()))
{
FatalErrorIn
(
"turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::"
"turbulentTemperatureCoupledBaffleMixedFvPatchScalarField\n"
"(\n"
" const fvPatch& p,\n"
" const DimensionedField<scalar, volMesh>& iF,\n"
" const dictionary& dict\n"
")\n"
) << "\n patch type '" << p.type()
<< "' not type '" << directMappedPatchBase::typeName << "'"
<< "\n for patch " << p.name()
<< " of field " << dimensionedInternalField().name()
<< " in file " << dimensionedInternalField().objectPath()
<< exit(FatalError);
}
fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
if (dict.found("refValue"))
{
// Full restart
refValue() = scalarField("refValue", dict, p.size());
refGrad() = scalarField("refGradient", dict, p.size());
valueFraction() = scalarField("valueFraction", dict, p.size());
fixesValue_ = readBool(dict.lookup("fixesValue"));
}
else
{
// Start from user entered data. Assume fixedValue.
refValue() = *this;
refGrad() = 0.0;
valueFraction() = 1.0;
fixesValue_ = true;
}
}
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField& wtcsf,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(wtcsf, iF),
neighbourFieldName_(wtcsf.neighbourFieldName_),
KName_(wtcsf.KName_),
fixesValue_(wtcsf.fixesValue_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<scalarField>
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::K() const
{
const fvMesh& mesh = patch().boundaryMesh().mesh();
if (KName_ == "none")
{
const compressible::RASModel& model =
db().lookupObject<compressible::RASModel>("RASProperties");
tmp<volScalarField> talpha = model.alphaEff();
const basicThermo& thermo =
db().lookupObject<basicThermo>("thermophysicalProperties");
return
talpha().boundaryField()[patch().index()]
*thermo.Cp()().boundaryField()[patch().index()];
}
else if (mesh.objectRegistry::foundObject<volScalarField>(KName_))
{
return patch().lookupPatchField<volScalarField, scalar>(KName_);
}
else if (mesh.objectRegistry::foundObject<volSymmTensorField>(KName_))
{
const symmTensorField& KWall =
patch().lookupPatchField<volSymmTensorField, scalar>(KName_);
vectorField n = patch().nf();
return n & KWall & n;
}
else
{
FatalErrorIn
(
"turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::K()"
" const"
) << "Did not find field " << KName_
<< " on mesh " << mesh.name() << " patch " << patch().name()
<< endl
<< "Please set 'K' to 'none', a valid volScalarField"
<< " or a valid volSymmTensorField." << exit(FatalError);
return scalarField(0);
}
}
void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
// Get the coupling information from the directMappedPatchBase
const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
(
patch().patch()
);
const polyMesh& nbrMesh = mpp.sampleMesh();
const fvPatch& nbrPatch = refCast<const fvMesh>
(
nbrMesh
).boundary()[mpp.samplePolyPatch().index()];
// Force recalculation of mapping and schedule
const mapDistribute& distMap = mpp.map();
(void)distMap.schedule();
tmp<scalarField> intFld = patchInternalField();
if (interfaceOwner(nbrMesh, nbrPatch.patch()))
{
// Note: other side information could be cached - it only needs
// to be updated the first time round the iteration (i.e. when
// switching regions) but unfortunately we don't have this information.
// Calculate the temperature by harmonic averaging
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&
nbrField =
refCast
<
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
>
(
nbrPatch.lookupPatchField<volScalarField, scalar>
(
neighbourFieldName_
)
);
// Swap to obtain full local values of neighbour internal field
scalarField nbrIntFld = nbrField.patchInternalField();
mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(), // what to send
distMap.constructMap(), // what to receive
nbrIntFld
);
// Swap to obtain full local values of neighbour K*delta
scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(), // what to send
distMap.constructMap(), // what to receive
nbrKDelta
);
tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
// Calculate common wall temperature. Reuse *this to store common value.
scalarField Twall
(
(myKDelta()*intFld() + nbrKDelta*nbrIntFld)
/ (myKDelta() + nbrKDelta)
);
// Assign to me
fvPatchScalarField::operator=(Twall);
// Distribute back and assign to neighbour
mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
nbrField.size(),
distMap.constructMap(), // reverse : what to send
distMap.subMap(),
Twall
);
const_cast<turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&>
(
nbrField
).fvPatchScalarField::operator=(Twall);
}
// Switch between fixed value (of harmonic avg) or gradient
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFixed = 0;
// Like snGrad but bypass switching on refValue/refGrad.
tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
if (debug)
{
scalar Q = gSum(K()*patch().magSf()*normalGradient());
Info<< patch().boundaryMesh().mesh().name() << ':'
<< patch().name() << ':'
<< this->dimensionedInternalField().name() << " -> "
<< nbrMesh.name() << ':'
<< nbrPatch.name() << ':'
<< this->dimensionedInternalField().name() << " :"
<< " heatFlux:" << Q
<< " walltemperature "
<< " min:" << gMin(*this)
<< " max:" << gMax(*this)
<< " avg:" << gAverage(*this)
<< endl;
}
forAll(*this, i)
{
// if outgoing flux use fixed value.
if (normalGradient()[i] < 0.0)
{
this->refValue()[i] = operator[](i);
this->refGrad()[i] = 0.0; // not used
this->valueFraction()[i] = 1.0;
nFixed++;
}
else
{
this->refValue()[i] = 0.0; // not used
this->refGrad()[i] = normalGradient()[i];
this->valueFraction()[i] = 0.0;
}
}
reduce(nFixed, sumOp<label>());
fixesValue_ = (nFixed > 0);
if (debug)
{
label nTotSize = returnReduce(this->size(), sumOp<label>());
Info<< patch().boundaryMesh().mesh().name() << ':'
<< patch().name() << ':'
<< this->dimensionedInternalField().name() << " -> "
<< nbrMesh.name() << ':'
<< nbrPatch.name() << ':'
<< this->dimensionedInternalField().name() << " :"
<< " patch:" << patch().name()
<< " out of:" << nTotSize
<< " fixedBC:" << nFixed
<< " gradient:" << nTotSize-nFixed << endl;
}
mixedFvPatchScalarField::updateCoeffs();
}
void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::write
(
Ostream& os
) const
{
mixedFvPatchScalarField::write(os);
os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
<< token::END_STATEMENT << nl;
os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace compressible
} // End namespace Foam
// ************************************************************************* //

204
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H Normal file
View File

@ -0,0 +1,204 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::compressible::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
Description
Mixed boundary condition for temperature, to be used for heat-transfer
on back-to-back baffles.
If my temperature is T1, neighbour is T2:
T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
Example usage:
myInterfacePatchName
{
type turbulentTemperatureCoupledBaffleMixed;
neighbourFieldName T;
K K; // or none
value uniform 300;
}
Needs to be on underlying directMapped(Wall)FvPatch.
Note: if K is "none" looks up RASModel and basicThermo, otherwise expects
the solver to calculate a 'K' field.
Note: runs in parallel with arbitrary decomposition. Uses directMapped
functionality to calculate exchange.
Note: lags interface data so both sides use same data.
- problem: schedule to calculate average would interfere
with standard processor swaps.
- so: updateCoeffs sets both to same Twall. Only need to do
this for last outer iteration but don't have access to this.
SourceFiles
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef turbulentTemperatureCoupledBaffleMixedFvPatchScalarField_H
#define turbulentTemperatureCoupledBaffleMixedFvPatchScalarField_H
//#include "fvPatchFields.H"
#include "mixedFvPatchFields.H"
//#include "fvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
/*---------------------------------------------------------------------------*\
Class turbulentTemperatureCoupledBaffleMixedFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
:
public mixedFvPatchScalarField
{
// Private data
//- Name of field on the neighbour region
const word neighbourFieldName_;
//- Name of thermal conductivity field
const word KName_;
bool fixesValue_;
// Private Member Functions
//- Am I or neighbour owner of interface
bool interfaceOwner(const polyMesh&, const polyPatch&) const;
public:
//- Runtime type information
TypeName("compressible::turbulentTemperatureCoupledBaffleMixed");
// Constructors
//- Construct from patch and internal field
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// turbulentTemperatureCoupledBaffleMixedFvPatchScalarField onto a
// new patch
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField>
(
new turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
*this
)
);
}
//- Construct as copy setting internal field reference
turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
(
*this,
iF
)
);
}
// Member functions
//- Get corresponding K field
tmp<scalarField> K() const;
//- Return true if this patch field fixes a value.
// Needed to check if a level has to be specified while solving
// Poissons equations.
virtual bool fixesValue() const
{
return fixesValue_;
}
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace compressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
View File

@ -239,7 +239,7 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
// Molecular Prandtl number
scalar Pr = muw[faceI]/alphaw[faceI];
// Molecular-to-turbulenbt Prandtl number ratio
// Molecular-to-turbulent Prandtl number ratio
scalar Prat = Pr/Prt_;
// Thermal sublayer thickness