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

Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry
2011-12-28 17:46:59 +00:00
parent a78a46e7c2 b51dbec8d0
commit b79aa9d27f
46 changed files with 1526 additions and 410 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -433,7 +433,7 @@ mtype {space}"MTYPE:"{space}
}
else
{
curGroupID = readLabel(groupStream);;
curGroupID = readLabel(groupStream);
}
BEGIN(cellStreams);

2
applications/utilities/mesh/manipulation/renumberMesh/Make/options
View File

@ -11,4 +11,4 @@ EXE_LIBS = \
-lfiniteVolume \
-lgenericPatchFields \
-lrenumberMethods \
-ldecompositionMethods
-ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp

127
applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
View File

@ -29,7 +29,7 @@ Description
renumbering all fields from all the time directories.
By default uses bandCompression (CuthillMcKee) but will
read system/renumberMeshDict if present and use the method from there.
read system/renumberMeshDict if -dict option is present
\*---------------------------------------------------------------------------*/
@ -45,6 +45,7 @@ Description
#include "renumberMethod.H"
#include "zeroGradientFvPatchFields.H"
#include "CuthillMcKeeRenumber.H"
#include "fvMeshSubset.H"
using namespace Foam;
@ -455,6 +456,69 @@ autoPtr<mapPolyMesh> reorderMesh
}
// Return new to old cell numbering
labelList regionRenumber
(
const renumberMethod& method,
const fvMesh& mesh,
const labelList& cellToRegion
)
{
Info<< "Determining cell order:" << endl;
labelList cellOrder(cellToRegion.size());
label nRegions = max(cellToRegion)+1;
labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
label cellI = 0;
forAll(regionToCells, regionI)
{
Info<< " region " << regionI << " starts at " << cellI << endl;
// Make sure no parallel comms
bool oldParRun = UPstream::parRun();
UPstream::parRun() = false;
// Per region do a reordering.
fvMeshSubset subsetter(mesh);
subsetter.setLargeCellSubset(cellToRegion, regionI);
const fvMesh& subMesh = subsetter.subMesh();
labelList subReverseCellOrder = method.renumber
(
subMesh,
subMesh.cellCentres()
);
labelList subCellOrder
(
invert
(
subMesh.nCells(),
subReverseCellOrder
)
);
// Restore state
UPstream::parRun() = oldParRun;
const labelList& cellMap = subsetter.cellMap();
forAll(subCellOrder, i)
{
cellOrder[cellI++] = cellMap[subCellOrder[i]];
}
}
Info<< endl;
return cellOrder;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@ -505,23 +569,28 @@ int main(int argc, char *argv[])
label blockSize = 0;
// Construct renumberMethod
autoPtr<IOdictionary> renumberDictPtr;
autoPtr<renumberMethod> renumberPtr;
if (readDict)
{
Info<< "Renumber according to renumberMeshDict." << nl << endl;
IOdictionary renumberDict
renumberDictPtr.reset
(
IOobject
new IOdictionary
(
"renumberMeshDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
IOobject
(
"renumberMeshDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
)
);
const IOdictionary renumberDict = renumberDictPtr();
renumberPtr = renumberMethod::New(renumberDict);
@ -562,7 +631,7 @@ int main(int argc, char *argv[])
<< endl;
}
writeMaps = readLabel(renumberDict.lookup("writeMaps"));
renumberDict.lookup("writeMaps") >> writeMaps;
if (writeMaps)
{
Info<< "Writing renumber maps (new to old) to polyMesh." << nl
@ -683,22 +752,15 @@ int main(int argc, char *argv[])
// fields is done correctly!
label nBlocks = mesh.nCells() / blockSize;
Info<< "nBlocks = " << nBlocks << endl;
Info<< "nBlocks = " << nBlocks << endl;
// Read decomposePar dictionary
IOdictionary decomposeDict
(
IOobject
(
"decomposeParDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// Read decompositionMethod dictionary
dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
decomposeDict.set("numberOfSubdomains", nBlocks);
bool oldParRun = UPstream::parRun();
UPstream::parRun() = false;
autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
(
decomposeDict
@ -713,6 +775,9 @@ int main(int argc, char *argv[])
)
);
// Restore state
UPstream::parRun() = oldParRun;
// For debugging: write out region
createScalarField
(
@ -726,23 +791,7 @@ int main(int argc, char *argv[])
<< nl << endl;
// Find point per region
pointField regionPoints(nBlocks, vector::zero);
forAll(cellToRegion, cellI)
{
regionPoints[cellToRegion[cellI]] = mesh.cellCentres()[cellI];
}
// Use block based renumbering.
// Detemines old to new cell ordering
labelList reverseCellOrder = renumberPtr().renumber
(
mesh,
cellToRegion,
regionPoints
);
cellOrder = invert(mesh.nCells(), reverseCellOrder);
cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
// Determine new to old face order with new cell numbering
faceOrder = getRegionFaceOrder

35
applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
View File

@ -22,10 +22,14 @@ writeMaps true;
// e.g. nonBlockingGaussSeidel.
sortCoupledFaceCells false;
// Optional entry: renumber on a block-by-block basis. This can be used on
// large cases to keep the blocks fitting in cache with all the the cache
// missed bunched at the end.
//blockSize 0;
// Optional entry: renumber on a block-by-block basis. It uses a
// blockCoeffs dictionary to construct a decompositionMethod to do
// a block subdivision) and then applies the renumberMethod to each
// block in turn. This can be used in large cases to keep the blocks
// fitting in cache with all the the cache misses bunched at the end.
// This number is the approximate size of the blocks - this gets converted
// to a number of blocks that is the input to the decomposition method.
//blockSize 1000;
// Optional entry: sort points into internal and boundary points
//orderPoints false;
@ -37,11 +41,11 @@ method CuthillMcKee;
//method random;
//method spring;
CuthillMcKeeCoeffs
{
// Reverse CuthillMcKee (RCM) or plain
reverse true;
}
//CuthillMcKeeCoeffs
//{
// // Reverse CuthillMcKee (RCM) or plain
// reverse true;
//}
manualCoeffs
@ -65,4 +69,17 @@ springCoeffs
}
blockCoeffs
{
method scotch;
//method hierarchical;
//hierarchicalCoeffs
//{
// n (1 2 1);
// delta 0.001;
// order xyz;
//}
}
// ************************************************************************* //

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

@ -528,7 +528,7 @@ int main(int argc, char *argv[])
<< "Cell number should be between 0 and "
<< mesh.nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< mesh.findCell(iter().position())
<< mesh.findCell(iter().position())
<< exit(FatalError);
}
@ -789,6 +789,14 @@ int main(int argc, char *argv[])
// Point fields
if
(
pointScalarFields.size()
|| pointVectorFields.size()
|| pointSphericalTensorFields.size()
|| pointSymmTensorFields.size()
|| pointTensorFields.size()
)
{
labelIOList pointProcAddressing
(

7
applications/utilities/parallelProcessing/redistributePar/redistributePar.C
View File

@ -795,9 +795,9 @@ int main(int argc, char *argv[])
);
fvMesh& mesh = meshPtr();
//Pout<< "Read mesh:" << endl;
//printMeshData(Pout, mesh);
//Pout<< endl;
// Print some statistics
Info<< "Before distribution:" << endl;
printMeshData(mesh);
@ -1022,7 +1022,6 @@ int main(int argc, char *argv[])
//map().distributeFaceData(faceCc);
// Print a bit
// Print some statistics
Info<< "After distribution:" << endl;
printMeshData(mesh);

4
applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
View File

@ -1,6 +1,6 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
#set -x
if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
then
@ -15,6 +15,8 @@ then
wmake libso vtkPV3Readers
PV3blockMeshReader/Allwmake
PV3FoamReader/Allwmake
else
echo "ERROR: ParaView not found in $ParaView_DIR"
fi
# ----------------------------------------------------------------- end-of-file

2
applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H
View File

@ -9,5 +9,5 @@ int USERD_get_number_of_files_in_dataset(void)
// use 1 insted of 0 which gives an un-necessary warning.
Num_dataset_files = 1;
return Num_dataset_files;;
return Num_dataset_files;
}

7
applications/utilities/postProcessing/wall/wallShearStress/Make/options
View File

@ -1,11 +1,14 @@
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleRASModels \
-lfiniteVolume \
-lgenericPatchFields

22
applications/utilities/postProcessing/wall/wallShearStress/createFields.H
View File

@ -1,22 +0,0 @@
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> RASModel
(
incompressible::RASModel::New(U, phi, laminarTransport)
);

153
applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C
View File

@ -25,23 +25,130 @@ Application
wallShearStress
Description
Calculates and writes the wall shear stress, for the specified times.
Calculates and reports wall shear stress for all patches, for the
specified times when using RAS turbulence models.
Default behaviour assumes operating in incompressible mode.
Use the -compressible option for compressible RAS cases.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
#include "RASModel.H"
#include "incompressible/RAS/RASModel/RASModel.H"
#include "basicPsiThermo.H"
#include "compressible/RAS/RASModel/RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void calcIncompressible
(
const fvMesh& mesh,
const Time& runTime,
const volVectorField& U,
volVectorField& wallShearStress
)
{
#include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> model
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
const volSymmTensorField Reff(model->devReff());
forAll(wallShearStress.boundaryField(), patchI)
{
wallShearStress.boundaryField()[patchI] =
(
-mesh.Sf().boundaryField()[patchI]
/mesh.magSf().boundaryField()[patchI]
) & Reff.boundaryField()[patchI];
}
}
void calcCompressible
(
const fvMesh& mesh,
const Time& runTime,
const volVectorField& U,
volVectorField& wallShearStress
)
{
IOobject rhoHeader
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
);
if (!rhoHeader.headerOk())
{
Info<< " no rho field" << endl;
return;
}
Info<< "Reading field rho\n" << endl;
volScalarField rho(rhoHeader, mesh);
#include "compressibleCreatePhi.H"
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
autoPtr<compressible::RASModel> model
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo
)
);
const volSymmTensorField Reff(model->devRhoReff());
forAll(wallShearStress.boundaryField(), patchI)
{
wallShearStress.boundaryField()[patchI] =
(
-mesh.Sf().boundaryField()[patchI]
/mesh.magSf().boundaryField()[patchI]
) & Reff.boundaryField()[patchI];
}
}
int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "addRegionOption.H"
argList::addBoolOption
(
"compressible",
"calculate compressible wall shear stress"
);
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMesh.H"
#include "createNamedMesh.H"
const bool compressible = args.optionFound("compressible");
forAll(timeDirs, timeI)
{
@ -49,10 +156,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();
#include "createFields.H"
volSymmTensorField Reff(RASModel->devReff());
volVectorField wallShearStress
(
IOobject
@ -67,19 +170,41 @@ int main(int argc, char *argv[])
dimensionedVector
(
"wallShearStress",
Reff.dimensions(),
sqr(dimLength)/sqr(dimTime),
vector::zero
)
);
forAll(wallShearStress.boundaryField(), patchi)
IOobject UHeader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
);
if (UHeader.headerOk())
{
wallShearStress.boundaryField()[patchi] =
(
-mesh.Sf().boundaryField()[patchi]
/mesh.magSf().boundaryField()[patchi]
) & Reff.boundaryField()[patchi];
Info<< "Reading field U\n" << endl;
volVectorField U(UHeader, mesh);
if (compressible)
{
calcCompressible(mesh, runTime, U, wallShearStress);
}
else
{
calcIncompressible(mesh, runTime, U, wallShearStress);
}
}
else
{
Info<< " no U field" << endl;
}
Info<< "Writing wall shear stress to field " << wallShearStress.name()
<< nl << endl;
wallShearStress.write();
}

