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' of ssh://hunt/~OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2008-09-10 13:58:03 +01:00
parent 913c912e87 58fa3c4074
commit f6d13e1a08
70 changed files with 2754 additions and 562 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/Lagrangian/kinematicParcelFoam/Make/files
View File

@ -1,3 +1,3 @@
kinematicParcelFoam.C
EXE = $(FOAM_USER_APPBIN)/kinematicParcelFoam
EXE = $(FOAM_APPBIN)/kinematicParcelFoam

3
applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
View File

@ -26,7 +26,8 @@ Application
kinematicParcelFoam
Description
Transient solver a single kinematicCloud.
Transient solver for a single kinematicCloud. Uses precalculated velocity
field to evolve a cloud.
\*---------------------------------------------------------------------------*/

2
applications/solvers/multiphase/interDyMFoam/Make/options
View File

@ -1,4 +1,6 @@
EXE_INC = \
-I../rasInterFoam \
-I../interFoam \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

34
applications/solvers/multiphase/interDyMFoam/UEqn.H
View File

@ -1,34 +0,0 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
UEqn.relax();
if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
)*mesh.magSf()
)
);
}

35
applications/solvers/multiphase/interDyMFoam/gammaEqn.H
View File

@ -1,35 +0,0 @@
{
word gammaScheme("div(phi,gamma)");
word gammarScheme("div(phirb,gamma)");
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cGamma()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
{
surfaceScalarField phiGamma =
fvc::flux
(
phi,
gamma,
gammaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
gamma,
gammarScheme
);
MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;
}
Info<< "Liquid phase volume fraction = "
<< gamma.weightedAverage(mesh.V()).value()
<< " Min(gamma) = " << min(gamma).value()
<< " Max(gamma) = " << max(gamma).value()
<< endl;
}

35
applications/solvers/multiphase/interDyMFoam/gammaEqnSubCycle.H
View File

@ -1,35 +0,0 @@
label nGammaCorr
(
readLabel(piso.lookup("nGammaCorr"))
);
label nGammaSubCycles
(
readLabel(piso.lookup("nGammaSubCycles"))
);
if (nGammaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
for
(
subCycle<volScalarField> gammaSubCycle(gamma, nGammaSubCycles);
!(++gammaSubCycle).end();
)
{
# include "gammaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
# include "gammaEqn.H"
}
interface.correct();
rho == gamma*rho1 + (scalar(1) - gamma)*rho2;

1
applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
View File

@ -41,7 +41,6 @@ Description
#include "twoPhaseMixture.H"
#include "incompressible/RASModel/RASModel.H"
#include "probes.H"
#include "EulerDdtScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

2
applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
View File

@ -697,7 +697,7 @@ endOfSection {space}")"{space}
/* ------ Ignore remaining space and \n s. ------ */
<*>{some_space}|\n {
<*>{some_space}|\n|\r {
}

4
applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C
View File

@ -28,7 +28,7 @@ License
#include "OFstream.H"
#include "floatScalar.H"
#include "writeFuns.H"
#include "emptyFvPatchFields.H"
#include "emptyFvsPatchFields.H"
#include "fvsPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -100,7 +100,7 @@ void writeSurfFields
const fvPatch& pp = mesh.boundary()[patchI];
if (isA<emptyFvPatchVectorField>(pf))
if (isA<emptyFvsPatchVectorField>(pf))
{
// Note: loop over polypatch size, not fvpatch size.
forAll(pp.patch(), i)

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

@ -248,7 +248,7 @@ int main(int argc, char *argv[])
indexedOctree<treeDataTriSurface> selectTree
(
treeDataTriSurface(selectSurf),
bb.extend(rndGen, 1E-3), // slightly randomize bb
bb.extend(rndGen, 1E-4), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity

4
etc/apps/paraview3/bashrc
View File

@ -45,9 +45,9 @@ export ParaView_INST_DIR=$WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
if [ "$PYTHONPATH" ]; then
export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping
export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
else
export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping
export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
fi

4
etc/apps/paraview3/cshrc
View File

@ -45,9 +45,9 @@ setenv ParaView_INST_DIR $WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
if ($?PYTHONPATH) then
setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping
setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
else
setenv PYTHONPATH $ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping
setenv PYTHONPATH $ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
endif
if ( -r $ParaView_INST_DIR ) then

90
src/OpenFOAM/containers/HashTables/HashSet/HashSet.C Normal file
View File

@ -0,0 +1,90 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
\*---------------------------------------------------------------------------*/
#ifndef HashSet_C
#define HashSet_C
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Key, class Hash>
bool HashSet<Key, Hash>::operator==(const HashSet<Key, Hash>& ht) const
{
const HashTable<empty, Key, Hash>& a = *this;
// Are all my elements in ht?
for
(
typename HashTable<empty, Key, Hash>::const_iterator iter = a.begin();
iter != a.end();
++iter
)
{
if (!ht.found(iter.key()))
{
return false;
}
}
// Are all ht elements in me?
for
(
typename HashTable<empty, Key, Hash>::const_iterator iter = ht.begin();
iter != ht.end();
++iter
)
{
if (!found(iter.key()))
{
return false;
}
}
return true;
}
template<class Key, class Hash>
bool HashSet<Key, Hash>::operator!=(const HashSet<Key, Hash>& ht) const
{
return !(operator==(ht));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

17
src/OpenFOAM/containers/HashTables/HashSet/HashSet.H
View File

@ -100,6 +100,17 @@ public:
{
return HashTable<empty, Key, Hash>::insert(key, empty());
}
// Member Operators
//- Equality. Two hashtables are equal if all contents of first are
// also in second and vice versa. So does not depend on table size or
// order!
bool operator==(const HashSet<Key, Hash>&) const;
//- The opposite of the equality operation.
bool operator!=(const HashSet<Key, Hash>&) const;
};
@ -112,6 +123,12 @@ typedef HashSet<> wordHashSet;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "HashSet.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

21
src/OpenFOAM/containers/Lists/List/List.C
View File

@ -321,7 +321,7 @@ List<T>::List(const BiIndirectList<T>& idl)
template<class T>
List<T>::~List()
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
}
@ -367,9 +367,8 @@ void List<T>::setSize(const label newSize)
register T* av = &nv[i];
while (i--) *--av = *--vv;
}
delete[] this->v_;
}
if (this->v_) delete[] this->v_;
this->size_ = newSize;
this->v_ = nv;
@ -400,7 +399,7 @@ void List<T>::setSize(const label newSize, const T& a)
template<class T>
void List<T>::clear()
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->size_ = 0;
this->v_ = 0;
}
@ -411,7 +410,7 @@ void List<T>::clear()
template<class T>
void List<T>::transfer(List<T>& a)
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->size_ = a.size_;
this->v_ = a.v_;
@ -457,7 +456,8 @@ void List<T>::operator=(const UList<T>& a)
{
if (a.size_ != this->size_)
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->v_ = 0;
this->size_ = a.size_;
if (this->size_) this->v_ = new T[this->size_];
}
@ -503,7 +503,8 @@ void List<T>::operator=(const SLList<T>& sll)
{
if (sll.size() != this->size_)
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->v_ = 0;
this->size_ = sll.size();
if (this->size_) this->v_ = new T[this->size_];
}
@ -530,7 +531,8 @@ void List<T>::operator=(const IndirectList<T>& idl)
{
if (idl.size() != this->size_)
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->v_ = 0;
this->size_ = idl.size();
if (this->size_) this->v_ = new T[this->size_];
}
@ -551,7 +553,8 @@ void List<T>::operator=(const BiIndirectList<T>& idl)
{
if (idl.size() != this->size_)
{
if (this->size_) delete[] this->v_;
if (this->v_) delete[] this->v_;
this->v_ = 0;
this->size_ = idl.size();
if (this->size_) this->v_ = new T[this->size_];
}

7
src/OpenFOAM/containers/Lists/List/List.H
View File

@ -137,6 +137,9 @@ public:
//- Return a null List
static const List<T>& null();
//- Return the number of elements in the UList.
inline label size() const;
// Edit
@ -156,6 +159,10 @@ public:
//- Return subscript-checked element of UList.
inline T& newElmt(const label);
//- Override size to be inconsistent with allocated storage.
// Use with care.
inline label& size();
// Member operators

14
src/OpenFOAM/containers/Lists/List/ListI.H
View File

@ -52,6 +52,20 @@ inline T& Foam::List<T>::newElmt(const label i)
}
template<class T>
inline Foam::label Foam::List<T>::size() const
{
return UList<T>::size_;
}
template<class T>
inline Foam::label& Foam::List<T>::size()
{
return UList<T>::size_;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class T>

7
src/OpenFOAM/db/Time/TimeIO.C
View File

@ -86,13 +86,6 @@ void Foam::Time::readDict()
purgeWrite_ = 0;
}
if (writeControl_ != wcTimeStep && purgeWrite_ > 0)
{
FatalIOErrorIn("Time::readDict()", controlDict_)
<< "writeControl must be set to timeStep for purgeWrite "
<< exit(FatalIOError);
}
}
if (controlDict_.found("timeFormat"))

42
src/OpenFOAM/meshes/boundBox/boundBox.C
View File

@ -87,9 +87,47 @@ boundBox::boundBox(Istream& is)
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const boundBox& b)
Ostream& operator<<(Ostream& os, const boundBox& bb)
{
return os << b.min() << token::SPACE << b.max();
if (os.format() == IOstream::ASCII)
{
os << bb.min_ << token::SPACE << bb.max_;
}
else
{
os.write
(
reinterpret_cast<const char*>(&bb.min_),
sizeof(boundBox)
);
}
// Check state of Ostream
os.check("Ostream& operator<<(Ostream&, const boundBox&)");
return os;
}
Istream& operator>>(Istream& is, boundBox& bb)
{
if (is.format() == IOstream::ASCII)
{
return is >> bb.min_ >> bb.max_;
}
else
{
is.read
(
reinterpret_cast<char*>(&bb.min_),
sizeof(boundBox)
);
}
// Check state of Istream
is.check("Istream& operator>>(Istream&, boundBox&)");
return is;
}

23
src/OpenFOAM/meshes/boundBox/boundBox.H
View File

@ -153,12 +153,31 @@ public:
}
// Ostream operator
// Friend Operators
friend Ostream& operator<<(Ostream& os, const boundBox& b);
friend bool operator==(const boundBox& a, const boundBox& b)
{
return (a.min_ == b.min_) && (a.max_ == b.max_);
}
friend bool operator!=(const boundBox& a, const boundBox& b)
{
return !(a == b);
}
// IOstream operator
friend Istream& operator>>(Istream& is, boundBox&);
friend Ostream& operator<<(Ostream& os, const boundBox&);
};
//- Specify data associated with boundBox type is contiguous
template<>
inline bool contiguous<boundBox>() {return contiguous<point>();}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

3
src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
View File

@ -282,8 +282,7 @@ const Foam::labelList& Foam::polyPatch::meshEdges() const
primitivePatch::meshEdges
(
boundaryMesh().mesh().edges(),
boundaryMesh().mesh().cellEdges(),
faceCells()
boundaryMesh().mesh().pointEdges()
)
);
}

34
src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
View File

@ -362,16 +362,27 @@ const Foam::labelList& Foam::faceZone::meshEdges() const
{
if (!mePtr_)
{
labelList faceCells(size());
const labelList& own = zoneMesh().mesh().faceOwner();
const labelList& faceLabels = *this;
forAll (faceCells, faceI)
{
faceCells[faceI] = own[faceLabels[faceI]];
}
//labelList faceCells(size());
//
//const labelList& own = zoneMesh().mesh().faceOwner();
//
//const labelList& faceLabels = *this;
//
//forAll (faceCells, faceI)
//{
// faceCells[faceI] = own[faceLabels[faceI]];
//}
//
//mePtr_ =
// new labelList
// (
// operator()().meshEdges
// (
// zoneMesh().mesh().edges(),
// zoneMesh().mesh().cellEdges(),
// faceCells
// )
// );
mePtr_ =
new labelList
@ -379,8 +390,7 @@ const Foam::labelList& Foam::faceZone::meshEdges() const
operator()().meshEdges
(
zoneMesh().mesh().edges(),
zoneMesh().mesh().cellEdges(),
faceCells
zoneMesh().mesh().pointEdges()
)
);
}

