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

Merge branch 'master' into master-copy

Browse Source
This commit is contained in:
sergio
2010-09-27 13:09:54 +01:00
parent 29cad82481 eb5f2954b9
commit a7b8be7a42
379 changed files with 14965 additions and 3392 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/combustion/dieselEngineFoam/Make/options
View File

@ -2,6 +2,7 @@ EXE_INC = \
-I../engineFoam \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/dieselSpray/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
@ -24,6 +25,7 @@ EXE_LIBS = \
-lreactionThermophysicalModels \
-lfiniteVolume \
-llagrangian \
-lmeshTools \
-ldieselSpray \
-lliquids \
-lliquidMixture \

2
applications/solvers/combustion/dieselFoam/Make/options
View File

@ -3,6 +3,7 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/dieselSpray/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
@ -21,6 +22,7 @@ EXE_LIBS = \
-lcompressibleLESModels \
-lreactionThermophysicalModels \
-llagrangian \
-lmeshTools \
-ldieselSpray \
-lliquids \
-lliquidMixture \

1
applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options
View File

@ -9,4 +9,3 @@ EXE_LIBS = \
-lfiniteVolume \
-llagrangian \
-ldsmc

3
applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/Make/files Normal file
View File

@ -0,0 +1,3 @@
incompressibleUncoupledKinematicParcelDyMFoam.C
EXE = $(FOAM_APPBIN)/incompressibleUncoupledKinematicParcelDyMFoam

33
applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/Make/options Normal file
View File

@ -0,0 +1,33 @@
EXE_INC = \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-llagrangian \
-llagrangianIntermediate \
-lthermophysicalFunctions \
-lbasicThermophysicalModels \
-lspecie \
-lradiation \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools \
-lsurfaceFilmModels \
-ldynamicMesh \
-ldynamicFvMesh \
-ltopoChangerFvMesh

147
applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/createFields.H Normal file
View File

@ -0,0 +1,147 @@
Info<< "\nReading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar rhoInfValue
(
transportProperties.lookup("rhoInf")
);
volScalarField rhoInf
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
rhoInfValue
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating turbulence model\n" << endl;
singlePhaseTransportModel laminarTransport(U, phi);
const volScalarField nu = laminarTransport.nu();
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
volScalarField mu
(
IOobject
(
"mu",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
nu*rhoInfValue
);
word kinematicCloudName("kinematicCloud");
args.optionReadIfPresent("cloudName", kinematicCloudName);
Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
basicKinematicCloud kinematicCloud
(
kinematicCloudName,
rhoInf,
U,
mu,
g
);
IOobject Hheader
(
"H",
runTime.timeName(),
mesh,
IOobject::NO_READ
);
autoPtr<volVectorField> HPtr_;
if (Hheader.headerOk())
{
Info<< "\nReading field H\n" << endl;
HPtr_.reset
(
new volVectorField
(
IOobject
(
"H",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
);
}
IOobject HdotGradHheader
(
"HdotGradH",
runTime.timeName(),
mesh,
IOobject::NO_READ
);
autoPtr<volVectorField> HdotGradHPtr_;
if (HdotGradHheader.headerOk())
{
Info<< "Reading field HdotGradH" << endl;
HdotGradHPtr_.reset
(
new volVectorField
(
IOobject
(
"HdotGradH",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
);
}

92
applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C Normal file
View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Application
uncoupledKinematicParcelFoam
Description
Transient solver for the passive transport of a single kinematic
particle could.
Uses a pre-calculated velocity field to evolve the cloud.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "basicKinematicCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addOption
(
"cloudName",
"name",
"specify alternative cloud name. default is 'kinematicCloud'"
);
#include "setRootCase.H"
#include "createTime.H"
# include "createDynamicFvMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
mesh.update();
U.correctBoundaryConditions();
Info<< "Evolving " << kinematicCloud.name() << endl;
laminarTransport.correct();
mu = nu*rhoInfValue;
kinematicCloud.evolve();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

1
applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
View File

@ -26,6 +26,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
- fvm::laplacian(turbulence->muEff(), Yi)
==
parcels.Srho(i)
+ surfaceFilm.Srho(i)
+ kappa*chemistry.RR(i)().dimensionedInternalField(),
mesh.solver("Yi")
);

7
applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
View File

@ -6,9 +6,10 @@
- fvm::laplacian(turbulence->alphaEff(), hs)
==
DpDt
+ parcels.Sh()
+ radiation->Shs(thermo)
+ chemistrySh
+ parcels.Sh()
+ surfaceFilm.Sh()
+ radiation->Shs(thermo)
+ chemistrySh
);
hsEqn.relax();

2
applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
View File

@ -24,6 +24,7 @@ if (transonic)
- fvm::laplacian(rho*rUA, p)
==
parcels.Srho()
+ surfaceFilm.Srho()
);
pEqn.solve();
@ -52,6 +53,7 @@ else
- fvm::laplacian(rho*rUA, p)
==
parcels.Srho()
+ surfaceFilm.Srho()
);
pEqn.solve();

1
applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H
View File

@ -36,6 +36,7 @@ Description
+ fvc::div(phi)
==
parcels.Srho()
+ surfaceFilm.Srho()
);
}

61
applications/test/findCell-octree/findCell-octree.C
View File

@ -29,6 +29,8 @@ License
#include "IStringStream.H"
#include "octree.H"
#include "octreeDataCell.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "OFstream.H"
using namespace Foam;
@ -44,11 +46,13 @@ int main(int argc, char *argv[])
# include "createTime.H"
# include "createMesh.H"
label nReps = 100000;
const point sample = args.argRead<point>(1);
treeBoundBox meshBb(mesh.points());
treeBoundBox meshBb(mesh.bounds());
// Calculate typical cell releated size to shift bb by.
// Calculate typical cell related size to shift bb by.
scalar typDim = meshBb.avgDim()/(2.0*Foam::cbrt(scalar(mesh.nCells())));
treeBoundBox shiftedBb
@ -57,16 +61,13 @@ int main(int argc, char *argv[])
meshBb.max() + vector(typDim, typDim, typDim)
);
Info<< "Mesh" << endl;
Info<< " bounding box : " << meshBb << endl;
Info<< " bounding box (shifted) : " << shiftedBb << endl;
Info<< " typical dimension : " << shiftedBb.typDim() << endl;
/*
* Now we have allBb and shiftedBb
*/
Info<< "Initialised mesh in "
<< runTime.cpuTimeIncrement() << " s" << endl;
// Wrap indices and mesh information into helper object
octreeDataCell shapes(mesh);
@ -80,14 +81,52 @@ int main(int argc, char *argv[])
10.0 // maximum ratio of cubes v.s. cells
);
Info<< "Point:" << sample << " is in shape "
<< oc.find(sample) << nl
<< "Point:" << sample << " is in cell "
<< mesh.findCell(sample) << endl;
for (label i = 0; i < nReps - 1 ; i++)
{
oc.find(sample);
}
Info<< "Point:" << sample << " is in shape "
<< oc.find(sample) << endl;
oc.printStats(Info);
Info<< "Found in octree " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
indexedOctree<treeDataCell> ioc
(
treeDataCell(true, mesh),
shiftedBb,
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
for (label i = 0; i < nReps - 1 ; i++)
{
ioc.findInside(sample);
}
Info<< "Point:" << sample << " is in shape "
<< ioc.findInside(sample)
<< ", where the possible cells were:" << nl
<< ioc.findIndices(sample)
<< endl;
Info<< "Found in indexedOctree " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
for (label i = 0; i < nReps - 1 ; i++)
{
mesh.findCell(sample);
}
Info<< "Point:" << sample << " is in cell "
<< mesh.findCell(sample) << endl;
Info<< "Found in mesh.findCell " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
Info<< "End\n" << endl;

227
applications/utilities/mesh/conversion/ansysToFoam/ansysToFoam.L
View File

@ -63,11 +63,16 @@ label maxNodei = 0;
SLPtrList<labelList> slCellLabels;
SLList<label> slCellMap;
SLList<label> slCellType;
label maxCelli = 0;
PtrList<SLList<label> > slPatchCells;
PtrList<SLList<label> > slPatchCellFaces;
// Cell types
Map<word> cellTypes;
label currentTypei = -1;
// Dummy yywrap to keep yylex happy at compile time.
// It is called by yylex but is not used as the mechanism to change file.
@ -108,6 +113,8 @@ value {floatNum}
node ^{space}"N"{cspace}
element ^{space}"EN"{cspace}
bface ^{space}"SFE"{cspace}
elementTypeName ^{space}"ET"{cspace}
elementType ^{space}"TYPE"{cspace}
%%
@ -160,6 +167,7 @@ bface ^{space}"SFE"{cspace}
slCellMap.append(celli);
slCellLabels.append(new labelList(labels));
slCellType.append(currentTypei);
}
@ -210,6 +218,37 @@ bface ^{space}"SFE"{cspace}
}
{elementTypeName}{label}{cspace}{identifier}{space}\n {
IStringStream elementStream(YYText());
char tag,c;
label cellTypei;
word cellTypeName;
elementStream
>> tag >> tag // skip 'ET'
>> c >> cellTypei
>> c >> cellTypeName;
Info<< "Read typeName " << cellTypeName
<< " for type " << cellTypei << endl;
cellTypes.insert(cellTypei, cellTypeName);
}
{elementType}{label}{space}\n {
IStringStream elementStream(YYText());
char tag,c;
label cellTypei;
elementStream
>> tag >> tag >> tag >> tag // skip 'TYPE'
>> c >> cellTypei;
currentTypei = cellTypei;
}
/* ------------------------------------------------------------------------- *\
------ Ignore remaining space and \n s. Any other characters are errors.
\* ------------------------------------------------------------------------- */
@ -231,6 +270,29 @@ bface ^{space}"SFE"{cspace}
#include <fstream>
using std::ifstream;
label findFace(const polyMesh& mesh, const face& f)
{
const labelList& pFaces = mesh.pointFaces()[f[0]];
forAll(pFaces, i)
{
label faceI = pFaces[i];
if (mesh.faces()[faceI] == f)
{
return faceI;
}
}
FatalErrorIn("findFace(const polyMesh&, const face&)")
<< "Cannot find a face matching " << f
<< exit(FatalError);
return -1;
}
int main(int argc, char *argv[])
{
argList::noParallel();
@ -253,7 +315,7 @@ int main(int argc, char *argv[])
# include "createTime.H"
const fileName ansysFile = args[1];
fileName ansysFile(args.additionalArgs()[0]);
ifstream ansysStream(ansysFile.c_str());
if (!ansysStream)
@ -377,6 +439,34 @@ int main(int argc, char *argv[])
}
}
const word defaultFacesName = "defaultFaces";
word defaultFacesType = emptyPolyPatch::typeName;
// Create dummy mesh just to find out what are internal/external
// faces
autoPtr<polyMesh> dummyMesh
(
new polyMesh
(
IOobject
(
"dummyMesh",
runTime.constant(),
runTime
),
xferCopy(points),
cellShapes,
faceListList(0),
wordList(0),
wordList(0),
defaultFacesName,
defaultFacesType,
wordList(0)
)
);
// Warning: tet face order has changed between version 1.9.6 and 2.0
//
label faceIndex[7][6] =
@ -423,10 +513,53 @@ int main(int argc, char *argv[])
boundary[patchI] = patchFaces;
patchNames[patchI] = word("patch") + name(patchI + 1);
}
//
// Lookup the face labels for all the boundary faces
//
labelListList boundaryFaceLabels(boundary.size());
forAll(boundary, patchI)
{
const faceList& bFaces = boundary[patchI];
labelList& bFaceLabels = boundaryFaceLabels[patchI];
bFaceLabels.setSize(bFaces.size());
forAll(bFaces, i)
{
bFaceLabels[i] = findFace(dummyMesh(), bFaces[i]);
}
}
// Now split the boundary faces into external and internal faces. All
// faces go into faceZones and external faces go into patches.
List<faceList> patchFaces(slPatchCells.size());
labelList patchNFaces(slPatchCells.size(), 0);
forAll(boundary, patchI)
{
const faceList& bFaces = boundary[patchI];
const labelList& bFaceLabels = boundaryFaceLabels[patchI];
patchFaces[patchI].setSize(bFaces.size());
forAll(bFaces, i)
{
if (!dummyMesh().isInternalFace(bFaceLabels[i]))
{
patchFaces[patchI][patchNFaces[patchI]++] = bFaces[i];
}
}
patchFaces[patchI].setSize(patchNFaces[patchI]);
Info<< "Patch " << patchI << " named " << patchNames[patchI]
<< ": " << boundary[patchI].size() << " faces" << endl;
}
// We no longer need the dummyMesh
dummyMesh.clear();
Info<< "ansysToFoam: " << endl
<< "Ansys file format does not provide information about the type of "
<< "the patch (eg. wall, symmetry plane, cyclic etc)." << endl
@ -434,10 +567,8 @@ int main(int argc, char *argv[])
<< "as type patch. Please reset after mesh conversion as necessary."
<< endl;
wordList patchTypes(boundary.size(), polyPatch::typeName);
word defaultFacesName = "defaultFaces";
word defaultFacesType = emptyPolyPatch::typeName;
wordList patchPhysicalTypes(boundary.size());
wordList patchTypes(patchFaces.size(), polyPatch::typeName);
wordList patchPhysicalTypes(patchFaces.size());
preservePatchTypes
(
@ -461,7 +592,7 @@ int main(int argc, char *argv[])
),
xferMove(points),
cellShapes,
boundary,
patchFaces,
patchNames,
patchTypes,
defaultFacesName,
@ -469,6 +600,90 @@ int main(int argc, char *argv[])
patchPhysicalTypes
);
if (cellTypes.size() > 0 || patchNames.size() > 0)
{
DynamicList<pointZone*> pz;
DynamicList<faceZone*> fz;
DynamicList<cellZone*> cz;
// FaceZones
forAll(boundaryFaceLabels, patchI)
{
if (boundaryFaceLabels[patchI].size())
{
// Re-do the boundaryFaceLabels since the boundary face
// labels will be different on the pShapeMesh.
const faceList& bFaces = boundary[patchI];
labelList& bFaceLabels = boundaryFaceLabels[patchI];
forAll(bFaceLabels, i)
{
bFaceLabels[i] = findFace(pShapeMesh, bFaces[i]);
}
Info<< "Creating faceZone " << patchNames[patchI]
<< " with " << bFaceLabels.size() << " faces" << endl;
fz.append
(
new faceZone
(
patchNames[patchI],
bFaceLabels,
boolList(bFaceLabels.size(), false),
fz.size(),
pShapeMesh.faceZones()
)
);
}
}
// CellZones
labelList types = cellTypes.sortedToc();
forAll(types, j)
{
label cellType = types[j];
// Pick up cells in zone
DynamicList<label> addr;
SLList<label>::iterator cellMapIter = slCellMap.begin();
SLList<label>::iterator typeIter = slCellType.begin();
for
(
;
typeIter != slCellType.end();
++typeIter, ++cellMapIter
)
{
if (typeIter() == cellType)
{
addr.append(cellMap[cellMapIter()]);
}
}
Info<< "Creating cellZone " << cellTypes[cellType]
<< " with " << addr.size() << " cells" << endl;
cz.append
(
new cellZone
(
cellTypes[cellType],
addr,
j,
pShapeMesh.cellZones()
)
);
}
pShapeMesh.addZones(pz, fz, cz);
}
// Set the precision of the points data to 10
IOstream::defaultPrecision(10);

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

@ -341,11 +341,14 @@ meshQualityControls
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum tet volume. Is absolute volume of the tet formed by the
// face-centre decomposition triangle and the cell centre.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minTetVol 1e-20;
//- Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
// the cell centre. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
// 1 = regular tet
minTetQuality 1e-9;
//- Minimum face area. Set to <0 to disable.
minArea -1;

52
applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
View File