20
applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
View File

@ -35,12 +35,13 @@ License
namespace Foam
{
template<class Type>
template<class Type, class CombineOp>
void MapConsistentVolFields
(
const IOobjectList& objects,
const meshToMesh& meshToMeshInterp,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const CombineOp& cop
)
{
const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@ -84,7 +85,13 @@ void MapConsistentVolFields
);
// Interpolate field
meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
// Write field
fieldTarget.write();
@ -97,7 +104,12 @@ void MapConsistentVolFields
GeometricField<Type, fvPatchField, volMesh> fieldTarget
(
fieldTargetIOobject,
meshToMeshInterp.interpolate(fieldSource, mapOrder)
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldSource,
mapOrder,
cop
)
);
// Write field

263
applications/utilities/preProcessing/mapFields/MapMeshes.H Normal file
View File

@ -0,0 +1,263 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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/>.
\*---------------------------------------------------------------------------*/
#ifndef MapMeshes_H
#define MapMeshes_H
#include "MapVolFields.H"
#include "MapConsistentVolFields.H"
#include "mapLagrangian.H"
#include "UnMapped.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<template<class> class CombineOp>
void MapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp(meshSource, meshTarget);
Info<< nl
<< "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
{
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapConsistentVolFields<scalar>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<scalar>()
);
MapConsistentVolFields<vector>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<vector>()
);
MapConsistentVolFields<sphericalTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapConsistentVolFields<symmTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<symmTensor>()
);
MapConsistentVolFields<tensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<tensor>()
);
}
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
}
mapLagrangian(meshToMeshInterp);
}
template<template<class> class CombineOp>
void MapSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::order& mapOrder
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
meshSource,
meshTarget,
patchMap,
cuttingPatches
);
Info<< nl
<< "Mapping fields for time " << meshSource.time().timeName()
<< nl << endl;
{
// Search for list of source objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<scalar>()
);
MapVolFields<vector>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
(
objects,
meshToMeshInterp,
mapOrder,
CombineOp<tensor>()
);
}
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
}
mapLagrangian(meshToMeshInterp);
}
template<template<class> class CombineOp>
void MapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
)
{
HashTable<word> patchMap;
HashTable<label> cuttingPatchTable;
forAll(meshTarget.boundary(), patchi)
{
if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
{
patchMap.insert
(
meshTarget.boundary()[patchi].name(),
meshTarget.boundary()[patchi].name()
);
}
else
{
cuttingPatchTable.insert
(
meshTarget.boundaryMesh()[patchi].name(),
-1
);
}
}
MapSubMesh<CombineOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatchTable.toc(),
mapOrder
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

13
applications/utilities/preProcessing/mapFields/MapVolFields.H
View File

@ -35,12 +35,13 @@ License
namespace Foam
{
template<class Type>
template<class Type, class CombineOp>
void MapVolFields
(
const IOobjectList& objects,
const meshToMesh& meshToMeshInterp,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const CombineOp& cop
)
{
const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@ -84,7 +85,13 @@ void MapVolFields
);
// Interpolate field
meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
meshToMeshInterp.interpolate//<Type, eqOp<Type> >
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
// Write field
fieldTarget.write();

217
applications/utilities/preProcessing/mapFields/mapFields.C
View File

@ -34,11 +34,8 @@ Description
#include "fvCFD.H"
#include "meshToMesh.H"
#include "MapVolFields.H"
#include "MapConsistentVolFields.H"
#include "UnMapped.H"
#include "processorFvPatch.H"
#include "mapLagrangian.H"
#include "MapMeshes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,56 +43,28 @@ void mapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp(meshSource, meshTarget);
Info<< nl
<< "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
if (subtract)
{
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapConsistentVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<vector>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<sphericalTensor>
MapConsistentMesh<minusEqOp>
(
objects,
meshToMeshInterp,
meshSource,
meshTarget,
mapOrder
);
MapConsistentVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
MapConsistentVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
}
else
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
MapConsistentMesh<eqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
mapLagrangian(meshToMeshInterp);
}
@ -105,57 +74,32 @@ void mapSubMesh
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
meshSource,
meshTarget,
patchMap,
cuttingPatches
);
Info<< nl
<< "Mapping fields for time " << meshSource.time().timeName()
<< nl << endl;
if (subtract)
{
// Search for list of source objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
MapVolFields<vector>(objects, meshToMeshInterp, mapOrder);
MapVolFields<sphericalTensor>(objects, meshToMeshInterp, mapOrder);
MapVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
MapVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
MapSubMesh<minusEqOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
else
{
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
MapSubMesh<eqOp>
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
mapLagrangian(meshToMeshInterp);
}
@ -163,40 +107,28 @@ void mapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::order& mapOrder
const meshToMesh::order& mapOrder,
const bool subtract
)
{
HashTable<word> patchMap;
HashTable<label> cuttingPatchTable;
forAll(meshTarget.boundary(), patchi)
if (subtract)
{
if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
{
patchMap.insert
(
meshTarget.boundary()[patchi].name(),
meshTarget.boundary()[patchi].name()
);
}
else
{
cuttingPatchTable.insert
(
meshTarget.boundaryMesh()[patchi].name(),
-1
);
}
MapConsistentSubMesh<minusEqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
else
{
MapConsistentSubMesh<eqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
cuttingPatchTable.toc(),
mapOrder
);
}
@ -288,6 +220,11 @@ int main(int argc, char *argv[])
"word",
"specify the mapping method"
);
argList::addBoolOption
(
"subtract",
"subtract mapped source from target"
);
argList args(argc, argv);
@ -350,6 +287,13 @@ int main(int argc, char *argv[])
Info<< "Mapping method: " << mapMethod << endl;
}
const bool subtract = args.optionFound("subtract");
if (subtract)
{
Info<< "Subtracting mapped source field from target" << endl;
}
#include "createTimes.H"
HashTable<word> patchMap;
@ -431,7 +375,13 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
@ -441,7 +391,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
cuttingPatches,
mapOrder
mapOrder,
subtract
);
}
}
@ -503,7 +454,13 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
@ -513,7 +470,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder
mapOrder,
subtract
);
}
}
@ -629,7 +587,8 @@ int main(int argc, char *argv[])
(
meshSource,
meshTarget,
mapOrder
mapOrder,
subtract
);
}
else
@ -640,7 +599,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder
mapOrder,
subtract
);
}
}
@ -679,7 +639,7 @@ int main(int argc, char *argv[])
if (consistent)
{
mapConsistentMesh(meshSource, meshTarget, mapOrder);
mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
}
else
{
@ -689,7 +649,8 @@ int main(int argc, char *argv[])
meshTarget,
patchMap,
cuttingPatches,
mapOrder
mapOrder,
subtract
);
}
}