11
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
View File

@ -348,7 +348,8 @@ public:
// index in the edge list. If the edge is not found, return -1
label whichEdge(const edge& e) const;
//- Return labels of patch edges in the global edge list
//- Return labels of patch edges in the global edge list using
// cell addressing
labelList meshEdges
(
const edgeList& allEdges,
@ -356,6 +357,14 @@ public:
const labelList& faceCells
) const;
//- Return labels of patch edges in the global edge list using
// basic edge addressing.
labelList meshEdges
(
const edgeList& allEdges,
const labelListList& pointEdges
) const;
//- Return face normals for patch
const Field<PointType>& faceNormals() const;

54
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
View File

@ -111,6 +111,60 @@ labelList PrimitivePatch<Face, FaceList, PointField, PointType>::meshEdges
}
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
labelList PrimitivePatch<Face, FaceList, PointField, PointType>::meshEdges
(
const edgeList& allEdges,
const labelListList& pointEdges
) const
{
if (debug)
{
Info<< "labelList PrimitivePatch<Face, FaceList, PointField, PointType>"
<< "::meshEdges() : "
<< "calculating labels of patch edges in mesh edge list"
<< endl;
}
// get reference to the list of edges on the patch
const edgeList& PatchEdges = edges();
// create the storage
labelList meshEdges(PatchEdges.size());
// get reference to the points on the patch
const labelList& pp = meshPoints();
// WARNING: Remember that local edges address into local point list;
// local-to-global point label translation is necessary
forAll (PatchEdges, edgeI)
{
const label globalPointI = pp[PatchEdges[edgeI].start()];
const edge curEdge(globalPointI, pp[PatchEdges[edgeI].end()]);
const labelList& pe = pointEdges[globalPointI];
forAll (pe, i)
{
if (allEdges[pe[i]] == curEdge)
{
meshEdges[edgeI] = pe[i];
break;
}
}
}
return meshEdges;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template

6
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
View File

@ -66,6 +66,9 @@ primitiveMesh::primitiveMesh()
ppPtr_(NULL),
cpPtr_(NULL),
allocSize_(0),
labels_(0),
cellCentresPtr_(NULL),
faceCentresPtr_(NULL),
cellVolumesPtr_(NULL),
@ -106,6 +109,9 @@ primitiveMesh::primitiveMesh
ppPtr_(NULL),
cpPtr_(NULL),
allocSize_(0),
labels_(0),
cellCentresPtr_(NULL),
faceCentresPtr_(NULL),
cellVolumesPtr_(NULL),

70
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
View File

@ -155,6 +155,17 @@ class primitiveMesh
mutable labelListList* cpPtr_;
// On-the-fly edge addresing storage
//- Temporary storage for addressing. allocSize is the real size
// of the labelList.
mutable label allocSize_;
mutable labelList labels_;
//- Temporary storage for addressing
mutable labelHashSet labelSet_;
// Geometric data
//- Cell centres
@ -209,8 +220,14 @@ class primitiveMesh
(
List<DynamicList<label> >&,
DynamicList<edge>&,
const label pA,
const label pB
const label,
const label
);
//- For on-the-fly addressing calculation
static label findFirstCommonElementFromSortedLists
(
const labelList&,
const labelList&
);
@ -667,6 +684,55 @@ public:
//- Print a list of all the currently allocated mesh data
void printAllocated() const;
// Per storage whether allocated
inline bool hasCellShapes() const;
inline bool hasEdges() const;
inline bool hasCellCells() const;
inline bool hasEdgeCells() const;
inline bool hasPointCells() const;
inline bool hasCells() const;
inline bool hasEdgeFaces() const;
inline bool hasPointFaces() const;
inline bool hasCellEdges() const;
inline bool hasFaceEdges() const;
inline bool hasPointEdges() const;
inline bool hasPointPoints() const;
inline bool hasCellPoints() const;
inline bool hasCellCentres() const;
inline bool hasFaceCentres() const;
inline bool hasCellVolumes() const;
inline bool hasFaceAreas() const;
// On-the-fly addressing calculation. These functions return either
// a reference to the full addressing (if already calculated) or
// a reference to member data labels_ so be careful when not storing
// result.
//- cellCells using cells
const labelList& cellCells(const label cellI) const;
//- cellPoints using cells
const labelList& cellPoints(const label cellI) const;
//- pointCells using pointFaces
const labelList& pointCells(const label pointI) const;
//- pointPoints using edges, pointEdges
const labelList& pointPoints(const label pointI) const;
//- faceEdges using pointFaces, edges, pointEdges
const labelList& faceEdges(const label faceI) const;
//- edgeFaces using pointFaces, edges, pointEdges
const labelList& edgeFaces(const label edgeI) const;
//- edgeCells using pointFaces, edges, pointEdges
const labelList& edgeCells(const label edgeI) const;
//- cellEdges using cells, pointFaces, edges, pointEdges
const labelList& cellEdges(const label cellI) const;
//- Clear geometry
void clearGeom();

47
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
View File

@ -105,6 +105,53 @@ const labelListList& primitiveMesh::cellCells() const
}
const labelList& primitiveMesh::cellCells(const label cellI) const
{
if (hasCellCells())
{
return cellCells()[cellI];
}
else
{
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cell& cFaces = cells()[cellI];
labels_.size() = allocSize_;
if (cFaces.size() > allocSize_)
{
labels_.clear();
allocSize_ = cFaces.size();
labels_.setSize(allocSize_);
}
label n = 0;
forAll(cFaces, i)
{
label faceI = cFaces[i];
if (faceI < nInternalFaces())
{
if (own[faceI] == cellI)
{
labels_[n++] = nei[faceI];
}
else
{
labels_[n++] = own[faceI];
}
}
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

46
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
View File

@ -53,6 +53,52 @@ const labelListList& primitiveMesh::cellPoints() const
}
const labelList& primitiveMesh::cellPoints(const label cellI) const
{
if (hasCellPoints())
{
return cellPoints()[cellI];
}
else
{
const faceList& fcs = faces();
const labelList& cFaces = cells()[cellI];
labelSet_.clear();
forAll(cFaces, i)
{
const labelList& f = fcs[cFaces[i]];
forAll(f, fp)
{
labelSet_.insert(f[fp]);
}
}
labels_.size() = allocSize_;
if (labelSet_.size() > allocSize_)
{
labels_.clear();
allocSize_ = labelSet_.size();
labels_.setSize(allocSize_);
}
label n = 0;
forAllConstIter(labelHashSet, labelSet_, iter)
{
labels_[n++] = iter.key();
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

7
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
View File

@ -771,11 +771,12 @@ bool primitiveMesh::checkPoints
}
}
const labelListList& pc = pointCells();
forAll (pc, pointI)
forAll (pf, pointI)
{
if (pc[pointI].size() == 0)
const labelList& pc = pointCells(pointI);
if (pc.size() == 0)
{
if (setPtr)
{

80
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C
View File

@ -51,6 +51,86 @@ const labelListList& primitiveMesh::edgeCells() const
}
const labelList& primitiveMesh::edgeCells(const label edgeI) const
{
if (hasEdgeCells())
{
return edgeCells()[edgeI];
}
else
{
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
// edge faces can either return labels_ or reference in edgeLabels.
labelList labelsCopy;
if (!hasEdgeFaces())
{
labelsCopy = edgeFaces(edgeI);
}
const labelList& eFaces =
(
hasEdgeFaces()
? edgeFaces()[edgeI]
: labelsCopy
);
labels_.size() = allocSize_;
// labels_ should certainly be big enough for edge cells.
label n = 0;
// Do quadratic insertion.
forAll(eFaces, i)
{
label faceI = eFaces[i];
{
label ownCellI = own[faceI];
// Check if not already in labels_
for (label j = 0; j < n; j++)
{
if (labels_[j] == ownCellI)
{
ownCellI = -1;
break;
}
}
if (ownCellI != -1)
{
labels_[n++] = ownCellI;
}
}
if (isInternalFace(faceI))
{
label neiCellI = nei[faceI];
for (label j = 0; j < n; j++)
{
if (labels_[j] == neiCellI)
{
neiCellI = -1;
break;
}
}
if (neiCellI != -1)
{
labels_[n++] = neiCellI;
}
}
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

53
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
View File

@ -51,6 +51,59 @@ const labelListList& primitiveMesh::edgeFaces() const
return *efPtr_;
}
const labelList& primitiveMesh::edgeFaces(const label edgeI) const
{
if (hasEdgeFaces())
{
return edgeFaces()[edgeI];
}
else
{
// Use the fact that pointEdges are sorted in incrementing edge order
const edge& e = edges()[edgeI];
const labelList& pFaces0 = pointFaces()[e[0]];
const labelList& pFaces1 = pointFaces()[e[1]];
label i0 = 0;
label i1 = 0;
label n = 0;
labels_.size() = allocSize_;
while (i0 < pFaces0.size() && i1 < pFaces1.size())
{
if (pFaces0[i0] < pFaces1[i1])
{
++i0;
}
else if (pFaces0[i0] > pFaces1[i1])
{
++i1;
}
else
{
// Equal. Append.
if (n == allocSize_)
{
// Have setSize copy contents so far
labels_.size() = n;
allocSize_ = allocSize_*2 + 1;
labels_.setSize(allocSize_);
}
labels_[n++] = pFaces0[i0];
++i0;
++i1;
}
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

125
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
View File

@ -461,6 +461,46 @@ void primitiveMesh::calcEdges(const bool doFaceEdges) const
}
label primitiveMesh::findFirstCommonElementFromSortedLists
(
const labelList& list1,
const labelList& list2
)
{
label result = -1;
labelList::const_iterator iter1 = list1.begin();
labelList::const_iterator iter2 = list2.begin();
while (iter1 != list1.end() && iter2 != list2.end())
{
if( *iter1 < *iter2)
{
++iter1;
}
else if (*iter1 > *iter2)
{
++iter2;
}
else
{
result = *iter1;
break;
}
}
if (result == -1)
{
FatalErrorIn
(
"primitiveMesh::findFirstCommonElementFromSortedLists"
"(const labelList&, const labelList&)"
) << "No common elements in lists " << list1 << " and " << list2
<< abort(FatalError);
}
return result;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const edgeList& primitiveMesh::edges() const
@ -542,6 +582,91 @@ void primitiveMesh::clearOutEdges()
deleteDemandDrivenData(edgesPtr_);
deleteDemandDrivenData(pePtr_);
deleteDemandDrivenData(fePtr_);
labels_.clear();
allocSize_ = 0;
}
const labelList& primitiveMesh::faceEdges(const label faceI) const
{
if (hasFaceEdges())
{
return faceEdges()[faceI];
}
else
{
const labelListList& pointEs = pointEdges();
const face& f = faces()[faceI];
labels_.size() = allocSize_;
if (f.size() > allocSize_)
{
labels_.clear();
allocSize_ = f.size();
labels_.setSize(allocSize_);
}
label n = 0;
forAll(f, fp)
{
labels_[n++] = findFirstCommonElementFromSortedLists
(
pointEs[f[fp]],
pointEs[f.nextLabel(fp)]
);
}
labels_.size() = n;
return labels_;
}
}
const labelList& primitiveMesh::cellEdges(const label cellI) const
{
if (hasCellEdges())
{
return cellEdges()[cellI];
}
else
{
const labelList& cFaces = cells()[cellI];
labelSet_.clear();
forAll(cFaces, i)
{
const labelList& fe = faceEdges(cFaces[i]);
forAll(fe, feI)
{
labelSet_.insert(fe[feI]);
}
}
labels_.size() = allocSize_;
if (labelSet_.size() > allocSize_)
{
labels_.clear();
allocSize_ = labelSet_.size();
labels_.setSize(allocSize_);
}
label n =0;
forAllConstIter(labelHashSet, labelSet_, iter)
{
labels_[n++] = iter.key();
}
labels_.size() = n;
return labels_;
}
}

102
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H
View File

@ -104,6 +104,108 @@ inline bool primitiveMesh::isInternalFace(const label faceIndex) const
}
inline bool primitiveMesh::hasCellShapes() const
{
return cellShapesPtr_;
}
inline bool primitiveMesh::hasEdges() const
{
return edgesPtr_;
}
inline bool primitiveMesh::hasCellCells() const
{
return ccPtr_;
}
inline bool primitiveMesh::hasEdgeCells() const
{
return ecPtr_;
}
inline bool primitiveMesh::hasPointCells() const
{
return pcPtr_;
}
inline bool primitiveMesh::hasCells() const
{
return cfPtr_;
}
inline bool primitiveMesh::hasEdgeFaces() const
{
return efPtr_;
}
inline bool primitiveMesh::hasPointFaces() const
{
return pfPtr_;
}
inline bool primitiveMesh::hasCellEdges() const
{
return cePtr_;
}
inline bool primitiveMesh::hasFaceEdges() const
{
return fePtr_;
}
inline bool primitiveMesh::hasPointEdges() const
{
return pePtr_;
}
inline bool primitiveMesh::hasPointPoints() const
{
return ppPtr_;
}
inline bool primitiveMesh::hasCellPoints() const
{
return cpPtr_;
}
inline bool primitiveMesh::hasCellCentres() const
{
return cellCentresPtr_;
}
inline bool primitiveMesh::hasFaceCentres() const
{
return faceCentresPtr_;
}
inline bool primitiveMesh::hasCellVolumes() const
{
return cellVolumesPtr_;
}
inline bool primitiveMesh::hasFaceAreas() const
{
return faceAreasPtr_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

64
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
View File

@ -114,6 +114,70 @@ const labelListList& primitiveMesh::pointCells() const
}
const labelList& primitiveMesh::pointCells(const label pointI) const
{
if (hasPointCells())
{
return pointCells()[pointI];
}
else
{
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const labelList& pFaces = pointFaces()[pointI];
labels_.size() = allocSize_;
label n = 0;
forAll(pFaces, i)
{
const label faceI = pFaces[i];
// Append owner
if (n == allocSize_)
{
labels_.size() = n;
allocSize_ = allocSize_*2 + 1;
labels_.setSize(allocSize_);
}
labels_[n++] = own[faceI];
// Append neighbour
if (faceI < nInternalFaces())
{
if (n == allocSize_)
{
labels_.size() = n;
allocSize_ = allocSize_*2 + 1;
labels_.setSize(allocSize_);
}
labels_[n++] = nei[faceI];
}
}
labels_.size() = n;
// Filter duplicates
sort(labels_);
n = 1;
for (label i = 1; i < labels_.size(); i++)
{
if (labels_[i] != labels_[i-1])
{
labels_[n++] = labels_[i];
}
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

34
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C
View File

@ -97,6 +97,40 @@ const labelListList& primitiveMesh::pointPoints() const
}
const labelList& primitiveMesh::pointPoints(const label pointI) const
{
if (hasPointPoints())
{
return pointPoints()[pointI];
}
else
{
const edgeList& edges = this->edges();
const labelList& pEdges = pointEdges()[pointI];
labels_.size() = allocSize_;
if (pEdges.size() > allocSize_)
{
// Set size() so memory allocation behaves as normal.
labels_.clear();
allocSize_ = pEdges.size();
labels_.setSize(allocSize_);
}
label n = 0;
forAll(pEdges, i)
{
labels_[n++] = edges[pEdges[i]].otherVertex(pointI);
}
labels_.size() = n;
return labels_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

23
src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
View File

@ -30,6 +30,7 @@ License
#include "triSurfaceMesh.H"
#include "refinementSurfaces.H"
#include "searchableSurfaces.H"
#include "orientedSurface.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -211,7 +212,27 @@ void Foam::shellSurfaces::orient()
refCast<const triSurfaceMesh>(s)
);
refinementSurfaces::orientSurface(outsidePt, shell);
// Flip surface so outsidePt is outside.
bool anyFlipped = orientedSurface::orient
(
shell,
outsidePt,
true
);
if (anyFlipped)
{
// orientedSurface will have done a clearOut of the surface.
// we could do a clearout of the triSurfaceMeshes::trees()
// but these aren't affected by orientation
// (except for cached
// sideness which should not be set at this point.
// !!Should check!)
Info<< "shellSurfaces : Flipped orientation of surface "
<< s.name()
<< " so point " << outsidePt << " is outside." << endl;
}
}
}
}

44
src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
View File

@ -29,8 +29,8 @@ License
#include "matchPoints.H"
#include "indirectPrimitivePatch.H"
#include "meshTools.H"
#include "octreeDataFace.H"
#include "octree.H"
#include "treeDataFace.H"
#include "indexedOctree.H"
#include "OFstream.H"
#include "IndirectList.H"
@ -1011,19 +1011,29 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
)
{
// Construct octree from all mesh0 boundary faces
octreeDataFace shapes(mesh0);
labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
forAll(bndFaces, i)
{
bndFaces[i] = mesh0.nInternalFaces() + i;
}
treeBoundBox overallBb(mesh0.points());
octree<octreeDataFace> tree
(
overallBb, // overall search domain
shapes, // all information needed to do checks on cells
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
10.0
);
Random rndGen(123456);
indexedOctree<treeDataFace> tree
(
treeDataFace // all information needed to search faces
(
false, // do not cache bb
mesh0,
bndFaces // boundary faces only
),
overallBb.extend(rndGen, 1E-4), // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
if (debug)
{
@ -1048,17 +1058,11 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
// Generate face centre (prevent cellCentres() reconstruction)
point fc(f1.centre(mesh1.points()));
// Search in bounding box of face only.
treeBoundBox tightest(static_cast<const pointField&>(f1.points(mesh1.points())));
pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
scalar tightestDist = GREAT;
label index = tree.findNearest(fc, tightest, tightestDist);
if (index != -1)
if (nearInfo.hit())
{
label mesh0FaceI = shapes.meshFaces()[index];
label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
// Check if points of f1 actually lie on top of mesh0 face
// This is the bit that might fail since if f0 is severely warped

354
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
View File

@ -62,6 +62,7 @@ namespace Foam
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::hexRef8::reorder
@ -80,7 +81,7 @@ void Foam::hexRef8::reorder
if (newI >= len)
{
FatalErrorIn("hexRef8::reorder") << abort(FatalError);
FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
}
if (newI >= 0)
@ -557,22 +558,11 @@ Foam::label Foam::hexRef8::getAnchorCell
}
// Problem.
// Pick up points of the cell
const labelList cPoints(cellPoints(cellI));
Perr<< "cell:" << cellI << " points:" << endl;
forAll(cPoints, i)
{
label pointI = cPoints[i];
Perr<< " " << pointI << " coord:" << mesh_.points()[pointI]
<< nl;
}
dumpCell(cellI);
Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI]
<< endl;
FatalErrorIn("hexRef8::getAnchorCell")
FatalErrorIn("hexRef8::getAnchorCell(..)")
<< "Could not find point " << pointI
<< " in the anchorPoints for cell " << cellI << endl
<< "Does your original mesh obey the 2:1 constraint and"
@ -690,9 +680,50 @@ Foam::label Foam::hexRef8::countAnchors
}
void Foam::hexRef8::dumpCell(const label cellI) const
{
OFstream str(mesh_.time().path()/"cell_" + Foam::name(cellI) + ".obj");
Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl;
const cell& cFaces = mesh_.cells()[cellI];
Map<label> pointToObjVert;
label objVertI = 0;
forAll(cFaces, i)
{
const face& f = mesh_.faces()[cFaces[i]];
forAll(f, fp)
{
if (pointToObjVert.insert(f[fp], objVertI))
{
meshTools::writeOBJ(str, mesh_.points()[f[fp]]);
objVertI++;
}
}
}
forAll(cFaces, i)
{
const face& f = mesh_.faces()[cFaces[i]];
forAll(f, fp)
{
label pointI = f[fp];
label nexPointI = f[f.fcIndex(fp)];
str << "l " << pointToObjVert[pointI]+1
<< ' ' << pointToObjVert[nexPointI]+1 << nl;
}
}
}
// Find point with certain pointLevel. Skip any higher levels.
Foam::label Foam::hexRef8::findLevel
(
const label faceI,
const face& f,
const label startFp,
const bool searchForward,
@ -707,7 +738,13 @@ Foam::label Foam::hexRef8::findLevel
if (pointLevel_[pointI] < wantedLevel)
{
FatalErrorIn("hexRef8::findLevel")
dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("hexRef8::findLevel(..)")
<< "face:" << f
<< " level:" << IndirectList<label>(pointLevel_, f)()
<< " startFp:" << startFp
@ -729,7 +766,13 @@ Foam::label Foam::hexRef8::findLevel
}
}
FatalErrorIn("hexRef8::findLevel")
dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("hexRef8::findLevel(..)")
<< "face:" << f
<< " level:" << IndirectList<label>(pointLevel_, f)()
<< " startFp:" << startFp
@ -763,15 +806,6 @@ Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
}
else
{
//FatalErrorIn("hexRef8::getAnchorLevel(const label) const")
// << "face:" << faceI
// << " centre:" << mesh_.faceCentres()[faceI]
// << " verts:" << f
// << " point levels:" << IndirectList<label>(pointLevel_, f)()
// << " own:" << mesh_.faceOwner()[faceI]
// << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]]
// << abort(FatalError);
return -1;
}
}
@ -797,7 +831,7 @@ void Foam::hexRef8::checkInternalOrientation
if ((dir & n) < 0)
{
FatalErrorIn("checkInternalOrientation")
FatalErrorIn("checkInternalOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << endl
<< " coords:" << compactPoints
@ -812,7 +846,7 @@ void Foam::hexRef8::checkInternalOrientation
if (s < 0.1 || s > 0.9)
{
FatalErrorIn("checkInternalOrientation")
FatalErrorIn("checkInternalOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << endl
<< " coords:" << compactPoints
@ -843,7 +877,7 @@ void Foam::hexRef8::checkBoundaryOrientation
if ((dir & n) < 0)
{
FatalErrorIn("checkBoundaryOrientation")
FatalErrorIn("checkBoundaryOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace
<< " coords:" << compactPoints
@ -858,7 +892,7 @@ void Foam::hexRef8::checkBoundaryOrientation
if (s < 0.7 || s > 1.3)
{
WarningIn("checkBoundaryOrientation")
WarningIn("checkBoundaryOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace
<< " coords:" << compactPoints
@ -1191,8 +1225,22 @@ void Foam::hexRef8::createInternalFaces
}
// Now the face mid point is the second cLevel+1 point
label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1);
label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1);
label edgeMid = findLevel
(
faceI,
f,
f.fcIndex(anchorFp),
true,
cLevel+1
);
label faceMid = findLevel
(
faceI,
f,
f.fcIndex(edgeMid),
true,
cLevel+1
);
faceMidPointI = f[faceMid];
}
@ -1205,7 +1253,13 @@ void Foam::hexRef8::createInternalFaces
}
else
{
FatalErrorIn("createInternalFaces")
dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("createInternalFaces(..)")
<< "nAnchors:" << nAnchors
<< " faceI:" << faceI
<< abort(FatalError);
@ -1243,9 +1297,11 @@ void Foam::hexRef8::createInternalFaces
if (edgeMidPointI == -1)
{
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI));
FatalErrorIn("createInternalFaces")
FatalErrorIn("createInternalFaces(..)")
<< "cell:" << cellI << " cLevel:" << cLevel
<< " cell points:" << cPoints
<< " pointLevel:"
@ -1264,7 +1320,7 @@ void Foam::hexRef8::createInternalFaces
else
{
// Search forward in face to clevel+1
label edgeMid = findLevel(f, fp1, true, cLevel+1);
label edgeMid = findLevel(faceI, f, fp1, true, cLevel+1);
edgeMidPointI = f[edgeMid];
}
@ -1314,9 +1370,11 @@ void Foam::hexRef8::createInternalFaces
if (edgeMidPointI == -1)
{
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI));
FatalErrorIn("createInternalFaces")
FatalErrorIn("createInternalFaces(..)")
<< "cell:" << cellI << " cLevel:" << cLevel
<< " cell points:" << cPoints
<< " pointLevel:"
@ -1335,7 +1393,14 @@ void Foam::hexRef8::createInternalFaces
else
{
// Search back to clevel+1
label edgeMid = findLevel(f, fpMin1, false, cLevel+1);
label edgeMid = findLevel
(
faceI,
f,
fpMin1,
false,
cLevel+1
);
edgeMidPointI = f[edgeMid];
}
@ -1591,9 +1656,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
if (mag(ownLevel-neiLevel) > 1)
{
dumpCell(own);
dumpCell(nei);
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const labelList&)"
"hexRef8::checkWantedRefinementLevels(const labelList&)"
) << "cell:" << own
<< " current level:" << cellLevel_[own]
<< " level after refinement:" << ownLevel
@ -1633,9 +1700,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
{
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
dumpCell(own);
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const labelList&)"
"hexRef8::checkWantedRefinementLevels(const labelList&)"
) << "Celllevel does not satisfy 2:1 constraint."
<< " On coupled face "
<< faceI
@ -1730,6 +1798,25 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
<< abort(FatalError);
}
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&)"
) << "Restarted from inconsistent cellLevel or pointLevel files."
<< endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError);
}
// Check refinement levels for consistency
checkRefinementLevels(-1, labelList(0));
@ -1806,7 +1893,24 @@ Foam::hexRef8::hexRef8
) << "History enabled but number of visible cells in it "
<< history_.visibleCells().size()
<< " is not equal to the number of cells in the mesh "
<< mesh_.nCells()
<< mesh_.nCells() << abort(FatalError);
}
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&, const labelList&"
", const labelList&, const refinementHistory&)"
) << "Incorrect cellLevel or pointLevel size." << endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError);
}
@ -1877,6 +1981,24 @@ Foam::hexRef8::hexRef8
savedPointLevel_(0),
savedCellLevel_(0)
{
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&, const labelList&"
", const labelList&)"
) << "Incorrect cellLevel or pointLevel size." << endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError);
}
// Check refinement levels for consistency
checkRefinementLevels(-1, labelList(0));
@ -2356,9 +2478,13 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
if (mag(ownLevel-neiLevel) > 1)
{
dumpCell(own);
dumpCell(nei);
FatalErrorIn
(
"hexRef8::consistentSlowRefinement"
"(const label, const labelList&, const labelList&"
", const label, const labelList&)"
) << "cell:" << own
<< " current level:" << cellLevel_[own]
<< " current refData:" << allCellInfo[own]
@ -2412,11 +2538,14 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
if (mag(ownLevel - nbrLevel) > 1)
{
dumpCell(own);
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const labelList&)"
"hexRef8::consistentSlowRefinement"
"(const label, const labelList&, const labelList&"
", const label, const labelList&)"
) << "Celllevel does not satisfy 2:1 constraint."
<< " On coupled face "
<< faceI
@ -2805,7 +2934,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
<< cellsToRefine.size() << " to " << newCellsToRefine.size()
<< " cells to refine." << endl;
//XXXX
// Check that newCellsToRefine obeys at least 2:1.
{
@ -2861,6 +2989,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
&& savedRefineCell.get(cellI) == 0
)
{
dumpCell(cellI);
FatalErrorIn
(
"hexRef8::consistentSlowRefinement2"
@ -2873,7 +3002,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
}
}
}
//XXXX
}
return newCellsToRefine;
@ -3342,6 +3470,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
{
if (nAnchorPoints[cellI] == 8)
{
dumpCell(cellI);
FatalErrorIn
(
"hexRef8::setRefinement(const labelList&"
@ -3365,6 +3494,8 @@ Foam::labelListList Foam::hexRef8::setRefinement
{
if (nAnchorPoints[cellI] != 8)
{
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI));
FatalErrorIn
@ -3841,7 +3972,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
if (minPointI != labelMax && minPointI != mesh_.nPoints())
{
FatalErrorIn("hexRef8::setRefinement")
FatalErrorIn("hexRef8::setRefinement(..)")
<< "Added point labels not consecutive to existing mesh points."
<< nl
<< "mesh_.nPoints():" << mesh_.nPoints()
@ -3995,13 +4126,6 @@ void Foam::hexRef8::updateMesh
if (oldCellI == -1)
{
//FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
// << "Problem : cell " << newCellI
// << " at " << mesh_.cellCentres()[newCellI]
// << " does not originate from another cell"
// << " (i.e. is inflated)." << nl
// << "Hence we cannot determine the new cellLevel"
// << " for it." << abort(FatalError);
newCellLevel[newCellI] = -1;
}
else
@ -4036,8 +4160,8 @@ void Foam::hexRef8::updateMesh
//{
// WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
// << "Problem : "
// << "cellLevel_ contains illegal value -1 after mapping at cell "
// << findIndex(cellLevel_, -1) << endl
// << "cellLevel_ contains illegal value -1 after mapping
// << " at cell " << findIndex(cellLevel_, -1) << endl
// << "This means that another program has inflated cells"
// << " (created cells out-of-nothing) and hence we don't know"
// << " their cell level. Continuing with illegal value."
@ -4174,7 +4298,7 @@ void Foam::hexRef8::subset
if (findIndex(cellLevel_, -1) != -1)
{
FatalErrorIn("hexRef8::subset")
FatalErrorIn("hexRef8::subset(..)")
<< "Problem : "
<< "cellLevel_ contains illegal value -1 after mapping:"
<< cellLevel_
@ -4195,7 +4319,7 @@ void Foam::hexRef8::subset
if (findIndex(pointLevel_, -1) != -1)
{
FatalErrorIn("hexRef8::subset")
FatalErrorIn("hexRef8::subset(..)")
<< "Problem : "
<< "pointLevel_ contains illegal value -1 after mapping:"
<< pointLevel_
@ -4292,6 +4416,7 @@ void Foam::hexRef8::checkMesh() const
if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
{
dumpCell(own);
FatalErrorIn("hexRef8::checkMesh()")
<< "Faces do not seem to be correct across coupled"
<< " boundaries" << endl
@ -4335,6 +4460,8 @@ void Foam::hexRef8::checkMesh() const
const face& f = mesh_.faces()[faceI];
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
dumpCell(mesh_.faceOwner()[faceI]);
FatalErrorIn("hexRef8::checkMesh()")
<< "Faces do not seem to be correct across coupled"
<< " boundaries" << endl
@ -4373,6 +4500,8 @@ void Foam::hexRef8::checkMesh() const
if (f.size() != nVerts[i])
{
dumpCell(mesh_.faceOwner()[faceI]);
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn("hexRef8::checkMesh()")
@ -4421,6 +4550,8 @@ void Foam::hexRef8::checkMesh() const
if (mag(anchorVec - anchorPoints[i]) > smallDim)
{
dumpCell(mesh_.faceOwner()[faceI]);
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn("hexRef8::checkMesh()")
@ -4487,6 +4618,9 @@ void Foam::hexRef8::checkRefinementLevels
if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
{
dumpCell(own);
dumpCell(nei);
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const label)"
@ -4520,6 +4654,8 @@ void Foam::hexRef8::checkRefinementLevels
if (mag(cellLevel_[own] - neiLevel[i]) > 1)
{
dumpCell(own);
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn
@ -4621,6 +4757,8 @@ void Foam::hexRef8::checkRefinementLevels
> maxPointDiff
)
{
dumpCell(cellI);
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const label)"
@ -4637,6 +4775,72 @@ void Foam::hexRef8::checkRefinementLevels
}
}
}
// Hanging points
{
// Any patches with points having only two edges.
boolList isHangingPoint(mesh_.nPoints(), false);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& meshPoints = pp.meshPoints();
forAll(meshPoints, i)
{
label pointI = meshPoints[i];
const labelList& pEdges = mesh_.pointEdges()[pointI];
if (pEdges.size() == 2)
{
isHangingPoint[pointI] = true;
}
}
}
syncTools::syncPointList
(
mesh_,
isHangingPoint,
andEqOp<bool>(), // only if all decide it is hanging point
true, // null
false // no separation
);
OFstream str(mesh_.time().path()/"hangingPoints.obj");
label nHanging = 0;
forAll(isHangingPoint, pointI)
{
if (isHangingPoint[pointI])
{
nHanging++;
Pout<< "Hanging boundary point " << pointI
<< " at " << mesh_.points()[pointI]
<< endl;
meshTools::writeOBJ(str, mesh_.points()[pointI]);
}
}
if (returnReduce(nHanging, sumOp<label>()) > 0)
{
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const label)"
) << "Detected a point used by two edges only (hanging point)"
<< nl << "This is not allowed"
<< abort(FatalError);
}
}
}
@ -4865,8 +5069,10 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
if (maxSet)
{
FatalErrorIn("hexRef8::consistentUnrefinement")
<< "maxSet not implemented yet."
FatalErrorIn
(
"hexRef8::consistentUnrefinement(const labelList&, const bool"
) << "maxSet not implemented yet."
<< abort(FatalError);
}
@ -4935,7 +5141,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{
if (unrefineCell.get(own) == 0)
{
FatalErrorIn("hexRef8::consistentUnrefinement")
FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError);
}
@ -4953,7 +5159,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{
if (unrefineCell.get(nei) == 0)
{
FatalErrorIn("hexRef8::consistentUnrefinement")
FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError);
}
@ -4989,7 +5195,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{
if (unrefineCell.get(own) == 0)
{
FatalErrorIn("hexRef8::consistentUnrefinement")
FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError);
}
@ -5003,7 +5209,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{
if (unrefineCell.get(own) == 1)
{
FatalErrorIn("hexRef8::consistentUnrefinement")
FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError);
}
@ -5086,8 +5292,10 @@ void Foam::hexRef8::setUnrefinement
{
if (!history_.active())
{
FatalErrorIn("hexRef8::setUnrefinement()")
<< "Only call if constructed with history capability"
FatalErrorIn
(
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Only call if constructed with history capability"
<< abort(FatalError);
}
@ -5102,8 +5310,10 @@ void Foam::hexRef8::setUnrefinement
{
if (cellLevel_[cellI] < 0)
{
FatalErrorIn("hexRef8::setUnrefinement()")
<< "Illegal cell level " << cellLevel_[cellI]
FatalErrorIn
(
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Illegal cell level " << cellLevel_[cellI]
<< " for cell " << cellI
<< abort(FatalError);
}
@ -5165,8 +5375,10 @@ void Foam::hexRef8::setUnrefinement
if (facesToRemove.size() != splitFaces.size())
{
FatalErrorIn("hexRef8::setUnrefinement")
<< "Ininitial set of split points to unrefine does not"
FatalErrorIn
(
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Ininitial set of split points to unrefine does not"
<< " seem to be consistent or not mid points of refined cells"
<< abort(FatalError);
}
@ -5187,8 +5399,10 @@ void Foam::hexRef8::setUnrefinement
// Check
if (pCells.size() != 8)
{
FatalErrorIn("hexRef8::setUnrefinement")
<< "splitPoint " << pointI
FatalErrorIn
(
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "splitPoint " << pointI
<< " should have 8 cells using it. It has " << pCells
<< abort(FatalError);
}
@ -5208,7 +5422,7 @@ void Foam::hexRef8::setUnrefinement
if (region == -1)
{
FatalErrorIn("hexRef8::setUnrefinement")
FatalErrorIn("hexRef8::setUnrefinement(..)")
<< "Ininitial set of split points to unrefine does not"
<< " seem to be consistent or not mid points"
<< " of refined cells" << nl
@ -5219,7 +5433,7 @@ void Foam::hexRef8::setUnrefinement
if (masterCellI != cellRegionMaster[region])
{
FatalErrorIn("hexRef8::setUnrefinement")
FatalErrorIn("hexRef8::setUnrefinement(..)")
<< "cell:" << cellI << " on splitPoint:" << pointI
<< " in region " << region
<< " has master:" << cellRegionMaster[region]

3
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
View File

@ -174,9 +174,12 @@ class hexRef8
label findMaxLevel(const labelList& f) const;
//- Count number of vertices <= anchorLevel
label countAnchors(const labelList&, const label) const;
//- Debugging: dump cell as .obj file
void dumpCell(const label cellI) const;
//- Find index of point with wantedLevel, starting from fp.
label findLevel
(
const label faceI,
const face& f,
const label startFp,
const bool searchForward,

4
src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
View File

@ -1295,7 +1295,7 @@ void Foam::polyTopoChange::calcFaceInflationMaps
selectFaces
(
mesh,
mesh.edgeFaces()[iter()],
mesh.edgeFaces(iter()),
true
)
);
@ -1309,7 +1309,7 @@ void Foam::polyTopoChange::calcFaceInflationMaps
selectFaces
(
mesh,
mesh.edgeFaces()[iter()],
mesh.edgeFaces(iter()),
false
)
);

20
src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
View File

@ -807,7 +807,7 @@ void Foam::removeFaces::setRefinement
// Edges to remove
labelHashSet edgesToRemove(faceLabels.size());
// Per face the region it is. -1 for removed faces, -2 for regions
// Per face the region it is in. -1 for removed faces, -2 for regions
// consisting of single face only.
labelList faceRegion(mesh_.nFaces(), -1);
@ -1258,10 +1258,15 @@ void Foam::removeFaces::setRefinement
// are only used by 2 unremoved edges.
{
// Usage of points by non-removed edges.
labelList nEdgesPerPoint(mesh_.nPoints(), labelMax);
labelList nEdgesPerPoint(mesh_.nPoints());
const labelListList& pointEdges = mesh_.pointEdges();
forAll(pointEdges, pointI)
{
nEdgesPerPoint[pointI] = pointEdges[pointI].size();
}
forAllConstIter(labelHashSet, edgesToRemove, iter)
{
// Edge will get removed.
@ -1269,16 +1274,7 @@ void Foam::removeFaces::setRefinement
forAll(e, i)
{
label pointI = e[i];
if (nEdgesPerPoint[pointI] == labelMax)
{
nEdgesPerPoint[pointI] = pointEdges[pointI].size()-1;
}
else
{
nEdgesPerPoint[pointI]--;
}
nEdgesPerPoint[e[i]]--;
}
}

6
src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
View File

@ -159,7 +159,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
// Get the scheduling information
const List<labelPair>& schedule = mpp.schedule();
const labelListList& sendCellLabels = mpp.sendCellLabels();
const labelListList& sendLabels = mpp.sendLabels();
const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
@ -177,7 +177,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
toProc<< IndirectList<Type>
(
this->internalField(),
sendCellLabels[recvProc]
sendLabels[recvProc]
)();
}
else
@ -204,7 +204,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
IndirectList<Type> fromFld
(
this->internalField(),
sendCellLabels[Pstream::myProcNo()]
sendLabels[Pstream::myProcNo()]
);
// Destination faces

4
src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
View File

@ -978,14 +978,14 @@ CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr
ddt0_<GeometricField<Type, fvPatchField, volMesh> >
(
"ddt0(" + U.name() + ')',
rho.dimensions()*U.dimensions()
U.dimensions()
);
DDt0Field<fluxFieldType>& dphidt0 =
ddt0_<fluxFieldType>
(
"ddt0(" + phi.name() + ')',
phi.dimensions()
U.dimensions()*dimArea
);
IOobject ddtIOobject

3
src/lagrangian/basic/Cloud/CloudIO.C
View File

@ -118,7 +118,8 @@ void Foam::Cloud<ParticleType>::checkFieldIOobject
"void Cloud<ParticleType>::checkFieldIOobject"
"(Cloud<ParticleType>, IOField<DataType>)"
) << "Size of " << data.name()
<< " field does not match the number of particles"
<< " field " << data.size()
<< " does not match the number of particles " << c.size()
<< abort(FatalError);
}
}

4
src/meshTools/Make/options
View File

@ -3,5 +3,5 @@ EXE_INC = \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
LIB_LIBS = \
-ltriSurface \
-llagrangian
-ltriSurface \
-llagrangian

327
src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
View File

@ -40,6 +40,18 @@ namespace Foam
addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word);
addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary);
template<>
const char* NamedEnum<directMappedPolyPatch::sampleMode, 3>::names[] =
{
"nearestCell",
"nearestPatchFace",
"nearestFace"
};
const NamedEnum<directMappedPolyPatch::sampleMode, 3>
directMappedPolyPatch::sampleModeNames_;
}
@ -113,107 +125,236 @@ void Foam::directMappedPolyPatch::collectSamples
void Foam::directMappedPolyPatch::findSamples
(
const pointField& samples,
labelList& sampleCellProcs,
labelList& sampleCells,
pointField& sampleCc
labelList& sampleProcs,
labelList& sampleIndices,
pointField& sampleLocations
) const
{
sampleCellProcs.setSize(samples.size());
sampleCells.setSize(samples.size());
sampleCc.setSize(samples.size());
sampleCc = point(-GREAT, -GREAT, -GREAT);
const polyMesh& mesh = boundaryMesh().mesh();
// All the info for nearest. Construct to miss
List<nearInfo> nearest(samples.size());
switch (mode_)
{
// Octree based search engine
meshSearch meshSearchEngine(boundaryMesh().mesh(), false);
forAll(samples, sampleI)
case NEARESTCELL:
{
sampleCells[sampleI] = meshSearchEngine.findCell(samples[sampleI]);
if (samplePatch_.size() != 0 && samplePatch_ != "none")
{
FatalErrorIn
(
"directMappedPolyPatch::findSamples(const pointField&,"
" labelList&, labelList&, pointField&) const"
) << "No need to supply a patch name when in "
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
if (sampleCells[sampleI] == -1)
// Octree based search engine
meshSearch meshSearchEngine(mesh, false);
forAll(samples, sampleI)
{
sampleCellProcs[sampleI] = -1;
const point& sample = samples[sampleI];
label cellI = meshSearchEngine.findCell(sample);
if (cellI == -1)
{
nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
{
const point& cc = mesh.cellCentres()[cellI];
nearest[sampleI].first() = pointIndexHit
(
true,
cc,
cellI
);
nearest[sampleI].second().first() = magSqr(cc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
}
else
break;
}
case NEARESTPATCHFACE:
{
label patchI = boundaryMesh().findPatchID(samplePatch_);
if (patchI == -1)
{
sampleCellProcs[sampleI] = Pstream::myProcNo();
sampleCc[sampleI] =
boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]];
FatalErrorIn("directMappedPolyPatch::findSamples(..)")
<< "Cannot find patch " << samplePatch_ << endl
<< "Valid patches are " << boundaryMesh().names()
<< exit(FatalError);
}
Random rndGen(123456);
const polyPatch& pp = boundaryMesh()[patchI];
// patch faces
const labelList patchFaces(identity(pp.size()) + pp.start());
const treeBoundBox patchBb
(
treeBoundBox(pp.points(), pp.meshPoints()).extend
(
rndGen,
1E-4
)
);
autoPtr<indexedOctree<treeDataFace> > boundaryTree
(
new indexedOctree<treeDataFace>
(
treeDataFace // all information needed to search faces
(
false, // do not cache bb
mesh,
patchFaces // boundary faces only
),
patchBb, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
pointIndexHit& nearInfo = nearest[sampleI].first();
nearInfo = boundaryTree().findNearest
(
sample,
magSqr(patchBb.max()-patchBb.min())
);
if (!nearInfo.hit())
{
nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
{
point fc(pp[nearInfo.index()].centre(pp.points()));
nearInfo.setPoint(fc);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
}
break;
}
case NEARESTFACE:
{
if (samplePatch_.size() != 0 && samplePatch_ != "none")
{
FatalErrorIn
(
"directMappedPolyPatch::findSamples(const pointField&,"
" labelList&, labelList&, pointField&) const"
) << "No need to supply a patch name when in "
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, false);
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
label faceI = meshSearchEngine.findNearestFace(sample);
if (faceI == -1)
{
nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
{
const point& fc = mesh.faceCentres()[faceI];
nearest[sampleI].first() = pointIndexHit
(
true,
fc,
faceI
);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
}
break;
}
default:
{
FatalErrorIn("directMappedPolyPatch::findSamples(..)")
<< "problem." << abort(FatalError);
}
}
// Use only that processor that contains the sample
Pstream::listCombineGather(sampleCells, maxEqOp<label>());
Pstream::listCombineScatter(sampleCells);
labelList minSampleCellProcs(sampleCellProcs);
Pstream::listCombineGather(sampleCellProcs, maxEqOp<label>());
Pstream::listCombineScatter(sampleCellProcs);
Pstream::listCombineGather(sampleCc, maxEqOp<point>());
Pstream::listCombineScatter(sampleCc);
// Find nearest.
Pstream::listCombineGather(nearest, nearestEqOp());
Pstream::listCombineScatter(nearest);
if (debug)
{
Info<< "directMappedPolyPatch::findSamples : " << endl;
forAll(sampleCells, sampleI)
forAll(nearest, sampleI)
{
Info<< " " << sampleI << " coord:" << samples[sampleI]
<< " found on processor:" << sampleCellProcs[sampleI]
<< " in cell:" << sampleCells[sampleI]
<< " with cc:" << sampleCc[sampleI] << endl;
label procI = nearest[sampleI].second().second();
label localI = nearest[sampleI].first().index();
Info<< " " << sampleI << " coord:"<< samples[sampleI]
<< " found on processor:" << procI
<< " in local cell/face:" << localI
<< " with cc:" << nearest[sampleI].first().rawPoint() << endl;
}
}
// Check for samples not being found
forAll(sampleCells, sampleI)
forAll(nearest, sampleI)
{
if (sampleCells[sampleI] == -1)
if (!nearest[sampleI].first().hit())
{
FatalErrorIn
(
"findSamples(const pointField&, labelList&, labelList&)"
"directMappedPolyPatch::findSamples"
"(const pointField&, labelList&"
", labelList&, pointField&)"
) << "Did not find sample " << samples[sampleI]
<< " on any processor." << exit(FatalError);
}
}
// Check for samples being found in multiple cells
{
forAll(minSampleCellProcs, sampleI)
{
if (minSampleCellProcs[sampleI] == -1)
{
minSampleCellProcs[sampleI] = labelMax;
}
}
Pstream::listCombineGather(minSampleCellProcs, minEqOp<label>());
Pstream::listCombineScatter(minSampleCellProcs);
// Convert back into proc+local index
sampleProcs.setSize(samples.size());
sampleIndices.setSize(samples.size());
sampleLocations.setSize(samples.size());
forAll(minSampleCellProcs, sampleI)
{
if (minSampleCellProcs[sampleI] != sampleCellProcs[sampleI])
{
FatalErrorIn
(
"directMappedPolyPatch::findSamples"
"(const pointField&, labelList&, labelList&)"
) << "Found sample " << samples[sampleI]
<< " on processor " << minSampleCellProcs[sampleI]
<< " and on processor " << sampleCellProcs[sampleI] << endl
<< "This is illegal. {lease move your sampling plane."
<< exit(FatalError);
}
}
forAll(nearest, sampleI)
{
sampleProcs[sampleI] = nearest[sampleI].second().second();
sampleIndices[sampleI] = nearest[sampleI].first().index();
sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
}
}
void Foam::directMappedPolyPatch::calcMapping() const
{
if (sendCellLabelsPtr_.valid())
if (sendLabelsPtr_.valid())
{
FatalErrorIn("directMappedPolyPatch::calcMapping() const")
<< "Mapping already calculated" << exit(FatalError);
@ -236,18 +377,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
pointField patchFc;
collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
// Find processor and cell samples are in
labelList sampleCellProcs;
labelList sampleCells;
pointField sampleCc;
findSamples(samples, sampleCellProcs, sampleCells, sampleCc);
// Find processor and cell/face samples are in and actual location.
labelList sampleProcs;
labelList sampleIndices;
pointField sampleLocations;
findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
// Now we have all the data we need:
// - where sample originates from (so destination when mapping):
// patchFaces, patchFaceProcs.
// - cell sample is in (so source when mapping)
// sampleCells, sampleCellProcs.
// - cell/face sample is in (so source when mapping)
// sampleIndices, sampleProcs.
if (debug && Pstream::master())
@ -267,18 +408,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
{
meshTools::writeOBJ(str, patchFc[i]);
vertI++;
meshTools::writeOBJ(str, sampleCc[i]);
meshTools::writeOBJ(str, sampleLocations[i]);
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
// Check that actual offset vector (sampleCc - patchFc) is more or less
// constant.
// Check that actual offset vector (sampleLocations - patchFc) is more or
// less constant.
if (Pstream::master())
{
const scalarField magOffset(mag(sampleCc - patchFc));
const scalarField magOffset(mag(sampleLocations - patchFc));
const scalar avgOffset(average(magOffset));
forAll(magOffset, sampleI)
@ -291,7 +432,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
<< " on a single plane."
<< " This might give numerical problems." << endl
<< " At patchface " << patchFc[sampleI]
<< " the sampled cell " << sampleCc[sampleI] << endl
<< " the sampled cell " << sampleLocations[sampleI] << endl
<< " is not on a plane " << avgOffset
<< " offset from the patch." << endl
<< " You might want to shift your plane offset."
@ -304,7 +445,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
// Determine schedule.
mapDistribute distMap(sampleCellProcs, patchFaceProcs);
mapDistribute distMap(sampleProcs, patchFaceProcs);
// Rework the schedule to cell data to send, face data to receive.
schedulePtr_.reset(new List<labelPair>(distMap.schedule()));
@ -313,17 +454,21 @@ void Foam::directMappedPolyPatch::calcMapping() const
const labelListList& constructMap = distMap.constructMap();
// Extract the particular data I need to send and receive.
sendCellLabelsPtr_.reset(new labelListList(subMap.size()));
labelListList& sendCellLabels = sendCellLabelsPtr_();
sendLabelsPtr_.reset(new labelListList(subMap.size()));
labelListList& sendLabels = sendLabelsPtr_();
forAll(subMap, procI)
{
sendCellLabels[procI] = IndirectList<label>(sampleCells, subMap[procI]);
sendLabels[procI] = IndirectList<label>
(
sampleIndices,
subMap[procI]
);
if (debug)
{
Pout<< "To proc:" << procI << " sending values of cells:"
<< sendCellLabels[procI] << endl;
Pout<< "To proc:" << procI << " sending values of cells/faces:"
<< sendLabels[procI] << endl;
}
}
@ -356,9 +501,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
)
:
polyPatch(name, size, start, index, bm),
mode_(NEARESTPATCHFACE),
samplePatch_("none"),
offset_(vector::zero),
schedulePtr_(NULL),
sendCellLabelsPtr_(NULL),
sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL)
{}
@ -372,9 +519,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
)
:
polyPatch(name, dict, index, bm),
mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
samplePatch_(dict.lookup("samplePatch")),
offset_(dict.lookup("offset")),
schedulePtr_(NULL),
sendCellLabelsPtr_(NULL),
sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL)
{}
@ -386,9 +535,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
)
:
polyPatch(pp, bm),
mode_(pp.mode_),
samplePatch_(pp.samplePatch_),
offset_(pp.offset_),
schedulePtr_(NULL),
sendCellLabelsPtr_(NULL),
sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL)
{}
@ -403,9 +554,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
)
:
polyPatch(pp, bm, index, newSize, newStart),
mode_(pp.mode_),
samplePatch_(pp.samplePatch_),
offset_(pp.offset_),
schedulePtr_(NULL),
sendCellLabelsPtr_(NULL),
sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL)
{}
@ -423,7 +576,7 @@ Foam::directMappedPolyPatch::~directMappedPolyPatch()
void Foam::directMappedPolyPatch::clearOut()
{
schedulePtr_.clear();
sendCellLabelsPtr_.clear();
sendLabelsPtr_.clear();
receiveFaceLabelsPtr_.clear();
}
@ -431,6 +584,10 @@ void Foam::directMappedPolyPatch::clearOut()
void Foam::directMappedPolyPatch::write(Ostream& os) const
{
polyPatch::write(os);
os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
<< token::END_STATEMENT << nl;
os.writeKeyword("samplePatch") << samplePatch_
<< token::END_STATEMENT << nl;
os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
}

84
src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
View File

@ -26,8 +26,8 @@ Class
Foam::directMappedPolyPatch
Description
Determines a mapping between patch face centres and mesh cell centres and
processors they're on.
Determines a mapping between patch face centres and mesh cell or face
centres and processors they're on.
Note
Storage is not optimal. It stores all face centres and cells on all
@ -43,6 +43,9 @@ SourceFiles
#include "polyPatch.H"
#include "labelPair.H"
#include "Tuple2.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,19 +60,41 @@ class directMappedPolyPatch
:
public polyPatch
{
public:
//- Mesh items to sample
enum sampleMode
{
NEARESTCELL,
NEARESTPATCHFACE,
NEARESTFACE
};
private:
// Private data
static const NamedEnum<sampleMode, 3> sampleModeNames_;
//- What to sample
sampleMode mode_;
//- Patch (only if NEARESTBOUNDARY)
const word samplePatch_;
//- Offset vector
const vector offset_;
// Derived information
//- Communication schedule.
mutable autoPtr<List<labelPair> > schedulePtr_;
//- Cells to sample per processor
mutable autoPtr<labelListList> sendCellLabelsPtr_;
//- Cells/faces to sample per processor
mutable autoPtr<labelListList> sendLabelsPtr_;
//- Patch faces to receive per processor
mutable autoPtr<labelListList> receiveFaceLabelsPtr_;
@ -86,19 +111,50 @@ class directMappedPolyPatch
pointField& patchFc
) const;
//- Find cells containing samples
//- Find cells/faces containing samples
void findSamples
(
const pointField&,
labelList& sampleCellProcs,
labelList& sampleCells,
pointField& sampleCc
labelList& sampleProcs, // processor containing sample
labelList& sampleIndices, // local index of cell/face
pointField& sampleLocations // actual representative location
) const;
//- Calculate matching
void calcMapping() const;
// Private class
//- Private class for finding nearest
// - point+local index
// - sqr(distance)
// - processor
typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
class nearestEqOp
{
public:
void operator()(nearInfo& x, const nearInfo& y) const
{
if (y.first().hit())
{
if (!x.first().hit())
{
x = y;
}
else if (y.second().first() < x.second().first())
{
x = y;
}
}
}
};
protected:
void clearOut();
@ -231,21 +287,25 @@ public:
return schedulePtr_();
}
//- Cells to sample per processor
const labelListList& sendCellLabels() const
//- Cells/faces to sample per processor
const labelListList& sendLabels() const
{
if (!sendCellLabelsPtr_.valid())
Pout<< "Asking for sendLabels." << endl;
if (!sendLabelsPtr_.valid())
{
Pout<< "Calculating mapping." << endl;
calcMapping();
}
return sendCellLabelsPtr_();
return sendLabelsPtr_();
}
//- Patch faces to receive per processor
const labelListList& receiveFaceLabels() const
{
Pout<< "Asking for receiveFaceLabels." << endl;
if (!receiveFaceLabelsPtr_.valid())
{
Pout<< "Calculating mapping." << endl;
calcMapping();
}
return receiveFaceLabelsPtr_();

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

@ -26,7 +26,7 @@ License
#include "indexedOctree.H"
#include "linePointRef.H"
#include "triSurface.H"
//#include "triSurface.H"
#include "meshTools.H"
#include "OFstream.H"

6
src/meshTools/indexedOctree/treeDataCell.C
View File

@ -26,7 +26,7 @@ License
#include "treeDataCell.H"
#include "indexedOctree.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -71,7 +71,7 @@ Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
Foam::treeDataCell::treeDataCell
(
const bool cacheBb,
const polyMesh& mesh,
const primitiveMesh& mesh,
const labelList& cellLabels
)
:
@ -94,7 +94,7 @@ Foam::treeDataCell::treeDataCell
Foam::treeDataCell::treeDataCell
(
const bool cacheBb,
const polyMesh& mesh
const primitiveMesh& mesh
)
:
mesh_(mesh),

10
src/meshTools/indexedOctree/treeDataCell.H
View File

@ -47,7 +47,7 @@ namespace Foam
{
// Forward declaration of classes
class polyMesh;
class primitiveMesh;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
@ -58,7 +58,7 @@ class treeDataCell
{
// Private data
const polyMesh& mesh_;
const primitiveMesh& mesh_;
//- Subset of cells to work on
const labelList cellLabels_;
@ -87,12 +87,12 @@ public:
treeDataCell
(
const bool cacheBb,
const polyMesh&,
const primitiveMesh&,
const labelList&
);
//- Construct from mesh. Uses all cells in mesh.
treeDataCell(const bool cacheBb, const polyMesh&);
treeDataCell(const bool cacheBb, const primitiveMesh&);
// Member Functions
@ -104,7 +104,7 @@ public:
return cellLabels_;
}
const polyMesh& mesh() const
const primitiveMesh& mesh() const
{
return mesh_;
}

409
src/meshTools/meshSearch/meshSearch.C
View File

@ -22,14 +22,11 @@ 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 "polyMesh.H"
#include "meshSearch.H"
#include "triPointRef.H"
#include "octree.H"
#include "polyMesh.H"
#include "indexedOctree.H"
#include "pointIndexHit.H"
#include "DynamicList.H"
#include "demandDrivenData.H"
@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::meshSearch::findNearer
(
const point& sample,
const pointField& points,
label& nearestI,
scalar& nearestDistSqr
)
{
bool nearer = false;
forAll(points, pointI)
{
scalar distSqr = magSqr(points[pointI] - sample);
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
nearestI = pointI;
nearer = true;
}
}
return nearer;
}
bool Foam::meshSearch::findNearer
(
const point& sample,
const pointField& points,
const labelList& indices,
label& nearestI,
scalar& nearestDistSqr
)
{
bool nearer = false;
forAll(indices, i)
{
label pointI = indices[i];
scalar distSqr = magSqr(points[pointI] - sample);
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
nearestI = pointI;
nearer = true;
}
}
return nearer;
}
// tree based searching
Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
{
treeBoundBox tightest(treeBoundBox::greatBox);
const indexedOctree<treeDataPoint>& tree = cellCentreTree();
scalar tightestDist = GREAT;
scalar span = mag(tree.bb().max() - tree.bb().min());
return cellCentreTree().findNearest(location, tightest, tightestDist);
pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
if (!info.hit())
{
info = tree.findNearest(location, Foam::sqr(GREAT));
}
return info.index();
}
@ -59,21 +117,18 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
{
const vectorField& centres = mesh_.cellCentres();
label nearestCelli = 0;
scalar minProximity = magSqr(centres[0] - location);
label nearestIndex = 0;
scalar minProximity = magSqr(centres[nearestIndex] - location);
forAll(centres, celli)
{
scalar proximity = magSqr(centres[celli] - location);
findNearer
(
location,
centres,
nearestIndex,
minProximity
);
if (proximity < minProximity)
{
nearestCelli = celli;
minProximity = proximity;
}
}
return nearestCelli;
return nearestIndex;
}
@ -92,40 +147,145 @@ Foam::label Foam::meshSearch::findNearestCellWalk
) << "illegal seedCell:" << seedCellI << exit(FatalError);
}
const vectorField& centres = mesh_.cellCentres();
const labelListList& cc = mesh_.cellCells();
// Walk in direction of face that decreases distance
label curCell = seedCellI;
scalar distanceSqr = magSqr(centres[curCell] - location);
label curCellI = seedCellI;
scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
bool closer;
do
{
closer = false;
// set the current list of neighbouring cells
const labelList& neighbours = cc[curCell];
forAll (neighbours, nI)
{
scalar curDistSqr = magSqr(centres[neighbours[nI]] - location);
// search through all the neighbours.
// If the cell is closer, reset current cell and distance
if (curDistSqr < distanceSqr)
{
distanceSqr = curDistSqr;
curCell = neighbours[nI];
closer = true; // a closer neighbour has been found
}
}
// Try neighbours of curCellI
closer = findNearer
(
location,
mesh_.cellCentres(),
mesh_.cellCells()[curCellI],
curCellI,
distanceSqr
);
} while (closer);
return curCell;
return curCellI;
}
// tree based searching
Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
{
// Search nearest cell centre.
const indexedOctree<treeDataPoint>& tree = cellCentreTree();
scalar span = mag(tree.bb().max() - tree.bb().min());
// Search with decent span
pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
if (!info.hit())
{
// Search with desparate span
info = tree.findNearest(location, Foam::sqr(GREAT));
}
// Now check any of the faces of the nearest cell
const vectorField& centres = mesh_.faceCentres();
const cell& ownFaces = mesh_.cells()[info.index()];
label nearestFaceI = ownFaces[0];
scalar minProximity = magSqr(centres[nearestFaceI] - location);
findNearer
(
location,
centres,
ownFaces,
nearestFaceI,
minProximity
);
return nearestFaceI;
}
// linear searching
Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
{
const vectorField& centres = mesh_.faceCentres();
label nearestFaceI = 0;
scalar minProximity = magSqr(centres[nearestFaceI] - location);
findNearer
(
location,
centres,
nearestFaceI,
minProximity
);
return nearestFaceI;
}
// walking from seed
Foam::label Foam::meshSearch::findNearestFaceWalk
(
const point& location,
const label seedFaceI
) const
{
if (seedFaceI < 0)
{
FatalErrorIn
(
"meshSearch::findNearestFaceWalk(const point&, const label)"
) << "illegal seedFace:" << seedFaceI << exit(FatalError);
}
const vectorField& centres = mesh_.faceCentres();
// Walk in direction of face that decreases distance
label curFaceI = seedFaceI;
scalar distanceSqr = magSqr(centres[curFaceI] - location);
while (true)
{
label betterFaceI = curFaceI;
findNearer
(
location,
centres,
mesh_.cells()[mesh_.faceOwner()[curFaceI]],
betterFaceI,
distanceSqr
);
if (mesh_.isInternalFace(curFaceI))
{
findNearer
(
location,
centres,
mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
betterFaceI,
distanceSqr
);
}
if (betterFaceI == curFaceI)
{
break;
}
curFaceI = betterFaceI;
}
return curFaceI;
}
@ -180,12 +340,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
const face& f = mesh_.faces()[curFaceI];
scalar minDist =
f.nearestPoint
(
location,
mesh_.points()
).distance();
scalar minDist = f.nearestPoint
(
location,
mesh_.points()
).distance();
bool closer;
@ -202,8 +361,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
forAll (myEdges, myEdgeI)
{
const labelList& neighbours =
mesh_.edgeFaces()[myEdges[myEdgeI]];
const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
// Check any face which uses edge, is boundary face and
// is not curFaceI itself.
@ -220,12 +378,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
{
const face& f = mesh_.faces()[faceI];
pointHit curHit =
f.nearestPoint
(
location,
mesh_.points()
);
pointHit curHit = f.nearestPoint
(
location,
mesh_.points()
);
// If the face is closer, reset current face and distance
if (curHit.distance() < minDist)
@ -283,9 +440,11 @@ Foam::meshSearch::~meshSearch()
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
const
{
if (!boundaryTreePtr_)
{
@ -294,17 +453,28 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
//
// all boundary faces (not just walls)
octreeDataFace shapes(mesh_);
labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(bndFaces, i)
{
bndFaces[i] = mesh_.nInternalFaces() + i;
}
treeBoundBox overallBb(mesh_.points());
boundaryTreePtr_ = new octree<octreeDataFace>
Random rndGen(123456);
boundaryTreePtr_ = new indexedOctree<treeDataFace>
(
overallBb, // overall search domain
shapes, // all information needed to do checks on cells
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
10.0
treeDataFace // all information needed to search faces
(
false, // do not cache bb
mesh_,
bndFaces // boundary faces only
),
overallBb.extend(rndGen, 1E-3), // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
}
@ -312,7 +482,8 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
}
const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
const
{
if (!cellTreePtr_)
{
@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
// Construct tree
//
octreeDataCell shapes(mesh_);
treeBoundBox overallBb(mesh_.points());
cellTreePtr_ = new octree<octreeDataCell>
cellTreePtr_ = new indexedOctree<treeDataCell>
(
treeDataCell
(
false, // not cache bb
mesh_
),
overallBb, // overall search domain
shapes, // all information needed to do checks on cells
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
2.0
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
}
@ -339,8 +512,8 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
}
const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
const
const Foam::indexedOctree<Foam::treeDataPoint>&
Foam::meshSearch::cellCentreTree() const
{
if (!cellCentreTreePtr_)
{
@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
// Construct tree
//
// shapes holds reference to cellCentres!
octreeDataPoint shapes(mesh_.cellCentres());
treeBoundBox overallBb(mesh_.cellCentres());
cellCentreTreePtr_ = new octree<octreeDataPoint>
cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
(
treeDataPoint(mesh_.cellCentres()),
overallBb, // overall search domain
shapes, // all information needed to do checks on cells
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
10, // max levels
10.0, // maximum ratio of cubes v.s. cells
100.0 // max. duplicity; n/a since no bounding boxes.
);
}
@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
forAll(f, fp)
{
pointHit inter =
f.ray
(
ctr,
dir,
mesh_.points(),
intersection::HALF_RAY,
intersection::VECTOR
);
pointHit inter = f.ray
(
ctr,
dir,
mesh_.points(),
intersection::HALF_RAY,
intersection::VECTOR
);
if (inter.hit())
{
@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell
}
Foam::label Foam::meshSearch::findNearestFace
(
const point& location,
const label seedFaceI,
const bool useTreeSearch
) const
{
if (seedFaceI == -1)
{
if (useTreeSearch)
{
return findNearestFaceTree(location);
}
else
{
return findNearestFaceLinear(location);
}
}
else
{
return findNearestFaceWalk(location, seedFaceI);
}
}
Foam::label Foam::meshSearch::findCell
(
const point& location,
@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
{
if (useTreeSearch)
{
treeBoundBox tightest(treeBoundBox::greatBox);
const indexedOctree<treeDataFace>& tree = boundaryTree();
scalar tightestDist = GREAT;
scalar span = mag(tree.bb().max() - tree.bb().min());
label index =
boundaryTree().findNearest
pointIndexHit info = boundaryTree().findNearest
(
location,
Foam::sqr(span)
);
if (!info.hit())
{
info = boundaryTree().findNearest
(
location,
tightest,
tightestDist
Foam::sqr(GREAT)
);
}
return boundaryTree().shapes().meshFaces()[index];
return tree.shapes().faceLabels()[info.index()];
}
else
{
@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
if (curHit.hit())
{
// Change index into octreeData into face label
curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]);
curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
}
return curHit;
}
@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
bool Foam::meshSearch::isInside(const point& p) const
{
return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE;
return
boundaryTree().getVolumeType(p)
== indexedOctree<treeDataFace>::INSIDE;
}

120
src/meshTools/meshSearch/meshSearch.H
View File

@ -26,7 +26,8 @@ Class
Foam::meshSearch
Description
Various searches on polyMesh; uses (demand driven) octree to search.
Various (local, not parallel) searches on polyMesh;
uses (demand driven) octree to search.
SourceFiles
meshSearch.C
@ -36,11 +37,10 @@ SourceFiles
#ifndef meshSearch_H
#define meshSearch_H
#include "octreeDataCell.H"
#include "octreeDataFace.H"
#include "octreeDataPoint.H"
#include "treeDataCell.H"
#include "treeDataFace.H"
#include "treeDataPoint.H"
#include "pointIndexHit.H"
#include "className.H"
#include "Cloud.H"
#include "passiveParticle.H"
@ -71,46 +71,77 @@ class meshSearch
//- demand driven octrees
mutable octree<octreeDataFace>* boundaryTreePtr_;
mutable octree<octreeDataCell>* cellTreePtr_;
mutable octree<octreeDataPoint>* cellCentreTreePtr_;
mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
mutable indexedOctree<treeDataCell>* cellTreePtr_;
mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
// Private Member Functions
//- nearest cell centre using octree
label findNearestCellTree(const point& location) const;
//- nearest cell centre going through all cells
label findNearestCellLinear(const point& location) const;
//- walk from seed. Does not 'go around' boundary, just returns
// last cell before boundary.
label findNearestCellWalk
//- Updates nearestI, nearestDistSqr from any closer ones.
static bool findNearer
(
const point& location,
const label seedCellI
) const;
const point& sample,
const pointField& points,
label& nearestI,
scalar& nearestDistSqr
);
//- cell containing location. Linear search.
label findCellLinear(const point& location) const;
//- walk from seed to find nearest boundary face. Gets stuck in
// local minimum.
label findNearestBoundaryFaceWalk
//- Updates nearestI, nearestDistSqr from any selected closer ones.
static bool findNearer
(
const point& location,
const label seedFaceI
) const;
const point& sample,
const pointField& points,
const labelList& indices,
label& nearestI,
scalar& nearestDistSqr
);
//- Calculate offset vector in direction dir with as length a fraction
// of the cell size (of the cell straddling boundary face)
vector offset
(
const point& bPoint,
const label bFaceI,
const vector& dir
) const;
// Cells
//- nearest cell centre using octree
label findNearestCellTree(const point&) const;
//- nearest cell centre going through all cells
label findNearestCellLinear(const point&) const;
//- walk from seed. Does not 'go around' boundary, just returns
// last cell before boundary.
label findNearestCellWalk(const point&, const label) const;
//- cell containing location. Linear search.
label findCellLinear(const point&) const;
// Cells
label findNearestFaceTree(const point&) const;
label findNearestFaceLinear(const point&) const;
label findNearestFaceWalk(const point&, const label) const;
// Boundary faces
//- walk from seed to find nearest boundary face. Gets stuck in
// local minimum.
label findNearestBoundaryFaceWalk
(
const point& location,
const label seedFaceI
) const;
//- Calculate offset vector in direction dir with as length a
// fraction of the cell size (of the cell straddling boundary face)
vector offset
(
const point& bPoint,
const label bFaceI,
const vector& dir
) const;
//- Disallow default bitwise copy construct
@ -154,13 +185,13 @@ public:
//- Get (demand driven) reference to octree holding all
// boundary faces
const octree<octreeDataFace>& boundaryTree() const;
const indexedOctree<treeDataFace>& boundaryTree() const;
//- Get (demand driven) reference to octree holding all cells
const octree<octreeDataCell>& cellTree() const;
const indexedOctree<treeDataCell>& cellTree() const;
//- Get (demand driven) reference to octree holding all cell centres
const octree<octreeDataPoint>& cellCentreTree() const;
const indexedOctree<treeDataPoint>& cellCentreTree() const;
// Queries
@ -181,6 +212,13 @@ public:
const bool useTreeSearch = true
) const;
label findNearestFace
(
const point& location,
const label seedFaceI = -1,
const bool useTreeSearch = true
) const;
//- Find cell containing (using pointInCell) location.
// If seed provided walks and falls back to linear/tree search.
// (so handles holes correctly)s
@ -206,7 +244,7 @@ public:
//- Find first intersection of boundary in segment [pStart, pEnd]
// (so inclusive of endpoints). Always octree for now
pointIndexHit intersection(const point& pStart, const point& pEnd)
const;
const;
//- Find all intersections of boundary within segment pStart .. pEnd
// Always octree for now

20
src/meshTools/meshTools/meshTools.C
View File

@ -182,7 +182,7 @@ Foam::vectorField Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)
"Foam::meshTools::calcBoxPointNormals"
"(const primitivePatch& pp)"
) << "No visible octant for point:" << pp.meshPoints()[pointI]
<< " cooord:" << pp.localPoints()[pointI] << nl
<< " cooord:" << pp.points()[pp.meshPoints()[pointI]] << nl
<< "Normal set to " << pn[pointI] << endl;
}
}
@ -299,7 +299,7 @@ bool Foam::meshTools::edgeOnCell
const label edgeI
)
{
return findIndex(mesh.edgeCells()[edgeI], cellI) != -1;
return findIndex(mesh.edgeCells(edgeI), cellI) != -1;
}
@ -310,7 +310,7 @@ bool Foam::meshTools::edgeOnFace
const label edgeI
)
{
return findIndex(mesh.faceEdges()[faceI], edgeI) != -1;
return findIndex(mesh.faceEdges(faceI), edgeI) != -1;
}
@ -403,7 +403,6 @@ Foam::label Foam::meshTools::getSharedEdge
const labelList& f0Edges = mesh.faceEdges()[f0];
const labelList& f1Edges = mesh.faceEdges()[f1];
forAll(f0Edges, f0EdgeI)
{
label edge0 = f0Edges[f0EdgeI];
@ -481,7 +480,7 @@ void Foam::meshTools::getEdgeFaces
label& face1
)
{
const labelList& eFaces = mesh.edgeFaces()[edgeI];
const labelList& eFaces = mesh.edgeFaces(edgeI);
face0 = -1;
face1 = -1;
@ -619,7 +618,7 @@ Foam::label Foam::meshTools::walkFace
const label nEdges
)
{
const labelList& fEdges = mesh.faceEdges()[faceI];
const labelList& fEdges = mesh.faceEdges(faceI);
label edgeI = startEdgeI;
@ -790,13 +789,4 @@ Foam::label Foam::meshTools::cutDirToEdge
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

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

@ -314,6 +314,7 @@ public:
//- Return slightly wider bounding box
// Extends all dimensions with s*span*Random::scalar01()
// and guarantees in any direction s*mag(span) minimum width
inline treeBoundBox extend(Random&, const scalar s) const;
// Friend Operators

10
src/meshTools/octree/treeBoundBoxI.H
View File

@ -437,7 +437,15 @@ inline treeBoundBox treeBoundBox::extend(Random& rndGen, const scalar s) const
{
treeBoundBox bb(*this);
const vector span(bb.max() - bb.min());
vector span(bb.max() - bb.min());
// Make 3D
scalar magSpan = Foam::mag(span);
for (direction dir = 0; dir < vector::nComponents; dir++)
{
span[dir] = Foam::max(s*magSpan, span[dir]);
}
bb.min() -= cmptMultiply(s*rndGen.vector01(), span);
bb.max() += cmptMultiply(s*rndGen.vector01(), span);

4
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -230,7 +230,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface(*this),
bb.extend(rndGen, 1E-3), // slightly randomize bb
bb.extend(rndGen, 1E-4), // slightly randomize bb
10, // maxLevel
10, // leafsize
3.0 // duplicity
@ -274,7 +274,7 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
localPoints(), // points
bEdges // selected edges
),
bb.extend(rndGen, 1E-3), // slightly randomize bb
bb.extend(rndGen, 1E-4), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity

2
src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
View File

@ -31,8 +31,6 @@ Description
#include "meshSearch.H"
#include "triSurface.H"
#include "triSurfaceSearch.H"
#include "octree.H"
#include "octreeDataTriSurface.H"
#include "cellClassification.H"
#include "addToRunTimeSelectionTable.H"

2
src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
View File

@ -30,8 +30,6 @@ Description
#include "polyMesh.H"
#include "meshSearch.H"
#include "triSurfaceSearch.H"
#include "octree.H"
#include "octreeDataTriSurface.H"
#include "addToRunTimeSelectionTable.H"

12
src/thermophysicalModels/combustion/mixtureThermos/mixtures/combustionMixture/combustionMixture.C
View File

@ -107,18 +107,6 @@ combustionMixture::combustionMixture
);
}
}
// Error check for single component mixtures - should be uniform 1
if (Y_.size() == 1)
{
if (average(Y_[0]).value() < 0.999)
{
FatalErrorIn("combustionMixture::combustionMixture")
<< "Mass fraction for single component mixture must be unity"
<< nl << "Please correct field: " << species_[0]
<< " (Ydefault)" << nl << abort(FatalError);
}
}
}