@ -6,6 +6,7 @@
#include "EdgeMap.H"
#include "wedgePolyPatch.H"
#include "unitConversion.H"
#include "polyMeshTetDecomposition.H"
// Find wedge with opposite orientation. Note: does not actually check that
@ -83,7 +84,7 @@ bool Foam::checkWedges
{
Info<< " Wedge " << pp.name() << " with angle "
<< radToDeg(wedgeAngle) << " degrees"
<< endl;
<< endl;
}
// Find opposite
@ -413,7 +414,6 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
{
faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFacePyramids(true, -SMALL, &faces))
@ -433,25 +433,6 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
{
faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100+1);
if (mesh.checkFaceTets(true, 0, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Info<< " <<Writing " << nFaces
<< " faces with incorrectly orientated triangles to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
}
{
faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
if (mesh.checkFaceSkewness(true, &faces))
@ -470,6 +451,35 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
if (allGeometry)
{
faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
if
(
polyMeshTetDecomposition::checkFaceTets
(
mesh,
polyMeshTetDecomposition::minTetQuality,
true,
&faces
)
)
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Info<< " <<Writing " << nFaces
<< " faces with low quality or negative volume "
<< "decomposition tets to set " << faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
}
if (allGeometry)
{
// Note use of nPoints since don't want edge construction.

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

@ -7,32 +7,6 @@
#include "pointSet.H"
#include "IOmanip.H"
bool Foam::checkSync(const wordList& names)
{
List<wordList> allNames(Pstream::nProcs());
allNames[Pstream::myProcNo()] = names;
Pstream::gatherList(allNames);
Pstream::scatterList(allNames);
bool hasError = false;
for (label procI = 1; procI < allNames.size(); procI++)
{
if (allNames[procI] != allNames[0])
{
hasError = true;
Info<< " ***Inconsistent zones across processors, "
"processor 0 has zones:" << allNames[0]
<< ", processor " << procI << " has zones:"
<< allNames[procI]
<< endl;
}
}
return hasError;
}
Foam::label Foam::checkTopology
(
const polyMesh& mesh,
@ -51,35 +25,22 @@ Foam::label Foam::checkTopology
mesh.boundaryMesh().checkParallelSync(true);
// Check names of zones are equal
if (checkSync(mesh.cellZones().names()))
mesh.cellZones().checkDefinition(true);
if (mesh.cellZones().checkParallelSync(true))
{
noFailedChecks++;
}
if (checkSync(mesh.faceZones().names()))
mesh.faceZones().checkDefinition(true);
if (mesh.faceZones().checkParallelSync(true))
{
noFailedChecks++;
}
if (checkSync(mesh.pointZones().names()))
mesh.pointZones().checkDefinition(true);
if (mesh.pointZones().checkParallelSync(true))
{
noFailedChecks++;
}
// Check contents of faceZones consistent
{
forAll(mesh.faceZones(), zoneI)
{
if (mesh.faceZones()[zoneI].checkParallelSync(false))
{
Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
<< " is not correctly synchronised"
<< " across coupled boundaries."
<< " (coupled faces are either not both "
<< " present in set or have same flipmap)" << endl;
noFailedChecks++;
}
}
}
{
pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
if (mesh.checkPoints(true, &points))

2
applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
View File

@ -5,7 +5,5 @@ namespace Foam
{
class polyMesh;
bool checkSync(const wordList& names);
label checkTopology(const polyMesh&, const bool, const bool);
}

28
applications/utilities/parallelProcessing/decomposePar/decomposePar.C
View File

@ -558,7 +558,6 @@ int main(int argc, char *argv[])
lagrangianScalarFieldFields
);
lagrangianFieldDecomposer::readFields
(
cloudI,
@ -687,6 +686,19 @@ int main(int argc, char *argv[])
)
);
labelIOList faceProcAddressing
(
IOobject
(
"faceProcAddressing",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
labelIOList cellProcAddressing
(
IOobject
@ -728,19 +740,6 @@ int main(int argc, char *argv[])
|| surfaceTensorFields.size()
)
{
labelIOList faceProcAddressing
(
IOobject
(
"faceProcAddressing",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
fvFieldDecomposer fieldDecomposer
(
mesh,
@ -814,6 +813,7 @@ int main(int argc, char *argv[])
(
mesh,
procMesh,
faceProcAddressing,
cellProcAddressing,
cloudDirs[cloudI],
lagrangianPositions[cloudI],

39
applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
View File

@ -35,6 +35,7 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
(
const polyMesh& mesh,
const polyMesh& procMesh,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const word& cloudName,
const Cloud<indexedParticle>& lagrangianPositions,
@ -47,6 +48,14 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
{
label pi = 0;
// faceProcAddressing not required currently
// labelList decodedProcFaceAddressing(faceProcAddressing.size());
// forAll(faceProcAddressing, i)
// {
// decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
// }
forAll(cellProcAddressing, procCelli)
{
label celli = cellProcAddressing[procCelli];
@ -60,13 +69,41 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
const indexedParticle& ppi = *iter();
particleIndices_[pi++] = ppi.index();
// label mappedTetFace = findIndex
// (
// decodedProcFaceAddressing,
// ppi.tetFace()
// );
// if (mappedTetFace == -1)
// {
// FatalErrorIn
// (
// "Foam::lagrangianFieldDecomposer"
// "::lagrangianFieldDecomposer"
// "("
// "const polyMesh& mesh, "
// "const polyMesh& procMesh, "
// "const labelList& faceProcAddressing, "
// "const labelList& cellProcAddressing, "
// "const word& cloudName, "
// "const Cloud<indexedParticle>& "
// "lagrangianPositions, "
// "const List<SLList<indexedParticle*>*>& "
// "cellParticles"
// ")"
// ) << "Face lookup failure." << nl
// << abort(FatalError);
// }
positions_.append
(
new passiveParticle
(
positions_,
ppi.position(),
procCelli
procCelli,
false
)
);
}

1
applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
View File

@ -84,6 +84,7 @@ public:
(
const polyMesh& mesh,
const polyMesh& procMesh,
const labelList& faceProcAddressing,
const labelList& cellProcAddressing,
const word& cloudName,
const Cloud<indexedParticle>& lagrangianPositions,

4
applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C
View File

@ -51,7 +51,7 @@ void Foam::lagrangianFieldDecomposer::readFields
)
);
label lagrangianFieldi=0;
label lagrangianFieldi = 0;
forAllIter(IOobjectList, lagrangianTypeObjects, iter)
{
lagrangianFields[cloudI].set
@ -91,7 +91,7 @@ void Foam::lagrangianFieldDecomposer::readFieldFields
)
);
label lagrangianFieldi=0;
label lagrangianFieldi = 0;
forAllIter(IOobjectList, lagrangianTypeObjectsA, iter)
{

2
applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options
View File

@ -1,10 +1,12 @@
EXE_INC = \
/* -DFULLDEBUG -g -O0 */ \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lgenericPatchFields \
-llagrangian

2
applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options
View File

@ -1,10 +1,12 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/conversion/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-llagrangian \
-lmeshTools \
-lgenericPatchFields \
-lconversion

2
applications/utilities/postProcessing/dataConversion/foamToFieldview9/Make/options
View File

@ -1,8 +1,10 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lgenericPatchFields \
-llagrangian

290
applications/utilities/postProcessing/dataConversion/foamToFieldview9/floatScalar.H
View File

@ -1,290 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Typedef
Foam::floatScalar
Description
single floating point number identical to float
SourceFiles
floatScalar.C
\*---------------------------------------------------------------------------*/
#ifndef floatScalar_H
#define floatScalar_H
#include "label.H"
#include "word.H"
#include <cmath>
#ifdef ibm
float lgamma(float);
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Define floatScalar as a float
namespace Foam
{
typedef float floatScalar;
}
/*
// Define floatScalar as a float
namespace Foam
{
typedef float floatScalar;
// Largest and smallest floatScalar values allowed in certain
// parts of the code (6 is the number of significant figures in an
// IEEE single precision number. See limits.h or float.h)
static const floatScalar GREAT = 1.0e+6;
static const floatScalar VGREAT = 1.0e+37;
static const floatScalar SMALL = 1.0e-6;
static const floatScalar VSMALL = 1.0e-37;
}
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "pTraits.H"
#include "products.H"
#include "direction.H"
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// template specialisation for pTraits<floatScalar>
template<>
class pTraits<floatScalar>
{
floatScalar p_;
public:
//- Component type
typedef floatScalar cmptType;
// Member constants
enum
{
dim = 3, // Dimensionality of space
rank = 0, // Rank od floatScalar is 0
nComponents = 1 // Number of components in floatScalar is 1
};
// Static data members
static const char* const typeName;
static const char* componentNames[];
static const floatScalar zero;
static const floatScalar one;
// Constructors
//- Construct from Istream
pTraits(Istream& is);
// Member Functions
operator floatScalar() const
{
return p_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//template<class Cmpt>
//class typeOfRank<Cmpt, 0>
//{
//public:
//
// typedef Cmpt type;
//};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Return a string representation of a floatScalar
word name(const floatScalar s)
{
return name(scalar(s));
}
#define MAXMINPOW(retType, type1, type2) \
\
MAXMIN(retType, type1, type2) \
\
inline float pow(const type1 s, const type2 e) \
{ \
return ::pow(float(s), float(e)); \
}
//MAXMINPOW(float, float, float)
//MAXMINPOW(float, float, int)
//MAXMINPOW(float, int, float)
//MAXMINPOW(float, float, long)
//MAXMINPOW(float, long, float)
//MAXMINPOW(float, float, int)
//MAXMINPOW(float, int, float)
//MAXMINPOW(float, float, long)
//MAXMINPOW(float, long, float)
#undef MAXMINPOW
inline floatScalar mag(const floatScalar s)
{
return ::fabs(s);
}
inline floatScalar sign(const floatScalar s)
{
return (s >= 0)? 1: -1;
}
inline floatScalar pos(const floatScalar s)
{
return (s >= 0)? 1: 0;
}
inline floatScalar neg(const floatScalar s)
{
return (s < 0)? 1: 0;
}
inline floatScalar limit(const floatScalar s1, const floatScalar s2)
{
return (mag(s1) < mag(s2))? s1: 0.0;
}
inline floatScalar magSqr(const floatScalar s)
{
return s*s;
}
inline floatScalar sqr(const floatScalar s)
{
return s*s;
}
inline floatScalar pow3(const floatScalar s)
{
return s*s*s;
}
inline floatScalar pow4(const floatScalar s)
{
return sqr(sqr(s));
}
inline floatScalar cmptAv(const floatScalar s)
{
return s;
}
inline floatScalar cmptMag(const floatScalar s)
{
return mag(s);
}
inline floatScalar scale(const floatScalar s, const floatScalar d)
{
return s*d;
}
#define transFunc(func) \
inline floatScalar func(const floatScalar s) \
{ \
return ::func(s); \
}
// Standard C++ transcendental functions
transFunc(sqrt)
transFunc(exp)
transFunc(log)
transFunc(log10)
transFunc(sin)
transFunc(cos)
transFunc(tan)
transFunc(asin)
transFunc(acos)
transFunc(atan)
transFunc(sinh)
transFunc(cosh)
transFunc(tanh)
transFunc(asinh)
transFunc(acosh)
transFunc(atanh)
// Standard ANSI-C (but not in <cmath>) transcendental functions
transFunc(erf)
transFunc(erfc)
transFunc(lgamma)
transFunc(j0)
transFunc(j1)
transFunc(y0)
transFunc(y1)
#undef transFunc
// Stabilisation around zero for division
inline floatScalar stabilise(const floatScalar s, const floatScalar small)
{
if (s >= 0)
{
return s + small;
}
else
{
return s - small;
}
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
//floatScalar readScalar(Istream& is);
Istream& operator>>(Istream&, floatScalar&);
Ostream& operator<<(Ostream&, const floatScalar);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
applications/utilities/postProcessing/dataConversion/foamToGMV/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/lagrangian/lnInclude \
@ -14,5 +15,6 @@ EXE_INC = \
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lgenericPatchFields \
-llagrangian

2
applications/utilities/postProcessing/graphics/ensightFoamReader/Make/options
View File

@ -1,5 +1,6 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/browser/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
@ -7,6 +8,7 @@ EXE_INC = \
LIB_LIBS = \
-lOpenFOAM \
-lfiniteVolume \
-lmeshTools \
-lgenericPatchFields \
-llagrangian \
$(PROJECT_LIBS)

2
applications/utilities/postProcessing/lagrangian/particleTracks/Make/options
View File

@ -1,8 +1,10 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lgenericPatchFields \
-llagrangian

1
applications/utilities/preProcessing/dsmcInitialise/Make/options
View File

@ -9,4 +9,3 @@ EXE_LIBS = \
-lfiniteVolume \
-llagrangian \
-ldsmc

136
applications/utilities/preProcessing/mapFields/MapLagrangianFields.H
View File

@ -37,6 +37,7 @@ Description
#include "GeometricField.H"
#include "meshToMesh.H"
#include "IOobjectList.H"
#include "IOFieldField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,37 +57,118 @@ void MapLagrangianFields
{
const fvMesh& meshTarget = meshToMeshInterp.toMesh();
IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
Info<< " mapping lagrangian field " << fieldIter()->name() << endl;
IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
// Read field (does not need mesh)
IOField<Type> fieldSource(*fieldIter());
// Map
IOField<Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
forAllIter(IOobjectList, fields, fieldIter)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
Info<< " mapping lagrangian field "
<< fieldIter()->name() << endl;
// Write field
fieldTarget.write();
// Read field (does not need mesh)
IOField<Type> fieldSource(*fieldIter());
// Map
IOField<Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(IOField<Field<Type> >::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
Info<< " mapping lagrangian fieldField "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
IOField<Field<Type> > fieldSource(*fieldIter());
// Map - use IOFieldField to automatically write in
// compact form for binary format.
IOFieldField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(IOFieldField<Field<Type>, Type>::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
Info<< " mapping lagrangian fieldField "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
IOFieldField<Field<Type>, Type> fieldSource(*fieldIter());
// Map
IOFieldField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
}

23
applications/utilities/preProcessing/mapFields/mapLagrangian.C
View File

@ -38,12 +38,13 @@ static const scalar perturbFactor = 1E-6;
// Special version of findCell that generates a cell guaranteed to be
// compatible with tracking.
static label findCell(const meshSearch& meshSearcher, const point& pt)
static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
{
const polyMesh& mesh = meshSearcher.mesh();
label cellI = -1;
label tetFaceI = -1;
label tetPtI = -1;
// Use tracking to find cell containing pt
label cellI = meshSearcher.findCell(pt);
cloud.findCellFacePt(pt, cellI, tetFaceI, tetPtI);
if (cellI >= 0)
{
@ -54,16 +55,24 @@ static label findCell(const meshSearch& meshSearcher, const point& pt)
// See if particle on face by finding nearest face and shifting
// particle.
const polyMesh& mesh = cloud.pMesh();
meshSearch meshSearcher(mesh, false);
label faceI = meshSearcher.findNearestBoundaryFace(pt);
if (faceI >= 0)
{
const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
return meshSearcher.findCell(perturbPt);
cloud.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI);
return cellI;
}
}
return -1;
}
@ -207,8 +216,6 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp)
if (unmappedSource.size())
{
meshSearch targetSearcher(meshTarget, false);
sourceParticleI = 0;
forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
@ -216,7 +223,7 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp)
if (unmappedSource.found(sourceParticleI))
{
label targetCell =
findCell(targetSearcher, iter().position());
findCell(targetParcels, iter().position());
if (targetCell >= 0)
{

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

@ -35,7 +35,6 @@ Description
#include "IOdictionary.H"
#include "boundBox.H"
#include "indexedOctree.H"
#include "octree.H"
#include "treeDataTriSurface.H"
#include "Random.H"

1
etc/controlDict
View File

@ -780,6 +780,7 @@ DebugSwitches
syringePressure 0;
tensorAverageField 0;
tensorField 0;
tetDecomposedPolyMesh 0;
thermoCloud 0;
thermophysicalFunction 0;
time 0;

3
src/Allwmake
View File

@ -20,8 +20,6 @@ Pstream/Allwmake
OSspecific/$WM_OSTYPE/Allwmake
wmake libso OpenFOAM
wmake libso lagrangian/basic
wmake libso fileFormats
wmake libso edgeMesh
wmake libso surfMesh
@ -33,6 +31,7 @@ parallel/decompose/AllwmakeLnInclude
dummyThirdParty/Allwmake
wmake libso meshTools
wmake libso lagrangian/basic
wmake libso finiteVolume
wmake libso genericPatchFields

16
src/ODE/ODESolvers/ODESolver/ODESolver.C
View File

@ -60,6 +60,8 @@ void Foam::ODESolver::solve
scalar x = xStart;
scalar h = hEst;
scalar hNext = 0;
scalar hPrev = 0;
for (label nStep=0; nStep<MAXSTP; nStep++)
{
@ -73,14 +75,24 @@ void Foam::ODESolver::solve
if ((x + h - xEnd)*(x + h - xStart) > 0.0)
{
h = xEnd - x;
hPrev = hNext;
}
scalar hNext, hDid;
hNext = 0;
scalar hDid;
solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext);
if ((x - xEnd)*(xEnd - xStart) >= 0.0)
{
hEst = hNext;
if (hPrev != 0)
{
hEst = hPrev;
}
else
{
hEst = hNext;
}
return;
}

2
src/OpenFOAM/Make/files
View File

@ -365,6 +365,8 @@ $(globalMeshData)/globalPoints.C
$(globalMeshData)/globalIndex.C
$(polyMesh)/syncTools/syncTools.C
$(polyMesh)/polyMeshTetDecomposition/polyMeshTetDecomposition.C
$(polyMesh)/polyMeshTetDecomposition/tetIndices.C
zone = $(polyMesh)/zones/zone
$(zone)/zone.C

26
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -33,6 +33,7 @@ License
#include "processorPolyPatch.H"
#include "OSspecific.H"
#include "demandDrivenData.H"
#include "polyMeshTetDecomposition.H"
#include "pointMesh.H"
@ -204,6 +205,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
bounds_(points_),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
pointZones_
(
IOobject
@ -392,6 +394,7 @@ Foam::polyMesh::polyMesh
bounds_(points_, syncPar),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
pointZones_
(
IOobject
@ -548,6 +551,7 @@ Foam::polyMesh::polyMesh
bounds_(points_, syncPar),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
pointZones_
(
IOobject
@ -847,6 +851,28 @@ Foam::label Foam::polyMesh::nSolutionD() const
}
const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
{
if (!tetBasePtIsPtr_)
{
if (debug)
{
WarningIn("const labelList& polyMesh::tetBasePtIs() const")
<< "Tet base point indices not available. "
<< "Forcing storage of base points."
<< endl;
}
tetBasePtIsPtr_ = new labelList
(
polyMeshTetDecomposition::findFaceBasePts(*this)
);
}
return *tetBasePtIsPtr_;
}
// Add boundary patches. Constructor helper
void Foam::polyMesh::addPatches
(

7
src/OpenFOAM/meshes/polyMesh/polyMesh.H
View File

@ -62,6 +62,7 @@ namespace Foam
class globalMeshData;
class mapPolyMesh;
class polyMeshTetDecomposition;
/*---------------------------------------------------------------------------*\
Class polyMesh Declaration
@ -127,6 +128,9 @@ private:
// defined according to the presence of empty patches
mutable Vector<label> solutionD_;
//- Base point for face decomposition into tets
mutable labelList* tetBasePtIsPtr_;
// Zoning information
@ -356,6 +360,9 @@ public:
//- Return the number of valid solved-for dimensions in the mesh
label nSolutionD() const;
//- Return the tetBasePtIs
const labelList& tetBasePtIs() const;
//- Return point zone mesh
const pointZoneMesh& pointZones() const
{

3
src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
View File

@ -71,6 +71,9 @@ void Foam::polyMesh::clearGeom()
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Remove the stored tet base points
deleteDemandDrivenData(tetBasePtIsPtr_);
pointMesh::Delete(*this);
}

2
src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
View File

@ -503,6 +503,7 @@ Foam::polyMesh::polyMesh
bounds_(points_, syncPar),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
pointZones_
(
IOobject
@ -775,6 +776,7 @@ Foam::polyMesh::polyMesh
bounds_(points_, syncPar),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
pointZones_
(
IOobject

627
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C Normal file
View File

@ -0,0 +1,627 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyMeshTetDecomposition.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = 1e-9;
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
(
const polyMesh& mesh,
label fI,
const point& nCc,
scalar tol,
bool report
)
{
const faceList& pFaces = mesh.faces();
const pointField& pPts = mesh.points();
const vectorField& pC = mesh.cellCentres();
const labelList& pOwner = mesh.faceOwner();
const face& f = pFaces[fI];
label oCI = pOwner[fI];
const point& oCc = pC[oCI];
forAll(f, faceBasePtI)
{
scalar thisBaseMinTetQuality = VGREAT;
const point& tetBasePt = pPts[f[faceBasePtI]];
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
List<scalar> tetQualities(2, 0.0);
{
// owner cell tet
label ptAI = f[facePtI];
label ptBI = f[otherFacePtI];
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(oCc, tetBasePt, pA, pB);
tetQualities[0] = tet.quality();
}
{
// neighbour cell tet
label ptAI = f[otherFacePtI];
label ptBI = f[facePtI];
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(nCc, tetBasePt, pA, pB);
tetQualities[1] = tet.quality();
}
if (min(tetQualities) < thisBaseMinTetQuality)
{
thisBaseMinTetQuality = min(tetQualities);
}
}
if (thisBaseMinTetQuality > tol)
{
return faceBasePtI;
}
}
// If a base point hasn't triggered a return by now, then there is
// non that can produce a good decomposition
return -1;
}
Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
(
const polyMesh& mesh,
label fI,
scalar tol,
bool report
)
{
return findSharedBasePoint
(
mesh,
fI,
mesh.cellCentres()[mesh.faceNeighbour()[fI]],
tol,
report
);
}
Foam::label Foam::polyMeshTetDecomposition::findBasePoint
(
const polyMesh& mesh,
label fI,
scalar tol,
bool report
)
{
const faceList& pFaces = mesh.faces();
const pointField& pPts = mesh.points();
const vectorField& pC = mesh.cellCentres();
const labelList& pOwner = mesh.faceOwner();
const face& f = pFaces[fI];
label cI = pOwner[fI];
bool own = (pOwner[fI] == cI);
const point& cC = pC[cI];
forAll(f, faceBasePtI)
{
scalar thisBaseMinTetQuality = VGREAT;
const point& tetBasePt = pPts[f[faceBasePtI]];
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label ptAI = -1;
label ptBI = -1;
if (own)
{
ptAI = f[facePtI];
ptBI = f[otherFacePtI];
}
else
{
ptAI = f[otherFacePtI];
ptBI = f[facePtI];
}
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(cC, tetBasePt, pA, pB);
scalar tetQuality = tet.quality();
if (tetQuality < thisBaseMinTetQuality)
{
thisBaseMinTetQuality = tetQuality;
}
}
if (thisBaseMinTetQuality > tol)
{
return faceBasePtI;
}
}
// If a base point hasn't triggered a return by now, then there is
// non that can produce a good decomposition
return -1;
}
Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
(
const polyMesh& mesh,
scalar tol,
bool report
)
{
const labelList& pOwner = mesh.faceOwner();
// Find a suitable base point for each face, considering both
// cells for interface faces or those on coupled patches
labelList tetBasePtIs(mesh.nFaces(), -1);
label nInternalFaces = mesh.nInternalFaces();
for (label fI = 0; fI < nInternalFaces; fI++)
{
tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
}
pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
for(label faceI = nInternalFaces; faceI < mesh.nFaces(); faceI++)
{
neighbourCellCentres[faceI - nInternalFaces] =
mesh.cellCentres()[pOwner[faceI]];
}
syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
SubList<label> boundaryFaceTetBasePtIs
(
tetBasePtIs,
mesh.nFaces() - nInternalFaces,
nInternalFaces
);
for
(
label fI = nInternalFaces, bFI = 0;
fI < mesh.nFaces();
fI++, bFI++
)
{
label patchI =
mesh.boundaryMesh().patchID()[bFI];
if (patches[patchI].coupled())
{
const coupledPolyPatch& pp =
refCast<const coupledPolyPatch>(patches[patchI]);
if (pp.owner())
{
boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
(
mesh,
fI,
neighbourCellCentres[bFI],
tol,
report
);
}
else
{
// Assign -2, to distinguish from a failed base point
// find, which returns -1.
boundaryFaceTetBasePtIs[bFI] = -2;
}
}
else
{
boundaryFaceTetBasePtIs[bFI] = findBasePoint
(
mesh,
fI,
tol,
report
);
}
}
// maxEqOp will replace the -2 values on the neighbour patches
// with the result from the owner base point find.
syncTools::syncBoundaryFaceList
(
mesh,
boundaryFaceTetBasePtIs,
maxEqOp<label>()
);
for
(
label fI = nInternalFaces, bFI = 0;
fI < mesh.nFaces();
fI++, bFI++
)
{
label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
if (bFTetBasePtI == -2)
{
FatalErrorIn
(
"Foam::labelList"
"Foam::polyMeshTetDecomposition::findFaceBasePts"
"("
"const polyMesh& mesh, "
"scalar tol, "
"bool report"
")"
)
<< "Coupled face base point exchange failure for face "
<< fI
<< abort(FatalError);
}
if (bFTetBasePtI < 1)
{
// If the base point is -1, it should be left as such to
// indicate a problem, if it is 0, then no action is required.
continue;
}
label patchI = mesh.boundaryMesh().patchID()[bFI];
if (patches[patchI].coupled())
{
const coupledPolyPatch& pp =
refCast<const coupledPolyPatch>(patches[patchI]);
// Calculated base points on coupled faces are those of
// the owner patch face. They need to be reindexed to for
// the non-owner face, which has the opposite order.
// So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
// i.e.:
// owner coupledPolyPatch face
// face (a b c d e f)
// fPtI 0 1 2 3 4 5
// +
// #
// neighbour coupledPolyPatch face
// face (a f e d c b)
// fPtI 0 1 2 3 4 5
// +
// #
// +: 6 - 1 = 5
// #: 6 - 2 = 4
if (!pp.owner())
{
bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
}
}
}
return tetBasePtIs;
}
bool Foam::polyMeshTetDecomposition::checkFaceTets
(
const polyMesh& mesh,
scalar tol,
const bool report,
labelHashSet* setPtr
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const vectorField& cc = mesh.cellCentres();
const vectorField& fc = mesh.faceCentres();
// Calculate coupled cell centre
pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
neiCc[faceI - mesh.nInternalFaces()] = cc[own[faceI]];
}
syncTools::swapBoundaryFacePositions(mesh, neiCc);
const faceList& fcs = mesh.faces();
const pointField& p = mesh.points();
label nErrorTets = 0;
forAll(fcs, faceI)
{
const face& f = fcs[faceI];
forAll(f, fPtI)
{
scalar tetQual = tetPointRef
(
p[f[fPtI]],
p[f.nextLabel(fPtI)],
fc[faceI],
cc[own[faceI]]
).quality();
if (tetQual > -tol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
break; // no need to check other tets
}
}
if (mesh.isInternalFace(faceI))
{
// Create the neighbour tet - it will have positive volume
const face& f = fcs[faceI];
forAll(f, fPtI)
{
scalar tetQual = tetPointRef
(
p[f[fPtI]],
p[f.nextLabel(fPtI)],
fc[faceI],
cc[nei[faceI]]
).quality();
if (tetQual < tol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
break;
}
}
if (findSharedBasePoint(mesh, faceI, tol, report) == -1)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
else
{
label patchI = patches.patchID()[faceI - mesh.nInternalFaces()];
if (patches[patchI].coupled())
{
if
(
findSharedBasePoint
(
mesh,
faceI,
neiCc[faceI - mesh.nInternalFaces()],
tol,
report
) == -1
)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
else
{
if (findBasePoint(mesh, faceI, tol, report) == -1)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
}
}
reduce(nErrorTets, sumOp<label>());
if (nErrorTets > 0)
{
if (report)
{
Info<< " ***Error in face tets: "
<< nErrorTets << " faces with low quality or negative volume "
<< "decomposition tets." << endl;
}
return true;
}
else
{
if (report)
{
Info<< " Face tets OK." << endl;
}
return false;
}
}
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
(
const polyMesh& mesh,
label fI,
label cI
)
{
const faceList& pFaces = mesh.faces();
const labelList& pOwner = mesh.faceOwner();
const face& f = pFaces[fI];
label nTets = f.size() - 2;
List<tetIndices> faceTets(nTets);
bool own = (pOwner[fI] == cI);
label tetBasePtI = mesh.tetBasePtIs()[fI];
if (tetBasePtI == -1)
{
FatalErrorIn
(
"Foam::List<Foam::FixedList<Foam::label, 4> >"
"Foam::Cloud<ParticleType>::"
"faceTetIndices(label fI, label cI) const"
)
<< "No base point for face " << fI << ", " << f
<< ", produces a valid tet decomposition."
<< abort(FatalError);
}
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
{
tetIndices& faceTetIs = faceTets[tetPtI - 1];
label facePtI = (tetPtI + tetBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
faceTetIs.cell() = cI;
faceTetIs.face() = fI;
faceTetIs.faceBasePt() = tetBasePtI;
if (own)
{
faceTetIs.facePtA() = facePtI;
faceTetIs.facePtB() = otherFacePtI;
}
else
{
faceTetIs.facePtA() = otherFacePtI;
faceTetIs.facePtB() = facePtI;
}
faceTetIs.tetPt() = tetPtI;
}
return faceTets;
}
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
(
const polyMesh& mesh,
label cI
)
{
const faceList& pFaces = mesh.faces();
const cellList& pCells = mesh.cells();
const cell& thisCell = pCells[cI];
label nTets = 0;
forAll(thisCell, cFI)
{
nTets += pFaces[thisCell[cFI]].size() - 2;
}
DynamicList<tetIndices> cellTets(nTets);
forAll(thisCell, cFI)
{
label fI = thisCell[cFI];
cellTets.append(faceTetIndices(mesh, fI, cI));
}
return cellTets;
}
// ************************************************************************* //

151
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H Normal file
View File

@ -0,0 +1,151 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::polyMeshTetDecomposition
Description
Tools for performing the minimum decomposition of faces of the
mesh into triangles so that the cells may be tet decomposed.
Includes functions for finding variable face starting (base)
points on each face to avoid the decomposition of cells into tets
that have negative or zero volume.
SourceFiles
polyMeshTetDecomposition.C
\*---------------------------------------------------------------------------*/
#ifndef polyMeshTetDecomposition_H
#define polyMeshTetDecomposition_H
#include "polyMesh.H"
#include "coupledPolyPatch.H"
#include "syncTools.H"
#include "tetPointRef.H"
#include "tetIndices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyMeshTetDecomposition Declaration
\*---------------------------------------------------------------------------*/
class polyMeshTetDecomposition
{
public:
// Static data members
//- Minimum tetrahedron quality
static const scalar minTetQuality;
// Member Functions
//- Find the first suitable base point to use for a minimum
// triangle decomposition of the face, suiting owner and
// neighbour cells. Finds the first base point on the face
// whose worst quality tet from either cell is better than
// tolerance. Neighbour cell centre supplied. For coupled
// patches.
static label findSharedBasePoint
(
const polyMesh& mesh,
label fI,
const point& nCc,
scalar tol,
bool report = false
);
//- As for findSharedBasePoint, but using neighbour cell
// centre from the mesh. For internal faces.
static label findSharedBasePoint
(
const polyMesh& mesh,
label fI,
scalar tol,
bool report = false
);
//- Find the base point to use for a minimum triangle
// decomposition of the face, using only the owner
// information. For non-coupled boundary faces.
static label findBasePoint
(
const polyMesh& mesh,
label fI,
scalar tol,
bool report = false
);
//- Find a suitable base point for each face for decomposition
// into tets
static labelList findFaceBasePts
(
const polyMesh& mesh,
scalar tol = minTetQuality,
bool report = false
);
//- Check face-decomposition tet volume
static bool checkFaceTets
(
const polyMesh& mesh,
scalar tol = minTetQuality,
const bool report = false,
labelHashSet* setPtr = NULL
);
//- Return the tet decomposition of the given face, with
// respect to the given cell
static List<tetIndices> faceTetIndices
(
const polyMesh& mesh,
label fI,
label cI
);
//- Return the tet decomposition of the given cell, see
// findFacePt for the meaning of the indices
static List<tetIndices> cellTetIndices
(
const polyMesh& mesh,
label cI
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

148
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C Normal file
View File

@ -0,0 +1,148 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "tetIndices.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::tetIndices::tetIndices()
:
cellI_(-1),
faceI_(-1),
faceBasePtI_(-1),
facePtAI_(-1),
facePtBI_(-1),
tetPtI_(-1)
{}
Foam::tetIndices::tetIndices
(
label cellI,
label faceI,
label faceBasePtI,
label facePtAI,
label facePtBI,
label tetPtI
)
:
cellI_(cellI),
faceI_(faceI),
faceBasePtI_(faceBasePtI),
facePtAI_(facePtAI),
facePtBI_(facePtBI),
tetPtI_(tetPtI)
{}
Foam::tetIndices::tetIndices
(
label cellI,
label faceI,
label tetPtI,
const polyMesh& mesh
)
:
cellI_(cellI),
faceI_(faceI),
faceBasePtI_(-1),
facePtAI_(-1),
facePtBI_(-1),
tetPtI_(tetPtI)
{
const faceList& pFaces = mesh.faces();
const labelList& pOwner = mesh.faceOwner();
const Foam::face& f = pFaces[faceI_];
bool own = (pOwner[faceI_] == cellI_);
faceBasePtI_ = mesh.tetBasePtIs()[faceI_];
label facePtI = (tetPtI_ + faceBasePtI_) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
if (own)
{
facePtAI_ = facePtI;
facePtBI_ = otherFacePtI;
}
else
{
facePtAI_ = otherFacePtI;
facePtBI_ = facePtI;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::tetIndices::~tetIndices()
{}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>(Istream& is, tetIndices& tI)
{
is >> tI.cell()
>> tI.face()
>> tI.faceBasePt()
>> tI.facePtA()
>> tI.facePtB()
>> tI.tetPt();
// Check state of Istream
is.check
(
"Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::tetIndices&)"
);
return is;
}
Foam::Ostream& Foam::operator<<(Ostream& os, const tetIndices& tI)
{
os << tI.cell() << token::SPACE
<< tI.face() << token::SPACE
<< tI.faceBasePt() << token::SPACE
<< tI.facePtA() << token::SPACE
<< tI.facePtB() << token::SPACE
<< tI.tetPt() << token::SPACE
<< endl;
// Check state of Ostream
os.check
(
"Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
"const Foam::tetIndices&)"
);
return os;
}
// ************************************************************************* //

201
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H Normal file
View File

@ -0,0 +1,201 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::tetIndices
Description
Storage and named access for the indices of a tet which is part of
the decomposition of a cell.
SourceFiles
tetIndicesI.H
tetIndices.C
\*---------------------------------------------------------------------------*/
#ifndef tetIndices_H
#define tetIndices_H
#include "label.H"
#include "tetPointRef.H"
#include "triPointRef.H"
#include "polyMesh.H"
#include "face.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class tetIndices Declaration
\*---------------------------------------------------------------------------*/
class tetIndices
{
// Private data
//- cell that this is a decomposed tet of
label cellI_;
//- face that holds this decomposed tet
label faceI_;
//- base point on the face
label faceBasePtI_;
//- point on the face such that the right-hand circulation
// {faceBasePtI_, facePtIA_, facePtBI_}
// forms a triangle that points out of the tet
label facePtAI_;
//- point on the face such that the right-hand circulation
// {faceBasePtI_, facePtIA_, facePtBI_}
// forms a triangle that points out of the tet
label facePtBI_;
//- point on the face, *relative to the base point*, which
// characterises this tet on the face.
label tetPtI_;
public:
// Constructors
//- Construct null
tetIndices();
//- Construct from components
tetIndices
(
label cellI,
label faceI,
label faceBasePtI,
label facePtAI,
label facePtBI,
label tetPtI
);
//- Construct from cell, face, pt description of tet
tetIndices
(
label cellI,
label faceI,
label tetPtI,
const polyMesh& mesh
);
//- Destructor
~tetIndices();
// Member Functions
// Access
//- Return the cell
inline label cell() const;
//- Return the face
inline label face() const;
//- Return the face base point
inline label faceBasePt() const;
//- Return face point A
inline label facePtA() const;
//- Return face point B
inline label facePtB() const;
//- Return the characterising tetPtI
inline label tetPt() const;
//- Return the geometry corresponding to this tet from the
// supplied mesh
inline tetPointRef tet(const polyMesh& mesh) const;
//- Return the geometry corresponding to this tet from the
// supplied mesh using the old positions
inline tetPointRef oldTet(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the
// mesh face for this tet from the supplied mesh
inline triPointRef faceTri(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the
// mesh face for this tet from the supplied mesh using
// the old position
inline triPointRef oldFaceTri(const polyMesh& mesh) const;
// Edit
//- Return non-const access to the cell
inline label& cell();
//- Return non-const access to the face
inline label& face();
//- Return non-const access to the face base point
inline label& faceBasePt();
//- Return non-const access to face point A
inline label& facePtA();
//- Return non-const access to face point B
inline label& facePtB();
//- Return non-const access to the characterising tetPtI
inline label& tetPt();
// Member Operators
inline bool operator==(const tetIndices&) const;
inline bool operator!=(const tetIndices&) const;
// IOstream Operators
friend Istream& operator>>(Istream&, tetIndices&);
friend Ostream& operator<<(Ostream&, const tetIndices&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "tetIndicesI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

202
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H Normal file
View File

@ -0,0 +1,202 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::tetIndices::cell() const
{
return cellI_;
}
Foam::label Foam::tetIndices::face() const
{
return faceI_;
}
Foam::label Foam::tetIndices::faceBasePt() const
{
return faceBasePtI_;
}
Foam::label Foam::tetIndices::facePtA() const
{
return facePtAI_;
}
Foam::label Foam::tetIndices::facePtB() const
{
return facePtBI_;
}
Foam::label Foam::tetIndices::tetPt() const
{
return tetPtI_;
}
Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const
{
const pointField& pPts = mesh.points();
const faceList& pFaces = mesh.faces();
const vectorField& pC = mesh.cellCentres();
const Foam::face& f = pFaces[faceI_];
return tetPointRef
(
pC[cellI_],
pPts[f[faceBasePtI_]],
pPts[f[facePtAI_]],
pPts[f[facePtBI_]]
);
}
Foam::tetPointRef Foam::tetIndices::oldTet(const polyMesh& mesh) const
{
const pointField& oldPPts = mesh.oldPoints();
const faceList& pFaces = mesh.faces();
// We need to reconstruct the old Cc from oldPoints (it isn't
// stored)
point oldC = mesh.cells()[cellI_].centre
(
oldPPts,
pFaces
);
const Foam::face& f = pFaces[faceI_];
return tetPointRef
(
oldC,
oldPPts[f[faceBasePtI_]],
oldPPts[f[facePtAI_]],
oldPPts[f[facePtBI_]]
);
}
Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
{
const pointField& pPts = mesh.points();
const faceList& pFaces = mesh.faces();
const Foam::face& f = pFaces[faceI_];
return triPointRef
(
pPts[f[faceBasePtI_]],
pPts[f[facePtAI_]],
pPts[f[facePtBI_]]
);
}
Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const
{
const pointField& oldPPts = mesh.oldPoints();
const faceList& pFaces = mesh.faces();
const Foam::face& f = pFaces[faceI_];
return triPointRef
(
oldPPts[f[faceBasePtI_]],
oldPPts[f[facePtAI_]],
oldPPts[f[facePtBI_]]
);
}
Foam::label& Foam::tetIndices::cell()
{
return cellI_;
}
Foam::label& Foam::tetIndices::face()
{
return faceI_;
}
Foam::label& Foam::tetIndices::faceBasePt()
{
return faceBasePtI_;
}
Foam::label& Foam::tetIndices::facePtA()
{
return facePtAI_;
}
Foam::label& Foam::tetIndices::facePtB()
{
return facePtBI_;
}
Foam::label& Foam::tetIndices::tetPt()
{
return tetPtI_;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline bool Foam::tetIndices::operator==(const Foam::tetIndices& rhs) const
{
return
(
cell() == rhs.cell()
&& face() == rhs.face()
&& faceBasePt() == rhs.faceBasePt()
&& facePtA() == rhs.facePtA()
&& facePtB() == rhs.facePtB()
&& tetPt() == rhs.tetPt()
);
}
inline bool Foam::tetIndices::operator!=(const Foam::tetIndices& rhs) const
{
return !(*this == rhs);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

13
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
View File

@ -229,6 +229,19 @@ public:
// Member functions
//- Return true only if this is a parallel run
virtual bool coupled() const
{
if (Pstream::parRun())
{
return true;
}
else
{
return false;
}
}
//- Return processor number
int myProcNo() const
{

6
src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
View File

@ -1223,7 +1223,7 @@ void Foam::syncTools::syncBoundaryFaceList
label patchStart = procPatch.start()-mesh.nInternalFaces();
UOPstream toNbr(procPatch.neighbProcNo(), pBufs);
toNbr <<
toNbr <<
SubField<T>
(
faceValues,
@ -1423,7 +1423,7 @@ void Foam::syncTools::syncFaceList
cop(t, val1);
faceValues[meshFace0] = t;
cop(val1, val0);
cop(val1, val0);
faceValues[meshFace1] = val1;
}
}
@ -1683,7 +1683,7 @@ void Foam::syncTools::syncEdgeList
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(patches[patchI]);
// Receive from neighbour.
// Receive from neighbour.
List<unsigned int> nbrPatchInfo(procPatch.nEdges());
{

79
src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
View File

@ -27,6 +27,7 @@ License
#include "entry.H"
#include "demandDrivenData.H"
#include "stringListOps.H"
#include "Pstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -372,6 +373,84 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkDefinition
}
template<class ZoneType, class MeshType>
bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
(
const bool report
) const
{
if (!Pstream::parRun())
{
return false;
}
const PtrList<ZoneType>& zones = *this;
bool hasError = false;
// Collect all names
List<wordList> allNames(Pstream::nProcs());
allNames[Pstream::myProcNo()] = this->names();
Pstream::gatherList(allNames);
Pstream::scatterList(allNames);
List<wordList> allTypes(Pstream::nProcs());
allTypes[Pstream::myProcNo()] = this->types();
Pstream::gatherList(allTypes);
Pstream::scatterList(allTypes);
// Have every processor check but only master print error.
for (label procI = 1; procI < allNames.size(); procI++)
{
if
(
(allNames[procI] != allNames[0])
|| (allTypes[procI] != allTypes[0])
)
{
hasError = true;
if (debug || (report && Pstream::master()))
{
Info<< " ***Inconsistent zones across processors, "
"processor 0 has zone names:" << allNames[0]
<< " zone types:" << allTypes[0]
<< " processor " << procI << " has zone names:"
<< allNames[procI]
<< " zone types:" << allTypes[procI]
<< endl;
}
}
}
// Check contents
if (!hasError)
{
forAll(zones, zoneI)
{
if (zones[zoneI].checkParallelSync(false))
{
hasError = true;
if (debug || (report && Pstream::master()))
{
Info<< " ***Zone " << zones[zoneI].name()
<< " of type " << zones[zoneI].type()
<< " is not correctly synchronised"
<< " across coupled boundaries."
<< " (coupled faces are either not both "
<< " present in set or have same flipmap)" << endl;
}
}
}
}
return hasError;
}
// Correct zone mesh after moving points
template<class ZoneType, class MeshType>
void Foam::ZoneMesh<ZoneType, MeshType>::movePoints(const pointField& p)

4
src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
View File

@ -149,6 +149,10 @@ public:
//- Check zone definition. Return true if in error.
bool checkDefinition(const bool report = false) const;
//- Check whether all procs have all zones and in same order. Return
// true if in error.
bool checkParallelSync(const bool report = false) const;
//- Correct zone mesh after moving points
void movePoints(const pointField&);

7
src/OpenFOAM/meshes/polyMesh/zones/cellZone/cellZone.H
View File

@ -209,6 +209,13 @@ public:
//- Check zone definition. Return true if in error.
virtual bool checkDefinition(const bool report = false) const;
//- Check whether zone is synchronised across coupled boundaries. Return
// true if in error.
virtual bool checkParallelSync(const bool report = false) const
{
return false;
}
//- Write dictionary
virtual void writeDict(Ostream&) const;

44
src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.C
View File

@ -29,6 +29,7 @@ License
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "demandDrivenData.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -134,6 +135,49 @@ bool Foam::pointZone::checkDefinition(const bool report) const
}
bool Foam::pointZone::checkParallelSync(const bool report) const
{
const polyMesh& mesh = zoneMesh().mesh();
labelList maxZone(mesh.nPoints(), -1);
labelList minZone(mesh.nPoints(), labelMax);
forAll(*this, i)
{
label pointI = operator[](i);
maxZone[pointI] = index();
minZone[pointI] = index();
}
syncTools::syncPointList(mesh, maxZone, maxEqOp<label>(), -1);
syncTools::syncPointList(mesh, minZone, minEqOp<label>(), labelMax);
bool error = false;
forAll(maxZone, pointI)
{
// Check point in zone on both sides
if (maxZone[pointI] != minZone[pointI])
{
if (report && !error)
{
Info<< " ***Problem with pointZone " << index()
<< " named " << name()
<< ". Point " << pointI
<< " at " << mesh.points()[pointI]
<< " is in zone "
<< (minZone[pointI] == labelMax ? -1 : minZone[pointI])
<< " on some processors and in zone "
<< maxZone[pointI]
<< " on some other processors."
<< endl;
}
error = true;
}
}
return error;
}
void Foam::pointZone::writeDict(Ostream& os) const
{
os << nl << name_ << nl << token::BEGIN_BLOCK << nl

4
src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.H
View File

@ -208,6 +208,10 @@ public:
//- Check zone definition. Return true if in error.
virtual bool checkDefinition(const bool report = false) const;
//- Check whether zone is synchronised across coupled boundaries. Return
// true if in error.
virtual bool checkParallelSync(const bool report = false) const;
//- Correct patch after moving points
virtual void movePoints(const pointField&)
{}

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

@ -66,6 +66,7 @@ SourceFiles
#include "HashSet.H"
#include "Map.H"
#include "EdgeMap.H"
#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -577,14 +578,6 @@ public:
labelHashSet* setPtr = NULL
) const;
//- Check face-decomposition tet volume
bool checkFaceTets
(
const bool report = false,
const scalar minTetVol = 0,
labelHashSet* setPtr = NULL
) const;
//- Check face skewness
bool checkFaceSkewness
(

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

@ -618,114 +618,6 @@ bool Foam::primitiveMesh::checkFacePyramids
}
bool Foam::primitiveMesh::checkFaceTets
(
const bool report,
const scalar minTetVol,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool primitiveMesh::checkFaceTets("
<< "const bool, const scalar, labelHashSet*) const: "
<< "checking face orientation" << endl;
}
// check whether face area vector points to the cell with higher label
const vectorField& cc = cellCentres();
const vectorField& fc = faceCentres();
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const faceList& fcs = faces();
const pointField& p = points();
label nErrorPyrs = 0;
forAll(fcs, faceI)
{
// Create the owner tets - they will have negative volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[own[faceI]]
).mag();
if (tetVol > -minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break; // no need to check other tets
}
}
if (isInternalFace(faceI))
{
// Create the neighbour tet - it will have positive volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[nei[faceI]]
).mag();
if (tetVol < minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break;
}
}
}
}
reduce(nErrorPyrs, sumOp<label>());
if (nErrorPyrs > 0)
{
if (debug || report)
{
Info<< " ***Error in face tets: "
<< nErrorPyrs << " faces have incorrectly oriented face"
<< " decomposition triangles." << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Face tets OK." << endl;
}
return false;
}
}
bool Foam::primitiveMesh::checkFaceSkewness
(
const bool report,

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

@ -32,32 +32,14 @@ License
// Is the point in the cell bounding box
bool Foam::primitiveMesh::pointInCellBB(const point& p, label celli) const
{
const pointField& points = this->points();
const faceList& f = faces();
const vectorField& centres = cellCentres();
const cellList& cf = cells();
labelList cellVertices = cf[celli].labels(f);
vector bbmax = -GREAT*vector::one;
vector bbmin = GREAT*vector::one;
forAll(cellVertices, vertexI)
{
bbmax = max(bbmax, points[cellVertices[vertexI]]);
bbmin = min(bbmin, points[cellVertices[vertexI]]);
}
scalar distance = mag(centres[celli] - p);
if ((distance - mag(bbmax - bbmin)) < SMALL)
{
return true;
}
else
{
return false;
}
return boundBox
(
cells()[celli].points
(
faces(),
points()
)
).contains(p);
}

1
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
View File

@ -30,7 +30,6 @@ Description
#include "triPointRef.H"
#include "scalarField.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// (Probably very inefficient) minimum containment sphere calculation.

29
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
View File

@ -42,6 +42,7 @@ SourceFiles
#include "point.H"
#include "primitiveFieldsFwd.H"
#include "pointHit.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -133,16 +134,42 @@ public:
inline vector Sd() const;
//- Return centre (centroid)
inline Point centre() const;
//- Return volume
inline scalar mag() const;
//- Return circum-centre
inline vector circumCentre() const;
inline Point circumCentre() const;
//- Return circum-radius
inline scalar circumRadius() const;
//- Return quality: Ratio of tetrahedron and circum-sphere
// volume, scaled so that a regular tetrahedron has a
// quality of 1
inline scalar quality() const;
//- Return a random point in the tetrahedron from a
// uniform distribution
inline Point randomPoint(Random& rndGen) const;
//- Calculate the barycentric coordinates of the given
// point, in the same order as a, b, c, d. Returns the
// determinant of the solution.
inline scalar barycentric
(
const point& pt,
List<scalar>& bary
) const;
//- Return nearest point to p on tetrahedron
inline pointHit nearestPoint
(
const point& p
) const;
//- Return (min)containment sphere, i.e. the smallest sphere with
// all points inside. Returns pointHit with:
// - hit : if sphere is equal to circumsphere

209
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
View File

@ -125,6 +125,13 @@ inline vector tetrahedron<Point, PointRef>::Sd() const
}
template<class Point, class PointRef>
inline Point tetrahedron<Point, PointRef>::centre() const
{
return 0.25*(a_ + b_ + c_ + d_);
}
template<class Point, class PointRef>
inline scalar tetrahedron<Point, PointRef>::mag() const
{
@ -133,19 +140,19 @@ inline scalar tetrahedron<Point, PointRef>::mag() const
template<class Point, class PointRef>
inline vector tetrahedron<Point, PointRef>::circumCentre() const
inline Point tetrahedron<Point, PointRef>::circumCentre() const
{
vector a = b_ - a_;
vector b = c_ - a_;
vector c = d_ - a_;
scalar lamda = magSqr(c) - (a&c);
scalar mu = magSqr(b) - (a&b);
scalar lambda = magSqr(c) - (a & c);
scalar mu = magSqr(b) - (a & b);
vector ba = b^a;
vector ca = c^a;
vector ba = b ^ a;
vector ca = c ^ a;
return a_ + 0.5*(a + (lamda*ba - mu*ca)/(c&ba));
return a_ + 0.5*(a + (lambda*ba - mu*ca)/((c & ba) + ROOTVSMALL));
}
@ -156,6 +163,196 @@ inline scalar tetrahedron<Point, PointRef>::circumRadius() const
}
template<class Point, class PointRef>
inline scalar tetrahedron<Point, PointRef>::quality() const
{
// Note: 8/(9*sqrt(3)) = 0.5132002393
return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
}
template<class Point, class PointRef>
inline Point tetrahedron<Point, PointRef>::randomPoint(Random& rndGen) const
{
// Adapted from
// http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
scalar s = rndGen.scalar01();
scalar t = rndGen.scalar01();
scalar u = rndGen.scalar01();
if (s + t > 1.0)
{
s = 1.0 - s;
t = 1.0 - t;
}
if (t + u > 1.0)
{
scalar tmp = u;
u = 1.0 - s - t;
t = 1.0 - tmp;
}
else if (s + t + u > 1.0)
{
scalar tmp = u;
u = s + t + u - 1.0;
s = 1.0 - t - tmp;
}
return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
}
template<class Point, class PointRef>
scalar tetrahedron<Point, PointRef>::barycentric
(
const point& pt,
List<scalar>& bary
) const
{
// From:
// http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
vector e0(a_ - d_);
vector e1(b_ - d_);
vector e2(c_ - d_);
tensor t
(
e0.x(), e1.x(), e2.x(),
e0.y(), e1.y(), e2.y(),
e0.z(), e1.z(), e2.z()
);
scalar detT = det(t);
if (Foam::mag(detT) < SMALL)
{
WarningIn
(
"List<scalar> tetrahedron<Point, PointRef>::barycentric"
"("
"const point& pt"
") const"
)
<< "Degenerate triangle - returning 1/4 barycentric coordinates."
<< endl;
bary = List<scalar>(4, 0.25);
return detT;
}
vector res = inv(t, detT) & (pt - d_);
bary.setSize(4);
bary[0] = res.x();
bary[1] = res.y();
bary[2] = res.z();
bary[3] = (1.0 - res.x() - res.y() - res.z());
return detT;
}
template<class Point, class PointRef>
inline pointHit tetrahedron<Point, PointRef>::nearestPoint
(
const point& p
) const
{
// Adapted from:
// Real-time collision detection, Christer Ericson, 2005, p142-144
// Assuming initially that the point is inside all of the
// halfspaces, so closest to itself.
point closestPt = p;
scalar minOutsideDistance = VGREAT;
bool inside = true;
if (((p - b_) & Sa()) >= 0)
{
// p is outside halfspace plane of tri
pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
inside = false;
if (info.distance() < minOutsideDistance)
{
closestPt = info.rawPoint();
minOutsideDistance = info.distance();
}
}
if (((p - a_) & Sb()) >= 0)
{
// p is outside halfspace plane of tri
pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
inside = false;
if (info.distance() < minOutsideDistance)
{
closestPt = info.rawPoint();
minOutsideDistance = info.distance();
}
}
if (((p - a_) & Sc()) >= 0)
{
// p is outside halfspace plane of tri
pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
inside = false;
if (info.distance() < minOutsideDistance)
{
closestPt = info.rawPoint();
minOutsideDistance = info.distance();
}
}
if (((p - a_) & Sd()) >= 0)
{
// p is outside halfspace plane of tri
pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
inside = false;
if (info.distance() < minOutsideDistance)
{
closestPt = info.rawPoint();
minOutsideDistance = info.distance();
}
}
// If the point is inside, then the distance to the closest point
// is zero
if (inside)
{
minOutsideDistance = 0;
}
return pointHit
(
inside,
closestPt,
minOutsideDistance,
!inside
);
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class point, class pointRef>

21
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
View File

@ -39,7 +39,7 @@ SourceFiles
#include "vector.H"
#include "tensor.H"
#include "pointHit.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -141,12 +141,14 @@ public:
inline vector normal() const;
//- Return circum-centre
inline vector circumCentre() const;
inline Point circumCentre() const;
//- Return circum-radius
inline scalar circumRadius() const;
//- Return quality: Ratio triangle and circum-circle area
//- Return quality: Ratio of triangle and circum-circle
// area, scaled so that an equilateral triangle has a
// quality of 1
inline scalar quality() const;
//- Return swept-volume
@ -160,6 +162,19 @@ public:
scalar density = 1.0
) const;
//- Return a random point on the triangle from a uniform
// distribution
inline Point randomPoint(Random& rndGen) const;
//- Calculate the barycentric coordinates of the given
// point, in the same order as a, b, c. Returns the
// determinant of the solution.
inline scalar barycentric
(
const point& pt,
List<scalar>& bary
) const;
//- Return point intersection with a ray.
// For a hit, the distance is signed. Positive number
// represents the point in front of triangle.

73
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View File

@ -259,7 +259,7 @@ inline vector triangle<Point, PointRef>::normal() const
template<class Point, class PointRef>
inline vector triangle<Point, PointRef>::circumCentre() const
inline Point triangle<Point, PointRef>::circumCentre() const
{
scalar d1 = (c_ - a_)&(b_ - a_);
scalar d2 = -(c_ - b_)&(b_ - a_);
@ -303,12 +303,14 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
template<class Point, class PointRef>
inline scalar triangle<Point, PointRef>::quality() const
{
// Note: 3*sqr(3)/(4*pi) = 0.4134966716
return
mag()
/ (
constant::mathematical::pi
*Foam::sqr(circumRadius())
*0.413497
*0.4134966716
+ VSMALL
);
}
@ -366,6 +368,72 @@ inline tensor triangle<Point, PointRef>::inertia
*density;
}
template<class Point, class PointRef>
inline Point triangle<Point, PointRef>::randomPoint(Random& rndGen) const
{
// Generating Random Points in Triangles
// by Greg Turk
// from "Graphics Gems", Academic Press, 1990
// http://tog.acm.org/GraphicsGems/gems/TriPoints.c
scalar s = rndGen.scalar01();
scalar t = sqrt(rndGen.scalar01());
return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
}
template<class Point, class PointRef>
scalar triangle<Point, PointRef>::barycentric
(
const point& pt,
List<scalar>& bary
) const
{
// From:
// Real-time collision detection, Christer Ericson, 2005, p47-48
vector v0 = b_ - a_;
vector v1 = c_ - a_;
vector v2 = pt - a_;
scalar d00 = v0 & v0;
scalar d01 = v0 & v1;
scalar d11 = v1 & v1;
scalar d20 = v2 & v0;
scalar d21 = v2 & v1;
scalar denom = d00*d11 - d01*d01;
if (Foam::mag(denom) < SMALL)
{
WarningIn
(
"List<scalar> triangle<Point, PointRef>::barycentric"
"("
"const point& pt"
") const"
)
<< "Degenerate triangle - returning 1/3 barycentric coordinates."
<< endl;
bary = List<scalar>(3, 1.0/3.0);
return denom;
}
bary.setSize(3);
bary[0] = (d11*d20 - d01*d21)/denom;
bary[1] = (d00*d21 - d01*d20)/denom;
bary[2] = 1.0 - bary[0] - bary[1];
return denom;
}
template<class Point, class PointRef>
inline pointHit triangle<Point, PointRef>::ray
(
@ -775,4 +843,3 @@ inline Ostream& operator<<(Ostream& os, const triangle<Point, PointRef>& t)
} // End namespace Foam
// ************************************************************************* //

2
src/OpenFOAM/primitives/Tensor/TensorI.H
View File

@ -481,7 +481,7 @@ inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
}
//- Return the inverse of a tensor give the determinant
//- Return the inverse of a tensor given the determinant
template <class Cmpt>
inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
{

24
src/dynamicMesh/motionSmoother/motionSmootherCheck.C
View File

@ -68,9 +68,9 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
const scalar minTetQuality
(
readScalar(dict.lookup("minTetVol", true))
readScalar(dict.lookup("minTetQuality", true))
);
const scalar maxConcave
(
@ -160,12 +160,12 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
if (minTetQuality > -GREAT)
{
polyMeshGeometry::checkFaceTets
(
report,
minTetVol,
minTetQuality,
mesh,
mesh.cellCentres(),
mesh.faceCentres(),
@ -177,8 +177,8 @@ bool Foam::motionSmoother::checkMesh
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
Info<< " faces with face-decomposition tet quality < "
<< setw(5) << minTetQuality << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
@ -443,9 +443,9 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
const scalar minTetQuality
(
readScalar(dict.lookup("minTetVol", true))
readScalar(dict.lookup("minTetQuality", true))
);
const scalar maxConcave
(
@ -531,12 +531,12 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
if (minTetQuality > -GREAT)
{
meshGeom.checkFaceTets
(
report,
minTetVol,
minTetQuality,
meshGeom.mesh().points(),
checkFaces,
baffles,
@ -545,8 +545,8 @@ bool Foam::motionSmoother::checkMesh
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
Info<< " faces with face-decomposition tet quality < "
<< setw(5) << minTetQuality << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;

158
src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "polyMeshGeometry.H"
#include "polyMeshTetDecomposition.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "syncTools.H"
@ -313,7 +314,7 @@ bool Foam::polyMeshGeometry::checkFaceTet
(
const polyMesh& mesh,
const bool report,
const scalar minTetVol,
const scalar minTetQuality,
const pointField& p,
const label faceI,
const point& fc, // face centre
@ -326,15 +327,15 @@ bool Foam::polyMeshGeometry::checkFaceTet
forAll(f, fp)
{
scalar tetVol = tetPointRef
scalar tetQual = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc,
cc
).mag();
).quality();
if (tetVol < minTetVol)
if (tetQual < minTetQuality)
{
if (report)
{
@ -345,7 +346,7 @@ bool Foam::polyMeshGeometry::checkFaceTet
<< "face " << faceI
<< " has a triangle that points the wrong way."
<< endl
<< "Tet volume: " << tetVol
<< "Tet quality: " << tetQual
<< " Face " << faceI
<< endl;
}
@ -423,16 +424,15 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
// Calculate coupled cell centre
pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
}
syncTools::swapBoundaryFacePositions(mesh, neiCc);
syncTools::swapBoundaryFacePositions(mesh, neiCc);
scalar minDDotS = GREAT;
@ -754,7 +754,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
"polyMeshGeometry::checkFacePyramids("
"const bool, const scalar, const pointField&"
", const labelList&, labelHashSet*)"
) << "Error in face pyramids: faces pointing the wrong way!"
) << "Error in face pyramids: faces pointing the wrong way."
<< endl;
}
@ -775,7 +775,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
bool Foam::polyMeshGeometry::checkFaceTets
(
const bool report,
const scalar minTetVol,
const scalar minTetQuality,
const polyMesh& mesh,
const vectorField& cellCentres,
const vectorField& faceCentres,
@ -785,11 +785,23 @@ bool Foam::polyMeshGeometry::checkFaceTets
labelHashSet* setPtr
)
{
// check whether face area vector points to the cell with higher label
// check whether decomposing each cell into tets results in
// positive volume, non-flat tets
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label nErrorPyrs = 0;
// Calculate coupled cell centre
pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
neiCc[faceI - mesh.nInternalFaces()] = cellCentres[own[faceI]];
}
syncTools::swapBoundaryFacePositions(mesh, neiCc);
label nErrorTets = 0;
forAll(checkFaces, i)
{
@ -801,36 +813,104 @@ bool Foam::polyMeshGeometry::checkFaceTets
(
mesh,
report,
minTetVol,
minTetQuality,
p,
faceI,
cellCentres[own[faceI]], // face centre
faceCentres[faceI], // cell centre
faceCentres[faceI], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
nErrorTets++;
}
if (mesh.isInternalFace(faceI))
{
// Create the neighbour pyramid - it will have positive volume
// Create the neighbour tets - they will have positive volume
bool tetError = checkFaceTet
(
mesh,
report,
minTetVol,
minTetQuality,
p,
faceI,
faceCentres[faceI], // face centre
faceCentres[faceI], // face centre
cellCentres[nei[faceI]], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
nErrorTets++;
}
if
(
polyMeshTetDecomposition::findSharedBasePoint
(
mesh,
faceI,
minTetQuality,
report
) == -1
)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
else
{
label patchI = patches.whichPatch(faceI);
if (patches[patchI].coupled())
{
if
(
polyMeshTetDecomposition::findSharedBasePoint
(
mesh,
faceI,
neiCc[faceI - mesh.nInternalFaces()],
minTetQuality,
report
) == -1
)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
else
{
if
(
polyMeshTetDecomposition::findBasePoint
(
mesh,
faceI,
minTetQuality,
report
) == -1
)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorTets++;
}
}
}
}
@ -844,7 +924,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
(
mesh,
report,
minTetVol,
minTetQuality,
p,
face0,
cellCentres[own[face0]], // face centre
@ -854,15 +934,15 @@ bool Foam::polyMeshGeometry::checkFaceTets
if (tetError)
{
nErrorPyrs++;
nErrorTets++;
}
// Create the neighbour pyramid - it will have positive volume
// Create the neighbour tets - they will have positive volume
tetError = checkFaceTet
(
mesh,
report,
minTetVol,
minTetQuality,
p,
face0,
faceCentres[face0], // face centre
@ -872,13 +952,33 @@ bool Foam::polyMeshGeometry::checkFaceTets
if (tetError)
{
nErrorPyrs++;
nErrorTets++;
}
if
(
polyMeshTetDecomposition::findSharedBasePoint
(
mesh,
face0,
cellCentres[own[face1]],
minTetQuality,
report
) == -1
)
{
if (setPtr)
{
setPtr->insert(face0);
}
nErrorTets++;
}
}
reduce(nErrorPyrs, sumOp<label>());
reduce(nErrorTets, sumOp<label>());
if (nErrorPyrs > 0)
if (nErrorTets > 0)
{
if (report)
{
@ -887,7 +987,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
"polyMeshGeometry::checkFaceTets("
"const bool, const scalar, const pointField&, const pointField&"
", const labelList&, labelHashSet*)"
) << "Error in face pyramids: faces pointing the wrong way!"
) << "Error in face decomposition: negative tets."
<< endl;
}
@ -2135,7 +2235,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
bool Foam::polyMeshGeometry::checkFaceTets
(
const bool report,
const scalar minTetVol,
const scalar minTetQuality,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
@ -2145,7 +2245,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
return checkFaceTets
(
report,
minTetVol,
minTetQuality,
mesh_,
cellCentres_,
faceCentres_,

5
src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
View File

@ -113,7 +113,7 @@ class polyMeshGeometry
(
const polyMesh&,
const bool report,
const scalar minTetVol,
const scalar minTetQuality,
const pointField& p,
const label faceI,
const point& fc, // face centre
@ -121,6 +121,7 @@ class polyMeshGeometry
labelHashSet* setPtr
);
public:
ClassName("polyMeshGeometry");
@ -350,7 +351,7 @@ public:
bool checkFaceTets
(
const bool report,
const scalar minTetVol,
const scalar minTetQuality,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,

4
src/finiteVolume/fields/fvsPatchFields/constraint/cyclic/cyclicFvsPatchField.C
View File

@ -56,7 +56,7 @@ cyclicFvsPatchField<Type>::cyclicFvsPatchField
coupledFvsPatchField<Type>(ptf, p, iF, mapper),
cyclicPatch_(refCast<const cyclicFvPatch>(p))
{
if (!isType<cyclicFvPatch>(this->patch()))
if (!isA<cyclicFvPatch>(this->patch()))
{
FatalErrorIn
(
@ -87,7 +87,7 @@ cyclicFvsPatchField<Type>::cyclicFvsPatchField
coupledFvsPatchField<Type>(p, iF, dict),
cyclicPatch_(refCast<const cyclicFvPatch>(p))
{
if (!isType<cyclicFvPatch>(p))
if (!isA<cyclicFvPatch>(p))
{
FatalIOErrorIn
(

19
src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H
View File

@ -38,6 +38,7 @@ Description
#include "typeInfo.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "tetIndices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -130,9 +131,23 @@ public:
virtual Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const = 0;
//- Interpolate field to the given point in the tetrahedron
// defined by the given indices. Calls interpolate function
// above here execpt where overridden by derived
// interpolation types.
virtual Type interpolate
(
const vector& position,
const tetIndices& tetIs,
const label faceI = -1
) const
{
return interpolate(position, tetIs.cell(), faceI);
}
};

4
src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.C
View File

@ -49,11 +49,11 @@ template<class Type>
Type interpolationCell<Type>::interpolate
(
const vector&,
const label celli,
const label cellI,
const label
) const
{
return this->psi_[celli];
return this->psi_[cellI];
}

4
src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H
View File

@ -72,8 +72,8 @@ public:
Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const;
};

359
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
View File

@ -25,10 +25,13 @@ License
#include "cellPointWeight.H"
#include "polyMesh.H"
#include "tetPointRef.H"
#include "polyMeshTetDecomposition.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::cellPointWeight::debug(debug::debugSwitch("cellPointWeight", 0));
Foam::scalar Foam::cellPointWeight::tol(SMALL);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -37,148 +40,101 @@ void Foam::cellPointWeight::findTetrahedron
(
const polyMesh& mesh,
const vector& position,
const label cellIndex
const label cellI
)
{
if (debug)
{
Pout<< "\nFoam::cellPointWeight::findTetrahedron" << nl
Pout<< nl << "Foam::cellPointWeight::findTetrahedron" << nl
<< "position = " << position << nl
<< "cellIndex = " << cellIndex << endl;
<< "cellI = " << cellI << endl;
}
// Initialise closest triangle variables
scalar minUVWClose = VGREAT;
label pointIClose = 0;
label faceClose = 0;
List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
(
mesh,
cellI
);
const vector& P0 = mesh.cellCentres()[cellIndex];
const labelList& cellFaces = mesh.cells()[cellIndex];
const scalar cellVolume = mesh.cellVolumes()[cellIndex];
const faceList& pFaces = mesh.faces();
const scalar cellVolume = mesh.cellVolumes()[cellI];
// Find the tet that the point occupies
forAll(cellFaces, faceI)
forAll(cellTets, tetI)
{
// Decompose each face into triangles, making a tet when
// augmented by the cell centre
const labelList& facePoints = mesh.faces()[cellFaces[faceI]];
const tetIndices& tetIs = cellTets[tetI];
label pointI = 1;
while ((pointI + 1) < facePoints.size())
const face& f = pFaces[tetIs.face()];
// Barycentric coordinates of the position
scalar det = tetIs.tet(mesh).barycentric(position, weights_);
if (mag(det/cellVolume) > tol)
{
// Cartesian co-ordinates of the triangle vertices
const vector& P1 = mesh.points()[facePoints[0]];
const vector& P2 = mesh.points()[facePoints[pointI]];
const vector& P3 = mesh.points()[facePoints[pointI + 1]];
const scalar& u = weights_[0];
const scalar& v = weights_[1];
const scalar& w = weights_[2];
// Edge vectors
const vector e1 = P1 - P0;
const vector e2 = P2 - P0;
const vector e3 = P3 - P0;
// Solve for interpolation weighting factors
// Source term
const vector rhs = position - P0;
// Determinant of coefficients matrix
// Note: if det(A) = 0 the tet is degenerate
const scalar detA =
e1.x()*e2.y()*e3.z() + e2.x()*e3.y()*e1.z()
+ e3.x()*e1.y()*e2.z() - e1.x()*e3.y()*e2.z()
- e2.x()*e1.y()*e3.z() - e3.x()*e2.y()*e1.z();
if (mag(detA/cellVolume) > tol)
if
(
(u + tol > 0)
&& (v + tol > 0)
&& (w + tol > 0)
&& (u + v + w < 1 + tol)
)
{
// Solve using Cramers' rule
const scalar u =
(
rhs.x()*e2.y()*e3.z() + e2.x()*e3.y()*rhs.z()
+ e3.x()*rhs.y()*e2.z() - rhs.x()*e3.y()*e2.z()
- e2.x()*rhs.y()*e3.z() - e3.x()*e2.y()*rhs.z()
)/detA;
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
const scalar v =
(
e1.x()*rhs.y()*e3.z() + rhs.x()*e3.y()*e1.z()
+ e3.x()*e1.y()*rhs.z() - e1.x()*e3.y()*rhs.z()
- rhs.x()*e1.y()*e3.z() - e3.x()*rhs.y()*e1.z()
)/detA;
const scalar w =
(
e1.x()*e2.y()*rhs.z() + e2.x()*rhs.y()*e1.z()
+ rhs.x()*e1.y()*e2.z() - e1.x()*rhs.y()*e2.z()
- e2.x()*e1.y()*rhs.z() - rhs.x()*e2.y()*e1.z()
)/detA;
// Check if point is in tet
// value = 0 indicates position lies on a tet face
if
(
(u + tol > 0) && (v + tol > 0) && (w + tol > 0)
&& (u + v + w < 1 + tol)
)
{
faceVertices_[0] = facePoints[0];
faceVertices_[1] = facePoints[pointI];
faceVertices_[2] = facePoints[pointI + 1];
weights_[0] = u;
weights_[1] = v;
weights_[2] = w;
weights_[3] = 1.0 - (u + v + w);
return;
}
else
{
scalar minU = mag(u);
scalar minV = mag(v);
scalar minW = mag(w);
if (minU > 1.0)
{
minU -= 1.0;
}
if (minV > 1.0)
{
minV -= 1.0;
}
if (minW > 1.0)
{
minW -= 1.0;
}
const scalar minUVW = mag(minU + minV + minW);
if (minUVW < minUVWClose)
{
minUVWClose = minUVW;
pointIClose = pointI;
faceClose = faceI;
}
}
return;
}
}
}
pointI++;
// A suitable point in a tetrahedron was not found, find the
// nearest.
scalar minNearDist = VGREAT;
label nearestTetI = -1;
forAll(cellTets, tetI)
{
const tetIndices& tetIs = cellTets[tetI];
scalar nearDist = tetIs.tet(mesh).nearestPoint(position).distance();
if (nearDist < minNearDist)
{
minNearDist = nearDist;
nearestTetI = tetI;
}
}
if (debug)
{
Pout<< "cellPointWeight::findTetrahedron" << nl
<< " Tetrahedron search failed; using closest tet values to "
<< "point " << nl << " cell: " << cellIndex << nl << endl;
<< " Tetrahedron search failed; using closest tet to point "
<< position << nl
<< " cell: "
<< cellI << nl
<< endl;
}
const labelList& facePointsClose = mesh.faces()[cellFaces[faceClose]];
faceVertices_[0] = facePointsClose[0];
faceVertices_[1] = facePointsClose[pointIClose];
faceVertices_[2] = facePointsClose[pointIClose + 1];
weights_[0] = 0.25;
weights_[1] = 0.25;
weights_[2] = 0.25;
weights_[3] = 0.25;
const tetIndices& tetIs = cellTets[nearestTetI];
const face& f = pFaces[tetIs.face()];
// Barycentric coordinates of the position, ignoring if the
// determinant is suitable. If not, the return from barycentric
// to weights_ is safe.
tetIs.tet(mesh).barycentric(position, weights_);
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
}
@ -186,119 +142,112 @@ void Foam::cellPointWeight::findTriangle
(
const polyMesh& mesh,
const vector& position,
const label faceIndex
const label faceI
)
{
if (debug)
{
Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
<< "position = " << position << nl
<< "faceIndex = " << faceIndex << endl;
<< "faceI = " << faceI << endl;
}
// Initialise closest triangle variables
scalar minUVClose = VGREAT;
label pointIClose = 0;
List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
(
mesh,
mesh.faceOwner()[faceI],
faceI
);
// Decompose each face into triangles, making a tet when
// augmented by the cell centre
const labelList& facePoints = mesh.faces()[faceIndex];
const scalar faceAreaSqr = magSqr(mesh.faceAreas()[faceI]);
const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);
const face& f = mesh.faces()[faceI];
label pointI = 1;
while ((pointI + 1) < facePoints.size())
forAll(faceTets, tetI)
{
// Cartesian co-ordinates of the triangle vertices
const vector& P1 = mesh.points()[facePoints[0]];
const vector& P2 = mesh.points()[facePoints[pointI]];
const vector& P3 = mesh.points()[facePoints[pointI + 1]];
const tetIndices& tetIs = faceTets[tetI];
// Direction vectors
vector v1 = position - P1;
const vector v2 = P2 - P1;
const vector v3 = P3 - P1;
List<scalar> triWeights(3);
// Plane normal
vector n = v2 ^ v3;
n /= mag(n);
// Barycentric coordinates of the position
scalar det = tetIs.faceTri(mesh).barycentric(position, triWeights);
// Remove any offset to plane
v1 -= (n & v1)*v1;
// Helper variables
const scalar d12 = v1 & v2;
const scalar d13 = v1 & v3;
const scalar d22 = v2 & v2;
const scalar d23 = v2 & v3;
const scalar d33 = v3 & v3;
// Determinant of coefficients matrix
// Note: if det(A) = 0 the triangle is degenerate
const scalar detA = d22*d33 - d23*d23;
if (0.25*detA/faceArea2 > tol)
if (0.25*mag(det)/faceAreaSqr > tol)
{
// Solve using Cramers' rule
const scalar u = (d12*d33 - d23*d13)/detA;
const scalar v = (d22*d13 - d12*d23)/detA;
const scalar& u = triWeights[0];
const scalar& v = triWeights[1];
// Check if point is in triangle
if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
if
(
(u + tol > 0)
&& (v + tol > 0)
&& (u + v < 1 + tol)
)
{
// Indices of the cell vertices making up the triangle
faceVertices_[0] = facePoints[0];
faceVertices_[1] = facePoints[pointI];
faceVertices_[2] = facePoints[pointI + 1];
// Weight[0] is for the cell centre.
weights_[0] = 0;
weights_[1] = triWeights[0];
weights_[2] = triWeights[1];
weights_[3] = triWeights[2];
weights_[0] = u;
weights_[1] = v;
weights_[2] = 1.0 - (u + v);
weights_[3] = 0.0;
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
return;
}
else
{
scalar minU = mag(u);
scalar minV = mag(v);
if (minU > 1.0)
{
minU -= 1.0;
}
if (minV > 1.0)
{
minV -= 1.0;
}
const scalar minUV = mag(minU + minV);
if (minUV < minUVClose)
{
minUVClose = minUV;
pointIClose = pointI;
}
}
}
}
pointI++;
// A suitable point in a triangle was not found, find the nearest.
scalar minNearDist = VGREAT;
label nearestTetI = -1;
forAll(faceTets, tetI)
{
const tetIndices& tetIs = faceTets[tetI];
scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();
if (nearDist < minNearDist)
{
minNearDist = nearDist;
nearestTetI = tetI;
}
}
if (debug)
{
Pout<< "Foam::cellPointWeight::findTriangle"
<< "Triangle search failed; using closest triangle to point" << nl
<< " cell face: " << faceIndex << nl << endl;
Pout<< "cellPointWeight::findTriangle" << nl
<< " Triangle search failed; using closest tri to point "
<< position << nl
<< " face: "
<< faceI << nl
<< endl;
}
// Indices of the cell vertices making up the triangle
faceVertices_[0] = facePoints[0];
faceVertices_[1] = facePoints[pointIClose];
faceVertices_[2] = facePoints[pointIClose + 1];
const tetIndices& tetIs = faceTets[nearestTetI];
weights_[0] = 1.0/3.0;
weights_[1] = 1.0/3.0;
weights_[2] = 1.0/3.0;
weights_[3] = 0.0;
// Barycentric coordinates of the position, ignoring if the
// determinant is suitable. If not, the return from barycentric
// to triWeights is safe.
List<scalar> triWeights(3);
tetIs.faceTri(mesh).barycentric(position, triWeights);
// Weight[0] is for the cell centre.
weights_[0] = 0;
weights_[1] = triWeights[0];
weights_[2] = triWeights[1];
weights_[3] = triWeights[2];
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
}
@ -308,21 +257,23 @@ Foam::cellPointWeight::cellPointWeight
(
const polyMesh& mesh,
const vector& position,
const label cellIndex,
const label faceIndex
const label cellI,
const label faceI
)
:
cellIndex_(cellIndex)
cellI_(cellI),
weights_(4),
faceVertices_(3)
{
if (faceIndex < 0)
if (faceI < 0)
{
// Face data not supplied
findTetrahedron(mesh, position, cellIndex);
findTetrahedron(mesh, position, cellI);
}
else
{
// Face data supplied
findTriangle(mesh, position, faceIndex);
findTriangle(mesh, position, faceI);
}
}

20
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.H
View File

@ -55,13 +55,13 @@ protected:
// Protected data
//- Cell index
const label cellIndex_;
const label cellI_;
//- Weights applied to tet vertices
FixedList<scalar, 4> weights_;
List<scalar> weights_;
//- Face vertex indices
FixedList<label, 3> faceVertices_;
List<label> faceVertices_;
// Protected Member Functions
@ -70,14 +70,14 @@ protected:
(
const polyMesh& mesh,
const vector& position,
const label cellIndex
const label cellI
);
void findTriangle
(
const polyMesh& mesh,
const vector& position,
const label faceIndex
const label faceI
);
@ -98,8 +98,8 @@ public:
(
const polyMesh& mesh,
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
);
@ -108,17 +108,17 @@ public:
//- Cell index
inline label cell() const
{
return cellIndex_;
return cellI_;
}
//- interpolation weights
inline const FixedList<scalar, 4>& weights() const
inline const List<scalar>& weights() const
{
return weights_;
}
//- interpolation addressing for points on face
inline const FixedList<label, 3>& faceVertices() const
inline const List<label>& faceVertices() const
{
return faceVertices_;
}

15
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H
View File

@ -25,7 +25,7 @@ Class
Foam::interpolationCellPoint
Description
Given cell centre values and point (vertex) values decompose into
Given cell centre values and point (vertex) values decompose into
tetrahedra and linear interpolate within them.
\*---------------------------------------------------------------------------*/
@ -82,8 +82,17 @@ public:
inline Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const;
//- Interpolate field to the given point in the tetrahedron
// defined by the given indices.
inline Type interpolate
(
const vector& position,
const tetIndices& tetIs,
const label faceI = -1
) const;
};

73
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H
View File

@ -31,13 +31,13 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
const cellPointWeight& cpw
) const
{
const FixedList<scalar, 4>& weights = cpw.weights();
const FixedList<label, 3>& faceVertices = cpw.faceVertices();
const List<scalar>& weights = cpw.weights();
const List<label>& faceVertices = cpw.faceVertices();
Type t = psip_[faceVertices[0]]*weights[0];
t += psip_[faceVertices[1]]*weights[1];
t += psip_[faceVertices[2]]*weights[2];
t += this->psi_[cpw.cell()]*weights[3];
Type t = this->psi_[cpw.cell()]*weights[0];
t += psip_[faceVertices[0]]*weights[1];
t += psip_[faceVertices[1]]*weights[2];
t += psip_[faceVertices[2]]*weights[3];
return t;
}
@ -47,11 +47,66 @@ template<class Type>
inline Type Foam::interpolationCellPoint<Type>::interpolate
(
const vector& position,
const label celli,
const label facei
const label cellI,
const label faceI
) const
{
return interpolate(cellPointWeight(this->pMesh_, position, celli, facei));
return interpolate(cellPointWeight(this->pMesh_, position, cellI, faceI));
}
template<class Type>
inline Type Foam::interpolationCellPoint<Type>::interpolate
(
const vector& position,
const tetIndices& tetIs,
const label faceI
) const
{
// Assumes that the position is consistent with the supplied
// tetIndices. Does not pay attention to whether or not faceI is
// supplied or not - the result will be essentially the same.
// Performs a consistency check, however.
if (faceI >= 0)
{
if (faceI != tetIs.face())
{
FatalErrorIn
(
"inline Type Foam::interpolationCellPoint<Type>::interpolate"
"("
"const vector& position, "
"const tetIndices& tetIs, "
"const label faceI"
") const"
)
<< "specified face " << faceI << " inconsistent with the face "
<< "stored by tetIndices: " << tetIs.face()
<< exit(FatalError);
}
}
List<scalar> weights;
tetIs.tet(this->pMesh_).barycentric(position, weights);
const faceList& pFaces = this->pMesh_.faces();
const face& f = pFaces[tetIs.face()];
// Order of weights is the same as that of the vertices of the tet, i.e.
// cellCentre, faceBasePt, facePtA, facePtB.
Type t = this->psi_[tetIs.cell()]*weights[0];
t += psip_[f[tetIs.faceBasePt()]]*weights[1];
t += psip_[f[tetIs.facePtA()]]*weights[2];
t += psip_[f[tetIs.facePtB()]]*weights[3];
return t;
}

78
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C
View File

@ -56,8 +56,8 @@ template<class Type>
Type interpolationCellPointFace<Type>::interpolate
(
const vector& position,
const label nCell,
const label facei
const label cellI,
const label faceI
) const
{
Type ts[4];
@ -68,10 +68,10 @@ Type interpolationCellPointFace<Type>::interpolate
Type t = pTraits<Type>::zero;
// only use face information when the position is on a face
if (facei < 0)
if (faceI < 0)
{
const vector& cellCentre = this->pMesh_.cellCentres()[nCell];
const labelList& cellFaces = this->pMesh_.cells()[nCell];
const vector& cellCentre = this->pMesh_.cellCentres()[cellI];
const labelList& cellFaces = this->pMesh_.cells()[cellI];
vector projection = position - cellCentre;
tetPoints[3] = cellCentre;
@ -85,9 +85,9 @@ Type interpolationCellPointFace<Type>::interpolate
label closestFace = -1;
scalar minDistance = GREAT;
forAll(cellFaces, facei)
forAll(cellFaces, faceI)
{
label nFace = cellFaces[facei];
label nFace = cellFaces[faceI];
vector normal = this->pMeshFaceAreas_[nFace];
normal /= mag(normal);
@ -160,10 +160,10 @@ Type interpolationCellPointFace<Type>::interpolate
{
minDistance = GREAT;
label facei = 0;
while (facei < cellFaces.size() && !foundTet)
label faceI = 0;
while (faceI < cellFaces.size() && !foundTet)
{
label nFace = cellFaces[facei];
label nFace = cellFaces[faceI];
if (nFace < this->pMeshFaceAreas_.size())
{
foundTet = findTet
@ -179,7 +179,7 @@ Type interpolationCellPointFace<Type>::interpolate
minDistance
);
}
facei++;
faceI++;
}
}
@ -217,16 +217,16 @@ Type interpolationCellPointFace<Type>::interpolate
}
else
{
label patchi =
label patchI =
this->pMesh_.boundaryMesh().whichPatch(closestFace);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
if (this->psi_.boundaryField()[patchI].size())
{
ts[2] = this->psi_.boundaryField()[patchi]
ts[2] = this->psi_.boundaryField()[patchI]
[
this->pMesh_.boundaryMesh()[patchi].whichFace
this->pMesh_.boundaryMesh()[patchI].whichFace
(
closestFace
)
@ -234,11 +234,11 @@ Type interpolationCellPointFace<Type>::interpolate
}
else
{
ts[2] = this->psi_[nCell];
ts[2] = this->psi_[cellI];
}
}
ts[3] = this->psi_[nCell];
ts[3] = this->psi_[cellI];
for (label n=0; n<4; n++)
{
@ -251,9 +251,9 @@ Type interpolationCellPointFace<Type>::interpolate
else
{
Info<< "interpolationCellPointFace<Type>::interpolate"
<< "(const vector&, const label nCell) const : "
<< "(const vector&, const label cellI) const : "
<< "search failed; using closest cellFace value" << endl
<< "cell number " << nCell << tab
<< "cell number " << cellI << tab
<< "position " << position << endl;
if (closestFace < psis_.size())
@ -262,16 +262,16 @@ Type interpolationCellPointFace<Type>::interpolate
}
else
{
label patchi =
label patchI =
this->pMesh_.boundaryMesh().whichPatch(closestFace);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
if (this->psi_.boundaryField()[patchI].size())
{
t = this->psi_.boundaryField()[patchi]
t = this->psi_.boundaryField()[patchI]
[
this->pMesh_.boundaryMesh()[patchi].whichFace
this->pMesh_.boundaryMesh()[patchI].whichFace
(
closestFace
)
@ -279,7 +279,7 @@ Type interpolationCellPointFace<Type>::interpolate
}
else
{
t = this->psi_[nCell];
t = this->psi_[cellI];
}
}
}
@ -289,7 +289,7 @@ Type interpolationCellPointFace<Type>::interpolate
bool foundTriangle = findTriangle
(
position,
facei,
faceI,
tetPointLabels,
phi
);
@ -304,48 +304,48 @@ Type interpolationCellPointFace<Type>::interpolate
}
// ... and the face value
if (facei < psis_.size())
if (faceI < psis_.size())
{
t += phi[2]*psis_[facei];
t += phi[2]*psis_[faceI];
}
else
{
label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
if (this->psi_.boundaryField()[patchI].size())
{
t += phi[2]*this->psi_.boundaryField()[patchi]
[this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
t += phi[2]*this->psi_.boundaryField()[patchI]
[this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
}
else
{
t += phi[2]*this->psi_[nCell];
t += phi[2]*this->psi_[cellI];
}
}
}
else
{
// use face value only
if (facei < psis_.size())
if (faceI < psis_.size())
{
t = psis_[facei];
t = psis_[faceI];
}
else
{
label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
if (this->psi_.boundaryField()[patchI].size())
{
t = this->psi_.boundaryField()[patchi]
[this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
t = this->psi_.boundaryField()[patchI]
[this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
}
else
{
t = this->psi_[nCell];
t = this->psi_[cellI];
}
}
}

4
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.H
View File

@ -99,8 +99,8 @@ public:
Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const;
};

29
src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.C
View File

@ -24,9 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "cellPointWeightWallModified.H"
#include "wallPolyPatch.H"
#include "polyMesh.H"
#include "polyBoundaryMesh.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -34,36 +31,30 @@ Foam::cellPointWeightWallModified::cellPointWeightWallModified
(
const polyMesh& mesh,
const vector& position,
const label cellIndex,
const label faceIndex
const label cellI,
const label faceI
)
:
cellPointWeight(mesh, position, cellIndex, faceIndex)
cellPointWeight(mesh, position, cellI, faceI)
{
if (faceIndex < 0)
{
findTetrahedron(mesh, position, cellIndex);
}
else
// findTetrahedron or findTriangle will already have been called
// by the cellPointWeight constructor
if (faceI >= 0)
{
const polyBoundaryMesh& bm = mesh.boundaryMesh();
label patchI = bm.whichPatch(faceIndex);
label patchI = bm.whichPatch(faceI);
if (patchI != -1)
{
if (isA<wallPolyPatch>(bm[patchI]))
{
// Apply cell centre value wall faces
weights_[0] = 0.0;
weights_[0] = 1.0;
weights_[1] = 0.0;
weights_[2] = 0.0;
weights_[3] = 1.0;
weights_[3] = 0.0;
}
}
else
{
// Interpolate
findTriangle(mesh, position, faceIndex);
}
}
}

7
src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.H
View File

@ -36,6 +36,9 @@ SourceFiles
#define cellPointWeightWallModified_H
#include "cellPointWeight.H"
#include "wallPolyPatch.H"
#include "polyMesh.H"
#include "polyBoundaryMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -61,8 +64,8 @@ public:
(
const polyMesh& mesh,
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
);
};

13
src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H
View File

@ -74,8 +74,17 @@ public:
inline Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const;
//- Interpolate field to the given point in the tetrahedron
// defined by the given indices.
inline Type interpolate
(
const vector& position,
const tetIndices& tetIs,
const label faceI = -1
) const;
};

84
src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H
View File

@ -31,13 +31,13 @@ inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
const cellPointWeightWallModified& cpw
) const
{
const FixedList<scalar, 4>& weights = cpw.weights();
const FixedList<label, 3>& faceVertices = cpw.faceVertices();
const List<scalar>& weights = cpw.weights();
const List<label>& faceVertices = cpw.faceVertices();
Type t = this->psip_[faceVertices[0]]*weights[0];
t += this->psip_[faceVertices[1]]*weights[1];
t += this->psip_[faceVertices[2]]*weights[2];
t += this->psi_[cpw.cell()]*weights[3];
Type t = this->psi_[cpw.cell()]*weights[0];
t += this->psip_[faceVertices[0]]*weights[1];
t += this->psip_[faceVertices[1]]*weights[2];
t += this->psip_[faceVertices[2]]*weights[3];
return t;
}
@ -47,21 +47,73 @@ template<class Type>
inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
(
const vector& position,
const label celli,
const label facei
const label cellI,
const label faceI
) const
{
return
interpolate
return interpolate
(
cellPointWeightWallModified
(
cellPointWeightWallModified
this->pMesh_,
position,
cellI,
faceI
)
);
}
template<class Type>
inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
(
const vector& position,
const tetIndices& tetIs,
const label faceI
) const
{
if (faceI >= 0)
{
if (faceI != tetIs.face())
{
FatalErrorIn
(
this->pMesh_,
position,
celli,
facei
"inline Type "
"Foam::interpolationCellPointWallModifie<Type>::interpolate"
"("
"const vector& position, "
"const tetIndices& tetIs, "
"const label faceI"
") const"
)
);
<< "specified face " << faceI << " inconsistent with the face "
<< "stored by tetIndices: " << tetIs.face()
<< exit(FatalError);
}
const polyBoundaryMesh& bm = this->pMesh_.boundaryMesh();
label patchI = bm.whichPatch(faceI);
if (patchI != -1)
{
if (isA<wallPolyPatch>(bm[patchI]))
{
Type t = this->psi_[tetIs.cell()];
return t;
}
}
}
// If the wall face selection did not return, then use the normal
// interpolate method
return interpolationCellPoint<Type>::interpolate
(
position,
tetIs,
faceI
);
}

4
src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H
View File

@ -81,8 +81,8 @@ public:
inline Type interpolate
(
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
) const;
};

6
src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H
View File

@ -39,13 +39,13 @@ template<class Type>
inline Type Foam::interpolationPoint<Type>::interpolate
(
const vector& position,
const label celli,
const label facei
const label cellI,
const label faceI
) const
{
return interpolate
(
pointMVCWeight(this->pMesh_, position, celli, facei)
pointMVCWeight(this->pMesh_, position, cellI, faceI)
);
}

4
src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H
View File

@ -119,8 +119,8 @@ public:
(
const polyMesh& mesh,
const vector& position,
const label nCell,
const label facei = -1
const label cellI,
const label faceI = -1
);

61
src/finiteVolume/interpolation/surfaceInterpolation/schemes/localBlended/localBlended.H
View File

@ -153,6 +153,67 @@ public:
blendingFactor*tScheme1_().interpolate(vf)
+ (scalar(1) - blendingFactor)*tScheme2_().interpolate(vf);
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return tScheme1_().corrected() || tScheme2_().corrected();
}
//- Return the explicit correction to the face-interpolate
// for the given field
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
correction
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
const surfaceScalarField& blendingFactor =
this->mesh().objectRegistry::
lookupObject<const surfaceScalarField>
(
word(vf.name() + "BlendingFactor")
);
if (tScheme1_().corrected())
{
if (tScheme2_().corrected())
{
return
(
blendingFactor
* tScheme1_().correction(vf)
+ (scalar(1.0) - blendingFactor)
* tScheme2_().correction(vf)
);
}
else
{
return
(
blendingFactor
* tScheme1_().correction(vf)
);
}
}
else if (tScheme2_().corrected())
{
return
(
(scalar(1.0) - blendingFactor)
* tScheme2_().correction(vf)
);
}
else
{
return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
(
NULL
);
}
}
};

340
src/lagrangian/basic/Cloud/Cloud.C
View File

@ -30,6 +30,41 @@ License
#include "mapPolyMesh.H"
#include "Time.H"
#include "OFstream.H"
#include "wallPolyPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::trackingCorrectionTol = 1e-5;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class ParticleType>
void Foam::Cloud<ParticleType>::calcCellWallFaces() const
{
cellWallFacesPtr_.reset(new PackedBoolList(pMesh().nCells(), false));
PackedBoolList& cellWallFaces = cellWallFacesPtr_();
const polyBoundaryMesh& patches = pMesh().boundaryMesh();
forAll(patches, patchI)
{
if (isA<wallPolyPatch>(patches[patchI]))
{
const polyPatch& patch = patches[patchI];
const labelList& pFaceCells = patch.faceCells();
forAll(pFaceCells, pFCI)
{
cellWallFaces[pFaceCells[pFCI]] = true;
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -43,7 +78,11 @@ Foam::Cloud<ParticleType>::Cloud
cloud(pMesh),
IDLList<ParticleType>(),
polyMesh_(pMesh),
particleCount_(0)
particleCount_(0),
labels_(),
cellTree_(),
nTrackingRescues_(),
cellWallFacesPtr_()
{
IDLList<ParticleType>::operator=(particles);
}
@ -60,7 +99,11 @@ Foam::Cloud<ParticleType>::Cloud
cloud(pMesh, cloudName),
IDLList<ParticleType>(),
polyMesh_(pMesh),
particleCount_(0)
particleCount_(0),
labels_(),
cellTree_(),
nTrackingRescues_(),
cellWallFacesPtr_()
{
IDLList<ParticleType>::operator=(particles);
}
@ -68,6 +111,250 @@ Foam::Cloud<ParticleType>::Cloud
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParticleType>
void Foam::Cloud<ParticleType>::findCellFacePt
(
const point& pt,
label& cellI,
label& tetFaceI,
label& tetPtI
) const
{
cellI = -1;
tetFaceI = -1;
tetPtI = -1;
const indexedOctree<treeDataCell>& tree = cellTree();
// Find nearest cell to the point
pointIndexHit info = tree.findNearest(pt, sqr(GREAT));
if (info.hit())
{
label nearestCellI = tree.shapes().cellLabels()[info.index()];
// Check the nearest cell to see if the point is inside.
findFacePt(nearestCellI, pt, tetFaceI, tetPtI);
if (tetFaceI != -1)
{
// Point was in the nearest cell
cellI = nearestCellI;
return;
}
else
{
// Check the other possible cells that the point may be in
labelList testCells = tree.findIndices(pt);
forAll(testCells, pCI)
{
label testCellI = tree.shapes().cellLabels()[testCells[pCI]];
if (testCellI == nearestCellI)
{
// Don't retest the nearest cell
continue;
}
// Check the test cell to see if the point is inside.
findFacePt(testCellI, pt, tetFaceI, tetPtI);
if (tetFaceI != -1)
{
// Point was in the test cell
cellI = testCellI;
return;
}
}
}
}
else
{
FatalErrorIn
(
"void Foam::Cloud<ParticleType>::findCellFacePt"
"("
"const point& pt, "
"label& cellI, "
"label& tetFaceI, "
"label& tetPtI"
") const"
) << "Did not find nearest cell in search tree."
<< abort(FatalError);
}
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::findFacePt
(
label cellI,
const point& pt,
label& tetFaceI,
label& tetPtI
) const
{
tetFaceI = -1;
tetPtI = -1;
List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
(
polyMesh_,
cellI
);
forAll(cellTets, tetI)
{
const tetIndices& cellTetIs = cellTets[tetI];
if (inTet(pt, cellTetIs.tet(polyMesh_)))
{
tetFaceI = cellTetIs.face();
tetPtI = cellTetIs.tetPt();
return;
}
}
}
template<class ParticleType>
bool Foam::Cloud<ParticleType>::inTet
(
const point& pt,
const tetPointRef& tet
) const
{
// For robustness, assuming that the point is in the tet unless
// "definitively" shown otherwise by obtaining a positive dot
// product greater than a tolerance of SMALL.
// The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
// vectors and base points for the half-space planes are:
// area[0] = tet.Sa();
// area[1] = tet.Sb();
// area[2] = tet.Sc();
// area[3] = tet.Sd();
// planeBase[0] = tetBasePt = tet.b()
// planeBase[1] = ptA = tet.c()
// planeBase[2] = tetBasePt = tet.b()
// planeBase[3] = tetBasePt = tet.b()
vector n = vector::zero;
{
// 0, a
const point& basePt = tet.b();
n = tet.Sa();
n /= (mag(n) + VSMALL);
if (((pt - basePt) & n) > SMALL)
{
return false;
}
}
{
// 1, b
const point& basePt = tet.c();
n = tet.Sb();
n /= (mag(n) + VSMALL);
if (((pt - basePt) & n) > SMALL)
{
return false;
}
}
{
// 2, c
const point& basePt = tet.b();
n = tet.Sc();
n /= (mag(n) + VSMALL);
if (((pt - basePt) & n) > SMALL)
{
return false;
}
}
{
// 3, d
const point& basePt = tet.b();
n = tet.Sd();
n /= (mag(n) + VSMALL);
if (((pt - basePt) & n) > SMALL)
{
return false;
}
}
return true;
}
template<class ParticleType>
const Foam::indexedOctree<Foam::treeDataCell>&
Foam::Cloud<ParticleType>::cellTree() const
{
if (cellTree_.empty())
{
treeBoundBox overallBb(polyMesh_.points());
Random rndGen(261782);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
cellTree_.reset
(
new indexedOctree<treeDataCell>
(
treeDataCell
(
false, // not cache bb
polyMesh_
),
overallBb,
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
return cellTree_();
}
template<class ParticleType>
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::cellHasWallFaces()
const
{
if (!cellWallFacesPtr_.valid())
{
calcCellWallFaces();
}
return cellWallFacesPtr_();
}
template<class ParticleType>
Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
{
@ -131,6 +418,9 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
pIter().stepFraction() = 0;
}
// Reset nTrackingRescues
nTrackingRescues_ = 0;
// While there are particles to transfer
while (true)
{
@ -162,29 +452,29 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
{
// If we are running in parallel and the particle is on a
// boundary face
if (Pstream::parRun() && p.facei_ >= pMesh().nInternalFaces())
if (Pstream::parRun() && p.faceI_ >= pMesh().nInternalFaces())
{
label patchi = pbm.whichPatch(p.facei_);
label patchI = pbm.whichPatch(p.faceI_);
// ... and the face is on a processor patch
// prepare it for transfer
if (procPatchIndices[patchi] != -1)
if (procPatchIndices[patchI] != -1)
{
label n = neighbourProcIndices
[
refCast<const processorPolyPatch>
(
pbm[patchi]
pbm[patchI]
).neighbProcNo()
];
p.prepareForParallelTransfer(patchi, td);
p.prepareForParallelTransfer(patchI, td);
particleTransferLists[n].append(this->remove(&p));
patchIndexTransferLists[n].append
(
procPatchNeighbours[patchi]
procPatchNeighbours[patchI]
);
}
}
@ -270,15 +560,22 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
{
ParticleType& newp = newpIter();
label patchi = procPatches[receivePatchIndex[pI++]];
label patchI = procPatches[receivePatchIndex[pI++]];
newp.correctAfterParallelTransfer(patchi, td);
newp.correctAfterParallelTransfer(patchI, td);
addParticle(newParticles.remove(&newp));
}
}
}
}
reduce(nTrackingRescues_, sumOp<label>());
if (nTrackingRescues_ > 0)
{
Info<< nTrackingRescues_ << " tracking rescue corrections" << endl;
}
}
@ -294,24 +591,30 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
const labelList& reverseCellMap = mapper.reverseCellMap();
const labelList& reverseFaceMap = mapper.reverseFaceMap();
// Reset stored data that relies on the mesh
cellTree_.clear();
cellWallFacesPtr_.clear();
forAllIter(typename Cloud<ParticleType>, *this, pIter)
{
if (reverseCellMap[pIter().celli_] >= 0)
if (reverseCellMap[pIter().cellI_] >= 0)
{
pIter().celli_ = reverseCellMap[pIter().celli_];
pIter().cellI_ = reverseCellMap[pIter().cellI_];
if (pIter().facei_ >= 0 && reverseFaceMap[pIter().facei_] >= 0)
if (pIter().faceI_ >= 0 && reverseFaceMap[pIter().faceI_] >= 0)
{
pIter().facei_ = reverseFaceMap[pIter().facei_];
pIter().faceI_ = reverseFaceMap[pIter().faceI_];
}
else
{
pIter().facei_ = -1;
pIter().faceI_ = -1;
}
pIter().initCellFacePt();
}
else
{
label trackStartCell = mapper.mergedCell(pIter().celli_);
label trackStartCell = mapper.mergedCell(pIter().cellI_);
if (trackStartCell < 0)
{
@ -319,9 +622,14 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
}
vector p = pIter().position();
const_cast<vector&>(pIter().position()) =
polyMesh_.cellCentres()[trackStartCell];
pIter().stepFraction() = 0;
pIter().initCellFacePt();
pIter().track(p);
}
}

95
src/lagrangian/basic/Cloud/Cloud.H
View File

@ -40,6 +40,11 @@ SourceFiles
#include "IOField.H"
#include "IOFieldField.H"
#include "polyMesh.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "tetPointRef.H"
#include "polyMeshTetDecomposition.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -78,15 +83,28 @@ class Cloud
//- Overall count of particles ever created. Never decreases.
mutable label particleCount_;
//- Temporary storage for addressing. Used in findFaces.
//- Temporary storage for addressing. Used in findTris.
mutable DynamicList<label> labels_;
//- Search tree to allow spatial tet searching
mutable autoPtr<indexedOctree<treeDataCell> > cellTree_;
//- Count of how many tracking rescue corrections have been
// applied
mutable label nTrackingRescues_;
//- Does the cell have wall faces
mutable autoPtr<PackedBoolList> cellWallFacesPtr_;
// Private Member Functions
//- Initialise cloud on IO constructor
void initCloud(const bool checkClass);
//- Find all cells which have wall faces
void calcCellWallFaces() const;
//- Read cloud properties dictionary
void readCloudUniformProperties();
@ -115,6 +133,10 @@ public:
//- Name of cloud properties dictionary
static word cloudPropertiesName;
//- Fraction of distance to tet centre to move a particle to
// 'rescue' it from a tracking problem
static const scalar trackingCorrectionTol;
// Constructors
@ -163,27 +185,27 @@ public:
}
//- Is this global face an internal face?
bool internalFace(const label facei) const
bool internalFace(const label faceI) const
{
return polyMesh_.isInternalFace(facei);
return polyMesh_.isInternalFace(faceI);
}
//- Is this global face a boundary face?
bool boundaryFace(const label facei) const
bool boundaryFace(const label faceI) const
{
return !internalFace(facei);
return !internalFace(faceI);
}
//- Which patch is this global face on
label facePatch(const label facei) const
label facePatch(const label faceI) const
{
return polyMesh_.boundaryMesh().whichPatch(facei);
return polyMesh_.boundaryMesh().whichPatch(faceI);
}
//- Which face of this patch is this global face
label patchFace(const label patchi, const label facei) const
label patchFace(const label patchI, const label faceI) const
{
return polyMesh_.boundaryMesh()[patchi].whichFace(facei);
return polyMesh_.boundaryMesh()[patchI].whichFace(faceI);
}
label size() const
@ -191,6 +213,61 @@ public:
return IDLList<ParticleType>::size();
};
//- Find the cell, tetFaceI and tetPtI for the given
// position
void findCellFacePt
(
const point& pt,
label& cellI,
label& tetFaceI,
label& tetPtI
) const;
//- Find the tetFaceI and tetPtI for the given position in
// the supplied cell, tetFaceI and tetPtI = -1 if not
// found
void findFacePt
(
label cellI,
const point& pt,
label& tetFaceI,
label& tetPtI
) const;
//- Test if the given position is inside the give tet
bool inTet
(
const point& pt,
const tetPointRef& tet
) const;
//- Build (if necessary) and return the cell search tree
const indexedOctree<treeDataCell>& cellTree() const;
//- Return nTrackingRescues
label nTrackingRescues() const
{
return nTrackingRescues_;
}
//- Increment the nTrackingRescues counter
void trackingRescue() const
{
nTrackingRescues_++;
}
//- Whether each cell has any wall faces (demand driven data)
const PackedBoolList& cellHasWallFaces() const;
//- Switch to specify if particles of the cloud can return
// non-zero wall distance values. By default, assume
// that they can't (default for wallImpactDistance in
// Particle is 0.0).
virtual bool hasWallImpactDistance() const
{
return false;
}
// Iterators

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

@ -126,6 +126,13 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
<< " " << ioP.path() << nl
<< " assuming the initial cloud contains 0 particles." << endl;
}
forAllIter(typename Cloud<ParticleType>, *this, pIter)
{
ParticleType& p = pIter();
p.initCellFacePt();
}
}
@ -140,7 +147,11 @@ Foam::Cloud<ParticleType>::Cloud
:
cloud(pMesh),
polyMesh_(pMesh),
particleCount_(0)
particleCount_(0),
labels_(),
cellTree_(),
nTrackingRescues_(),
cellWallFacesPtr_()
{
initCloud(checkClass);
}
@ -156,7 +167,11 @@ Foam::Cloud<ParticleType>::Cloud
:
cloud(pMesh, cloudName),
polyMesh_(pMesh),
particleCount_(0)
particleCount_(0),
labels_(),
cellTree_(),
nTrackingRescues_(),
cellWallFacesPtr_()
{
initCloud(checkClass);
}

5
src/lagrangian/basic/Make/options
View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
-lmeshTools

665
src/lagrangian/basic/Particle/Particle.C
View File

@ -29,73 +29,20 @@ License
#include "symmetryPolyPatch.H"
#include "cyclicPolyPatch.H"
#include "processorPolyPatch.H"
#include "wallPolyPatch.H"
#include "transform.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class ParticleType>
void Foam::Particle<ParticleType>::findFaces
(
const vector& position,
DynamicList<label>& faceList
) const
{
const polyMesh& mesh = cloud_.polyMesh_;
const labelList& faces = mesh.cells()[celli_];
const vector& C = mesh.cellCentres()[celli_];
faceList.clear();
forAll(faces, i)
{
label facei = faces[i];
scalar lam = lambda(C, position, facei);
if ((lam > 0) && (lam < 1.0))
{
faceList.append(facei);
}
}
}
template<class ParticleType>
void Foam::Particle<ParticleType>::findFaces
(
const vector& position,
const label celli,
const scalar stepFraction,
DynamicList<label>& faceList
) const
{
const polyMesh& mesh = cloud_.pMesh();
const labelList& faces = mesh.cells()[celli];
const vector& C = mesh.cellCentres()[celli];
faceList.clear();
forAll(faces, i)
{
label facei = faces[i];
scalar lam = lambda(C, position, facei, stepFraction);
if ((lam > 0) && (lam < 1.0))
{
faceList.append(facei);
}
}
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::prepareForParallelTransfer
(
const label patchi,
const label patchI,
TrackData& td
)
{
// Convert the face index to be local to the processor patch
facei_ = patchFace(patchi, facei_);
faceI_ = patchFace(patchI, faceI_);
}
@ -103,15 +50,15 @@ template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::correctAfterParallelTransfer
(
const label patchi,
const label patchI,
TrackData& td
)
{
const processorPolyPatch& ppp =
refCast<const processorPolyPatch>
(cloud_.pMesh().boundaryMesh()[patchi]);
(cloud_.pMesh().boundaryMesh()[patchI]);
celli_ = ppp.faceCells()[facei_];
cellI_ = ppp.faceCells()[faceI_];
if (!ppp.parallel())
{
@ -123,7 +70,7 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
}
else
{
const tensor& T = ppp.forwardT()[facei_];
const tensor& T = ppp.forwardT()[faceI_];
transformPosition(T);
static_cast<ParticleType&>(*this).transformProperties(T);
}
@ -140,23 +87,48 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
}
else
{
position_ -= ppp.separation()[facei_];
position_ -= ppp.separation()[faceI_];
static_cast<ParticleType&>(*this).transformProperties
(
-ppp.separation()[facei_]
-ppp.separation()[faceI_]
);
}
}
tetFaceI_ = faceI_ + ppp.start();
// Faces either side of a coupled patch have matched base indices,
// tetPtI is specified relative to the base point, already and
// opposite circulation directions by design, so if the vertices
// are:
// source:
// face (a b c d e f)
// fPtI 0 1 2 3 4 5
// +
// destination:
// face (a f e d c b)
// fPtI 0 1 2 3 4 5
// +
// where a is the base point of the face are matching , and we
// have fPtI = 1 on the source processor face, i.e. vertex b, then
// this because of the face circulation direction change, vertex c
// is the characterising point on the destination processor face,
// giving the destination fPtI as:
// fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
// This relationship can be verified for other points and sizes of
// face.
tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
// Reset the face index for the next tracking operation
if (stepFraction_ > (1.0 - SMALL))
{
stepFraction_ = 1.0;
facei_ = -1;
faceI_ = -1;
}
else
{
facei_ += ppp.start();
faceI_ += ppp.start();
}
}
@ -168,27 +140,59 @@ Foam::Particle<ParticleType>::Particle
(
const Cloud<ParticleType>& cloud,
const vector& position,
const label celli
const label cellI,
const label tetFaceI,
const label tetPtI
)
:
cloud_(cloud),
position_(position),
celli_(celli),
facei_(-1),
cellI_(cellI),
faceI_(-1),
stepFraction_(0.0),
tetFaceI_(tetFaceI),
tetPtI_(tetPtI),
origProc_(Pstream::myProcNo()),
origId_(cloud_.getNewParticleID())
{}
template<class ParticleType>
Foam::Particle<ParticleType>::Particle
(
const Cloud<ParticleType>& cloud,
const vector& position,
const label cellI,
bool doCellFacePt
)
:
cloud_(cloud),
position_(position),
cellI_(cellI),
faceI_(-1),
stepFraction_(0.0),
tetFaceI_(-1),
tetPtI_(-1),
origProc_(Pstream::myProcNo()),
origId_(cloud_.getNewParticleID())
{
if (doCellFacePt)
{
initCellFacePt();
}
}
template<class ParticleType>
Foam::Particle<ParticleType>::Particle(const Particle<ParticleType>& p)
:
cloud_(p.cloud_),
position_(p.position_),
celli_(p.celli_),
facei_(p.facei_),
cellI_(p.cellI_),
faceI_(p.faceI_),
stepFraction_(p.stepFraction_),
tetFaceI_(p.tetFaceI_),
tetPtI_(p.tetPtI_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@ -204,7 +208,7 @@ Foam::label Foam::Particle<ParticleType>::track
TrackData& td
)
{
facei_ = -1;
faceI_ = -1;
// Tracks to endPosition or stop on boundary
while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
@ -212,7 +216,7 @@ Foam::label Foam::Particle<ParticleType>::track
stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
}
return facei_;
return faceI_;
}
@ -224,6 +228,7 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
return track(endPosition, dummyTd);
}
template<class ParticleType>
template<class TrackData>
Foam::scalar Foam::Particle<ParticleType>::trackToFace
@ -234,175 +239,414 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
{
const polyMesh& mesh = cloud_.polyMesh_;
DynamicList<label>& faces = cloud_.labels_;
findFaces(endPosition, faces);
const faceList& pFaces = mesh.faces();
const pointField& pPts = mesh.points();
const vectorField& pC = mesh.cellCentres();
faceI_ = -1;
facei_ = -1;
scalar trackFraction = 0.0;
if (faces.empty()) // inside cell
{
trackFraction = 1.0;
position_ = endPosition;
}
else // hit face
{
scalar lambdaMin = GREAT;
// Minimum tetrahedron decomposition of each cell of the mesh into
// using the cell centre, base point on face, and further two
// points on the face. For each face of n points, there are n - 2
// tets generated.
if (faces.size() == 1)
// The points for each tet are organised to match those used in the
// tetrahedron class, supplying them in the order:
// Cc, basePt, pA, pB
// where:
// + Cc is the cell centre;
// + basePt is the base point on the face;
// + pA and pB are the remaining points on the face, such that
// the circulation, {basePt, pA, pB} produces a positive
// normal by the right-hand rule. pA and pB are chosen from
// tetPtI_ do accomplish this depending if the cell owns the
// face, tetPtI_ is the vertex that characterises the tet, and
// is the first vertex on the tet when circulating around the
// face. Therefore, the same tetPtI represents the same face
// triangle for both the owner and neighbour cell.
//
// Each tet has its four triangles represented in the same order:
// 0) tri joining a tet to the tet across the face in next cell.
// This is the triangle opposite Cc.
// 1) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the edge of the current face
// that doesn't contain basePt. This is the triangle opposite
// basePt.
// 2) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the tet-edge (basePt - pB).
// This may be on the same face, or a different one. This is
// the triangle opposite basePt. This is the triangle opposite
// pA.
// 4) tri joining a tet to the tet that is in the same cell, but
// belongs to the face that shares the tet-edge (basePt - pA).
// This may be on the same face, or a different one. This is
// the triangle opposite basePt. This is the triangle opposite
// pA.
// Which tri (0..3) of the tet has been crossed
label triI = -1;
// Determine which face was actually crossed. lambdaMin < SMALL
// is considered a trigger for a tracking correction towards the
// current tet centre.
scalar lambdaMin = VGREAT;
DynamicList<label>& tris = cloud_.labels_;
// Tet indices that will be set by hitWallFaces if a wall face is
// to be hit, or are set when any wall tri of a tet is hit.
// Carries the description of the tet on which the cell face has
// been hit. For the case of being set in hitWallFaces, this may
// be a different tet to the one that the particle occupies.
tetIndices faceHitTetIs;
do
{
if (triI != -1)
{
lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
facei_ = faces[0];
// Change tet ownership because a tri face has been crossed
tetNeighbour(triI);
}
const Foam::face& f = pFaces[tetFaceI_];
bool own = (mesh.faceOwner()[tetFaceI_] == cellI_);
label tetBasePtI = mesh.tetBasePtIs()[tetFaceI_];
label basePtI = f[tetBasePtI];
label facePtI = (tetPtI_ + tetBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label fPtAI = -1;
label fPtBI = -1;
if (own)
{
fPtAI = facePtI;
fPtBI = otherFacePtI;
}
else
{
// If the particle has to cross more than one cell to reach the
// endPosition, we check which way to go.
// If one of the faces is a boundary face and the particle is
// outside, we choose the boundary face.
// The particle is outside if one of the lambda's is > 1 or < 0
forAll(faces, i)
fPtAI = otherFacePtI;
fPtBI = facePtI;
}
tetPointRef tet
(
pC[cellI_],
pPts[basePtI],
pPts[f[fPtAI]],
pPts[f[fPtBI]]
);
if (lambdaMin < SMALL)
{
// Apply tracking correction towards tet centre
position_ +=
Cloud<ParticleType>::trackingCorrectionTol
*(tet.centre() - position_);
cloud_.trackingRescue();
return trackFraction;
}
if (triI != -1 && mesh.moving())
{
// Mesh motion requires stepFraction to be correct for
// each tracking portion, so trackToFace must return after
// every lambda calculation.
return trackFraction;
}
FixedList<vector, 4> tetAreas;
tetAreas[0] = tet.Sa();
tetAreas[1] = tet.Sb();
tetAreas[2] = tet.Sc();
tetAreas[3] = tet.Sd();
FixedList<label, 4> tetPlaneBasePtIs;
tetPlaneBasePtIs[0] = basePtI;
tetPlaneBasePtIs[1] = f[fPtAI];
tetPlaneBasePtIs[2] = basePtI;
tetPlaneBasePtIs[3] = basePtI;
findTris(endPosition, tris, tet, tetAreas, tetPlaneBasePtIs);
// Reset variables for new track
triI = -1;
lambdaMin = VGREAT;
// Sets a value for lambdaMin and faceI_ if a wall face is hit
// by the track.
hitWallFaces(position_, endPosition, lambdaMin, faceHitTetIs);
// Did not hit any tet tri faces, and no wall face has been
// found to hit.
if (tris.empty() && faceI_ < 0)
{
position_ = endPosition;
return 1.0;
}
else
{
// Loop over all found tris and see if any of them find a
// lambda value smaller than that found for a wall face.
forAll(tris, i)
{
scalar lam =
lambda(position_, endPosition, faces[i], stepFraction_);
label tI = tris[i];
scalar lam = tetLambda
(
position_,
endPosition,
triI,
tetAreas[tI],
tetPlaneBasePtIs[tI],
cellI_,
tetFaceI_,
tetPtI_
);
if (lam < lambdaMin)
{
lambdaMin = lam;
facei_ = faces[i];
triI = tI;
}
}
}
bool internalFace = cloud_.internalFace(facei_);
if (triI == 0)
{
// This must be a cell face crossing
faceI_ = tetFaceI_;
// For warped faces the particle can be 'outside' the cell.
// This will yield a lambda larger than 1, or smaller than 0
// For values < 0, the particle travels away from the cell
// and we don't move the particle, only change cell.
// For values larger than 1, we move the particle to endPosition only.
if (lambdaMin > 0.0)
// Set the faceHitTetIs to those for the current tet in case a
// wall interaction is required with the cell face
faceHitTetIs = tetIndices
(
cellI_,
tetFaceI_,
tetBasePtI,
fPtAI,
fPtBI,
tetPtI_
);
}
else if (triI > 0)
{
// A tri was found to be crossed before a wall face was hit (if any)
faceI_ = -1;
}
// The particle can be 'outside' the tet. This will yield a
// lambda larger than 1, or smaller than 0. For values < 0,
// the particle travels away from the tet and we don't move
// the particle, only change tet/cell. For values larger than
// 1, we move the particle to endPosition before the tet/cell
// change.
if (lambdaMin > SMALL)
{
if (lambdaMin <= 1.0)
{
trackFraction = lambdaMin;
position_ += trackFraction*(endPosition - position_);
trackFraction += lambdaMin*(1 - trackFraction);
position_ += lambdaMin*(endPosition - position_);
}
else
{
trackFraction = 1.0;
position_ = endPosition;
}
}
else if (static_cast<ParticleType&>(*this).softImpact())
{
// Soft-sphere particles can travel outside the domain
// but we don't use lambda since this the particle
// is going away from face
trackFraction = 1.0;
position_ = endPosition;
}
// change cell
if (internalFace) // Internal face
{
if (celli_ == mesh.faceOwner()[facei_])
{
celli_ = mesh.faceNeighbour()[facei_];
}
else if (celli_ == mesh.faceNeighbour()[facei_])
{
celli_ = mesh.faceOwner()[facei_];
}
else
{
FatalErrorIn
(
"Particle::trackToFace(const vector&, TrackData&)"
)<< "addressing failure" << nl
<< abort(FatalError);
position_ = endPosition;
}
}
else
{
ParticleType& p = static_cast<ParticleType&>(*this);
// Set lambdaMin to zero to force a towards-tet-centre
// correction.
lambdaMin = 0.0;
}
// Soft-sphere algorithm ignores the boundary
if (p.softImpact())
} while (faceI_ < 0);
if (cloud_.internalFace(faceI_))
{
// Change tet ownership because a tri face has been crossed,
// in general this is:
// tetNeighbour(triI);
// but triI must be 0;
// No modifications are required for triI = 0, no call required to
// tetNeighbour(0);
if (cellI_ == mesh.faceOwner()[faceI_])
{
cellI_ = mesh.faceNeighbour()[faceI_];
}
else if (cellI_ == mesh.faceNeighbour()[faceI_])
{
cellI_ = mesh.faceOwner()[faceI_];
}
else
{
FatalErrorIn
(
"Particle::trackToFace(const vector&, TrackData&)"
) << "addressing failure" << nl
<< abort(FatalError);
}
}
else
{
ParticleType& p = static_cast<ParticleType&>(*this);
label origFaceI = faceI_;
label patchI = patch(faceI_);
// No action taken for tetPtI_ for tetFaceI_ here, handled by
// patch interaction call or later during processor transfer.
if
(
!p.hitPatch
(
mesh.boundaryMesh()[patchI],
td,
patchI,
trackFraction,
faceHitTetIs
)
)
{
// Did patch interaction model switch patches?
if (faceI_ != origFaceI)
{
trackFraction = 1.0;
position_ = endPosition;
patchI = patch(faceI_);
}
label origFacei = facei_;
label patchi = patch(facei_);
const polyPatch& patch = mesh.boundaryMesh()[patchI];
if (!p.hitPatch(mesh.boundaryMesh()[patchi], td, patchi))
if (isA<wedgePolyPatch>(patch))
{
// Did patch interaction model switch patches?
if (facei_ != origFacei)
{
patchi = patch(facei_);
}
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (isA<wedgePolyPatch>(patch))
{
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td
);
}
else
{
p.hitPatch(patch, td);
}
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
);
}
else
{
p.hitPatch(patch, td);
}
}
}
// If the trackFraction = 0 something went wrong.
// Either the particle is flipping back and forth across a face perhaps
// due to velocity interpolation errors or it is in a "hole" in the mesh
// caused by face warpage.
// In both cases resolve the positional ambiguity by moving the particle
// slightly towards the cell-centre.
if (trackFraction < SMALL)
if (lambdaMin < SMALL)
{
position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
// Apply tracking correction towards tet centre.
// Generate current tet to find centre to apply correction.
tetPointRef tet = currentTet();
position_ +=
Cloud<ParticleType>::trackingCorrectionTol
*(tet.centre() - position_);
if
(
cloud_.hasWallImpactDistance()
&& !cloud_.internalFace(faceHitTetIs.face())
&& cloud_.cellHasWallFaces()[faceHitTetIs.cell()]
)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label fI = faceHitTetIs.face();
label patchI = patches.patchID()[fI - mesh.nInternalFaces()];
if (isA<wallPolyPatch>(patches[patchI]))
{
// In the case of collision with a wall where there is
// a non-zero wallImpactDistance, it is possible for
// there to be a tracking correction required to bring
// the particle into the domain, but the position of
// the particle is further from the wall than the tet
// centre, in which case the normal correction can be
// counter-productive, i.e. pushes the particle
// further out of the domain. In this case it is the
// position that hit the wall that is in need of a
// rescue correction.
triPointRef wallTri = faceHitTetIs.faceTri(mesh);
tetPointRef wallTet = faceHitTetIs.tet(mesh);
vector nHat = wallTri.normal();
nHat /= mag(nHat);
const ParticleType& p = static_cast<const ParticleType&>(*this);
scalar r = p.wallImpactDistance(nHat);
// Removing (approximately) the wallTri normal
// component of the existing correction, to avoid the
// situation where the existing correction in the wall
// normal direction is larger towards the wall than
// the new correction is away from it.
position_ +=
Cloud<ParticleType>::trackingCorrectionTol
*(
(wallTet.centre() - (position_ + r*nHat))
- (nHat & (tet.centre() - position_))*nHat
);
}
}
cloud_.trackingRescue();
}
return trackFraction;
}
template<class ParticleType>
Foam::scalar Foam::Particle<ParticleType>::trackToFace
(
@ -413,6 +657,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
return trackToFace(endPosition, dummyTd);
}
template<class ParticleType>
void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
{
@ -436,7 +681,9 @@ bool Foam::Particle<ParticleType>::hitPatch
(
const polyPatch&,
TrackData&,
const label
const label,
const scalar,
const tetIndices&
)
{
return false;
@ -451,7 +698,17 @@ void Foam::Particle<ParticleType>::hitWedgePatch
TrackData&
)
{
vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
FatalErrorIn
(
"void Foam::Particle<ParticleType>::hitWedgePatch"
"("
"const wedgePolyPatch& wpp, "
"TrackData&"
")"
) << "Hitting a wedge patch should not be possible."
<< abort(FatalError);
vector nf = normal();
nf /= mag(nf);
static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
@ -466,7 +723,7 @@ void Foam::Particle<ParticleType>::hitSymmetryPatch
TrackData&
)
{
vector nf = spp.faceAreas()[spp.whichFace(facei_)];
vector nf = normal();
nf /= mag(nf);
static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
@ -481,11 +738,16 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
TrackData&
)
{
// label patchFacei_ = cpp.whichFace(facei_);
// label patchFaceI_ = cpp.whichFace(faceI_);
facei_ = cpp.transformGlobalFace(facei_);
faceI_ = cpp.transformGlobalFace(faceI_);
celli_ = cloud_.polyMesh_.faceOwner()[facei_];
cellI_ = cloud_.polyMesh_.faceOwner()[faceI_];
tetFaceI_ = faceI_;
// See note in correctAfterParallelTransfer for tetPtI_ addressing.
tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
// Now the particle is on the receiving side
@ -522,7 +784,8 @@ template<class TrackData>
void Foam::Particle<ParticleType>::hitWallPatch
(
const wallPolyPatch& spp,
TrackData&
TrackData&,
const tetIndices&
)
{}
@ -549,7 +812,7 @@ bool Foam::operator==
return
(
pA.origProc() == pB.origProc()
&& pA.origId() == pB.origId()
&& pA.origId() == pB.origId()
);
}

180
src/lagrangian/basic/Particle/Particle.H
View File

@ -38,6 +38,10 @@ Description
#include "faceList.H"
#include "typeInfo.H"
#include "OFstream.H"
#include "tetPointRef.H"
#include "FixedList.H"
#include "polyMeshTetDecomposition.H"
#include "wallPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -125,14 +129,23 @@ protected:
vector position_;
//- Index of the cell it is in
label celli_;
label cellI_;
//- Face index if the particle is on a face otherwise -1
label facei_;
label faceI_;
//- Fraction of time-step completed
scalar stepFraction_;
//- Index of the face that owns the decomposed tet that the
// particle is in
label tetFaceI_;
//- Index of the point on the face that defines the decomposed
// tet that the particle is in. Relative to the face base
// point.
label tetPtI_;
//- Originating processor id
label origProc_;
@ -142,54 +155,82 @@ protected:
// Private Member Functions
//- Return the 'lambda' value for the position, p, on the face,
// where, p = from + lamda*(to - from)
// for non-static meshes
inline scalar lambda
//- Find the tet tri faces between position and tet centre
inline void findTris
(
const vector& position,
DynamicList<label>& faceList,
const tetPointRef& tet,
const FixedList<vector, 4>& tetAreas,
const FixedList<label, 4>& tetPlaneBasePtIs
) const;
//- Find the lambda value for the line to-from across the
// given tri face, where p = from + lambda*(to - from)
inline scalar tetLambda
(
const vector& from,
const vector& to,
const label facei,
const scalar stepFraction
const label triI,
const vector& tetArea,
const label tetPlaneBasePtI,
const label cellI,
const label tetFaceI,
const label tetPtI
) const;
//- Return the 'lambda' value for the position, p, on the face,
// where, p = from + lamda*(to - from)
// for static meshes
inline scalar lambda
//- Find the lambda value for a moving tri face
inline scalar movingTetLambda
(
const vector& from,
const vector& to,
const label triI,
const vector& tetArea,
const label tetPlaneBasePtI,
const label cellI,
const label tetFaceI,
const label tetPtI
) const;
//- Modify the tet owner data by crossing triI
inline void tetNeighbour(label triI);
//- Cross the from the given face across the given edge of the
// given cell to find the resulting face and tetPtI
inline void crossEdgeConnectedFace
(
const label& cellI,
label& tetFaceI,
label& tetPtI,
const edge& e
);
//- Hit wall faces in the current cell if the
//- wallImpactDistance is non-zero. They may not be in
//- different tets to the current.
inline void hitWallFaces
(
const vector& from,
const vector& to,
const label facei
) const;
//- Find the faces between position and cell centre
void findFaces
(
const vector& position,
DynamicList<label>& faceList
) const;
//- Find the faces between position and cell centre
void findFaces
(
const vector& position,
const label celli,
const scalar stepFraction,
DynamicList<label>& faceList
) const;
scalar& lambdaMin,
tetIndices& closestTetIs
);
// Patch interactions
//- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions
//- Overridable function to handle the particle hitting a
// patch. Executed before other patch-hitting functions.
// trackFraction is passed in to allow mesh motion to
// interpolate in time to the correct face state.
template<class TrackData>
bool hitPatch
(
const polyPatch&,
TrackData& td,
const label patchI
const label patchI,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a wedgePatch
@ -231,7 +272,8 @@ protected:
void hitWallPatch
(
const wallPolyPatch&,
TrackData& td
TrackData& td,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a
@ -264,12 +306,12 @@ protected:
//- Convert global addressing to the processor patch
// local equivalents
template<class TrackData>
void prepareForParallelTransfer(const label patchi, TrackData& td);
void prepareForParallelTransfer(const label patchI, TrackData& td);
//- Convert processor patch addressing to the global equivalents
// and set the celli to the face-neighbour
// and set the cellI to the face-neighbour
template<class TrackData>
void correctAfterParallelTransfer(const label patchi, TrackData& td);
void correctAfterParallelTransfer(const label patchI, TrackData& td);
public:
@ -293,7 +335,19 @@ public:
(
const Cloud<ParticleType>&,
const vector& position,
const label celli
const label cellI,
const label tetFaceI,
const label tetPtI
);
//- Construct from components, tetFaceI_ and tetPtI_ are not
// supplied so they will be deduced by a search
Particle
(
const Cloud<ParticleType>&,
const vector& position,
const label cellI,
bool doCellFacePt = true
);
//- Construct from Istream
@ -355,17 +409,6 @@ public:
// Access
//- Return true if particle is in cell
inline bool inCell() const;
//- Return true if position is in cell i
inline bool inCell
(
const vector& position,
const label celli,
const scalar stepFraction
) const;
//- Return current particle position
inline const vector& position() const;
@ -378,6 +421,35 @@ public:
//- Return current cell particle is in
inline label cell() const;
//- Return current tet face particle is in
inline label& tetFace();
//- Return current tet face particle is in
inline label tetFace() const;
//- Return current tet face particle is in
inline label& tetPt();
//- Return current tet face particle is in
inline label tetPt() const;
//- Return the indices of the current tet that the
// particle occupies.
inline tetIndices currentTetIndices() const;
//- Return the geometry of the current tet that the
// particle occupies.
inline tetPointRef currentTet() const;
//- Return the normal of the tri on tetFaceI_ for the
// current tet.
inline vector normal() const;
//- Return the normal of the tri on tetFaceI_ for the
// current tet at the start of the timestep, i.e. based
// on oldPoints
inline vector oldNormal() const;
//- Return current face particle is on otherwise -1
inline label& face();
@ -396,17 +468,21 @@ public:
// Check
//- Check the stored cell value (setting if necessary) and
// initialise the tetFace and tetPt values
inline void initCellFacePt();
//- Is the particle on the boundary/(or outside the domain)?
inline bool onBoundary() const;
//- Which patch is particle on
inline label patch(const label facei) const;
inline label patch(const label faceI) const;
//- Which face of this patch is this particle on
inline label patchFace
(
const label patchi,
const label facei
const label patchI,
const label faceI
) const;
//- The nearest distance to a wall that

1124
src/lagrangian/basic/Particle/ParticleI.H
View File

File diff suppressed because it is too large Load Diff

56
src/lagrangian/basic/Particle/ParticleIO.C
View File

@ -31,7 +31,8 @@ License
template<class ParticleType>
Foam::string Foam::Particle<ParticleType>::propHeader =
"(Px Py Pz) cellI origProc origId";
"(Px Py Pz) cellI tetFaceI tetPtI origProc origId";
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -44,20 +45,24 @@ Foam::Particle<ParticleType>::Particle
)
:
cloud_(cloud),
facei_(-1),
position_(),
cellI_(-1),
faceI_(-1),
stepFraction_(0.0),
tetFaceI_(-1),
tetPtI_(-1),
origProc_(Pstream::myProcNo()),
origId_(-1)
{
// readFields : read additional data. Should be consistent with writeFields.
if (is.format() == IOstream::ASCII)
{
is >> position_ >> celli_;
is >> position_ >> cellI_;
if (readFields)
{
is >> origProc_ >> origId_;
is >> tetFaceI_ >> tetPtI_ >> origProc_ >> origId_;
}
}
else
@ -69,9 +74,11 @@ Foam::Particle<ParticleType>::Particle
(
reinterpret_cast<char*>(&position_),
sizeof(position_)
+ sizeof(celli_)
+ sizeof(facei_)
+ sizeof(cellI_)
+ sizeof(faceI_)
+ sizeof(stepFraction_)
+ sizeof(tetFaceI_)
+ sizeof(tetPtI_)
+ sizeof(origProc_)
+ sizeof(origId_)
);
@ -82,18 +89,13 @@ Foam::Particle<ParticleType>::Particle
(
reinterpret_cast<char*>(&position_),
sizeof(position_)
+ sizeof(celli_)
+ sizeof(facei_)
+ sizeof(cellI_)
+ sizeof(faceI_)
+ sizeof(stepFraction_)
);
}
}
if (celli_ == -1)
{
celli_ = cloud_.pMesh().findCell(position_);
}
// Check state of Istream
is.check("Particle<ParticleType>::Particle(Istream&)");
}
@ -176,29 +178,33 @@ void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
if (writeFields)
{
// Write the additional entries
os << position_
<< token::SPACE << celli_
<< token::SPACE << origProc_
<< token::SPACE << origId_;
os << position_
<< token::SPACE << cellI_
<< token::SPACE << tetFaceI_
<< token::SPACE << tetPtI_
<< token::SPACE << origProc_
<< token::SPACE << origId_;
}
else
{
os << position_
<< token::SPACE << celli_;
os << position_
<< token::SPACE << cellI_;
}
}
else
{
// In binary write both celli_ and facei_, needed for parallel transfer
// In binary write both cellI_ and faceI_, needed for parallel transfer
if (writeFields)
{
os.write
(
reinterpret_cast<const char*>(&position_),
sizeof(position_)
+ sizeof(celli_)
+ sizeof(facei_)
+ sizeof(cellI_)
+ sizeof(faceI_)
+ sizeof(stepFraction_)
+ sizeof(tetFaceI_)
+ sizeof(tetPtI_)
+ sizeof(origProc_)
+ sizeof(origId_)
);
@ -209,8 +215,8 @@ void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
(
reinterpret_cast<const char*>(&position_),
sizeof(position_)
+ sizeof(celli_)
+ sizeof(facei_)
+ sizeof(cellI_)
+ sizeof(faceI_)
+ sizeof(stepFraction_)
);
}

19
src/lagrangian/basic/indexedParticle/indexedParticle.H
View File

@ -67,11 +67,26 @@ public:
(
const Cloud<indexedParticle>& c,
const vector& position,
const label celli,
const label cellI,
const label tetFaceI,
const label tetPtI,
const label index = 0
)
:
Particle<indexedParticle>(c, position, celli),
Particle<indexedParticle>(c, position, cellI, tetFaceI, tetPtI),
index_(index)
{}
//- Construct from components, with searching for tetFace and tetPt
indexedParticle
(
const Cloud<indexedParticle>& c,
const vector& position,
const label cellI,
const label index = 0
)
:
Particle<indexedParticle>(c, position, cellI),
index_(index)
{}

19
src/lagrangian/basic/passiveParticle/passiveParticle.H
View File

@ -60,10 +60,25 @@ public:
(
const Cloud<passiveParticle>& c,
const vector& position,
const label celli
const label cellI,
const label tetFaceI,
const label tetPtI
)
:
Particle<passiveParticle>(c, position, celli)
Particle<passiveParticle>(c, position, cellI, tetFaceI, tetPtI)
{}
//- Construct from components, with searching for tetFace and
// tetPt unless disabled by doCellFacePt = false.
passiveParticle
(
const Cloud<passiveParticle>& c,
const vector& position,
const label cellI,
bool doCellFacePt = true
)
:
Particle<passiveParticle>(c, position, cellI, doCellFacePt)
{}
//- Construct from Istream

19
src/lagrangian/coalCombustion/coalParcel/coalParcel.C
View File

@ -31,10 +31,19 @@ Foam::coalParcel::coalParcel
(
ReactingMultiphaseCloud<coalParcel>& owner,
const vector& position,
const label cellI
const label cellI,
const label tetFaceI,
const label tetPtI
)
:
ReactingMultiphaseParcel<coalParcel>(owner, position, cellI)
ReactingMultiphaseParcel<coalParcel>
(
owner,
position,
cellI,
tetFaceI,
tetPtI
)
{}
@ -43,9 +52,12 @@ Foam::coalParcel::coalParcel
ReactingMultiphaseCloud<coalParcel>& owner,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const label typeId,
const scalar nParticle0,
const scalar d0,
const scalar dTarget0,
const vector& U0,
const vector& f0,
const vector& pi0,
@ -62,9 +74,12 @@ Foam::coalParcel::coalParcel
owner,
position,
cellI,
tetFaceI,
tetPtI,
typeId,
nParticle0,
d0,
dTarget0,
U0,
f0,
pi0,

7
src/lagrangian/coalCombustion/coalParcel/coalParcel.H
View File

@ -63,7 +63,9 @@ public:
(
ReactingMultiphaseCloud<coalParcel>& owner,
const vector& position,
const label cellI
const label cellI,
const label tetFaceI,
const label tetPtI
);
//- Construct from components
@ -72,9 +74,12 @@ public:
ReactingMultiphaseCloud<coalParcel>& owner,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const label typeId,
const scalar nParticle0,
const scalar d0,
const scalar dTarget0,
const vector& U0,
const vector& f0,
const vector& pi0,

2
src/lagrangian/dieselSpray/Make/options
View File

@ -1,5 +1,6 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
@ -16,6 +17,7 @@ EXE_INC = \
LIB_LIBS = \
-llagrangian \
-lmeshTools \
-lfiniteVolume \
-lcompressibleRASModels \
-lcompressibleLESModels \

14
src/lagrangian/dieselSpray/parcel/boundaryTreatment.H
View File

@ -15,12 +15,11 @@ if (isA<wallPolyPatch>(pbMesh[patch(face())]))
else if (isA<wedgePolyPatch>(pbMesh[patch(face())]))
{
// check if parcel is trying to move out of the domain
label patchi = patch(face());
label patchFacei = patchFace(patchi, face());
const polyPatch& patch = mesh.boundaryMesh()[patchi];
vector nf = patch.faceAreas()[patchFacei];
vector nf = normal();
scalar Un = U() & nf;
if (Un > 0)
{
scalar Un2 = U() & n();
@ -30,12 +29,11 @@ else if (isA<wedgePolyPatch>(pbMesh[patch(face())]))
else if (isA<symmetryPolyPatch>(pbMesh[patch(face())]))
{
// check if parcel is trying to move out of the domain
label patchi = patch(face());
label patchFacei = patchFace(patchi, face());
const polyPatch& patch = mesh.boundaryMesh()[patchi];
vector nf = patch.faceAreas()[patchFacei];
vector nf = normal();
scalar Un = U() & nf;
if (Un > 0)
{
if (sDB.twoD())

73
src/lagrangian/dieselSpray/parcel/parcel.C
View File

@ -50,6 +50,8 @@ Foam::parcel::parcel
const Cloud<parcel>& cloud,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const vector& n,
const scalar d,
const scalar T,
@ -67,7 +69,7 @@ Foam::parcel::parcel
const List<word>& liquidNames
)
:
Particle<parcel>(cloud, position, cellI),
Particle<parcel>(cloud, position, cellI, tetFaceI, tetPtI),
liquidComponents_
(
liquidNames
@ -103,12 +105,13 @@ bool Foam::parcel::move(spray& sDB)
label Nf = fuels.components().size();
label Ns = sDB.composition().Y().size();
tetIndices tetIs = this->currentTetIndices();
// Calculate the interpolated gas properties at the position of the parcel
vector Up = sDB.UInterpolator().interpolate(position(), cell())
+ Uturb();
scalar rhog = sDB.rhoInterpolator().interpolate(position(), cell());
scalar pg = sDB.pInterpolator().interpolate(position(), cell());
scalar Tg = sDB.TInterpolator().interpolate(position(), cell());
vector Up = sDB.UInterpolator().interpolate(position(), tetIs) + Uturb();
scalar rhog = sDB.rhoInterpolator().interpolate(position(), tetIs);
scalar pg = sDB.pInterpolator().interpolate(position(), tetIs);
scalar Tg = sDB.TInterpolator().interpolate(position(), tetIs);
scalarField Yfg(Nf, 0.0);
@ -119,8 +122,8 @@ bool Foam::parcel::move(spray& sDB)
if (sDB.isLiquidFuel()[i])
{
label j = sDB.gasToLiquidIndex()[i];
scalar Yicelli = Yi[cell()];
Yfg[j] = Yicelli;
scalar YicellI = Yi[cell()];
Yfg[j] = YicellI;
}
cpMixture += Yi[cell()]*sDB.gasProperties()[i].Cp(Tg);
}
@ -199,8 +202,8 @@ bool Foam::parcel::move(spray& sDB)
// remember which cell the parcel is in
// since this will change if a face is hit
label celli = cell();
scalar p = sDB.p()[celli];
label cellI = cell();
scalar p = sDB.p()[cellI];
// track parcel to face, or end of trajectory
if (keepParcel)
@ -259,7 +262,7 @@ bool Foam::parcel::move(spray& sDB)
(
dt,
sDB,
celli,
cellI,
face()
);
@ -283,11 +286,11 @@ bool Foam::parcel::move(spray& sDB)
// Update the Spray Source Terms
forAll(nMass, i)
{
sDB.srhos()[i][celli] += oMass[i] - nMass[i];
sDB.srhos()[i][cellI] += oMass[i] - nMass[i];
}
sDB.sms()[celli] += oMom - nMom;
sDB.sms()[cellI] += oMom - nMom;
sDB.shs()[celli] += oTotMass*(oH + oPE) - m()*(nH + nPE);
sDB.shs()[cellI] += oTotMass*(oH + oPE) - m()*(nH + nPE);
// Remove evaporated mass from stripped mass
ms() -= ms()*(oTotMass-m())/oTotMass;
@ -300,11 +303,11 @@ bool Foam::parcel::move(spray& sDB)
// ... and add the removed 'stuff' to the gas
forAll(nMass, i)
{
sDB.srhos()[i][celli] += nMass[i];
sDB.srhos()[i][cellI] += nMass[i];
}
sDB.sms()[celli] += nMom;
sDB.shs()[celli] += m()*(nH + nPE);
sDB.sms()[cellI] += nMom;
sDB.shs()[cellI] += m()*(nH + nPE);
}
if (onBoundary() && keepParcel)
@ -327,8 +330,8 @@ void Foam::parcel::updateParcelProperties
(
const scalar dt,
spray& sDB,
const label celli,
const label facei
const label cellI,
const label faceI
)
{
const liquidMixture& fuels = sDB.fuels();
@ -340,34 +343,34 @@ void Foam::parcel::updateParcelProperties
scalar W = 0.0;
for (label i=0; i<Ns; i++)
{
W += sDB.composition().Y()[i][celli]/sDB.gasProperties()[i].W();
W += sDB.composition().Y()[i][cellI]/sDB.gasProperties()[i].W();
}
W = 1.0/W;
// Calculate the interpolated gas properties at the position of the parcel
vector Up = sDB.UInterpolator().interpolate(position(), celli, facei)
vector Up = sDB.UInterpolator().interpolate(position(), cellI, faceI)
+ Uturb();
scalar rhog = sDB.rhoInterpolator().interpolate(position(), celli, facei);
scalar pg = sDB.pInterpolator().interpolate(position(), celli, facei);
scalar Tg0 = sDB.TInterpolator().interpolate(position(), celli, facei);
scalar rhog = sDB.rhoInterpolator().interpolate(position(), cellI, faceI);
scalar pg = sDB.pInterpolator().interpolate(position(), cellI, faceI);
scalar Tg0 = sDB.TInterpolator().interpolate(position(), cellI, faceI);
// correct the gaseous temperature for evaporated fuel
scalar cpMix = 0.0;
for (label i=0; i<Ns; i++)
{
cpMix += sDB.composition().Y()[i][celli]
cpMix += sDB.composition().Y()[i][cellI]
*sDB.gasProperties()[i].Cp(T());
}
scalar cellV = sDB.mesh().V()[celli];
scalar rho = sDB.rho()[celli];
scalar cellV = sDB.mesh().V()[cellI];
scalar rho = sDB.rho()[cellI];
scalar cellMass = rho*cellV;
scalar dh = sDB.shs()[celli];
scalar dh = sDB.shs()[cellI];
scalarField addedMass(Nf, 0.0);
forAll(addedMass, i)
{
addedMass[i] += sDB.srhos()[i][celli]*cellV;
addedMass[i] += sDB.srhos()[i][cellI]*cellV;
}
scalar Tg = Tg0 + dh/(cpMix*cellMass);
@ -378,7 +381,7 @@ void Foam::parcel::updateParcelProperties
{
label j = sDB.liquidToGasIndex()[i];
const volScalarField& Yj = sDB.composition().Y()[j];
scalar Yfg0 = Yj[celli];
scalar Yfg0 = Yj[cellI];
Yfg[i] = (Yfg0*cellMass + addedMass[i])/(addedMass[i] + cellMass);
}
@ -389,7 +392,7 @@ void Foam::parcel::updateParcelProperties
setRelaxationTimes
(
celli,
cellI,
tauMomentum,
tauEvaporation,
tauHeatTransfer,
@ -492,14 +495,14 @@ void Foam::parcel::updateParcelProperties
{
label j = sDB.liquidToGasIndex()[i];
const volScalarField& Yj = sDB.composition().Y()[j];
scalar Yfg0 = Yj[celli];
scalar Yfg0 = Yj[cellI];
Yfg[i] = (Yfg0*cellMass + addedMass[i] + dm)
/(addedMass[i] + cellMass + dm);
}
setRelaxationTimes
(
celli,
cellI,
tauMomentum,
tauEvaporation,
tauHeatTransfer,
@ -530,7 +533,7 @@ void Foam::parcel::updateParcelProperties
}
else
{
scalar Y = sDB.composition().Y()[i][celli];
scalar Y = sDB.composition().Y()[i][cellI];
cpMix += Y*sDB.gasProperties()[i].Cp(Taverage);
}
}
@ -579,7 +582,7 @@ void Foam::parcel::updateParcelProperties
{
// can not go above boiling temperature
scalar Terr = 1.0e-3;
label n=0;
label n = 0;
scalar dT = 1.0;
scalar pOld = pAtSurface;
while (dT > Terr)

10
src/lagrangian/dieselSpray/parcel/parcel.H
View File

@ -113,7 +113,7 @@ class parcel
//- Set the relaxation times
void setRelaxationTimes
(
label celli,
label cellI,
scalar& tauMomentum,
scalarField& tauEvaporation,
scalar& tauHeatTransfer,
@ -133,8 +133,8 @@ class parcel
(
const scalar dt,
spray& sprayData,
const label celli,
const label facei
const label cellI,
const label faceI
);
@ -150,7 +150,9 @@ public:
(
const Cloud<parcel>& cloud,
const vector& position,
const label celli,
const label cellI,
const label tetFaceI,
const label tetPtI,
const vector& n,
const scalar d,
const scalar T,

10
src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C
View File

@ -36,7 +36,7 @@ License
void Foam::parcel::setRelaxationTimes
(
label celli,
label cellI,
scalar& tauMomentum,
scalarField& tauEvaporation,
scalar& tauHeatTransfer,
@ -70,7 +70,7 @@ void Foam::parcel::setRelaxationTimes
for (label i=0; i<Ns; i++)
{
scalar Y = sDB.composition().Y()[i][celli];
scalar Y = sDB.composition().Y()[i][cellI];
W += Y/sDB.gasProperties()[i].W();
// Using mass-fractions to average...
kMixture += Y*sDB.gasProperties()[i].kappa(Tf);
@ -87,7 +87,7 @@ void Foam::parcel::setRelaxationTimes
for (label i=0; i<Nf; i++)
{
label j = sDB.liquidToGasIndex()[i];
scalar Y = sDB.composition().Y()[j][celli];
scalar Y = sDB.composition().Y()[j][cellI];
scalar Wi = sDB.gasProperties()[j].W();
Yf[i] = Y;
Xf[i] = Y*W/Wi;
@ -264,10 +264,10 @@ void Foam::parcel::setRelaxationTimes
forAll(sDB.gasProperties(), k)
{
vapourSurfaceEnthalpy +=
sDB.composition().Y()[k][celli]
sDB.composition().Y()[k][cellI]
*sDB.gasProperties()[k].H(tBoilingSurface);
vapourFarEnthalpy +=
sDB.composition().Y()[k][celli]
sDB.composition().Y()[k][cellI]
*sDB.gasProperties()[k].H(temperature);
}

9
src/lagrangian/dieselSpray/spray/findInjectorCell.H
View File

@ -22,7 +22,14 @@ reduce(foundCell, orOp<bool>());
if (!foundCell)
{
injectionPosition = it->position(n);
injectorCell = mesh_.findCell(injectionPosition);
findCellFacePt
(
injectionPosition,
injectorCell,
injectorTetFaceI,
injectorTetPtI
);
if (injectorCell >= 0)
{

Some files were not shown because too many files have changed in this diff Show More