2
src/Allwmake
View File

@ -46,7 +46,7 @@ wmake $makeType genericPatchFields
# Build the proper scotchDecomp, metisDecomp etc.
parallel/Allwmake $*
wmake $makeType renumberMethods
renumber/Allwmake $*
wmake $makeType conversion

6
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View File

@ -144,11 +144,13 @@ public:
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
// mesh face for this tet from the supplied mesh. Normal of
// the tri points out of the cell.
inline triPointRef faceTri(const polyMesh& mesh) const;
//- Return the point indices corresponding to the tri on the mesh
// face for this tet from the supplied mesh
// face for this tet from the supplied mesh. Normal of
// the tri points out of the cell.
inline triFace faceTriIs(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the

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

@ -277,9 +277,9 @@ Foam::scalar Foam::triangle<Point, PointRef>::barycentric
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];
bary[1] = (d11*d20 - d01*d21)/denom;
bary[2] = (d00*d21 - d01*d20)/denom;
bary[0] = 1.0 - bary[1] - bary[2];
return denom;
}

3
src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
View File

@ -196,9 +196,10 @@ bool Foam::pimpleControl::loop()
bool completed = false;
if (criteriaSatisfied())
{
Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
Info<< algorithmName_ << ": converged in " << corr_ - 1 << " iterations"
<< endl;
completed = true;
corr_ = 0;
}
else
{

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

@ -83,8 +83,8 @@ void Foam::cellPointWeight::findTetrahedron
)
{
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
return;
}
@ -191,8 +191,8 @@ void Foam::cellPointWeight::findTriangle
weights_[3] = triWeights[2];
faceVertices_[0] = f[tetIs.faceBasePt()];
faceVertices_[1] = f[tetIs.facePtA()];;
faceVertices_[2] = f[tetIs.facePtB()];;
faceVertices_[1] = f[tetIs.facePtA()];
faceVertices_[2] = f[tetIs.facePtB()];
return;
}