1
src/turbulenceModels/RAS/incompressible/Make/files
View File

@ -5,6 +5,7 @@ laminar/laminar.C
kEpsilon/kEpsilon.C
RNGkEpsilon/RNGkEpsilon.C
realizableKE/realizableKE.C
kOmega/kOmega.C
kOmegaSST/kOmegaSST.C
SpalartAllmaras/SpalartAllmaras.C
LRR/LRR.C

276
src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C Normal file
View File

@ -0,0 +1,276 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "kOmega.H"
#include "addToRunTimeSelectionTable.H"
#include "wallFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kOmega, 0);
addToRunTimeSelectionTable(RASModel, kOmega, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
kOmega::kOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& lamTransportModel
)
:
RASModel(typeName, U, phi, lamTransportModel),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"betaStar",
coeffDict_,
0.09
)
),
beta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta",
coeffDict_,
0.072
)
),
alpha_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alpha",
coeffDict_,
0.52
)
),
alphaK_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK",
coeffDict_,
0.5
)
),
alphaOmega_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega",
coeffDict_,
0.5
)
),
omega0_("omega0", dimless/dimTime, SMALL),
omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
k_
(
IOobject
(
"k",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
nut_(k_/(omega_ + omegaSmall_))
{
# include "kOmegaWallViscosityI.H"
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> kOmega::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
k_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> kOmega::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
{
return
(
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))
);
}
bool kOmega::read()
{
if (RASModel::read())
{
Cmu_.readIfPresent(coeffDict_);
beta_.readIfPresent(coeffDict_);
alphaK_.readIfPresent(coeffDict_);
alphaOmega_.readIfPresent(coeffDict_);
return true;
}
else
{
return false;
}
}
void kOmega::correct()
{
transportModel_.correct();
if (!turbulence_)
{
return;
}
RASModel::correct();
volScalarField G = nut_*2*magSqr(symm(fvc::grad(U_)));
# include "kOmegaWallFunctionsI.H"
// Turbulence specific dissipation rate equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(omega_)
+ fvm::div(phi_, omega_)
- fvm::Sp(fvc::div(phi_), omega_)
- fvm::laplacian(DomegaEff(), omega_)
==
alpha_*G*omega_/k_
- fvm::Sp(beta_*omega_, omega_)
);
omegaEqn().relax();
# include "wallOmegaI.H"
solve(omegaEqn);
bound(omega_, omega0_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::Sp(fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(Cmu_*omega_, k_)
);
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
// Re-calculate viscosity
nut_ = k_/omega_;
# include "kOmegaWallViscosityI.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

210
src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H Normal file
View File

@ -0,0 +1,210 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::incompressible::RASModels::kOmega
Description
Standard high Reynolds-number k-omega turbulence model for
incompressible flows.
References:
@verbatim
"Turbulence Modeling for CFD"
D. C. Wilcox,
DCW Industries, Inc., La Canada,
California, 1998.
See also:
http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model
@endverbatim
The default model coefficients correspond to the following:
@verbatim
kOmegaCoeffs
{
Cmu 0.09; // Equivalent to betaStar
alpha 0.52;
beta 0.072;
alphak 0.5;
alphaOmega 0.5;
}
@endverbatim
SourceFiles
kOmega.C
\*---------------------------------------------------------------------------*/
#ifndef kOmega_H
#define kOmega_H
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kOmega Declaration
\*---------------------------------------------------------------------------*/
class kOmega
:
public RASModel
{
// Private data
// Model coefficients
dimensionedScalar Cmu_;
dimensionedScalar beta_;
dimensionedScalar alpha_;
dimensionedScalar alphaK_;
dimensionedScalar alphaOmega_;
dimensionedScalar omega0_;
dimensionedScalar omegaSmall_;
// Fields
volScalarField k_;
volScalarField omega_;
volScalarField nut_;
public:
//- Runtime type information
TypeName("kOmega");
// Constructors
//- Construct from components
kOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport
);
// Destructor
~kOmega()
{}
// Member Functions
//- Return the turbulence viscosity
tmp<volScalarField> nut() const
{
return nut_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff() const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK_*nut_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff() const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaOmega_*nut_ + nu())
);
}
//- Return the turbulence kinetic energy
tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence specific dissipation rate
tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the turbulence kinetic energy dissipation rate
tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
mesh_.time().timeName(),
mesh_
),
Cmu_*k_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the Reynolds stress tensor
tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar stress
tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
//- Solve the turbulence equations and correct the turbulence viscosity
void correct();
//- Read turbulenceProperties dictionary
bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

