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 /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2011-12-22 12:27:34 +00:00
parent ddfd2bfb38 b752c60cd1
commit c2b0d66f3f
37 changed files with 1474 additions and 337 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -629,8 +629,9 @@ Foam::multiphaseSystem::dragCoeffs() const
( (
max max
( (
fvc::average(dm.phase1())*fvc::average(dm.phase2()), //fvc::average(dm.phase1()*dm.phase2()),
//dm.phase1()*dm.phase2(), //fvc::average(dm.phase1())*fvc::average(dm.phase2()),
dm.phase1()*dm.phase2(),
dm.residualPhaseFraction() dm.residualPhaseFraction()
) )
*dm.K *dm.K

8
applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C
View File

@ -212,10 +212,10 @@ bool Foam::phaseModel::read(const dictionary& phaseDict)
//if (nuModel_->read(phaseDict_)) //if (nuModel_->read(phaseDict_))
{ {
phaseDict_.lookup("nu") >> nu_; phaseDict_.lookup("nu") >> nu_.value();
phaseDict_.lookup("kappa") >> kappa_; phaseDict_.lookup("kappa") >> kappa_.value();
phaseDict_.lookup("Cp") >> Cp_; phaseDict_.lookup("Cp") >> Cp_.value();
phaseDict_.lookup("rho") >> rho_; phaseDict_.lookup("rho") >> rho_.value();
return true; return true;
} }

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

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

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

@ -29,7 +29,7 @@ Description
renumbering all fields from all the time directories. renumbering all fields from all the time directories.
By default uses bandCompression (CuthillMcKee) but will 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 "renumberMethod.H"
#include "zeroGradientFvPatchFields.H" #include "zeroGradientFvPatchFields.H"
#include "CuthillMcKeeRenumber.H" #include "CuthillMcKeeRenumber.H"
#include "fvMeshSubset.H"
using namespace Foam; 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[]) int main(int argc, char *argv[])
@ -505,13 +569,16 @@ int main(int argc, char *argv[])
label blockSize = 0; label blockSize = 0;
// Construct renumberMethod // Construct renumberMethod
autoPtr<IOdictionary> renumberDictPtr;
autoPtr<renumberMethod> renumberPtr; autoPtr<renumberMethod> renumberPtr;
if (readDict) if (readDict)
{ {
Info<< "Renumber according to renumberMeshDict." << nl << endl; Info<< "Renumber according to renumberMeshDict." << nl << endl;
IOdictionary renumberDict renumberDictPtr.reset
(
new IOdictionary
( (
IOobject IOobject
( (
@ -521,7 +588,9 @@ int main(int argc, char *argv[])
IOobject::MUST_READ_IF_MODIFIED, IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE IOobject::NO_WRITE
) )
)
); );
const IOdictionary renumberDict = renumberDictPtr();
renumberPtr = renumberMethod::New(renumberDict); renumberPtr = renumberMethod::New(renumberDict);
@ -562,7 +631,7 @@ int main(int argc, char *argv[])
<< endl; << endl;
} }
writeMaps = readLabel(renumberDict.lookup("writeMaps")); renumberDict.lookup("writeMaps") >> writeMaps;
if (writeMaps) if (writeMaps)
{ {
Info<< "Writing renumber maps (new to old) to polyMesh." << nl Info<< "Writing renumber maps (new to old) to polyMesh." << nl
@ -685,20 +754,13 @@ int main(int argc, char *argv[])
label nBlocks = mesh.nCells() / blockSize; label nBlocks = mesh.nCells() / blockSize;
Info<< "nBlocks = " << nBlocks << endl; Info<< "nBlocks = " << nBlocks << endl;
// Read decomposePar dictionary // Read decompositionMethod dictionary
IOdictionary decomposeDict dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
(
IOobject
(
"decomposeParDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
decomposeDict.set("numberOfSubdomains", nBlocks); decomposeDict.set("numberOfSubdomains", nBlocks);
bool oldParRun = UPstream::parRun();
UPstream::parRun() = false;
autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
( (
decomposeDict decomposeDict
@ -713,6 +775,9 @@ int main(int argc, char *argv[])
) )
); );
// Restore state
UPstream::parRun() = oldParRun;
// For debugging: write out region // For debugging: write out region
createScalarField createScalarField
( (
@ -726,23 +791,7 @@ int main(int argc, char *argv[])
<< nl << endl; << nl << endl;
// Find point per region cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
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);
// Determine new to old face order with new cell numbering // Determine new to old face order with new cell numbering
faceOrder = getRegionFaceOrder faceOrder = getRegionFaceOrder

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

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

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

@ -789,6 +789,14 @@ int main(int argc, char *argv[])
// Point fields // Point fields
if
(
pointScalarFields.size()
|| pointVectorFields.size()
|| pointSphericalTensorFields.size()
|| pointSymmTensorFields.size()
|| pointTensorFields.size()
)
{ {
labelIOList pointProcAddressing labelIOList pointProcAddressing
( (

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

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

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

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

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

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

28
src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
View File

@ -107,16 +107,29 @@ void Foam::GaussSeidelSmoother::smooth
// To compensate for this, it is necessary to turn the // To compensate for this, it is necessary to turn the
// sign of the contribution. // sign of the contribution.
FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size()); //FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size());
//
//forAll(mBouCoeffs, patchi)
//{
// if (interfaces_.set(patchi))
// {
// mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]);
// }
//}
FieldField<Field, scalar>& mBouCoeffs =
const_cast<FieldField<Field, scalar>&>
(
interfaceBouCoeffs_
);
forAll(mBouCoeffs, patchi) forAll(mBouCoeffs, patchi)
{ {
if (interfaces_.set(patchi)) if (interfaces_.set(patchi))
{ {
mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]); mBouCoeffs[patchi].negate();
} }
} }
for (label sweep=0; sweep<nSweeps; sweep++) for (label sweep=0; sweep<nSweeps; sweep++)
{ {
bPrime = source; bPrime = source;
@ -170,6 +183,15 @@ void Foam::GaussSeidelSmoother::smooth
psiPtr[cellI] = curPsi; psiPtr[cellI] = curPsi;
} }
} }
// Restore interfaceBouCoeffs_
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
} }

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

@ -144,11 +144,13 @@ public:
inline tetPointRef oldTet(const polyMesh& mesh) const; inline tetPointRef oldTet(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the //- 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; inline triPointRef faceTri(const polyMesh& mesh) const;
//- Return the point indices corresponding to the tri on the mesh //- 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; inline triFace faceTriIs(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the //- Return the geometry corresponding to the tri on the

4
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
View File

@ -266,7 +266,7 @@ void Foam::cyclicPolyPatch::calcTransforms
// Calculate using the given rotation axis and centre. Do not // Calculate using the given rotation axis and centre. Do not
// use calculated normals. // use calculated normals.
vector n0 = findFaceMaxRadius(half0Ctrs); vector n0 = findFaceMaxRadius(half0Ctrs);
vector n1 = findFaceMaxRadius(half1Ctrs); vector n1 = -findFaceMaxRadius(half1Ctrs);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
@ -424,7 +424,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
case ROTATIONAL: case ROTATIONAL:
{ {
vector n0 = findFaceMaxRadius(half0Ctrs); vector n0 = findFaceMaxRadius(half0Ctrs);
vector n1 = findFaceMaxRadius(half1Ctrs); vector n1 = -findFaceMaxRadius(half1Ctrs);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;

12
src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
View File

@ -44,11 +44,15 @@ void Foam::MRFZone::relativeRhoFlux
const vector& origin = origin_.value(); const vector& origin = origin_.value();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();
scalarField& phii = phi.internalField();
// Internal faces // Internal faces
forAll(internalFaces_, i) forAll(internalFaces_, i)
{ {
label facei = internalFaces_[i]; label facei = internalFaces_[i];
phi[facei] -= rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei]; phii[facei] -= rho[facei]*(Omega ^ (Cfi[facei] - origin)) & Sfi[facei];
} }
// Included patches // Included patches
@ -91,11 +95,15 @@ void Foam::MRFZone::absoluteRhoFlux
const vector& origin = origin_.value(); const vector& origin = origin_.value();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();
scalarField& phii = phi.internalField();
// Internal faces // Internal faces
forAll(internalFaces_, i) forAll(internalFaces_, i)
{ {
label facei = internalFaces_[i]; label facei = internalFaces_[i];
phi[facei] += rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei]; phii[facei] += rho[facei]*(Omega ^ (Cfi[facei] - origin)) & Sfi[facei];
} }
// Included patches // Included patches

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

@ -34,6 +34,7 @@ License
#include "globalIndex.H" #include "globalIndex.H"
#include "mapDistribute.H" #include "mapDistribute.H"
#include "interpolationCellPoint.H" #include "interpolationCellPoint.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,6 +42,57 @@ defineTypeNameAndDebug(Foam::streamLine, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * 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() void Foam::streamLine::track()
{ {
const Time& runTime = obr_.time(); const Time& runTime = obr_.time();
@ -72,7 +124,7 @@ void Foam::streamLine::track()
label nSeeds = returnReduce(particles.size(), sumOp<label>()); label nSeeds = returnReduce(particles.size(), sumOp<label>());
Info<< "Seeded " << nSeeds << " particles." << endl; Info<< type() << " : seeded " << nSeeds << " particles." << endl;
// Read or lookup fields // Read or lookup fields
PtrList<volScalarField> vsFlds; PtrList<volScalarField> vsFlds;
@ -163,7 +215,6 @@ void Foam::streamLine::track()
vsInterp.set vsInterp.set
( (
nScalar++, nScalar++,
//new interpolationCellPoint<scalar>(f)
interpolation<scalar>::New interpolation<scalar>::New
( (
interpolationScheme_, interpolationScheme_,
@ -186,7 +237,6 @@ void Foam::streamLine::track()
vvInterp.set vvInterp.set
( (
nVector++, nVector++,
//new interpolationCellPoint<vector>(f)
interpolation<vector>::New interpolation<vector>::New
( (
interpolationScheme_, interpolationScheme_,
@ -241,6 +291,120 @@ void Foam::streamLine::track()
allVectors_[i].setCapacity(nSeeds); 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 // additional particle info
streamLineParticle::trackingData td streamLineParticle::trackingData td
( (
@ -248,8 +412,16 @@ void Foam::streamLine::track()
vsInterp, vsInterp,
vvInterp, vvInterp,
UIndex, // index of U in vvInterp UIndex, // index of U in vvInterp
trackForward_, // track in +u direction trackForward_, // track in +u direction?
nSubCycle_, // step through cells in steps? 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_, allTracks_,
allScalars_, allScalars_,
allVectors_ allVectors_
@ -279,7 +451,8 @@ Foam::streamLine::streamLine
name_(name), name_(name),
obr_(obr), obr_(obr),
loadFromFiles_(loadFromFiles), loadFromFiles_(loadFromFiles),
active_(true) active_(true),
nSubCycle_(0)
{ {
// Only active if a fvMesh is available // Only active if a fvMesh is available
if (isA<fvMesh>(obr_)) if (isA<fvMesh>(obr_))
@ -334,6 +507,16 @@ void Foam::streamLine::read(const dictionary& dict)
dict.lookup("U") >> UName_; 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("trackForward") >> trackForward_;
dict.lookup("lifeTime") >> lifeTime_; dict.lookup("lifeTime") >> lifeTime_;
if (lifeTime_ < 1) if (lifeTime_ < 1)
@ -343,11 +526,49 @@ void Foam::streamLine::read(const dictionary& dict)
<< exit(FatalError); << exit(FatalError);
} }
dict.lookup("nSubCycle") >> nSubCycle_;
bool subCycling = dict.found("nSubCycle");
bool fixedLength = dict.found("trackLength");
if (subCycling && fixedLength)
{
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) if (nSubCycle_ < 1)
{ {
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 interpolationScheme_ = dict.lookupOrDefault
( (

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

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

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

@ -43,7 +43,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
const vector& U const vector& U
) const ) const
{ {
streamLineParticle testParticle(*this); particle testParticle(*this);
bool oldKeepParticle = td.keepParticle; bool oldKeepParticle = td.keepParticle;
bool oldSwitchProcessor = td.switchProcessor; bool oldSwitchProcessor = td.switchProcessor;
@ -59,7 +59,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
( (
const trackingData& td, const trackingData& td,
const point& position, const point& position,
const label cellI const label cellI,
const label faceI
) )
{ {
if (cellI == -1) if (cellI == -1)
@ -76,7 +77,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vsInterp_[scalarI].interpolate td.vsInterp_[scalarI].interpolate
( (
position, position,
cellI cellI,
faceI
) )
); );
} }
@ -89,7 +91,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vvInterp_[vectorI].interpolate td.vvInterp_[vectorI].interpolate
( (
position, position,
cellI cellI,
faceI
) )
); );
} }
@ -201,13 +204,24 @@ bool Foam::streamLineParticle::move
// Store current position and sampled velocity. // Store current position and sampled velocity.
sampledPositions_.append(position()); sampledPositions_.append(position());
vector U = interpolateFields(td, position(), cell()); vector U = interpolateFields(td, position(), cell(), face());
if (!td.trackForward_) if (!td.trackForward_)
{ {
U = -U; 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); scalar magU = mag(U);
if (magU < SMALL) if (magU < SMALL)
@ -219,7 +233,13 @@ bool Foam::streamLineParticle::move
U /= magU; 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 // Adapt dt to cross cell in a few steps
dt = calcSubCycleDeltaT(td, dt, U); dt = calcSubCycleDeltaT(td, dt, U);
@ -264,7 +284,8 @@ bool Foam::streamLineParticle::move
if (debug) if (debug)
{ {
Pout<< "streamLineParticle : Removing stagnant particle:" Pout<< "streamLineParticle : Removing stagnant particle:"
<< p << " sampled positions:" << sampledPositions_.size() << p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl; << endl;
} }
td.keepParticle = false; td.keepParticle = false;
@ -273,12 +294,13 @@ bool Foam::streamLineParticle::move
{ {
// Normal exit. Store last position and fields // Normal exit. Store last position and fields
sampledPositions_.append(position()); sampledPositions_.append(position());
interpolateFields(td, position(), cell()); interpolateFields(td, position(), cell(), face());
if (debug) if (debug)
{ {
Pout<< "streamLineParticle : Removing particle:" Pout<< "streamLineParticle : Removing particle:"
<< p << " sampled positions:" << sampledPositions_.size() << p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl; << endl;
} }
} }
@ -373,10 +395,44 @@ void Foam::streamLineParticle::hitWallPatch
trackingData& td, trackingData& td,
const tetIndices& const tetIndices&
) )
{
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 // Remove particle
td.keepParticle = false; td.keepParticle = false;
} }
}
void Foam::streamLineParticle::hitPatch void Foam::streamLineParticle::hitPatch

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

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

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

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

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

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

1
src/renumberMethods/Make/files
View File

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

201
src/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/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
View File

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@ -28,9 +28,16 @@ streamLines
// Steps particles can travel before being removed // Steps particles can travel before being removed
lifeTime 10000; lifeTime 10000;
//- 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. // Number of steps per cell (estimate). Set to 1 to disable subcycling.
nSubCycle 5; nSubCycle 5;
// Cloud name to use // Cloud name to use
cloudName particleTracks; cloudName particleTracks;

119
tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/constant/transportProperties
View File

@ -69,7 +69,7 @@ phases
diameterModel constant; diameterModel constant;
constantCoeffs constantCoeffs
{ {
d 1e-3; d 3e-3;
} }
} }
); );
@ -96,9 +96,9 @@ interfaceCompression
virtualMass virtualMass
( (
(air water) 0.5 (air water) 0
(air oil) 0.5 (air oil) 0
(air mercury) 0.5 (air mercury) 0
(water oil) 0.5 (water oil) 0.5
(water mercury) 0.5 (water mercury) 0.5
(oil mercury) 0.5 (oil mercury) 0.5
@ -107,22 +107,133 @@ virtualMass
drag drag
( (
(air water) (air water)
{
type blended;
air
{ {
type SchillerNaumann; type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
water
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3; residualPhaseFraction 1e-3;
residualSlip 1e-3; residualSlip 1e-3;
} }
(air oil) (air oil)
{
type blended;
air
{ {
type SchillerNaumann; type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
oil
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3; residualPhaseFraction 1e-3;
residualSlip 1e-3; residualSlip 1e-3;
} }
(air mercury) (air mercury)
{
type blended;
air
{ {
type SchillerNaumann; type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
mercury
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3;
residualSlip 1e-3;
}
(water oil)
{
type blended;
water
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
oil
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3;
residualSlip 1e-3;
}
(water mercury)
{
type blended;
water
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
mercury
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3;
residualSlip 1e-3;
}
(oil mercury)
{
type blended;
oil
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
mercury
{
type SchillerNaumann;
residualPhaseFraction 0;
residualSlip 0;
}
residualPhaseFraction 1e-3; residualPhaseFraction 1e-3;
residualSlip 1e-3; residualSlip 1e-3;
} }

2
tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/controlDict
View File

@ -17,7 +17,7 @@ FoamFile
application multiphaseInterFoam; application multiphaseInterFoam;
startFrom latestTime; startFrom startTime;
startTime 0; startTime 0;