26
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -286,6 +286,10 @@ bool Foam::meshRefinement::isCollapsedFace
const label faceI
) const
{
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(maxNonOrtho));
vector s = mesh_.faces()[faceI].normal(points);
scalar magS = mag(s);
@ -301,11 +305,11 @@ bool Foam::meshRefinement::isCollapsedFace
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
vector d = ownCc - mesh_.cellCentres()[nei];
vector d = mesh_.cellCentres()[nei] - ownCc;
scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
if (dDotS < maxNonOrtho)
if (dDotS < severeNonorthogonalityThreshold)
{
return true;
}
@ -320,11 +324,11 @@ bool Foam::meshRefinement::isCollapsedFace
if (mesh_.boundaryMesh()[patchI].coupled())
{
vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
vector d = neiCc[faceI-mesh_.nInternalFaces()] - ownCc;
scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
if (dDotS < maxNonOrtho)
if (dDotS < severeNonorthogonalityThreshold)
{
return true;
}
@ -337,7 +341,7 @@ bool Foam::meshRefinement::isCollapsedFace
{
// Collapsing normal boundary face does not cause problems with
// non-orthogonality
return true;
return false;
}
}
}
@ -663,7 +667,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
)
{
nPrevented++;
//Pout<< "Preventing collapse of 8 anchor point cell "
//Pout<< "Preventing baffling/removal of 8 anchor point"
// << " cell "
// << cellI << " at " << mesh_.cellCentres()[cellI]
// << " since new volume "
// << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
@ -747,7 +752,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
)
{
nPrevented++;
//Pout<< "Preventing collapse of 7 anchor cell "
//Pout<< "Preventing baffling of 7 anchor cell "
// << cellI
// << " at " << mesh_.cellCentres()[cellI]
// << " since new volume "
@ -853,7 +858,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
)
{
nPrevented++;
//Pout<< "Preventing collapse of face " << faceI
//Pout<< "Preventing baffling (to avoid collapse) of face "
// << faceI
// << " with all boundary edges "
// << " at " << mesh_.faceCentres()[faceI]
// << endl;
@ -909,7 +915,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
)
{
nPrevented++;
//Pout<< "Preventing collapse of coupled face "
//Pout<< "Preventing baffling of coupled face "
// << faceI
// << " with all boundary edges "
// << " at " << mesh_.faceCentres()[faceI]
@ -943,7 +949,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
{
Info<< "markFacesOnProblemCells : prevented "
<< returnReduce(nPrevented, sumOp<label>())
<< " internal faces fom getting converted into baffles."
<< " internal faces from getting converted into baffles."
<< endl;
}

19
src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
View File

@ -236,8 +236,8 @@ void Foam::decompositionMethod::calcCellCells
label globalNei = globalNeighbour[bFaceI];
if
(
globalAgglom.isLocal(globalNei)
&& globalAgglom.toLocal(globalNei) != own
!globalAgglom.isLocal(globalNei)
|| globalAgglom.toLocal(globalNei) != own
)
{
nFacesPerCell[own]++;
@ -285,10 +285,11 @@ void Foam::decompositionMethod::calcCellCells
label own = agglom[faceOwner[faceI]];
label globalNei = globalNeighbour[bFaceI];
if
(
globalAgglom.isLocal(globalNei)
&& globalAgglom.toLocal(globalNei) != own
!globalAgglom.isLocal(globalNei)
|| globalAgglom.toLocal(globalNei) != own
)
{
m[offsets[own] + nFacesPerCell[own]++] = globalNei;
@ -304,15 +305,15 @@ void Foam::decompositionMethod::calcCellCells
// Check for duplicates connections between cells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Done as postprocessing step since we now have cellCells.
label startIndex = 0;
label newIndex = 0;
labelHashSet nbrCells;
forAll(cellCells, cellI)
{
nbrCells.clear();
nbrCells.insert(cellI);
nbrCells.insert(globalAgglom.toGlobal(cellI));
label& endIndex = cellCells.offsets()[cellI+1];
label startIndex = cellCells.offsets()[cellI];
label endIndex = cellCells.offsets()[cellI+1];
for (label i = startIndex; i < endIndex; i++)
{
@ -321,9 +322,7 @@ void Foam::decompositionMethod::calcCellCells
cellCells.m()[newIndex++] = cellCells.m()[i];
}
}
startIndex = endIndex;
endIndex = newIndex;
cellCells.offsets()[cellI+1] = newIndex;
}
cellCells.m().setSize(newIndex);

23
src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
View File

@ -371,29 +371,6 @@ Foam::label Foam::ptscotchDecomp::decompose
Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl;
}
if (debug)
{
Pout<< "nProcessors_:" << nProcessors_ << endl;
globalIndex globalCells(xadj.size()-1);
Pout<< "Xadj:" << endl;
for (label cellI = 0; cellI < xadj.size()-1; cellI++)
{
Pout<< "cell:" << cellI
<< " global:" << globalCells.toGlobal(cellI)
<< " connected to:" << endl;
label start = xadj[cellI];
label end = xadj[cellI+1];
for (label i = start; i < end; i++)
{
Pout<< " cell:" << adjncy[i] << endl;
}
}
Pout<< endl;
}
// Dump graph
if (decompositionDict_.found("ptscotchCoeffs"))
{

239
src/postProcessing/functionObjects/field/streamLine/streamLine.C
View File

@ -34,6 +34,7 @@ License
#include "globalIndex.H"
#include "mapDistribute.H"
#include "interpolationCellPoint.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,6 +42,57 @@ defineTypeNameAndDebug(Foam::streamLine, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::streamLine::wallPatch() const
{
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
nFaces += patches[patchI].size();
}
}
labelList addressing(nFaces);
nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
addressing[nFaces++] = pp.start()+i;
}
}
}
return autoPtr<indirectPrimitivePatch>
(
new indirectPrimitivePatch
(
IndirectList<face>
(
mesh.faces(),
addressing
),
mesh.points()
)
);
}
void Foam::streamLine::track()
{
const Time& runTime = obr_.time();
@ -72,7 +124,7 @@ void Foam::streamLine::track()
label nSeeds = returnReduce(particles.size(), sumOp<label>());
Info<< "Seeded " << nSeeds << " particles." << endl;
Info<< type() << " : seeded " << nSeeds << " particles." << endl;
// Read or lookup fields
PtrList<volScalarField> vsFlds;
@ -163,7 +215,6 @@ void Foam::streamLine::track()
vsInterp.set
(
nScalar++,
//new interpolationCellPoint<scalar>(f)
interpolation<scalar>::New
(
interpolationScheme_,
@ -186,7 +237,6 @@ void Foam::streamLine::track()
vvInterp.set
(
nVector++,
//new interpolationCellPoint<vector>(f)
interpolation<vector>::New
(
interpolationScheme_,
@ -241,6 +291,120 @@ void Foam::streamLine::track()
allVectors_[i].setCapacity(nSeeds);
}
//- For surface following: per face the normal to constrain the velocity
// with.
pointField patchNormals;
if (followSurface_)
{
// Determine geometric data on the non-constraint patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - point normals
// - neighbouring cells
autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
// Calculate parallel consistent normal (on boundaryPatch points only)
pointField boundaryNormals
(
PatchTools::pointNormals
(
mesh,
boundaryPatch(),
boundaryPatch().addressing()
)
);
// Make sure all points know about point normal. There might be points
// which are not on boundaryPatch() but coupled to points that are.
pointField pointNormals(mesh.nPoints(), vector::zero);
UIndirectList<point>
(
pointNormals,
boundaryPatch().meshPoints()
) = boundaryNormals;
boundaryNormals.clear();
mesh.globalData().syncPointData
(
pointNormals,
maxMagSqrEqOp<point>(),
mapDistribute::transform()
);
// Calculate per-face a single constraint normal:
// - faces on patch: face-normal
// - faces on cell with point on patch: point-normal
// - rest: vector::zero
//
// Note: the disadvantage is that we can only handle faces with one
// point on the patch - otherwise we might pick up the wrong
// velocity.
patchNormals.setSize(mesh.nFaces(), vector::zero);
label nPointNormals = 0;
// 1. faces on patch
UIndirectList<point>(patchNormals, boundaryPatch().addressing()) =
boundaryPatch().faceNormals();
// 2. faces with a point on the patch
forAll(mesh.faces(), faceI)
{
if (patchNormals[faceI] == vector::zero)
{
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
label pointI = f[fp];
if (pointNormals[pointI] != vector::zero)
{
patchNormals[faceI] = pointNormals[pointI];
nPointNormals++;
}
}
}
}
// 3. cells on faces with a point on the patch
forAll(mesh.points(), pointI)
{
if (pointNormals[pointI] != vector::zero)
{
const labelList& pCells = mesh.pointCells(pointI);
forAll(pCells, i)
{
const cell& cFaces = mesh.cells()[pCells[i]];
forAll(cFaces, j)
{
label faceI = cFaces[j];
if (patchNormals[faceI] == vector::zero)
{
patchNormals[faceI] = pointNormals[pointI];
nPointNormals++;
}
}
}
}
}
Info<< type() << " : calculated constrained-tracking normals" << nl
<< "Number of faces in mesh : "
<< returnReduce(mesh.nFaces(), sumOp<label>()) << nl
<< " of which on walls : "
<< returnReduce
(
boundaryPatch().addressing().size(),
sumOp<label>()
) << nl
<< " of which on cell with point on boundary : "
<< returnReduce(nPointNormals, sumOp<label>()) << nl
<< endl;
}
// additional particle info
streamLineParticle::trackingData td
(
@ -248,8 +412,16 @@ void Foam::streamLine::track()
vsInterp,
vvInterp,
UIndex, // index of U in vvInterp
trackForward_, // track in +u direction
nSubCycle_, // step through cells in steps?
trackForward_, // track in +u direction?
nSubCycle_, // automatic track control:step through cells in steps?
trackLength_, // fixed track length
//- For surface-constrained streamlines
followSurface_, // stay on surface?
patchNormals, // per face normal to constrain tracking to
//pointNormals, // per point ,,
trackHeight_, // track height
allTracks_,
allScalars_,
allVectors_
@ -279,7 +451,8 @@ Foam::streamLine::streamLine
name_(name),
obr_(obr),
loadFromFiles_(loadFromFiles),
active_(true)
active_(true),
nSubCycle_(0)
{
// Only active if a fvMesh is available
if (isA<fvMesh>(obr_))
@ -334,6 +507,16 @@ void Foam::streamLine::read(const dictionary& dict)
dict.lookup("U") >> UName_;
}
}
if (findIndex(fields_, UName_) == -1)
{
FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
<< "Velocity field for tracking " << UName_
<< " should be present in the list of fields " << fields_
<< exit(FatalIOError);
}
dict.lookup("trackForward") >> trackForward_;
dict.lookup("lifeTime") >> lifeTime_;
if (lifeTime_ < 1)
@ -343,10 +526,48 @@ void Foam::streamLine::read(const dictionary& dict)
<< exit(FatalError);
}
dict.lookup("nSubCycle") >> nSubCycle_;
if (nSubCycle_ < 1)
bool subCycling = dict.found("nSubCycle");
bool fixedLength = dict.found("trackLength");
if (subCycling && fixedLength)
{
nSubCycle_ = 1;
FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
<< "Cannot both specify automatic time stepping (through '"
<< "nSubCycle' specification) and fixed track length (through '"
<< "trackLength')"
<< exit(FatalIOError);
}
nSubCycle_ = 1;
if (dict.readIfPresent("nSubCycle", nSubCycle_))
{
trackLength_ = VGREAT;
if (nSubCycle_ < 1)
{
nSubCycle_ = 1;
}
Info<< type() << " : automatic track length specified through"
<< " number of sub cycles : " << nSubCycle_ << nl << endl;
}
else
{
dict.lookup("trackLength") >> trackLength_;
Info<< type() << " : fixed track length specified : "
<< trackLength_ << nl << endl;
}
followSurface_ = false;
trackHeight_ = VGREAT;
if (dict.readIfPresent("followSurface", followSurface_))
{
dict.lookup("trackHeight") >> trackHeight_;
Info<< type() << " : surface-constrained streamline at "
<< trackHeight_ << " distance above the surface" << nl << endl;
}
interpolationScheme_ = dict.lookupOrDefault

11
src/postProcessing/functionObjects/field/streamLine/streamLine.H
View File

@ -43,6 +43,7 @@ SourceFiles
#include "vectorList.H"
#include "polyMesh.H"
#include "writer.H"
#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -95,9 +96,17 @@ class streamLine
//- Maximum lifetime (= number of cells) of particle
label lifeTime_;
//- Surface-constrained?
bool followSurface_;
scalar trackHeight_;
//- Number of subcycling steps
label nSubCycle_;
//- Track length
scalar trackLength_;
//- Optional specified name of particles
word cloudName_;
@ -141,6 +150,8 @@ class streamLine
List<DynamicList<vectorList> > allVectors_;
//- Construct patch out of all wall patch faces
autoPtr<indirectPrimitivePatch> wallPatch() const;
//- Do all seeding and tracking
void track();

78
src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
View File

@ -43,7 +43,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
const vector& U
) const
{
streamLineParticle testParticle(*this);
particle testParticle(*this);
bool oldKeepParticle = td.keepParticle;
bool oldSwitchProcessor = td.switchProcessor;
@ -59,7 +59,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
(
const trackingData& td,
const point& position,
const label cellI
const label cellI,
const label faceI
)
{
if (cellI == -1)
@ -76,7 +77,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vsInterp_[scalarI].interpolate
(
position,
cellI
cellI,
faceI
)
);
}
@ -89,7 +91,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vvInterp_[vectorI].interpolate
(
position,
cellI
cellI,
faceI
)
);
}
@ -201,13 +204,24 @@ bool Foam::streamLineParticle::move
// Store current position and sampled velocity.
sampledPositions_.append(position());
vector U = interpolateFields(td, position(), cell());
vector U = interpolateFields(td, position(), cell(), face());
if (!td.trackForward_)
{
U = -U;
}
if (td.followSurface_)
{
//- Simple: remove wall-normal component. Should keep particle
// at same height. Does not work well on warped faces.
const vector& n = td.patchNormals_[tetFaceI_];
if (n != vector::zero)
{
U -= (U&n)*n;
}
}
scalar magU = mag(U);
if (magU < SMALL)
@ -219,7 +233,13 @@ bool Foam::streamLineParticle::move
U /= magU;
if (subIter == 0 && td.nSubCycle_ > 1)
if (td.trackLength_ < GREAT)
{
dt = td.trackLength_;
//Pout<< " subiteration " << subIter
// << " : fixed length: updated dt:" << dt << endl;
}
else if (subIter == 0 && td.nSubCycle_ > 1)
{
// Adapt dt to cross cell in a few steps
dt = calcSubCycleDeltaT(td, dt, U);
@ -264,7 +284,8 @@ bool Foam::streamLineParticle::move
if (debug)
{
Pout<< "streamLineParticle : Removing stagnant particle:"
<< p << " sampled positions:" << sampledPositions_.size()
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
td.keepParticle = false;
@ -273,12 +294,13 @@ bool Foam::streamLineParticle::move
{
// Normal exit. Store last position and fields
sampledPositions_.append(position());
interpolateFields(td, position(), cell());
interpolateFields(td, position(), cell(), face());
if (debug)
{
Pout<< "streamLineParticle : Removing particle:"
<< p << " sampled positions:" << sampledPositions_.size()
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
}
@ -374,8 +396,42 @@ void Foam::streamLineParticle::hitWallPatch
const tetIndices&
)
{
// Remove particle
td.keepParticle = false;
if (td.followSurface_)
{
// Push slightly in. Since we're pushing to tet centre we should still
// be in the same tet.
tetPointRef tet = currentTet();
const point oldPosition(position());
position() += trackingCorrectionTol*(tet.centre() - position());
if (debug)
{
label patchID = patch(face());
Pout<< "hitWallPatch : on wall "
<< mesh().boundaryMesh()[patchID].name()
<< " moved from " << oldPosition
<< " to " << position()
<< " due to centre " << tet.centre() << endl;
}
// Check that still in tet (should always be the case
if (debug && !tet.inside(position()))
{
FatalErrorIn("streamLineParticle::hitWallPatch()")
<< "Pushing particle " << *this
<< " away from wall " << wpp.name()
<< " moved it outside of tet."
<< exit(FatalError);
}
}
else
{
// Remove particle
td.keepParticle = false;
}
}

29
src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
View File

@ -73,6 +73,13 @@ public:
const label UIndex_;
const bool trackForward_;
const label nSubCycle_;
const scalar trackLength_;
//- For surface-constrained tracking
const bool followSurface_;
const pointField& patchNormals_;
const scalar trackHeight_;
DynamicList<vectorList>& allPositions_;
List<DynamicList<scalarList> >& allScalars_;
@ -89,6 +96,12 @@ public:
const label UIndex,
const bool trackForward,
const label nSubCycle,
const scalar trackLength,
const bool followSurface,
const pointField& patchNormals,
const scalar trackHeight,
DynamicList<List<point> >& allPositions,
List<DynamicList<scalarList> >& allScalars,
List<DynamicList<vectorList> >& allVectors
@ -100,6 +113,12 @@ public:
UIndex_(UIndex),
trackForward_(trackForward),
nSubCycle_(nSubCycle),
trackLength_(trackLength),
followSurface_(followSurface),
patchNormals_(patchNormals),
trackHeight_(trackHeight),
allPositions_(allPositions),
allScalars_(allScalars),
allVectors_(allVectors)
@ -135,12 +154,20 @@ private:
const vector& U
) const;
void constrainVelocity
(
trackingData& td,
const scalar dt,
vector& U
);
//- Interpolate all quantities; return interpolated velocity.
vector interpolateFields
(
const trackingData&,
const point&,
const label cellI
const label cellI,
const label faceI
);

30
src/renumber/Allwmake Executable file
View File

@ -0,0 +1,30 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
makeType=${1:-libso}
## get ZOLTAN_ARCH_PATH
#if settings=`$WM_PROJECT_DIR/bin/foamEtcFile config/zoltan.sh`
#then
# . $settings
# echo "using ZOLTAN_ARCH_PATH=$ZOLTAN_ARCH_PATH"
#else
# echo
# echo "Error: no config/zoltan.sh settings"
# echo
#fi
set -x
wmake $makeType renumberMethods
#if [ -n "$ZOLTAN_ARCH_PATH" ]
#then
# wmake $makeType zoltanRenumber
#else
# echo
# echo "Skipping zoltanRenumber"
# echo
#fi
# ----------------------------------------------------------------- end-of-file

14
src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C → src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
View File

@ -63,7 +63,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
)
) const
{
CompactListList<label> cellCells;
decompositionMethod::calcCellCells
@ -79,10 +79,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
if (reverse_)
{
for (label i = 0; i < orderedToOld.size()/2; i++)
{
Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
}
reverse(orderedToOld);
}
return invert(orderedToOld.size(), orderedToOld);
@ -93,16 +90,13 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
(
const labelListList& cellCells,
const pointField& points
)
) const
{
labelList orderedToOld = bandCompression(cellCells);
if (reverse_)
{
for (label i = 0; i < orderedToOld.size()/2; i++)
{
Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
}
reverse(orderedToOld);
}
return invert(orderedToOld.size(), orderedToOld);

6
src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H → src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
View File

@ -82,7 +82,7 @@ public:
//- Return for every coordinate the wanted processor number.
// We need a polyMesh (to be able to load the file)
virtual labelList renumber(const pointField&)
virtual labelList renumber(const pointField&) const
{
notImplemented("CuthillMcKeeRenumber::renumber(const pointField&)");
return labelList(0);
@ -94,7 +94,7 @@ public:
(
const polyMesh& mesh,
const pointField& cc
);
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
@ -103,7 +103,7 @@ public:
(
const labelListList& cellCells,
const pointField& cc
);
) const;
};