125
src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallFunctionsI.H Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
kOmegaWallFunctions
Description
Calculate wall generation and frequency omega from wall-functions.
\*---------------------------------------------------------------------------*/
{
labelList cellBoundaryFaceCount(omega_.size(), 0);
scalar Cmu25 = pow(Cmu_.value(), 0.25);
const fvPatchList& patches = mesh_.boundary();
//- Initialise the near-wall omega and G fields to zero
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
omega_[faceCelli] = 0.0;
G[faceCelli] = 0.0;
}
}
}
//- Accumulate the wall face contributions to omega and G
// Increment cellBoundaryFaceCount for each face for averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
# include "checkkOmegaPatchFieldTypes.H"
const scalarField& nuw = nu().boundaryField()[patchi];
const scalarField& nutw = nut_.boundaryField()[patchi];
scalarField magFaceGradU =
mag(U_.boundaryField()[patchi].snGrad());
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y_[patchi][facei]
*sqrt(k_[faceCelli])
/nuw[facei];
// For corner cells (with two boundary or more faces),
// omega and G in the near-wall cell are calculated
// as an average
cellBoundaryFaceCount[faceCelli]++;
omega_[faceCelli] +=
sqrt(k_[faceCelli])
/(Cmu25*kappa_.value()*y_[patchi][facei]);
if (yPlus > yPlusLam_)
{
G[faceCelli] +=
(nutw[facei] + nuw[facei])
*magFaceGradU[facei]
*Cmu25*sqrt(k_[faceCelli])
/(kappa_.value()*y_[patchi][facei]);
}
}
}
}
// Perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
omega_[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
}
}
}
}
// ************************************************************************* //

