mapFields: Reinstated mapFields from OpenFOAM-2.2.x and renamed the current mapFields -> mapFieldsPar

This required the addition of the meshToMesh class in the sampling
library from OpenFOAM-2.2.x which is now named meshToMesh0.
This commit is contained in:
Henry
2015-05-26 11:32:46 +01:00
parent 526d7efeb0
commit 71c6af8fe9
45 changed files with 3804 additions and 278 deletions

View File

@ -0,0 +1,130 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 MapConsistentVolFields_H
#define MapConsistentVolFields_H
#include "GeometricField.H"
#include "meshToMesh0.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type, class CombineOp>
void MapConsistentVolFields
(
const IOobjectList& objects,
const meshToMesh0& meshToMesh0Interp,
const meshToMesh0::order& mapOrder,
const CombineOp& cop
)
{
const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
word fieldClassName
(
GeometricField<Type, fvPatchField, volMesh>::typeName
);
IOobjectList fields = objects.lookupClass(fieldClassName);
forAllIter(IOobjectList, fields, fieldIter)
{
Info<< " interpolating " << fieldIter()->name()
<< endl;
// Read field
GeometricField<Type, fvPatchField, volMesh> fieldSource
(
*fieldIter(),
meshSource
);
IOobject fieldTargetIOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
meshTarget,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (fieldTargetIOobject.headerOk())
{
// Read fieldTarget
GeometricField<Type, fvPatchField, volMesh> fieldTarget
(
fieldTargetIOobject,
meshTarget
);
// Interpolate field
meshToMesh0Interp.interpolate
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
// Write field
fieldTarget.write();
}
else
{
fieldTargetIOobject.readOpt() = IOobject::NO_READ;
// Interpolate field
GeometricField<Type, fvPatchField, volMesh> fieldTarget
(
fieldTargetIOobject,
meshToMesh0Interp.interpolate
(
fieldSource,
mapOrder,
cop
)
);
// Write field
fieldTarget.write();
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,7 +35,7 @@ Description
#include "cloud.H"
#include "GeometricField.H"
#include "meshToMesh.H"
#include "meshToMesh0.H"
#include "IOobjectList.H"
#include "CompactIOField.H"
@ -51,18 +51,19 @@ void MapLagrangianFields
(
const string& cloudName,
const IOobjectList& objects,
const polyMesh& meshTarget,
const meshToMesh0& meshToMesh0Interp,
const labelList& addParticles
)
{
const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
{
IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian field " << fieldName << endl;
Info<< " mapping lagrangian field "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
IOField<Type> fieldSource(*fieldIter());
@ -72,7 +73,7 @@ void MapLagrangianFields
(
IOobject
(
fieldName,
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
@ -99,9 +100,8 @@ void MapLagrangianFields
forAllIter(IOobjectList, fieldFields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian fieldField " << fieldName << endl;
Info<< " mapping lagrangian fieldField "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
IOField<Field<Type> > fieldSource(*fieldIter());
@ -112,7 +112,7 @@ void MapLagrangianFields
(
IOobject
(
fieldName,
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,9 +27,9 @@ License
#define MapMeshes_H
#include "MapVolFields.H"
#include "MapConsistentVolFields.H"
#include "mapLagrangian.H"
#include "UnMapped.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -37,62 +37,64 @@ namespace Foam
{
template<template<class> class CombineOp>
void MapMesh
void MapConsistentMesh
(
const meshToMesh& interp,
const HashSet<word>& selectedFields,
const bool noLagrangian
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh0::order& mapOrder
)
{
{
const polyMesh& meshSource = interp.srcRegion();
// Create the interpolation scheme
meshToMesh0 meshToMesh0Interp(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
// ~~~~~~~~~~~~~
MapVolFields<scalar>
MapConsistentVolFields<scalar>
(
objects,
selectedFields,
interp,
meshToMesh0Interp,
mapOrder,
CombineOp<scalar>()
);
MapVolFields<vector>
MapConsistentVolFields<vector>
(
objects,
selectedFields,
interp,
meshToMesh0Interp,
mapOrder,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
MapConsistentVolFields<sphericalTensor>
(
objects,
selectedFields,
interp,
meshToMesh0Interp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
MapConsistentVolFields<symmTensor>
(
objects,
selectedFields,
interp,
meshToMesh0Interp,
mapOrder,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
MapConsistentVolFields<tensor>
(
objects,
selectedFields,
interp,
meshToMesh0Interp,
mapOrder,
CombineOp<tensor>()
);
}
{
const polyMesh& meshTarget = interp.tgtRegion();
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
@ -113,10 +115,140 @@ void MapMesh
UnMapped<pointTensorField>(objects);
}
if (!noLagrangian)
mapLagrangian(meshToMesh0Interp);
}
template<template<class> class CombineOp>
void MapSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh0::order& mapOrder
)
{
// Create the interpolation scheme
meshToMesh0 meshToMesh0Interp
(
meshSource,
meshTarget,
patchMap,
cuttingPatches
);
Info<< nl
<< "Mapping fields for time " << meshSource.time().timeName()
<< nl << endl;
{
mapLagrangian(interp);
// Search for list of source objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>
(
objects,
meshToMesh0Interp,
mapOrder,
CombineOp<scalar>()
);
MapVolFields<vector>
(
objects,
meshToMesh0Interp,
mapOrder,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
(
objects,
meshToMesh0Interp,
mapOrder,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
(
objects,
meshToMesh0Interp,
mapOrder,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
(
objects,
meshToMesh0Interp,
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(meshToMesh0Interp);
}
template<template<class> class CombineOp>
void MapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh0::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
);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,11 +23,11 @@ License
\*---------------------------------------------------------------------------*/
#ifndef MapConsistentVolFields_H
#define MapConsistentVolFields_H
#ifndef MapVolFields_H
#define MapVolFields_H
#include "GeometricField.H"
#include "meshToMesh.H"
#include "meshToMesh0.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -39,55 +39,62 @@ template<class Type, class CombineOp>
void MapVolFields
(
const IOobjectList& objects,
const HashSet<word>& selectedFields,
const meshToMesh& interp,
const meshToMesh0& meshToMesh0Interp,
const meshToMesh0::order& mapOrder,
const CombineOp& cop
)
{
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
const fvMesh& meshSource = static_cast<const fvMesh&>(interp.srcRegion());
const fvMesh& meshTarget = static_cast<const fvMesh&>(interp.tgtRegion());
word fieldClassName
(
GeometricField<Type, fvPatchField, volMesh>::typeName
);
IOobjectList fields = objects.lookupClass(fieldType::typeName);
IOobjectList fields = objects.lookupClass(fieldClassName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
IOobject fieldTargetIOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
meshTarget,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (selectedFields.empty() || selectedFields.found(fieldName))
if (fieldTargetIOobject.headerOk())
{
Info<< " interpolating " << fieldName << endl;
Info<< " interpolating " << fieldIter()->name()
<< endl;
const fieldType fieldSource(*fieldIter(), meshSource);
IOobject targetIO
// Read field fieldSource
GeometricField<Type, fvPatchField, volMesh> fieldSource
(
fieldName,
meshTarget.time().timeName(),
meshTarget,
IOobject::MUST_READ
*fieldIter(),
meshSource
);
if (targetIO.headerOk())
{
fieldType fieldTarget(targetIO, meshTarget);
// Read fieldTarget
GeometricField<Type, fvPatchField, volMesh> fieldTarget
(
fieldTargetIOobject,
meshTarget
);
interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
// Interpolate field
meshToMesh0Interp.interpolate
(
fieldTarget,
fieldSource,
mapOrder,
cop
);
fieldTarget.write();
}
else
{
targetIO.readOpt() = IOobject::NO_READ;
tmp<fieldType>
tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
fieldType fieldTarget(targetIO, tfieldTarget);
fieldTarget.write();
}
// Write field
fieldTarget.write();
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -1,11 +1,15 @@
Info<< "\nCreate databases as time" << endl;
HashTable<string> srcOptions(args.options());
srcOptions.erase("case");
srcOptions.insert("case", fileName(rootDirSource/caseDirSource));
Time runTimeSource
(
Time::controlDictName,
rootDirSource,
caseDirSource
);
argList argsSrc(args, srcOptions, false, false, false);
Time runTimeSource(Time::controlDictName, argsSrc);
Time runTimeTarget(Time::controlDictName, args);
Time runTimeTarget
(
Time::controlDictName,
rootDirTarget,
caseDirTarget
);

View File

@ -27,12 +27,14 @@ Application
Description
Maps volume fields from one mesh to another, reading and
interpolating all fields present in the time directory of both cases.
Parallel and non-parallel cases are handled without the need to reconstruct
them first.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "meshToMesh.H"
#include "processorPolyPatch.H"
#include "meshToMesh0.H"
#include "processorFvPatch.H"
#include "MapMeshes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -41,33 +43,26 @@ void mapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
const meshToMesh0::order& mapOrder,
const bool subtract
)
{
Info<< nl << "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMesh interp(meshSource, meshTarget, mapMethod);
if (subtract)
{
MapMesh<minusEqOp>
MapConsistentMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
meshSource,
meshTarget,
mapOrder
);
}
else
{
MapMesh<plusEqOp>
MapConsistentMesh<eqOp>
(
interp,
selectedFields,
noLagrangian
meshSource,
meshTarget,
mapOrder
);
}
}
@ -79,40 +74,59 @@ void mapSubMesh
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
const meshToMesh0::order& mapOrder,
const bool subtract
)
{
Info<< nl << "Creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMesh interp
(
meshSource,
meshTarget,
mapMethod,
patchMap,
cuttingPatches
);
if (subtract)
{
MapMesh<minusEqOp>
MapSubMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
else
{
MapMesh<plusEqOp>
MapSubMesh<eqOp>
(
interp,
selectedFields,
noLagrangian
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder
);
}
}
void mapConsistentSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh0::order& mapOrder,
const bool subtract
)
{
if (subtract)
{
MapConsistentSubMesh<minusEqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
else
{
MapConsistentSubMesh<eqOp>
(
meshSource,
meshTarget,
mapOrder
);
}
}
@ -125,20 +139,30 @@ wordList addProcessorPatches
)
{
// Add the processor patches to the cutting list
HashSet<word> cuttingPatchTable;
HashTable<label> cuttingPatchTable;
forAll(cuttingPatches, i)
{
cuttingPatchTable.insert(cuttingPatches[i]);
cuttingPatchTable.insert(cuttingPatches[i], i);
}
const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();
forAll(pbm, patchI)
forAll(meshTarget.boundary(), patchi)
{
if (isA<processorPolyPatch>(pbm[patchI]))
if (isA<processorFvPatch>(meshTarget.boundary()[patchi]))
{
const word& patchName = pbm[patchI].name();
cuttingPatchTable.insert(patchName);
if
(
!cuttingPatchTable.found
(
meshTarget.boundaryMesh()[patchi].name()
)
)
{
cuttingPatchTable.insert
(
meshTarget.boundaryMesh()[patchi].name(),
-1
);
}
}
}
@ -154,7 +178,7 @@ int main(int argc, char *argv[])
(
"map volume fields from one mesh to another"
);
argList::noParallel();
argList::validArgs.append("sourceCase");
argList::addOption
@ -176,6 +200,16 @@ int main(int argc, char *argv[])
"specify the target region"
);
argList::addBoolOption
(
"parallelSource",
"the source is decomposed"
);
argList::addBoolOption
(
"parallelTarget",
"the target is decomposed"
);
argList::addBoolOption
(
"consistent",
"source and target geometry and boundary conditions identical"
@ -191,21 +225,14 @@ int main(int argc, char *argv[])
"subtract",
"subtract mapped source from target"
);
argList::addOption
(
"fields",
"list",
"specify a list of fields to be mapped. Eg, '(U T p)' - "
"regular expressions not currently supported"
);
argList::addBoolOption
(
"noLagrangian",
"skip mapping lagrangian positions and fields"
);
argList args(argc, argv);
if (!args.check())
{
FatalError.exit();
}
fileName rootDirTarget(args.rootPath());
fileName caseDirTarget(args.globalCaseName());
@ -229,17 +256,35 @@ int main(int argc, char *argv[])
Info<< "Target region: " << targetRegion << endl;
}
const bool parallelSource = args.optionFound("parallelSource");
const bool parallelTarget = args.optionFound("parallelTarget");
const bool consistent = args.optionFound("consistent");
meshToMesh::interpolationMethod mapMethod =
meshToMesh::imCellVolumeWeight;
meshToMesh0::order mapOrder = meshToMesh0::INTERPOLATE;
if (args.optionFound("mapMethod"))
{
mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]];
const word mapMethod(args["mapMethod"]);
if (mapMethod == "mapNearest")
{
mapOrder = meshToMesh0::MAP;
}
else if (mapMethod == "interpolate")
{
mapOrder = meshToMesh0::INTERPOLATE;
}
else if (mapMethod == "cellPointInterpolate")
{
mapOrder = meshToMesh0::CELL_POINT_INTERPOLATE;
}
else
{
FatalErrorIn(args.executable())
<< "Unknown mapMethod " << mapMethod << ". Valid options are: "
<< "mapNearest, interpolate and cellPointInterpolate"
<< exit(FatalError);
}
Info<< "Mapping method: "
<< meshToMesh::interpolationMethodNames_[mapMethod] << endl;
Info<< "Mapping method: " << mapMethod << endl;
}
const bool subtract = args.optionFound("subtract");
@ -248,13 +293,6 @@ int main(int argc, char *argv[])
Info<< "Subtracting mapped source field from target" << endl;
}
HashSet<word> selectedFields;
if (args.optionFound("fields"))
{
args.optionLookup("fields")() >> selectedFields;
}
const bool noLagrangian = args.optionFound("noLagrangian");
#include "createTimes.H"
@ -280,58 +318,341 @@ int main(int argc, char *argv[])
mapFieldsDict.lookup("cuttingPatches") >> cuttingPatches;
}
#include "setTimeIndex.H"
Info<< "\nCreate meshes\n" << endl;
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "Source mesh size: " << meshSource.nCells() << tab
<< "Target mesh size: " << meshTarget.nCells() << nl << endl;
if (consistent)
if (parallelSource && !parallelTarget)
{
mapConsistentMesh
IOdictionary decompositionDict
(
meshSource,
meshTarget,
mapMethod,
subtract,
selectedFields,
noLagrangian
IOobject
(
"decomposeParDict",
runTimeSource.system(),
runTimeSource,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
int nProcs(readInt(decompositionDict.lookup("numberOfSubdomains")));
Info<< "Create target mesh\n" << endl;
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "Target mesh size: " << meshTarget.nCells() << endl;
for (int procI=0; procI<nProcs; procI++)
{
Info<< nl << "Source processor " << procI << endl;
Time runTimeSource
(
Time::controlDictName,
rootDirSource,
caseDirSource/fileName(word("processor") + name(procI))
);
#include "setTimeIndex.H"
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
Info<< "mesh size: " << meshSource.nCells() << endl;
if (consistent)
{
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder,
subtract
);
}
}
}
else if (!parallelSource && parallelTarget)
{
IOdictionary decompositionDict
(
IOobject
(
"decomposeParDict",
runTimeTarget.system(),
runTimeTarget,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
int nProcs(readInt(decompositionDict.lookup("numberOfSubdomains")));
Info<< "Create source mesh\n" << endl;
#include "setTimeIndex.H"
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
Info<< "Source mesh size: " << meshSource.nCells() << endl;
for (int procI=0; procI<nProcs; procI++)
{
Info<< nl << "Target processor " << procI << endl;
Time runTimeTarget
(
Time::controlDictName,
rootDirTarget,
caseDirTarget/fileName(word("processor") + name(procI))
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "mesh size: " << meshTarget.nCells() << endl;
if (consistent)
{
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder,
subtract
);
}
}
}
else if (parallelSource && parallelTarget)
{
IOdictionary decompositionDictSource
(
IOobject
(
"decomposeParDict",
runTimeSource.system(),
runTimeSource,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
int nProcsSource
(
readInt(decompositionDictSource.lookup("numberOfSubdomains"))
);
IOdictionary decompositionDictTarget
(
IOobject
(
"decomposeParDict",
runTimeTarget.system(),
runTimeTarget,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
int nProcsTarget
(
readInt(decompositionDictTarget.lookup("numberOfSubdomains"))
);
List<boundBox> bbsTarget(nProcsTarget);
List<bool> bbsTargetSet(nProcsTarget, false);
for (int procISource=0; procISource<nProcsSource; procISource++)
{
Info<< nl << "Source processor " << procISource << endl;
Time runTimeSource
(
Time::controlDictName,
rootDirSource,
caseDirSource/fileName(word("processor") + name(procISource))
);
#include "setTimeIndex.H"
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
Info<< "mesh size: " << meshSource.nCells() << endl;
boundBox bbSource(meshSource.bounds());
for (int procITarget=0; procITarget<nProcsTarget; procITarget++)
{
if
(
!bbsTargetSet[procITarget]
|| (
bbsTargetSet[procITarget]
&& bbsTarget[procITarget].overlaps(bbSource)
)
)
{
Info<< nl << "Target processor " << procITarget << endl;
Time runTimeTarget
(
Time::controlDictName,
rootDirTarget,
caseDirTarget/fileName(word("processor")
+ name(procITarget))
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "mesh size: " << meshTarget.nCells() << endl;
bbsTarget[procITarget] = meshTarget.bounds();
bbsTargetSet[procITarget] = true;
if (bbsTarget[procITarget].overlaps(bbSource))
{
if (consistent)
{
mapConsistentSubMesh
(
meshSource,
meshTarget,
mapOrder,
subtract
);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapOrder,
subtract
);
}
}
}
}
}
}
else
{
mapSubMesh
#include "setTimeIndex.H"
Info<< "Create meshes\n" << endl;
fvMesh meshSource
(
meshSource,
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapMethod,
subtract,
selectedFields,
noLagrangian
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "Source mesh size: " << meshSource.nCells() << tab
<< "Target mesh size: " << meshTarget.nCells() << nl << endl;
if (consistent)
{
mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
cuttingPatches,
mapOrder,
subtract
);
}
}
Info<< "\nEnd\n" << endl;

View File

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / O peration | Version: 2.2.2 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
@ -14,12 +14,6 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Specify how to map patches. There are three different options:
// - patch exists in the source case: specify mapping (patchMap)
// - patch should be interpolated from internal values in source case
// (cuttingPatches)
// - patch should not be mapped. Default if not in patchMap or cuttingPatches
// List of pairs of target/source patches for mapping
patchMap
(

View File

@ -80,17 +80,27 @@ static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
}
void mapLagrangian(const meshToMesh& interp)
void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
{
// Determine which particles are in meshTarget
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const polyMesh& meshSource = interp.srcRegion();
const polyMesh& meshTarget = interp.tgtRegion();
const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
// target to source cell map
const labelList& cellAddressing = meshToMesh0Interp.cellAddressing();
// Invert celladdressing to get source to target(s).
// Note: could use sparse addressing but that is too storage inefficient
// (Map<labelList>)
labelListList sourceToTargets
(
invertOneToMany(meshToMesh0Interp.fromMesh().nCells(), cellAddressing)
);
const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs
(
readDir
@ -110,7 +120,7 @@ void mapLagrangian(const meshToMesh& interp)
cloud::prefix/cloudDirs[cloudI]
);
IOobject* positionsPtr = objects.lookup(word("positions"));
IOobject* positionsPtr = objects.lookup("positions");
if (positionsPtr)
{
@ -158,7 +168,7 @@ void mapLagrangian(const meshToMesh& interp)
if (iter().cell() >= 0)
{
const labelList& targetCells =
sourceToTarget[iter().cell()];
sourceToTargets[iter().cell()];
// Particle probably in one of the targetcells. Try
// all by tracking from their cell centre to the parcel
@ -200,7 +210,7 @@ void mapLagrangian(const meshToMesh& interp)
sourceParticleI++;
}
Info<< " after meshToMesh addressing found "
Info<< " after meshToMesh0 addressing found "
<< targetParcels.size()
<< " parcels in target mesh." << endl;
@ -249,47 +259,17 @@ void mapLagrangian(const meshToMesh& interp)
// ~~~~~~~~~~~~~~~~~~~~~
MapLagrangianFields<label>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
MapLagrangianFields<scalar>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
MapLagrangianFields<vector>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
MapLagrangianFields<sphericalTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
MapLagrangianFields<symmTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
MapLagrangianFields<tensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
(cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,7 +35,7 @@ SourceFiles
#ifndef mapLagrangian_H
#define mapLagrangian_H
#include "meshToMesh.H"
#include "meshToMesh0.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -43,7 +43,7 @@ namespace Foam
{
//- Maps lagrangian positions and fields
void mapLagrangian(const meshToMesh& interp);
void mapLagrangian(const meshToMesh0& meshToMesh0Interp);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -1,4 +1,3 @@
{
instantList sourceTimes = runTimeSource.times();
label sourceTimeIndex = runTimeSource.timeIndex();
if (args.optionFound("sourceTime"))
@ -30,4 +29,3 @@
Info<< "\nSource time: " << runTimeSource.value()
<< "\nTarget time: " << runTimeTarget.value()
<< endl;
}

View File

@ -0,0 +1,4 @@
mapLagrangian.C
mapFieldsPar.C
EXE = $(FOAM_APPBIN)/mapFieldsPar

View File

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

View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
InNamespace
Foam
Description
Gets the indices of (source)particles that have been appended to the
target cloud and maps the lagrangian fields accordingly.
\*---------------------------------------------------------------------------*/
#ifndef MapLagrangianFields_H
#define MapLagrangianFields_H
#include "cloud.H"
#include "GeometricField.H"
#include "meshToMesh.H"
#include "IOobjectList.H"
#include "CompactIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Gets the indices of (source)particles that have been appended to the
// target cloud and maps the lagrangian fields accordingly.
template<class Type>
void MapLagrangianFields
(
const string& cloudName,
const IOobjectList& objects,
const polyMesh& meshTarget,
const labelList& addParticles
)
{
{
IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian field " << fieldName << endl;
// Read field (does not need mesh)
IOField<Type> fieldSource(*fieldIter());
// Map
IOField<Type> fieldTarget
(
IOobject
(
fieldName,
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(IOField<Field<Type> >::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian fieldField " << fieldName << endl;
// Read field (does not need mesh)
IOField<Field<Type> > fieldSource(*fieldIter());
// Map - use CompactIOField to automatically write in
// compact form for binary format.
CompactIOField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldName,
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(CompactIOField<Field<Type>, Type>::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
Info<< " mapping lagrangian fieldField "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
CompactIOField<Field<Type>, Type> fieldSource(*fieldIter());
// Map
CompactIOField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 "mapLagrangian.H"
#include "UnMapped.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<template<class> class CombineOp>
void MapMesh
(
const meshToMesh& interp,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
{
const polyMesh& meshSource = interp.srcRegion();
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>
(
objects,
selectedFields,
interp,
CombineOp<scalar>()
);
MapVolFields<vector>
(
objects,
selectedFields,
interp,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
(
objects,
selectedFields,
interp,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
(
objects,
selectedFields,
interp,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
(
objects,
selectedFields,
interp,
CombineOp<tensor>()
);
}
{
const polyMesh& meshTarget = interp.tgtRegion();
// 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);
}
if (!noLagrangian)
{
mapLagrangian(interp);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,104 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 MapConsistentVolFields_H
#define MapConsistentVolFields_H
#include "GeometricField.H"
#include "meshToMesh.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type, class CombineOp>
void MapVolFields
(
const IOobjectList& objects,
const HashSet<word>& selectedFields,
const meshToMesh& interp,
const CombineOp& cop
)
{
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const fvMesh& meshSource = static_cast<const fvMesh&>(interp.srcRegion());
const fvMesh& meshTarget = static_cast<const fvMesh&>(interp.tgtRegion());
IOobjectList fields = objects.lookupClass(fieldType::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
if (selectedFields.empty() || selectedFields.found(fieldName))
{
Info<< " interpolating " << fieldName << endl;
const fieldType fieldSource(*fieldIter(), meshSource);
IOobject targetIO
(
fieldName,
meshTarget.time().timeName(),
meshTarget,
IOobject::MUST_READ
);
if (targetIO.headerOk())
{
fieldType fieldTarget(targetIO, meshTarget);
interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
fieldTarget.write();
}
else
{
targetIO.readOpt() = IOobject::NO_READ;
tmp<fieldType>
tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
fieldType fieldTarget(targetIO, tfieldTarget);
fieldTarget.write();
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 UnMapped_H
#define UnMapped_H
#include "IOobjectList.H"
#include "OSspecific.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type>
void UnMapped(const IOobjectList& objects)
{
IOobjectList fields = objects.lookupClass(Type::typeName);
forAllConstIter(IOobjectList, fields, fieldIter)
{
mvBak(fieldIter()->objectPath(), "unmapped");
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,11 @@
Info<< "\nCreate databases as time" << endl;
HashTable<string> srcOptions(args.options());
srcOptions.erase("case");
srcOptions.insert("case", fileName(rootDirSource/caseDirSource));
argList argsSrc(args, srcOptions, false, false, false);
Time runTimeSource(Time::controlDictName, argsSrc);
Time runTimeTarget(Time::controlDictName, args);

View File

@ -0,0 +1,36 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object mapFieldsParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Specify how to map patches. There are three different options:
// - patch exists in the source case: specify mapping (patchMap)
// - patch should be interpolated from internal values in source case
// (cuttingPatches)
// - patch should not be mapped. Default if not in patchMap or cuttingPatches
// List of pairs of target/source patches for mapping
patchMap
(
lid movingWall
);
// List of target patches cutting the source domain (these need to be
// handled specially e.g. interpolated from internal values)
cuttingPatches
(
fixedWalls
);
// ************************************************************************* //

View File

@ -0,0 +1,343 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
Application
mapFieldsPar
Description
Maps volume fields from one mesh to another, reading and
interpolating all fields present in the time directory of both cases.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "meshToMesh.H"
#include "processorPolyPatch.H"
#include "MapMeshes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void mapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMesh::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
Info<< nl << "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMesh interp(meshSource, meshTarget, mapMethod);
if (subtract)
{
MapMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
else
{
MapMesh<plusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
}
void mapSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMesh::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
Info<< nl << "Creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMesh interp
(
meshSource,
meshTarget,
mapMethod,
patchMap,
cuttingPatches
);
if (subtract)
{
MapMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
else
{
MapMesh<plusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
}
wordList addProcessorPatches
(
const fvMesh& meshTarget,
const wordList& cuttingPatches
)
{
// Add the processor patches to the cutting list
HashSet<word> cuttingPatchTable;
forAll(cuttingPatches, i)
{
cuttingPatchTable.insert(cuttingPatches[i]);
}
const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();
forAll(pbm, patchI)
{
if (isA<processorPolyPatch>(pbm[patchI]))
{
const word& patchName = pbm[patchI].name();
cuttingPatchTable.insert(patchName);
}
}
return cuttingPatchTable.toc();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"map volume fields from one mesh to another"
);
argList::validArgs.append("sourceCase");
argList::addOption
(
"sourceTime",
"scalar|'latestTime'",
"specify the source time"
);
argList::addOption
(
"sourceRegion",
"word",
"specify the source region"
);
argList::addOption
(
"targetRegion",
"word",
"specify the target region"
);
argList::addBoolOption
(
"consistent",
"source and target geometry and boundary conditions identical"
);
argList::addOption
(
"mapMethod",
"word",
"specify the mapping method"
);
argList::addBoolOption
(
"subtract",
"subtract mapped source from target"
);
argList::addOption
(
"fields",
"list",
"specify a list of fields to be mapped. Eg, '(U T p)' - "
"regular expressions not currently supported"
);
argList::addBoolOption
(
"noLagrangian",
"skip mapping lagrangian positions and fields"
);
argList args(argc, argv);
fileName rootDirTarget(args.rootPath());
fileName caseDirTarget(args.globalCaseName());
const fileName casePath = args[1];
const fileName rootDirSource = casePath.path();
const fileName caseDirSource = casePath.name();
Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;
word sourceRegion = fvMesh::defaultRegion;
if (args.optionFound("sourceRegion"))
{
sourceRegion = args["sourceRegion"];
Info<< "Source region: " << sourceRegion << endl;
}
Info<< "Target: " << rootDirTarget << " " << caseDirTarget << endl;
word targetRegion = fvMesh::defaultRegion;
if (args.optionFound("targetRegion"))
{
targetRegion = args["targetRegion"];
Info<< "Target region: " << targetRegion << endl;
}
const bool consistent = args.optionFound("consistent");
meshToMesh::interpolationMethod mapMethod =
meshToMesh::imCellVolumeWeight;
if (args.optionFound("mapMethod"))
{
mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]];
Info<< "Mapping method: "
<< meshToMesh::interpolationMethodNames_[mapMethod] << endl;
}
const bool subtract = args.optionFound("subtract");
if (subtract)
{
Info<< "Subtracting mapped source field from target" << endl;
}
HashSet<word> selectedFields;
if (args.optionFound("fields"))
{
args.optionLookup("fields")() >> selectedFields;
}
const bool noLagrangian = args.optionFound("noLagrangian");
#include "createTimes.H"
HashTable<word> patchMap;
wordList cuttingPatches;
if (!consistent)
{
IOdictionary mapFieldsParDict
(
IOobject
(
"mapFieldsParDict",
runTimeTarget.system(),
runTimeTarget,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
);
mapFieldsParDict.lookup("patchMap") >> patchMap;
mapFieldsParDict.lookup("cuttingPatches") >> cuttingPatches;
}
#include "setTimeIndex.H"
Info<< "\nCreate meshes\n" << endl;
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "Source mesh size: " << meshSource.nCells() << tab
<< "Target mesh size: " << meshTarget.nCells() << nl << endl;
if (consistent)
{
mapConsistentMesh
(
meshSource,
meshTarget,
mapMethod,
subtract,
selectedFields,
noLagrangian
);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapMethod,
subtract,
selectedFields,
noLagrangian
);
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,303 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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 "MapLagrangianFields.H"
#include "passiveParticleCloud.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
static const scalar perturbFactor = 1e-6;
// Special version of findCell that generates a cell guaranteed to be
// compatible with tracking.
static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
{
label cellI = -1;
label tetFaceI = -1;
label tetPtI = -1;
const polyMesh& mesh = cloud.pMesh();
mesh.findCellFacePt(pt, cellI, tetFaceI, tetPtI);
if (cellI >= 0)
{
return cellI;
}
else
{
// See if particle on face by finding nearest face and shifting
// particle.
meshSearch meshSearcher
(
mesh,
polyMesh::FACE_PLANES // no decomposition needed
);
label faceI = meshSearcher.findNearestBoundaryFace(pt);
if (faceI >= 0)
{
const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
mesh.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI);
return cellI;
}
}
return -1;
}
void mapLagrangian(const meshToMesh& interp)
{
// Determine which particles are in meshTarget
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const polyMesh& meshSource = interp.srcRegion();
const polyMesh& meshTarget = interp.tgtRegion();
const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs
(
readDir
(
meshSource.time().timePath()/cloud::prefix,
fileName::DIRECTORY
)
);
forAll(cloudDirs, cloudI)
{
// Search for list of lagrangian objects for this time
IOobjectList objects
(
meshSource,
meshSource.time().timeName(),
cloud::prefix/cloudDirs[cloudI]
);
IOobject* positionsPtr = objects.lookup(word("positions"));
if (positionsPtr)
{
Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
// Read positions & cell
passiveParticleCloud sourceParcels
(
meshSource,
cloudDirs[cloudI],
false
);
Info<< " read " << sourceParcels.size()
<< " parcels from source mesh." << endl;
// Construct empty target cloud
passiveParticleCloud targetParcels
(
meshTarget,
cloudDirs[cloudI],
IDLList<passiveParticle>()
);
particle::TrackingData<passiveParticleCloud> td(targetParcels);
label sourceParticleI = 0;
// Indices of source particles that get added to targetParcels
DynamicList<label> addParticles(sourceParcels.size());
// Unmapped particles
labelHashSet unmappedSource(sourceParcels.size());
// Initial: track from fine-mesh cell centre to particle position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This requires there to be no boundary in the way.
forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
{
bool foundCell = false;
// Assume that cell from read parcel is the correct one...
if (iter().cell() >= 0)
{
const labelList& targetCells =
sourceToTarget[iter().cell()];
// Particle probably in one of the targetcells. Try
// all by tracking from their cell centre to the parcel
// position.
forAll(targetCells, i)
{
// Track from its cellcentre to position to make sure.
autoPtr<passiveParticle> newPtr
(
new passiveParticle
(
meshTarget,
targetCc[targetCells[i]],
targetCells[i]
)
);
passiveParticle& newP = newPtr();
label faceI = newP.track(iter().position(), td);
if (faceI < 0 && newP.cell() >= 0)
{
// Hit position.
foundCell = true;
addParticles.append(sourceParticleI);
targetParcels.addParticle(newPtr.ptr());
break;
}
}
}
if (!foundCell)
{
// Store for closer analysis
unmappedSource.insert(sourceParticleI);
}
sourceParticleI++;
}
Info<< " after meshToMesh addressing found "
<< targetParcels.size()
<< " parcels in target mesh." << endl;
// Do closer inspection for unmapped particles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (unmappedSource.size())
{
sourceParticleI = 0;
forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
{
if (unmappedSource.found(sourceParticleI))
{
label targetCell =
findCell(targetParcels, iter().position());
if (targetCell >= 0)
{
unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI);
iter().cell() = targetCell;
targetParcels.addParticle
(
sourceParcels.remove(&iter())
);
}
}
sourceParticleI++;
}
}
addParticles.shrink();
Info<< " after additional mesh searching found "
<< targetParcels.size() << " parcels in target mesh." << endl;
if (addParticles.size())
{
IOPosition<passiveParticleCloud>(targetParcels).write();
// addParticles now contains the indices of the sourceMesh
// particles that were appended to the target mesh.
// Map lagrangian fields
// ~~~~~~~~~~~~~~~~~~~~~
MapLagrangianFields<label>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<scalar>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<vector>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<sphericalTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<symmTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<tensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
InNamespace
Foam
Description
Maps lagrangian positions and fields
SourceFiles
mapLagrangian.C
\*---------------------------------------------------------------------------*/
#ifndef mapLagrangian_H
#define mapLagrangian_H
#include "meshToMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Maps lagrangian positions and fields
void mapLagrangian(const meshToMesh& interp);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,33 @@
{
instantList sourceTimes = runTimeSource.times();
label sourceTimeIndex = runTimeSource.timeIndex();
if (args.optionFound("sourceTime"))
{
if (args["sourceTime"] == "latestTime")
{
sourceTimeIndex = sourceTimes.size() - 1;
}
else
{
sourceTimeIndex = Time::findClosestTimeIndex
(
sourceTimes,
args.optionRead<scalar>("sourceTime")
);
}
}
else
{
sourceTimeIndex = Time::findClosestTimeIndex
(
sourceTimes,
runTimeTarget.time().value()
);
}
runTimeSource.setTime(sourceTimes[sourceTimeIndex], sourceTimeIndex);
Info<< "\nSource time: " << runTimeSource.value()
<< "\nTarget time: " << runTimeTarget.value()
<< endl;
}