1
src/renumberMethods/Make/files → src/renumber/renumberMethods/Make/files
View File

@ -3,5 +3,6 @@ manualRenumber/manualRenumber.C
CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
randomRenumber/randomRenumber.C
springRenumber/springRenumber.C
boundaryFirstRenumber/boundaryFirstRenumber.C
LIB = $(FOAM_LIBBIN)/librenumberMethods

0
src/renumberMethods/Make/options → src/renumber/renumberMethods/Make/options
View File

201
src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.C Normal file
View File

@ -0,0 +1,201 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 "boundaryFirstRenumber.H"
#include "addToRunTimeSelectionTable.H"
#include "bandCompression.H"
#include "decompositionMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(boundaryFirstRenumber, 0);
addToRunTimeSelectionTable
(
renumberMethod,
boundaryFirstRenumber,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::boundaryFirstRenumber::boundaryFirstRenumber
(
const dictionary& renumberDict
)
:
renumberMethod(renumberDict),
reverse_
(
renumberDict.found(typeName + "Coeffs")
? Switch(renumberDict.subDict(typeName + "Coeffs").lookup("reverse"))
: Switch(true)
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::boundaryFirstRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
) const
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
// Distance to boundary. Minimum of neighbours.
labelList distance(mesh.nCells(), -1);
// New order
labelList newOrder(mesh.nCells());
label cellInOrder = 0;
// Starting faces for walk. These are the zero distance faces.
DynamicList<label> frontFaces(mesh.nFaces()-mesh.nInternalFaces());
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
//- Note: cannot check for empty since these are introduced by
// fvMeshSubset. Also most loops don't care about patch type.
//if (!isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
frontFaces.append(pp.start()+i);
}
}
}
// Loop over all frontFaces
label currentDistance = 0;
while (frontFaces.size() > 0)
{
DynamicList<label> frontCells(frontFaces.size());
// Set all of frontFaces' neighbours to current distance
forAll(frontFaces, i)
{
label faceI = frontFaces[i];
if (mesh.isInternalFace(faceI))
{
label ownCellI = own[faceI];
if (distance[ownCellI] == -1)
{
distance[ownCellI] = currentDistance;
frontCells.append(ownCellI);
}
label neiCellI = nei[faceI];
if (distance[neiCellI] == -1)
{
distance[neiCellI] = currentDistance;
frontCells.append(neiCellI);
}
}
else
{
label ownCellI = own[faceI];
if (distance[ownCellI] == -1)
{
distance[ownCellI] = currentDistance;
frontCells.append(ownCellI);
}
}
}
//Pout<< "For distance:" << currentDistance
// << " from " << frontFaces.size() << " faces to "
// << frontCells.size() << " cells." << endl;
// TBD. Determine order within current shell (frontCells). For now
// just add them.
forAll(frontCells, i)
{
newOrder[cellInOrder] = frontCells[i];
cellInOrder++;
}
// From cells to faces
frontFaces.clear();
forAll(frontCells, i)
{
label cellI = frontCells[i];
const cell& cFaces = mesh.cells()[cellI];
forAll(cFaces, i)
{
label faceI = cFaces[i];
if (mesh.isInternalFace(faceI))
{
label nbrCellI =
(
mesh.faceOwner()[faceI] == cellI
? mesh.faceNeighbour()[faceI]
: mesh.faceOwner()[faceI]
);
if (distance[nbrCellI] == -1)
{
frontFaces.append(faceI);
}
}
}
}
currentDistance++;
}
// Return furthest away cell first
if (reverse_)
{
reverse(newOrder);
}
//forAll(newOrder, i)
//{
// label cellI = newOrder[i];
//
// Pout<< "cell:" << cellI << endl;
// Pout<< " at distance:" << distance[cellI]
// << endl;
//}
return invert(newOrder.size(), newOrder);
}
// ************************************************************************* //