70
src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallViscosityI.H Normal file
View File

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
kOmegaWallViscosity
Description
Calculate wall viscosity from wall-functions.
\*---------------------------------------------------------------------------*/
{
scalar Cmu25 = pow(Cmu_.value(), 0.25);
const fvPatchList& patches = mesh_.boundary();
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
const scalarField& nuw = nu().boundaryField()[patchi];
scalarField& nutw = nut_.boundaryField()[patchi];
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei];
if (yPlus > yPlusLam_)
{
nutw[facei] =
nuw[facei]
*(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1);
}
else
{
nutw[facei] = 0.0;
}
}
}
}
}
// ************************************************************************* //

51
src/turbulenceModels/RAS/incompressible/kOmega/wallOmegaI.H Normal file
View File

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
wallOmega
Description
Set wall dissipation in the omega matrix
\*---------------------------------------------------------------------------*/
{
const fvPatchList& patches = mesh_.boundary();
forAll(patches, patchi)
{
const fvPatch& p = patches[patchi];
if (isType<wallFvPatch>(p))
{
omegaEqn().setValues
(
p.faceCells(),
omega_.boundaryField()[patchi].patchInternalField()
);
}
}
}
// ************************************************************************* //

1
src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.H
View File

@ -235,6 +235,7 @@ public:
return k_;
}
//- Return the turbulence specific dissipation rate
tmp<volScalarField> omega() const
{
return omega_;

12
tutorials/boundaryFoam/boundaryWallFunctions/constant/polyMesh/blockMeshDict
View File

@ -16,7 +16,7 @@ FoamFile
convertToMeters 0.05;
vertices
vertices
(
(0 -1 0)
(0 0 0)
@ -32,24 +32,24 @@ vertices
(0.1 1 0.1)
);
blocks
blocks
(
hex (0 3 4 1 6 9 10 7) (1 40 1) simpleGrading (1 1 1)
hex (1 4 5 2 7 10 11 8) (1 40 1) simpleGrading (1 1 1)
);
edges
edges
(
);
patches
patches
(
wall lowerWall
wall lowerWall
(
(0 3 9 6)
)
wall upperWall
wall upperWall
(
(2 8 11 5)
)

4
tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
View File

@ -8,7 +8,7 @@
FoamFile
{
version 2.0;
format ascii;
`format' ascii;
class dictionary;
object blockMeshDict;
}
@ -16,7 +16,7 @@ FoamFile
// General m4 macros
changecom(//)changequote([,])
define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; print ($1)')])
define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; printf ($1)')])
define(VCOUNT, 0)
define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])

5
tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary
View File

@ -20,8 +20,9 @@ inlet
{
nFaces 30;
startFace 27238;
type directMappedPatch;
offset ( 0.05 0 0 );
offset ( 0.0495 0 0 );
sampleMode nearestCell;
samplePatch none;
}
outlet

2
tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
View File

@ -23,6 +23,8 @@ dictionaryReplacement
{
type directMappedPatch;
offset (0.0495 0 0);
sampleMode nearestCell;
samplePatch none;
}
}
}