129
src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.H Normal file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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::boundaryFirstRenumber
Description
Cuthill-McKee renumbering
SourceFiles
boundaryFirstRenumber.C
\*---------------------------------------------------------------------------*/
#ifndef boundaryFirstRenumber_H
#define boundaryFirstRenumber_H
#include "renumberMethod.H"
#include "Switch.H"
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class boundaryFirstRenumber Declaration
\*---------------------------------------------------------------------------*/
class boundaryFirstRenumber
:
public renumberMethod
{
// Private data
const Switch reverse_;
// Private Member Functions
//- Disallow default bitwise copy construct and assignment
void operator=(const boundaryFirstRenumber&);
boundaryFirstRenumber(const boundaryFirstRenumber&);
public:
//- Runtime type information
TypeName("boundaryFirst");
// Constructors
//- Construct given the renumber dictionary
boundaryFirstRenumber(const dictionary& renumberDict);
//- Destructor
virtual ~boundaryFirstRenumber()
{}
// Member Functions
//- Return for every coordinate the wanted processor number.
// We need a polyMesh (to be able to load the file)
virtual labelList renumber(const pointField&) const
{
notImplemented
(
"boundaryFirstRenumber::renumber(const pointField&)"
);
return labelList(0);
}
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
virtual labelList renumber
(
const polyMesh& mesh,
const pointField& cc
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
// - the connections are across coupled patches
virtual labelList renumber
(
const labelListList& cellCells,
const pointField& cc
) const
{
notImplemented
(
"boundaryFirstRenumber::renumber"
"(const labelListList&, const pointField&)"
);
return labelList(0);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/renumberMethods/manualRenumber/manualRenumber.C → src/renumber/renumberMethods/manualRenumber/manualRenumber.C
View File

@ -61,7 +61,7 @@ Foam::labelList Foam::manualRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
)
) const
{
labelIOList oldToNew
(

6
src/renumberMethods/manualRenumber/manualRenumber.H → src/renumber/renumberMethods/manualRenumber/manualRenumber.H
View File

@ -81,7 +81,7 @@ public:
//- Return for every coordinate the wanted processor number.
// We need a polyMesh (to be able to load the file)
virtual labelList renumber(const pointField&)
virtual labelList renumber(const pointField&) const
{
notImplemented("manualRenumber::renumber(const pointField&)");
return labelList(0);
@ -93,7 +93,7 @@ public:
(
const polyMesh& mesh,
const pointField& cc
);
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
@ -102,7 +102,7 @@ public:
(
const labelListList& cellCells,
const pointField& cc
)
) const
{
notImplemented
(

6
src/renumberMethods/randomRenumber/randomRenumber.C → src/renumber/renumberMethods/randomRenumber/randomRenumber.C
View File

@ -55,7 +55,7 @@ Foam::randomRenumber::randomRenumber(const dictionary& renumberDict)
Foam::labelList Foam::randomRenumber::renumber
(
const pointField& points
)
) const
{
Random rndGen(0);
@ -77,7 +77,7 @@ Foam::labelList Foam::randomRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
)
) const
{
return renumber(points);
}
@ -87,7 +87,7 @@ Foam::labelList Foam::randomRenumber::renumber
(
const labelListList& cellCells,
const pointField& points
)
) const
{
return renumber(points);
}

6
src/renumberMethods/randomRenumber/randomRenumber.H → src/renumber/renumberMethods/randomRenumber/randomRenumber.H
View File

@ -76,7 +76,7 @@ public:
//- Return for every coordinate the wanted processor number.
// We need a polyMesh (to be able to load the file)
virtual labelList renumber(const pointField&);
virtual labelList renumber(const pointField&) const;
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
@ -84,7 +84,7 @@ public:
(
const polyMesh& mesh,
const pointField& cc
);
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
@ -93,7 +93,7 @@ public:
(
const labelListList& cellCells,
const pointField& cc
);
) const;
};

4
src/renumberMethods/renumberMethod/renumberMethod.C → src/renumber/renumberMethods/renumberMethod/renumberMethod.C
View File

@ -72,7 +72,7 @@ Foam::labelList Foam::renumberMethod::renumber
(
const polyMesh& mesh,
const pointField& points
)
) const
{
CompactListList<label> cellCells;
decompositionMethod::calcCellCells
@ -94,7 +94,7 @@ Foam::labelList Foam::renumberMethod::renumber
const polyMesh& mesh,
const labelList& fineToCoarse,
const pointField& coarsePoints
)
) const
{
CompactListList<label> coarseCellCells;
decompositionMethod::calcCellCells

8
src/renumberMethods/renumberMethod/renumberMethod.H → src/renumber/renumberMethods/renumberMethod/renumberMethod.H
View File

@ -111,7 +111,7 @@ public:
//- Return for every cell the new cell label.
// This is only defined for geometric renumberMethods.
virtual labelList renumber(const pointField&)
virtual labelList renumber(const pointField&) const
{
notImplemented
(
@ -122,7 +122,7 @@ public:
//- Return for every cell the new cell label. Use the
// mesh connectivity (if needed)
virtual labelList renumber(const polyMesh&, const pointField&);
virtual labelList renumber(const polyMesh&, const pointField&) const;
//- Return for every cell the new cell label. Gets
// passed agglomeration map (from fine to coarse cells) and coarse
@ -136,7 +136,7 @@ public:
const polyMesh& mesh,
const labelList& cellToRegion,
const pointField& regionPoints
);
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
@ -145,7 +145,7 @@ public:
(
const labelListList& cellCells,
const pointField& cc
) = 0;
) const = 0;
};

4
src/renumberMethods/springRenumber/springRenumber.C → src/renumber/renumberMethods/springRenumber/springRenumber.C
View File

@ -60,7 +60,7 @@ Foam::labelList Foam::springRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
)
) const
{
CompactListList<label> cellCells;
decompositionMethod::calcCellCells
@ -80,7 +80,7 @@ Foam::labelList Foam::springRenumber::renumber
(
const labelListList& cellCells,
const pointField& points
)
) const
{
// Look at cell index as a 1D position parameter.
// Move cells to the average 'position' of their neighbour.

6
src/renumberMethods/springRenumber/springRenumber.H → src/renumber/renumberMethods/springRenumber/springRenumber.H
View File

@ -96,7 +96,7 @@ public:
// Member Functions
//- Return for every coordinate the wanted processor number.
virtual labelList renumber(const pointField&)
virtual labelList renumber(const pointField&) const
{
notImplemented("springRenumber::renumber(const pointField&)");
return labelList(0);
@ -108,7 +108,7 @@ public:
(
const polyMesh& mesh,
const pointField& cc
);
) const;
//- Return for every cell the new cell label.
// The connectivity is equal to mesh.cellCells() except
@ -117,7 +117,7 @@ public:
(
const labelListList& cellCells,
const pointField& cc
);
) const;
};

45
src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
View File

@ -219,84 +219,93 @@ public:
// Interpolation
//- Map field
template<class Type>
template<class Type, class CombineOp>
void mapField
(
Field<Type>&,
const Field<Type>&,
const labelList& adr
const labelList& adr,
const CombineOp& cop
) const;
//- Interpolate field using inverse-distance weights
template<class Type>
template<class Type, class CombineOp>
void interpolateField
(
Field<Type>&,
const GeometricField<Type, fvPatchField, volMesh>&,
const labelList& adr,
const scalarListList& weights
const scalarListList& weights,
const CombineOp& cop
) const;
//- Interpolate field using cell-point interpolation
template<class Type>
template<class Type, class CombineOp>
void interpolateField
(
Field<Type>&,
const GeometricField<Type, fvPatchField, volMesh>&,
const labelList& adr,
const vectorField& centres
const vectorField& centres,
const CombineOp& cop
) const;
//- Interpolate internal volume field
template<class Type>
template<class Type, class CombineOp>
void interpolateInternalField
(
Field<Type>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
template<class Type>
template<class Type, class CombineOp>
void interpolateInternalField
(
Field<Type>&,
const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
//- Interpolate volume field
template<class Type>
template<class Type, class CombineOp>
void interpolate
(
GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
template<class Type>
template<class Type, class CombineOp>
void interpolate
(
GeometricField<Type, fvPatchField, volMesh>&,
const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
//- Interpolate volume field
template<class Type>
template<class Type, class CombineOp>
tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
template<class Type>
template<class Type, class CombineOp>
tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
(
const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
order=INTERPOLATE
order=INTERPOLATE,
const CombineOp& cop = eqOp<Type>()
) const;
};

121
src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
View File

@ -31,12 +31,13 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::mapField
(
Field<Type>& toF,
const Field<Type>& fromVf,
const labelList& adr
const labelList& adr,
const CombineOp& cop
) const
{
// Direct mapping of nearest-cell values
@ -45,7 +46,7 @@ void Foam::meshToMesh::mapField
{
if (adr[celli] != -1)
{
toF[celli] = fromVf[adr[celli]];
cop(toF[celli], fromVf[adr[celli]]);
}
}
@ -53,13 +54,14 @@ void Foam::meshToMesh::mapField
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolateField
(
Field<Type>& toF,
const GeometricField<Type, fvPatchField, volMesh>& fromVf,
const labelList& adr,
const scalarListList& weights
const scalarListList& weights,
const CombineOp& cop
) const
{
// Inverse distance weighted interpolation
@ -74,24 +76,27 @@ void Foam::meshToMesh::interpolateField
const labelList& neighbours = cc[adr[celli]];
const scalarList& w = weights[celli];
toF[celli] = fromVf[adr[celli]]*w[0];
Type f = fromVf[adr[celli]]*w[0];
for (label ni = 1; ni < w.size(); ni++)
{
toF[celli] += fromVf[neighbours[ni - 1]]*w[ni];
f += fromVf[neighbours[ni - 1]]*w[ni];
}
cop(toF[celli], f);
}
}
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolateField
(
Field<Type>& toF,
const GeometricField<Type, fvPatchField, volMesh>& fromVf,
const labelList& adr,
const vectorField& centres
const vectorField& centres,
const CombineOp& cop
) const
{
// Cell-Point interpolation
@ -101,31 +106,36 @@ void Foam::meshToMesh::interpolateField
{
if (adr[celli] != -1)
{
toF[celli] = interpolator.interpolate
cop
(
centres[celli],
adr[celli]
toF[celli],
interpolator.interpolate
(
centres[celli],
adr[celli]
)
);
}
}
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolateInternalField
(
Field<Type>& toF,
const GeometricField<Type, fvPatchField, volMesh>& fromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
if (fromVf.mesh() != fromMesh_)
{
FatalErrorIn
(
"meshToMesh::interpolateInternalField(Field<Type>& toF, "
"const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
"meshToMesh::order ord) const"
"meshToMesh::interpolateInternalField(Field<Type>&, "
"const GeometricField<Type, fvPatchField, volMesh>&, "
"meshToMesh::order, const CombineOp&) const"
) << "the argument field does not correspond to the right mesh. "
<< "Field size: " << fromVf.size()
<< " mesh size: " << fromMesh_.nCells()
@ -136,9 +146,9 @@ void Foam::meshToMesh::interpolateInternalField
{
FatalErrorIn
(
"meshToMesh::interpolateInternalField(Field<Type>& toF, "
"const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
"meshToMesh::order ord) const"
"meshToMesh::interpolateInternalField(Field<Type>&, "
"const GeometricField<Type, fvPatchField, volMesh>&, "
"meshToMesh::order, const CombineOp&) const"
) << "the argument field does not correspond to the right mesh. "
<< "Field size: " << toF.size()
<< " mesh size: " << toMesh_.nCells()
@ -148,7 +158,7 @@ void Foam::meshToMesh::interpolateInternalField
switch(ord)
{
case MAP:
mapField(toF, fromVf, cellAddressing_);
mapField(toF, fromVf, cellAddressing_, cop);
break;
case INTERPOLATE:
@ -157,7 +167,8 @@ void Foam::meshToMesh::interpolateInternalField
toF,
fromVf,
cellAddressing_,
inverseDistanceWeights()
inverseDistanceWeights(),
cop
);
break;
@ -167,44 +178,47 @@ void Foam::meshToMesh::interpolateInternalField
toF,
fromVf,
cellAddressing_,
toMesh_.cellCentres()
toMesh_.cellCentres(),
cop
);
break;
default:
FatalErrorIn
(
"meshToMesh::interpolateInternalField(Field<Type>& toF, "
"const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
"meshToMesh::order ord) const"
"meshToMesh::interpolateInternalField(Field<Type>&, "
"const GeometricField<Type, fvPatchField, volMesh>&, "
"meshToMesh::order, const CombineOp&) const"
) << "unknown interpolation scheme " << ord
<< exit(FatalError);
}
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolateInternalField
(
Field<Type>& toF,
const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
interpolateInternalField(toF, tfromVf(), ord);
interpolateInternalField(toF, tfromVf(), ord, cop);
tfromVf.clear();
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolate
(
GeometricField<Type, fvPatchField, volMesh>& toVf,
const GeometricField<Type, fvPatchField, volMesh>& fromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
interpolateInternalField(toVf, fromVf, ord);
interpolateInternalField(toVf, fromVf, ord, cop);
forAll(toMesh_.boundaryMesh(), patchi)
{
@ -219,7 +233,8 @@ void Foam::meshToMesh::interpolate
(
toVf.boundaryField()[patchi],
fromVf,
boundaryAddressing_[patchi]
boundaryAddressing_[patchi],
cop
);
break;
@ -229,7 +244,8 @@ void Foam::meshToMesh::interpolate
toVf.boundaryField()[patchi],
fromVf,
boundaryAddressing_[patchi],
toPatch.Cf()
toPatch.Cf(),
cop
);
break;
@ -239,7 +255,8 @@ void Foam::meshToMesh::interpolate
toVf.boundaryField()[patchi],
fromVf,
boundaryAddressing_[patchi],
toPatch.Cf()
toPatch.Cf(),
cop
);
break;
@ -247,9 +264,9 @@ void Foam::meshToMesh::interpolate
FatalErrorIn
(
"meshToMesh::interpolate("
"GeometricField<Type, fvPatchField, volMesh>& toVf, "
"const GeometricField<Type, fvPatchField, volMesh>& "
"fromVf, meshToMesh::order ord) const"
"GeometricField<Type, fvPatchField, volMesh>&, "
"const GeometricField<Type, fvPatchField, volMesh>&, "
"meshToMesh::order, const CombineOp&) const"
) << "unknown interpolation scheme " << ord
<< exit(FatalError);
}
@ -286,37 +303,40 @@ void Foam::meshToMesh::interpolate
[
fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
],
boundaryAddressing_[patchi]
boundaryAddressing_[patchi],
cop
);
}
}
}
template<class Type>
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolate
(
GeometricField<Type, fvPatchField, volMesh>& toVf,
const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
interpolate(toVf, tfromVf(), ord);
interpolate(toVf, tfromVf(), ord, cop);
tfromVf.clear();
}
template<class Type>
template<class Type, class CombineOp>
Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& fromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
// Create and map the internal-field values
Field<Type> internalField(toMesh_.nCells());
interpolateInternalField(internalField, fromVf, ord);
interpolateInternalField(internalField, fromVf, ord, cop);
// check whether both meshes have got the same number
// of boundary patches
@ -325,8 +345,8 @@ Foam::meshToMesh::interpolate
FatalErrorIn
(
"meshToMesh::interpolate"
"(const GeometricField<Type, fvPatchField, volMesh>& fromVf,"
"meshToMesh::order ord) const"
"(const GeometricField<Type, fvPatchField, volMesh>&,"
"meshToMesh::order, const CombineOp&) const"
) << "Incompatible meshes: different number of boundaries, "
"only internal field may be interpolated"
<< exit(FatalError);
@ -381,16 +401,17 @@ Foam::meshToMesh::interpolate
}
template<class Type>
template<class Type, class CombineOp>
Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::interpolate
(
const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
meshToMesh::order ord
meshToMesh::order ord,
const CombineOp& cop
) const
{
tmp<GeometricField<Type, fvPatchField, volMesh> > tint =
interpolate(tfromVf(), ord);
interpolate(tfromVf(), ord, cop);
tfromVf.clear();
return tint;

11
tutorials/incompressible/simpleFoam/motorBike/system/streamLines
View File

@ -28,8 +28,15 @@ streamLines
// Steps particles can travel before being removed
lifeTime 10000;
// Number of steps per cell (estimate). Set to 1 to disable subcycling.
nSubCycle 5;
//- Specify either absolute length of steps (trackLength) or a number
// of subcycling steps per cell (nSubCycle)
// Size of single track segment [m]
//trackLength 1e-3;
// Number of steps per cell (estimate). Set to 1 to disable subcycling.
nSubCycle 5;
// Cloud name to use
cloudName particleTracks;

2
tutorials/incompressible/simpleFoam/turbineSiting/system/snappyHexMeshDict
View File

@ -173,7 +173,7 @@ castellatedMeshControls
}
}
resolveFeatureAngle 2;
resolveFeatureAngle 30;
// Region-wise refinement