From 631a47a37810d7d5904e8dc7330126852c845b0f Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 9 Jan 2013 11:50:23 +0000
Subject: [PATCH] ENH: Added new mesh-to-mesh interpolation class
---
src/sampling/Make/files | 3 +
.../meshToMeshNew/meshToMeshNew.C | 955 ++++++++++++++++++
.../meshToMeshNew/meshToMeshNew.H | 494 +++++++++
.../meshToMeshNew/meshToMeshNewI.H | 64 ++
.../meshToMeshNew/meshToMeshNewParallelOps.C | 884 ++++++++++++++++
.../meshToMeshNew/meshToMeshNewTemplates.C | 537 ++++++++++
6 files changed, 2937 insertions(+)
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 5dc8a0ac8fb..747097a781c 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -59,5 +59,8 @@ $(meshToMesh)/meshToMesh.C
$(meshToMesh)/calculateMeshToMeshAddressing.C
$(meshToMesh)/calculateMeshToMeshWeights.C
+meshToMeshNew = meshToMeshInterpolation/meshToMeshNew
+$(meshToMeshNew)/meshToMeshNew.C
+$(meshToMeshNew)/meshToMeshNewParallelOps.C
LIB = $(FOAM_LIBBIN)/libsampling
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
new file mode 100644
index 00000000000..b71e768f428
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
@@ -0,0 +1,955 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H"
+#include "OFstream.H"
+#include "Time.H"
+#include "globalIndex.H"
+#include "mergePoints.H"
+#include "treeBoundBox.H"
+#include "tetOverlapVolume.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(meshToMeshNew, 0);
+}
+
+Foam::scalar Foam::meshToMeshNew::tolerance_ = 1e-6;
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::meshToMeshNew::writeConnectivity
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelListList& srcToTargetAddr
+) const
+{
+ Pout<< "Source size = " << src.nCells() << endl;
+ Pout<< "Target size = " << tgt.nCells() << endl;
+
+ word fName("addressing_" + src.name() + "_to_" + tgt.name());
+
+ if (Pstream::parRun())
+ {
+ fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
+ }
+
+ OFstream os(src.time().path()/fName + ".obj");
+
+ label vertI = 0;
+ forAll(srcToTargetAddr, i)
+ {
+ const labelList& tgtAddress = srcToTargetAddr[i];
+ forAll(tgtAddress, j)
+ {
+ label tgtI = tgtAddress[j];
+ const vector& c0 = src.cellCentres()[i];
+
+ const cell& c = tgt.cells()[tgtI];
+ const pointField pts(c.points(tgt.faces(), tgt.points()));
+ forAll(pts, j)
+ {
+ const point& p = pts[j];
+ os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+ vertI++;
+ os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
+ << nl;
+ vertI++;
+ os << "l " << vertI - 1 << ' ' << vertI << nl;
+ }
+ }
+ }
+}
+
+
+Foam::labelList Foam::meshToMeshNew::maskCells
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+) const
+{
+ boundBox intersectBb
+ (
+ max(src.bounds().min(), tgt.bounds().min()),
+ min(src.bounds().max(), tgt.bounds().max())
+ );
+
+ intersectBb.inflate(0.01);
+
+ const cellList& srcCells = src.cells();
+ const faceList& srcFaces = src.faces();
+ const pointField& srcPts = src.points();
+
+ DynamicList<label> cells(src.size());
+ forAll(srcCells, srcI)
+ {
+ boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
+ if (intersectBb.overlaps(cellBb))
+ {
+ cells.append(srcI);
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "participating source mesh cells: " << cells.size() << endl;
+ }
+
+ return cells;
+}
+
+
+bool Foam::meshToMeshNew::findInitialSeeds
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+) const
+{
+ const cellList& srcCells = src.cells();
+ const faceList& srcFaces = src.faces();
+ const pointField& srcPts = src.points();
+
+ for (label i = startSeedI; i < srcCellIDs.size(); i++)
+ {
+ label srcI = srcCellIDs[i];
+
+ if (mapFlag[srcI])
+ {
+ const pointField
+ pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
+
+ forAll(pts, ptI)
+ {
+ const point& pt = pts[ptI];
+ label tgtI = tgt.cellTree().findInside(pt);
+
+ if (tgtI != -1 && intersect(src, tgt, srcI, tgtI))
+ {
+ srcSeedI = srcI;
+ tgtSeedI = tgtI;
+
+ return true;
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "could not find starting seed" << endl;
+ }
+
+ return false;
+}
+
+
+void Foam::meshToMeshNew::appendToDirectSeeds
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ boolList& mapFlag,
+ labelList& srcTgtSeed,
+ DynamicList<label>& srcSeeds,
+ label& srcSeedI,
+ label& tgtSeedI
+) const
+{
+ const labelList& srcNbr = src.cellCells()[srcSeedI];
+ const labelList& tgtNbr = tgt.cellCells()[tgtSeedI];
+
+ const vectorField& srcCentre = src.cellCentres();
+
+ forAll(srcNbr, i)
+ {
+ label srcI = srcNbr[i];
+
+ if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1))
+ {
+ // source cell srcI not yet mapped
+
+ // identfy if target cell exists for source cell srcI
+ bool found = false;
+ forAll(tgtNbr, j)
+ {
+ label tgtI = tgtNbr[j];
+
+ if (tgt.pointInCell(srcCentre[srcI], tgtI))
+ {
+ // new match - append to lists
+ found = true;
+
+ srcTgtSeed[srcI] = tgtI;
+ srcSeeds.append(srcI);
+
+ break;
+ }
+ }
+
+ if (!found)
+ {
+ // no match available for source cell srcI
+ mapFlag[srcI] = false;
+ }
+ }
+ }
+
+ if (srcSeeds.size())
+ {
+ srcSeedI = srcSeeds.remove();
+ tgtSeedI = srcTgtSeed[srcSeedI];
+ }
+ else
+ {
+ srcSeedI = -1;
+ tgtSeedI = -1;
+ }
+}
+
+
+void Foam::meshToMeshNew::calcDirect
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcSeedI,
+ const label tgtSeedI
+)
+{
+ // store a list of src cells already mapped
+ boolList srcSeedFlag(src.nCells(), true);
+ labelList srcTgtSeed(src.nCells(), -1);
+
+ List<DynamicList<label> > srcToTgt(src.nCells());
+ List<DynamicList<label> > tgtToSrc(tgt.nCells());
+
+ DynamicList<label> srcSeeds;
+
+ const scalarField& srcVc = src.cellVolumes();
+
+ label srcCellI = srcSeedI;
+ label tgtCellI = tgtSeedI;
+
+ do
+ {
+ // store src/tgt cell pair
+ srcToTgt[srcCellI].append(tgtCellI);
+ tgtToSrc[tgtCellI].append(srcCellI);
+
+ // mark source cell srcSeedI as matched
+ srcSeedFlag[srcCellI] = false;
+
+ // accumulate intersection volume
+ V_ += srcVc[srcCellI];
+
+ // find new source seed cell
+ appendToDirectSeeds
+ (
+ src,
+ tgt,
+ srcSeedFlag,
+ srcTgtSeed,
+ srcSeeds,
+ srcCellI,
+ tgtCellI
+ );
+
+ }
+ while (srcCellI >= 0);
+
+ // transfer addressing into persistent storage
+ forAll(srcToTgtCellAddr_, i)
+ {
+ srcToTgtCellAddr_[i].transfer(srcToTgt[i]);
+ srcToTgtCellWght_[i] = scalarList(srcToTgtCellAddr_[i].size(), 1.0);
+ }
+
+ forAll(tgtToSrcCellAddr_, i)
+ {
+ tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]);
+ tgtToSrcCellWght_[i] = scalarList(tgtToSrcCellAddr_[i].size(), 1.0);
+ }
+}
+
+
+void Foam::meshToMeshNew::normaliseWeights
+(
+ const word& descriptor,
+ const scalarField& cellVolumes,
+ const labelListList& addr,
+ scalarListList& wght
+) const
+{
+ const label nCell = returnReduce(wght.size(), sumOp<label>());
+
+ if (nCell > 0)
+ {
+ scalar minW = GREAT;
+ scalar maxW = -GREAT;
+
+ forAll(wght, cellI)
+ {
+ scalarList& w = wght[cellI];
+ scalar s = sum(w);
+ scalar Vc = cellVolumes[cellI];
+
+ forAll(w, i)
+ {
+ w[i] /= Vc;
+ }
+
+ minW = min(minW, s/Vc);
+ maxW = max(maxW, s/Vc);
+ }
+
+ Info<< type() << ": " << descriptor << " weights min/max = "
+ << returnReduce(minW, minOp<scalar>()) << ", "
+ << returnReduce(maxW, maxOp<scalar>()) << endl;
+ }
+}
+
+
+void Foam::meshToMeshNew::appendNbrTgtCells
+(
+ const label tgtCellI,
+ const polyMesh& tgt,
+ const DynamicList<label>& visitedTgtCells,
+ DynamicList<label>& nbrTgtCellIDs
+) const
+{
+ const labelList& nbrCells = tgt.cellCells()[tgtCellI];
+
+ // filter out cells already visited from cell neighbours
+ forAll(nbrCells, i)
+ {
+ label nbrCellI = nbrCells[i];
+
+ if
+ (
+ (findIndex(visitedTgtCells, nbrCellI) == -1)
+ && (findIndex(nbrTgtCellIDs, nbrCellI) == -1)
+ )
+ {
+ nbrTgtCellIDs.append(nbrCellI);
+ }
+ }
+}
+
+
+void Foam::meshToMeshNew::setNextCells
+(
+ label& startSeedI,
+ label& srcCellI,
+ label& tgtCellI,
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const DynamicList<label>& visitedCells,
+ labelList& seedCells
+) const
+{
+ const labelList& srcNbrCells = src.cellCells()[srcCellI];
+
+ // set possible seeds for later use by querying all src cell neighbours
+ // with all visited target cells
+ bool valuesSet = false;
+ forAll(srcNbrCells, i)
+ {
+ label cellS = srcNbrCells[i];
+
+ if (mapFlag[cellS] && seedCells[cellS] == -1)
+ {
+ forAll(visitedCells, j)
+ {
+ label cellT = visitedCells[j];
+
+ if (intersect(src, tgt, cellS, cellT))
+ {
+ seedCells[cellS] = cellT;
+
+ if (!valuesSet)
+ {
+ srcCellI = cellS;
+ tgtCellI = cellT;
+ valuesSet = true;
+ }
+ }
+ }
+ }
+ }
+
+ // set next src and tgt cells if not set above
+ if (valuesSet)
+ {
+ return;
+ }
+ else
+ {
+ // try to use existing seed
+ bool foundNextSeed = false;
+ for (label i = startSeedI; i < srcCellIDs.size(); i++)
+ {
+ label cellS = srcCellIDs[i];
+
+ if (mapFlag[cellS])
+ {
+ if (!foundNextSeed)
+ {
+ startSeedI = i;
+ foundNextSeed = true;
+ }
+
+ if (seedCells[cellS] != -1)
+ {
+ srcCellI = cellS;
+ tgtCellI = seedCells[cellS];
+
+ return;
+ }
+ }
+ }
+
+ // perform new search to find match
+ if (debug)
+ {
+ Pout<< "Advancing front stalled: searching for new "
+ << "target cell" << endl;
+ }
+
+ bool restart =
+ findInitialSeeds
+ (
+ src,
+ tgt,
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcCellI,
+ tgtCellI
+ );
+
+ if (restart)
+ {
+ // successfully found new starting seed-pair
+ return;
+ }
+ }
+
+ // if we have got to here, there are no more src/tgt cell intersections
+ srcCellI = -1;
+ tgtCellI = -1;
+}
+
+
+bool Foam::meshToMeshNew::intersect
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcCellI,
+ const label tgtCellI
+) const
+{
+ scalar threshold = tolerance_*src.cellVolumes()[srcCellI];
+
+ tetOverlapVolume overlapEngine;
+
+ treeBoundBox bbTgtCell
+ (
+ pointField
+ (
+ tgt.points(),
+ tgt.cellPoints()[tgtCellI]
+ )
+ );
+
+ return overlapEngine.cellCellOverlapMinDecomp
+ (
+ src,
+ srcCellI,
+ tgt,
+ tgtCellI,
+ bbTgtCell,
+ threshold
+ );
+}
+
+
+Foam::scalar Foam::meshToMeshNew::interVol
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcCellI,
+ const label tgtCellI
+) const
+{
+ tetOverlapVolume overlapEngine;
+
+ treeBoundBox bbTgtCell
+ (
+ pointField
+ (
+ tgt.points(),
+ tgt.cellPoints()[tgtCellI]
+ )
+ );
+
+ scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
+ (
+ src,
+ srcCellI,
+ tgt,
+ tgtCellI,
+ bbTgtCell
+ );
+
+ return vol;
+}
+
+
+void Foam::meshToMeshNew::calcIndirect
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs,
+ boolList& mapFlag,
+ label& startSeedI
+)
+{
+ label srcCellI = srcSeedI;
+ label tgtCellI = tgtSeedI;
+
+ List<DynamicList<label> > srcToTgtAddr(src.nCells());
+ List<DynamicList<scalar> > srcToTgtWght(src.nCells());
+
+ List<DynamicList<label> > tgtToSrcAddr(tgt.nCells());
+ List<DynamicList<scalar> > tgtToSrcWght(tgt.nCells());
+
+ // list of tgt cell neighbour cells
+ DynamicList<label> nbrTgtCells(10);
+
+ // list of tgt cells currently visited for srcCellI to avoid multiple hits
+ DynamicList<label> visitedTgtCells(10);
+
+ // list to keep track of tgt cells used to seed src cells
+ labelList seedCells(src.nCells(), -1);
+ seedCells[srcCellI] = tgtCellI;
+
+ const scalarField& srcVol = src.cellVolumes();
+
+ do
+ {
+ nbrTgtCells.clear();
+ visitedTgtCells.clear();
+
+ // append initial target cell and neighbours
+ nbrTgtCells.append(tgtCellI);
+ appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+
+ do
+ {
+ tgtCellI = nbrTgtCells.remove();
+ visitedTgtCells.append(tgtCellI);
+
+ scalar vol = interVol(src, tgt, srcCellI, tgtCellI);
+
+ // accumulate addressing and weights for valid intersection
+ if (vol/srcVol[srcCellI] > tolerance_)
+ {
+ // store src/tgt cell pair
+ srcToTgtAddr[srcCellI].append(tgtCellI);
+ srcToTgtWght[srcCellI].append(vol);
+
+ tgtToSrcAddr[tgtCellI].append(srcCellI);
+ tgtToSrcWght[tgtCellI].append(vol);
+
+ appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+
+ // accumulate intersection volume
+ V_ += vol;
+ }
+ }
+ while (!nbrTgtCells.empty());
+
+ mapFlag[srcCellI] = false;
+
+ // find new source seed cell
+ setNextCells
+ (
+ startSeedI,
+ srcCellI,
+ tgtCellI,
+ src,
+ tgt,
+ srcCellIDs,
+ mapFlag,
+ visitedTgtCells,
+ seedCells
+ );
+ }
+ while (srcCellI != -1);
+
+ // transfer addressing into persistent storage
+ forAll(srcToTgtCellAddr_, i)
+ {
+ srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]);
+ srcToTgtCellWght_[i].transfer(srcToTgtWght[i]);
+ }
+
+ forAll(tgtToSrcCellAddr_, i)
+ {
+ tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]);
+ tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]);
+ }
+}
+
+
+void Foam::meshToMeshNew::calcAddressing
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+)
+{
+ srcToTgtCellAddr_.setSize(src.nCells());
+ srcToTgtCellWght_.setSize(src.nCells());
+
+ tgtToSrcCellAddr_.setSize(tgt.nCells());
+ tgtToSrcCellWght_.setSize(tgt.nCells());
+
+ if (!src.nCells() || !tgt.nCells())
+ {
+ if (debug)
+ {
+ Pout<< "mesh interpolation: cells not on processor: Source cells = "
+ << src.nCells() << ", target cells = " << tgt.nCells()
+ << endl;
+ }
+ }
+
+ if (!src.nCells())
+ {
+ return;
+ }
+ else if (!tgt.nCells())
+ {
+ if (debug)
+ {
+ Pout<< "mesh interpolation: hhave " << src.nCells() << " source "
+ << " cells but no target cells" << endl;
+ }
+
+ return;
+ }
+
+ // (potentially) participating source mesh cells
+ const labelList srcCellIDs = maskCells(src, tgt);
+
+ // list to keep track of whether src cell can be mapped
+ boolList mapFlag(src.nCells(), false);
+ UIndirectList<bool>(mapFlag, srcCellIDs) = true;
+
+ // find initial point in tgt mesh
+ label srcSeedI = -1;
+ label tgtSeedI = -1;
+ label startSeedI = 0;
+
+ bool startWalk =
+ findInitialSeeds
+ (
+ src,
+ tgt,
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcSeedI,
+ tgtSeedI
+ );
+
+ if (!startWalk)
+ {
+ // if meshes are collocated, after inflating the source mesh bounding
+ // box tgt mesh cells may be transferred, but may still not overlap
+ // with the source mesh
+ return;
+ }
+
+
+ if (directMapping_)
+ {
+ calcDirect(src, tgt, srcSeedI, tgtSeedI);
+ }
+ else
+ {
+ calcIndirect
+ (
+ src,
+ tgt,
+ srcSeedI,
+ tgtSeedI,
+ srcCellIDs,
+ mapFlag,
+ startSeedI
+ );
+ }
+
+
+ if (debug)
+ {
+ writeConnectivity(src, tgt, srcToTgtCellAddr_);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::meshToMeshNew::meshToMeshNew
+(
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const bool directMapping
+)
+:
+ srcRegionName_(src.name()),
+ tgtRegionName_(tgt.name()),
+ srcToTgtCellAddr_(),
+ tgtToSrcCellAddr_(),
+ srcToTgtCellWght_(),
+ tgtToSrcCellWght_(),
+ V_(0.0),
+ singleMeshProc_(-1),
+ srcMapPtr_(NULL),
+ tgtMapPtr_(NULL),
+ directMapping_(directMapping)
+{
+ Info<< "Creating mesh-to-mesh addressing for " << src.name()
+ << " and " << tgt.name() << " regions" << endl;
+
+ singleMeshProc_ = calcDistribution(src, tgt);
+
+ if (singleMeshProc_ == -1)
+ {
+ // create global indexing for src and tgt meshes
+ globalIndex globalSrcCells(src.nCells());
+ globalIndex globalTgtCells(tgt.nCells());
+
+ // Create processor map of overlapping cells. This map gets
+ // (possibly remote) cells from the tgt mesh such that they (together)
+ // cover all of the src mesh
+ autoPtr<mapDistribute> mapPtr = calcProcMap(src, tgt);
+ const mapDistribute& map = mapPtr();
+
+ pointField newTgtPoints;
+ faceList newTgtFaces;
+ labelList newTgtFaceOwners;
+ labelList newTgtFaceNeighbours;
+ labelList newTgtCellIDs;
+
+ distributeAndMergeCells
+ (
+ map,
+ tgt,
+ globalTgtCells,
+ newTgtPoints,
+ newTgtFaces,
+ newTgtFaceOwners,
+ newTgtFaceNeighbours,
+ newTgtCellIDs
+ );
+
+
+ // create a new target mesh
+ polyMesh newTgt
+ (
+ IOobject
+ (
+ "newTgt::" + Foam::name(Pstream::myProcNo()),
+ tgt.time().timeName(),
+ tgt.time(),
+ IOobject::NO_READ
+ ),
+ xferMove(newTgtPoints),
+ xferMove(newTgtFaces),
+ xferMove(newTgtFaceOwners),
+ xferMove(newTgtFaceNeighbours),
+ false // no parallel comms
+ );
+
+ // create some dummy patch info
+ List<polyPatch*> patches(1);
+ patches[0] = new polyPatch
+ (
+ "defaultFaces",
+ newTgt.nFaces() - newTgt.nInternalFaces(),
+ newTgt.nInternalFaces(),
+ 0,
+ newTgt.boundaryMesh(),
+ word::null
+ );
+
+ newTgt.addPatches(patches);
+
+ // force calculation of tet-base points used for point-in-cell
+ (void)newTgt.tetBasePtIs();
+
+ // force construction of cell tree
+// (void)newTgt.cellTree();
+
+ if (debug)
+ {
+ Pout<< "Created newTgt mesh:" << nl
+ << " old cells = " << tgt.nCells()
+ << ", new cells = " << newTgt.nCells() << nl
+ << " old faces = " << tgt.nFaces()
+ << ", new faces = " << newTgt.nFaces() << endl;
+
+ if (debug > 1)
+ {
+ Pout<< "Writing newTgt mesh: " << newTgt.name() << endl;
+ newTgt.write();
+ }
+ }
+
+ calcAddressing(src, newTgt);
+
+ // per source cell the target cell address in newTgt mesh
+ forAll(srcToTgtCellAddr_, i)
+ {
+ labelList& addressing = srcToTgtCellAddr_[i];
+ forAll(addressing, addrI)
+ {
+ addressing[addrI] = newTgtCellIDs[addressing[addrI]];
+ }
+ }
+
+ // convert target addresses in newTgtMesh into global cell numbering
+ forAll(tgtToSrcCellAddr_, i)
+ {
+ labelList& addressing = tgtToSrcCellAddr_[i];
+ forAll(addressing, addrI)
+ {
+ addressing[addrI] = globalSrcCells.toGlobal(addressing[addrI]);
+ }
+ }
+
+ // set up as a reverse distribute
+ mapDistribute::distribute
+ (
+ Pstream::nonBlocking,
+ List<labelPair>(),
+ tgt.nCells(),
+ map.constructMap(),
+ map.subMap(),
+ tgtToSrcCellAddr_,
+ ListPlusEqOp<label>(),
+ labelList()
+ );
+
+ // set up as a reverse distribute
+ mapDistribute::distribute
+ (
+ Pstream::nonBlocking,
+ List<labelPair>(),
+ tgt.nCells(),
+ map.constructMap(),
+ map.subMap(),
+ tgtToSrcCellWght_,
+ ListPlusEqOp<scalar>(),
+ scalarList()
+ );
+
+ // weights normalisation
+ normaliseWeights
+ (
+ "source",
+ src.cellVolumes(),
+ srcToTgtCellAddr_,
+ srcToTgtCellWght_
+ );
+
+ normaliseWeights
+ (
+ "target",
+ tgt.cellVolumes(),
+ tgtToSrcCellAddr_,
+ tgtToSrcCellWght_
+ );
+
+ // cache maps and reset addresses
+ List<Map<label> > cMap;
+ srcMapPtr_.reset
+ (
+ new mapDistribute(globalSrcCells, tgtToSrcCellAddr_, cMap)
+ );
+ tgtMapPtr_.reset
+ (
+ new mapDistribute(globalTgtCells, srcToTgtCellAddr_, cMap)
+ );
+
+ // collect volume intersection contributions
+ reduce(V_, sumOp<scalar>());
+ }
+ else
+ {
+ calcAddressing(src, tgt);
+
+ normaliseWeights
+ (
+ "source",
+ src.cellVolumes(),
+ srcToTgtCellAddr_,
+ srcToTgtCellWght_
+ );
+
+ normaliseWeights
+ (
+ "target",
+ tgt.cellVolumes(),
+ tgtToSrcCellAddr_,
+ tgtToSrcCellWght_
+ );
+ }
+
+ Info<< " Overlap volume: " << V_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::meshToMeshNew::~meshToMeshNew()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
new file mode 100644
index 00000000000..6b5211ef43d
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
@@ -0,0 +1,494 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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::meshToMeshNew
+
+Description
+ Class to calculate the cell-addressing between two overlapping meshes
+
+SourceFiles
+ meshToMeshNew.C
+ meshToMeshNewTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshToMeshNew_H
+#define meshToMeshNew_H
+
+#include "polyMesh.H"
+#include "boundBox.H"
+#include "mapDistribute.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class meshToMeshNew Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshToMeshNew
+{
+ // Private data
+
+ //- Name of source mesh region
+ const word srcRegionName_;
+
+ //- Name of target mesh region
+ const word tgtRegionName_;
+
+ //- Source to target cell addressing
+ labelListList srcToTgtCellAddr_;
+
+ //- Target to source cell addressing
+ labelListList tgtToSrcCellAddr_;
+
+ //- Source to target cell interplation weights
+ scalarListList srcToTgtCellWght_;
+
+ //- Target to source cell interpolation weights
+ scalarListList tgtToSrcCellWght_;
+
+ //- Cell total volume in overlap region [m3]
+ scalar V_;
+
+ //- Index of processor that holds all of both sides. -1 in all other
+ // cases
+ label singleMeshProc_;
+
+ //- Source map pointer - parallel running only
+ autoPtr<mapDistribute> srcMapPtr_;
+
+ //- Target map pointer - parallel running only
+ autoPtr<mapDistribute> tgtMapPtr_;
+
+ //- Flag to indicate that direct (one-to-one) mapping should be applied
+ bool directMapping_;
+
+ //- Tolerance used in volume overlap calculations
+ static scalar tolerance_;
+
+
+ // Private Member Functions
+
+ //- Helper function to add a constant offset to a list
+ template<class Type>
+ void add(UList<Type>& fld, const label offset) const;
+
+ //- Write the connectivity (debugging)
+ void writeConnectivity
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelListList& srcToTargetAddr
+ ) const;
+
+ //- Return src cell IDs for the overlap region
+ labelList maskCells(const polyMesh& src, const polyMesh& tgt) const;
+
+ //- Find indices of overlapping cells in src and tgt meshes - returns
+ // true if found a matching pair
+ bool findInitialSeeds
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+
+ // Direct (one-to-one) mapping
+
+ //- Append to list of src mesh seed indices
+ void appendToDirectSeeds
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ boolList& mapFlag,
+ labelList& srcTgtSeed,
+ DynamicList<label>& srcSeeds,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+ //- Main driver routine for direct mapping
+ void calcDirect
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcSeedI,
+ const label tgtSeedI
+ );
+
+
+ // Indirect (non-conformal) mapping
+
+ //- Normalise the interpolation weights
+ void normaliseWeights
+ (
+ const word& descriptor,
+ const scalarField& cellVolumes,
+ const labelListList& addr,
+ scalarListList& wght
+ ) const;
+
+ //- Append target cell neihgbour cells to cellIDs list
+ void appendNbrTgtCells
+ (
+ const label tgtCellI,
+ const polyMesh& tgt,
+ const DynamicList<label>& visitedTgtCells,
+ DynamicList<label>& nbrTgtCellIDs
+ ) const;
+
+ //- Set the next cells in the advancing front algorithm
+ void setNextCells
+ (
+ label& startSeedI,
+ label& srcCellI,
+ label& tgtCellI,
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const DynamicList<label>& visitedCells,
+ labelList& seedCells
+ ) const;
+
+ //- Return the true if cells intersect
+ bool intersect
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcCellI,
+ const label tgtCellI
+ ) const;
+
+ //- Return the intersection volume between two cells
+ scalar interVol
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcCellI,
+ const label tgtCellI
+ ) const;
+
+ //- Main driver routine for indirect mapping
+ void calcIndirect
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs,
+ boolList& mapFlag,
+ label& startSeedI
+ );
+
+
+ //- Calculate the addressing between overalping regions of src and tgt
+ // meshes - main driver function
+ void calcAddressing(const polyMesh& src, const polyMesh& tgt);
+
+
+ // Parallel operations
+
+ //- Determine whether the meshes are split across multiple pocessors
+ label calcDistribution
+ (
+ const polyMesh& src,
+ const polyMesh& tgt
+ ) const;
+
+ //- Determine which processor bounding-boxes overlap
+ label calcOverlappingProcs
+ (
+ const List<boundBox>& procBb,
+ const boundBox& bb,
+ boolList& overlaps
+ ) const;
+
+ //- Calculate the mapping between processors
+ autoPtr<mapDistribute> calcProcMap
+ (
+ const polyMesh& src,
+ const polyMesh& tgt
+ ) const;
+
+ //- Distribute mesh info from 'my' processor to others
+ void distributeCells
+ (
+ const mapDistribute& map,
+ const polyMesh& tgtMesh,
+ const globalIndex& globalI,
+ List<pointField>& points,
+ List<label>& nInternalFaces,
+ List<faceList>& faces,
+ List<labelList>& faceOwner,
+ List<labelList>& faceNeighbour,
+ List<labelList>& cellIDs,
+ List<labelList>& nbrProcIDs,
+ List<labelList>& procLocalFaceIDs
+ ) const;
+
+ //- Collect pieces of tgt mesh from other procssors and restructure
+ void distributeAndMergeCells
+ (
+ const mapDistribute& map,
+ const polyMesh& tgt,
+ const globalIndex& globalI,
+ pointField& tgtPoints,
+ faceList& tgtFaces,
+ labelList& tgtFaceOwners,
+ labelList& tgtFaceNeighbours,
+ labelList& tgtCellIDs
+ ) const;
+
+
+ //- Disallow default bitwise copy construct
+ meshToMeshNew(const meshToMeshNew&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const meshToMeshNew&);
+
+
+public:
+
+ //- Run-time type information
+ TypeName("meshToMeshNew");
+
+
+ //- Construct from source and target meshes
+ meshToMeshNew
+ (
+ const polyMesh& src,
+ const polyMesh& tgt,
+ const bool directMapping
+ );
+
+
+ //- Destructor
+ virtual ~meshToMeshNew();
+
+
+ // Member Functions
+
+ // Access
+
+ //- Return const access to the source to target cell addressing
+ inline const labelListList& srcToTgtCellAddr() const;
+
+ //- Return const access to the target to source cell addressing
+ inline const labelListList& tgtToSrcCellAddr() const;
+
+ //- Return const access to the source to target cell weights
+ inline const scalarListList& srcToTgtCellWght() const;
+
+ //- Return const access to the target to source cell weights
+ inline const scalarListList& tgtToSrcCellWght() const;
+
+ //- Return const access to the overlap volume
+ inline scalar V() const;
+
+
+ // Evaluation
+
+ // Source-to-target field mapping
+
+ //- Map field from src to tgt mesh with defined operation
+ // Values passed in via 'result' are used to initialise the
+ // return value
+ template<class Type, class CombineOp>
+ void mapSrcToTgt
+ (
+ const UList<Type>& srcFld,
+ const CombineOp& cop,
+ List<Type>& result
+ ) const;
+
+ //- Return the src field mapped to the tgt mesh with a defined
+ // operation. Initial values of the result are set to zero
+ template<class Type, class CombineOp>
+ tmp<Field<Type> > mapSrcToTgt
+ (
+ const Field<Type>& srcFld,
+ const CombineOp& cop
+ ) const;
+
+ //- Convenience function to map a tmp field to the tgt mesh
+ // with a defined operation
+ template<class Type, class CombineOp>
+ tmp<Field<Type> > mapSrcToTgt
+ (
+ const tmp<Field<Type> >& tsrcFld,
+ const CombineOp& cop
+ ) const;
+
+ //- Convenience function to map a field to the tgt mesh with a
+ // default operation (plusEqOp)
+ template<class Type>
+ tmp<Field<Type> > mapSrcToTgt
+ (
+ const Field<Type>& srcFld
+ ) const;
+
+ //- Convenience function to map a tmp field to the tgt mesh
+ // with a default operation (plusEqOp)
+ template<class Type>
+ tmp<Field<Type> > mapSrcToTgt
+ (
+ const tmp<Field<Type> >& tsrcFld
+ ) const;
+
+
+ // Target-to-source field mapping
+
+ //- Map field from tgt to src mesh with defined operation
+ // Values passed in via 'result' are used to initialise the
+ // return value
+ template<class Type, class CombineOp>
+ void mapTgtToSrc
+ (
+ const UList<Type>& tgtFld,
+ const CombineOp& cop,
+ List<Type>& result
+ ) const;
+
+ //- Return the tgt field mapped to the src mesh with a defined
+ // operation. Initial values of the result are set to zero
+ template<class Type, class CombineOp>
+ tmp<Field<Type> > mapTgtToSrc
+ (
+ const Field<Type>& tgtFld,
+ const CombineOp& cop
+ ) const;
+
+ //- Convenience function to map a tmp field to the src mesh
+ // with a defined operation
+ template<class Type, class CombineOp>
+ tmp<Field<Type> > mapTgtToSrc
+ (
+ const tmp<Field<Type> >& ttgtFld,
+ const CombineOp& cop
+ ) const;
+
+ //- Convenience function to map a field to the src mesh with a
+ // default operation (plusEqOp)
+ template<class Type>
+ tmp<Field<Type> > mapTgtToSrc
+ (
+ const Field<Type>& tgtFld
+ ) const;
+
+ //- Convenience function to map a tmp field to the src mesh
+ // with a default operation (plusEqOp)
+ template<class Type>
+ tmp<Field<Type> > mapTgtToSrc
+ (
+ const tmp<Field<Type> >& ttgtFld
+ ) const;
+
+
+ // Volume field mapping
+
+ //- Interpolare a field with a defined operation. Values
+ // passed in via 'result' are used to initialise the return
+ // value. Optionally interpolate patch values
+ template<class Type, class CombineOp>
+ void interpolate
+ (
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const CombineOp& cop,
+ GeometricField<Type, fvPatchField, volMesh>& result,
+ const bool interpPatches = true
+ ) const;
+
+ //- Interpolare a field with a defined operation. The initial
+ // values of the result are set to zero. Optionally
+ // interpolate patch values
+ template<class Type, class CombineOp>
+ tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+ (
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const CombineOp& cop,
+ const bool interpPatches = true
+ ) const;
+
+ //- Interpolare a tmp field with a defined operation. The
+ // initial values of the result are set to zero. Optionally
+ // interpolate patch values
+ template<class Type, class CombineOp>
+ tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+ (
+ const tmp<GeometricField<Type, fvPatchField, volMesh> >&
+ tfield,
+ const CombineOp& cop,
+ const bool interpPatches = true
+ ) const;
+
+ //- Convenience function to map a field with a default
+ // operation (plusEqOp). Optionally interpolate patch values
+ template<class Type>
+ tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+ (
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const bool interpPatches = true
+ ) const;
+
+ //- Convenience function to map a tmp field with a default
+ // operation (plusEqOp). Optionally interpolate patch values
+ template<class Type>
+ tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+ (
+ const tmp<GeometricField<Type, fvPatchField, volMesh> >&
+ tfield,
+ const bool interpPatches = true
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "meshToMeshNewI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "meshToMeshNewTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H
new file mode 100644
index 00000000000..41068e8c455
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H"
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+inline const Foam::labelListList&
+Foam::meshToMeshNew::srcToTgtCellAddr() const
+{
+ return srcToTgtCellAddr_;
+}
+
+
+inline const Foam::labelListList&
+Foam::meshToMeshNew::tgtToSrcCellAddr() const
+{
+ return tgtToSrcCellAddr_;
+}
+
+
+inline const Foam::scalarListList&
+Foam::meshToMeshNew::srcToTgtCellWght() const
+{
+ return srcToTgtCellWght_;
+}
+
+
+inline const Foam::scalarListList&
+Foam::meshToMeshNew::tgtToSrcCellWght() const
+{
+ return tgtToSrcCellWght_;
+}
+
+
+inline Foam::scalar Foam::meshToMeshNew::V() const
+{
+ return V_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C
new file mode 100644
index 00000000000..3cb9d0597d8
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C
@@ -0,0 +1,884 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H"
+#include "OFstream.H"
+#include "Time.H"
+#include "globalIndex.H"
+#include "mergePoints.H"
+#include "processorPolyPatch.H"
+#include "SubField.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::label Foam::meshToMeshNew::calcDistribution
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+) const
+{
+ label procI = 0;
+
+ if (Pstream::parRun())
+ {
+ List<label> cellsPresentOnProc(Pstream::nProcs(), 0);
+ if ((src.nCells() > 0) || (tgt.nCells() > 0))
+ {
+ cellsPresentOnProc[Pstream::myProcNo()] = 1;
+ }
+ else
+ {
+ cellsPresentOnProc[Pstream::myProcNo()] = 0;
+ }
+
+ Pstream::gatherList(cellsPresentOnProc);
+ Pstream::scatterList(cellsPresentOnProc);
+
+ label nHaveCells = sum(cellsPresentOnProc);
+
+ if (nHaveCells > 1)
+ {
+ procI = -1;
+ if (debug)
+ {
+ Info<< "meshToMeshNew::calcDistribution: "
+ << "Meshes split across multiple processors" << endl;
+ }
+ }
+ else if (nHaveCells == 1)
+ {
+ procI = findIndex(cellsPresentOnProc, 1);
+ if (debug)
+ {
+ Info<< "meshToMeshNew::calcDistribution: "
+ << "Meshes local to processor" << procI << endl;
+ }
+ }
+ }
+
+ return procI;
+}
+
+
+Foam::label Foam::meshToMeshNew::calcOverlappingProcs
+(
+ const List<boundBox>& procBb,
+ const boundBox& bb,
+ boolList& overlaps
+) const
+{
+ overlaps = false;
+
+ label nOverlaps = 0;
+
+ forAll(procBb, procI)
+ {
+ const boundBox& bbp = procBb[procI];
+
+ if (bbp.overlaps(bb))
+ {
+ overlaps[procI] = true;
+ nOverlaps++;
+ }
+ }
+
+ return nOverlaps;
+}
+
+
+Foam::autoPtr<Foam::mapDistribute> Foam::meshToMeshNew::calcProcMap
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+) const
+{
+ // get decomposition of cells on src mesh
+ List<boundBox> procBb(Pstream::nProcs());
+
+ if (src.nCells() > 0)
+ {
+ // bounding box for my mesh - do not parallel reduce
+ procBb[Pstream::myProcNo()] = boundBox(src.points(), false);
+
+ // slightly increase size of bounding boxes to allow for cases where
+ // bounding boxes are perfectly alligned
+ procBb[Pstream::myProcNo()].inflate(0.01);
+ }
+ else
+ {
+ procBb[Pstream::myProcNo()] = boundBox();
+ }
+
+
+ Pstream::gatherList(procBb);
+ Pstream::scatterList(procBb);
+
+
+ if (debug)
+ {
+ Info<< "Determining extent of src mesh per processor:" << nl
+ << "\tproc\tbb" << endl;
+ forAll(procBb, procI)
+ {
+ Info<< '\t' << procI << '\t' << procBb[procI] << endl;
+ }
+ }
+
+
+ // determine which cells of tgt mesh overlaps src mesh per proc
+ const cellList& cells = tgt.cells();
+ const faceList& faces = tgt.faces();
+ const pointField& points = tgt.points();
+
+ labelListList sendMap;
+
+ {
+ // per processor indices into all segments to send
+ List<DynamicList<label> > dynSendMap(Pstream::nProcs());
+ label iniSize = floor(tgt.nCells()/Pstream::nProcs());
+
+ forAll(dynSendMap, procI)
+ {
+ dynSendMap[procI].setCapacity(iniSize);
+ }
+
+ // work array - whether src processor bb overlaps the tgt cell bounds
+ boolList procBbOverlaps(Pstream::nProcs());
+ forAll(cells, cellI)
+ {
+ const cell& c = cells[cellI];
+
+ // determine bounding box of tgt cell
+ boundBox cellBb(point::max, point::min);
+ forAll(c, faceI)
+ {
+ const face& f = faces[c[faceI]];
+ forAll(f, fp)
+ {
+ cellBb.min() = min(cellBb.min(), points[f[fp]]);
+ cellBb.max() = max(cellBb.max(), points[f[fp]]);
+ }
+ }
+
+ // find the overlapping tgt cells on each src processor
+ (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
+
+ forAll(procBbOverlaps, procI)
+ {
+ if (procBbOverlaps[procI])
+ {
+ dynSendMap[procI].append(cellI);
+ }
+ }
+ }
+
+ // convert dynamicList to labelList
+ sendMap.setSize(Pstream::nProcs());
+ forAll(sendMap, procI)
+ {
+ sendMap[procI].transfer(dynSendMap[procI]);
+ }
+ }
+
+ // debug printing
+ if (debug)
+ {
+ Pout<< "Of my " << cells.size() << " target cells I need to send to:"
+ << nl << "\tproc\tcells" << endl;
+ forAll(sendMap, procI)
+ {
+ Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
+ }
+ }
+
+
+ // send over how many tgt cells I need to receive from each processor
+ labelListList sendSizes(Pstream::nProcs());
+ sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+ forAll(sendMap, procI)
+ {
+ sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+ }
+ Pstream::gatherList(sendSizes);
+ Pstream::scatterList(sendSizes);
+
+
+ // determine order of receiving
+ labelListList constructMap(Pstream::nProcs());
+
+ label segmentI = 0;
+ forAll(constructMap, procI)
+ {
+ // what I need to receive is what other processor is sending to me
+ label nRecv = sendSizes[procI][Pstream::myProcNo()];
+ constructMap[procI].setSize(nRecv);
+
+ for (label i = 0; i < nRecv; i++)
+ {
+ constructMap[procI][i] = segmentI++;
+ }
+ }
+
+ autoPtr<mapDistribute> mapPtr
+ (
+ new mapDistribute
+ (
+ segmentI, // size after construction
+ sendMap.xfer(),
+ constructMap.xfer()
+ )
+ );
+
+ return mapPtr;
+}
+
+
+void Foam::meshToMeshNew::distributeCells
+(
+ const mapDistribute& map,
+ const polyMesh& tgtMesh,
+ const globalIndex& globalI,
+ List<pointField>& points,
+ List<label>& nInternalFaces,
+ List<faceList>& faces,
+ List<labelList>& faceOwner,
+ List<labelList>& faceNeighbour,
+ List<labelList>& cellIDs,
+ List<labelList>& nbrProcIDs,
+ List<labelList>& procLocalFaceIDs
+) const
+{
+ PstreamBuffers pBufs(Pstream::nonBlocking);
+
+ points.setSize(Pstream::nProcs());
+ nInternalFaces.setSize(Pstream::nProcs(), 0);
+ faces.setSize(Pstream::nProcs());
+ faceOwner.setSize(Pstream::nProcs());
+ faceNeighbour.setSize(Pstream::nProcs());
+ cellIDs.setSize(Pstream::nProcs());
+
+ nbrProcIDs.setSize(Pstream::nProcs());;
+ procLocalFaceIDs.setSize(Pstream::nProcs());;
+
+
+ for (label domain = 0; domain < Pstream::nProcs(); domain++)
+ {
+ const labelList& sendElems = map.subMap()[domain];
+
+ if (sendElems.size())
+ {
+ // reverse cell map
+ labelList reverseCellMap(tgtMesh.nCells(), -1);
+ forAll(sendElems, subCellI)
+ {
+ reverseCellMap[sendElems[subCellI]] = subCellI;
+ }
+
+ DynamicList<face> subFaces(tgtMesh.nFaces());
+ DynamicList<label> subFaceOwner(tgtMesh.nFaces());
+ DynamicList<label> subFaceNeighbour(tgtMesh.nFaces());
+
+ DynamicList<label> subNbrProcIDs(tgtMesh.nFaces());
+ DynamicList<label> subProcLocalFaceIDs(tgtMesh.nFaces());
+
+ label nInternal = 0;
+
+ // internal faces
+ forAll(tgtMesh.faceNeighbour(), faceI)
+ {
+ label own = tgtMesh.faceOwner()[faceI];
+ label nbr = tgtMesh.faceNeighbour()[faceI];
+ label subOwn = reverseCellMap[own];
+ label subNbr = reverseCellMap[nbr];
+
+ if (subOwn != -1 && subNbr != -1)
+ {
+ nInternal++;
+
+ if (subOwn < subNbr)
+ {
+ subFaces.append(tgtMesh.faces()[faceI]);
+ subFaceOwner.append(subOwn);
+ subFaceNeighbour.append(subNbr);
+ subNbrProcIDs.append(-1);
+ subProcLocalFaceIDs.append(-1);
+ }
+ else
+ {
+ subFaces.append(tgtMesh.faces()[faceI].reverseFace());
+ subFaceOwner.append(subNbr);
+ subFaceNeighbour.append(subOwn);
+ subNbrProcIDs.append(-1);
+ subProcLocalFaceIDs.append(-1);
+ }
+ }
+ }
+
+ // boundary faces for new region
+ forAll(tgtMesh.faceNeighbour(), faceI)
+ {
+ label own = tgtMesh.faceOwner()[faceI];
+ label nbr = tgtMesh.faceNeighbour()[faceI];
+ label subOwn = reverseCellMap[own];
+ label subNbr = reverseCellMap[nbr];
+
+ if (subOwn != -1 && subNbr == -1)
+ {
+ subFaces.append(tgtMesh.faces()[faceI]);
+ subFaceOwner.append(subOwn);
+ subFaceNeighbour.append(subNbr);
+ subNbrProcIDs.append(-1);
+ subProcLocalFaceIDs.append(-1);
+ }
+ else if (subOwn == -1 && subNbr != -1)
+ {
+ subFaces.append(tgtMesh.faces()[faceI].reverseFace());
+ subFaceOwner.append(subNbr);
+ subFaceNeighbour.append(subOwn);
+ subNbrProcIDs.append(-1);
+ subProcLocalFaceIDs.append(-1);
+ }
+ }
+
+ // boundary faces of existing region
+ forAll(tgtMesh.boundaryMesh(), patchI)
+ {
+ const polyPatch& pp = tgtMesh.boundaryMesh()[patchI];
+
+ label nbrProcI = -1;
+
+ // store info for faces on processor patches
+ if (isA<processorPolyPatch>(pp))
+ {
+ const processorPolyPatch& ppp =
+ dynamic_cast<const processorPolyPatch&>(pp);
+
+ nbrProcI = ppp.neighbProcNo();
+ }
+
+ forAll(pp, i)
+ {
+ label faceI = pp.start() + i;
+ label own = tgtMesh.faceOwner()[faceI];
+
+ if (reverseCellMap[own] != -1)
+ {
+ subFaces.append(tgtMesh.faces()[faceI]);
+ subFaceOwner.append(reverseCellMap[own]);
+ subFaceNeighbour.append(-1);
+ subNbrProcIDs.append(nbrProcI);
+ subProcLocalFaceIDs.append(i);
+ }
+ }
+ }
+
+ // reverse point map
+ labelList reversePointMap(tgtMesh.nPoints(), -1);
+ DynamicList<point> subPoints(tgtMesh.nPoints());
+ forAll(subFaces, subFaceI)
+ {
+ face& f = subFaces[subFaceI];
+ forAll(f, fp)
+ {
+ label pointI = f[fp];
+ if (reversePointMap[pointI] == -1)
+ {
+ reversePointMap[pointI] = subPoints.size();
+ subPoints.append(tgtMesh.points()[pointI]);
+ }
+
+ f[fp] = reversePointMap[pointI];
+ }
+ }
+
+ // tgt cells into global numbering
+ labelList globalElems(sendElems.size());
+ forAll(sendElems, i)
+ {
+ if (debug)
+ {
+ Pout<< "tgtProc:" << Pstream::myProcNo()
+ << " sending tgt cell " << sendElems[i]
+ << "[" << globalI.toGlobal(sendElems[i]) << "]"
+ << " to srcProc " << domain << endl;
+ }
+
+ globalElems[i] = globalI.toGlobal(sendElems[i]);
+ }
+
+ // pass data
+ if (domain == Pstream::myProcNo())
+ {
+ // allocate my own data
+ points[Pstream::myProcNo()] = subPoints;
+ nInternalFaces[Pstream::myProcNo()] = nInternal;
+ faces[Pstream::myProcNo()] = subFaces;
+ faceOwner[Pstream::myProcNo()] = subFaceOwner;
+ faceNeighbour[Pstream::myProcNo()] = subFaceNeighbour;
+ cellIDs[Pstream::myProcNo()] = globalElems;
+ nbrProcIDs[Pstream::myProcNo()] = subNbrProcIDs;
+ procLocalFaceIDs[Pstream::myProcNo()] = subProcLocalFaceIDs;
+ }
+ else
+ {
+ // send data to other processor domains
+ UOPstream toDomain(domain, pBufs);
+
+ toDomain
+ << subPoints
+ << nInternal
+ << subFaces
+ << subFaceOwner
+ << subFaceNeighbour
+ << globalElems
+ << subNbrProcIDs
+ << subProcLocalFaceIDs;
+ }
+ }
+ }
+
+ // Start receiving
+ pBufs.finishedSends();
+
+ // Consume
+ for (label domain = 0; domain < Pstream::nProcs(); domain++)
+ {
+ const labelList& recvElems = map.constructMap()[domain];
+
+ if (domain != Pstream::myProcNo() && recvElems.size())
+ {
+ UIPstream str(domain, pBufs);
+
+ str >> points[domain]
+ >> nInternalFaces[domain]
+ >> faces[domain]
+ >> faceOwner[domain]
+ >> faceNeighbour[domain]
+ >> cellIDs[domain]
+ >> nbrProcIDs[domain]
+ >> procLocalFaceIDs[domain];
+ }
+
+ if (debug)
+ {
+ Pout<< "Target mesh send sizes[" << domain << "]"
+ << ": points="<< points[domain].size()
+ << ", faces=" << faces[domain].size()
+ << ", nInternalFaces=" << nInternalFaces[domain]
+ << ", faceOwn=" << faceOwner[domain].size()
+ << ", faceNbr=" << faceNeighbour[domain].size()
+ << ", cellIDs=" << cellIDs[domain].size() << endl;
+ }
+ }
+}
+
+
+void Foam::meshToMeshNew::distributeAndMergeCells
+(
+ const mapDistribute& map,
+ const polyMesh& tgt,
+ const globalIndex& globalI,
+ pointField& tgtPoints,
+ faceList& tgtFaces,
+ labelList& tgtFaceOwners,
+ labelList& tgtFaceNeighbours,
+ labelList& tgtCellIDs
+) const
+{
+ // Exchange per-processor data
+ List<pointField> allPoints;
+ List<label> allNInternalFaces;
+ List<faceList> allFaces;
+ List<labelList> allFaceOwners;
+ List<labelList> allFaceNeighbours;
+ List<labelList> allTgtCellIDs;
+
+ // Per target mesh face the neighbouring proc and index in
+ // processor patch (all -1 for normal boundary face)
+ List<labelList> allNbrProcIDs;
+ List<labelList> allProcLocalFaceIDs;
+
+ distributeCells
+ (
+ map,
+ tgt,
+ globalI,
+ allPoints,
+ allNInternalFaces,
+ allFaces,
+ allFaceOwners,
+ allFaceNeighbours,
+ allTgtCellIDs,
+ allNbrProcIDs,
+ allProcLocalFaceIDs
+ );
+
+ // Convert lists into format that can be used to generate a valid polyMesh
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ //
+ // Points and cells are collected into single flat lists:
+ // - i.e. proc0, proc1 ... procN
+ //
+ // Faces need to be sorted after collection to that internal faces are
+ // contiguous, followed by all boundary faces
+ //
+ // Processor patch faces between included cells on neighbouring processors
+ // are converted into internal faces
+ //
+ // Face list structure:
+ // - Per processor:
+ // - internal faces
+ // - processor faces that have been converted into internal faces
+ // - Followed by all boundary faces
+ // - from 'normal' boundary faces
+ // - from singularly-sided processor patch faces
+
+
+ // Number of internal+coupled faces
+ labelList allNIntCoupledFaces(allNInternalFaces);
+
+ // Starting offset for points
+ label nPoints = 0;
+ labelList pointOffset(Pstream::nProcs(), 0);
+ forAll(allPoints, procI)
+ {
+ pointOffset[procI] = nPoints;
+ nPoints += allPoints[procI].size();
+ }
+
+ // Starting offset for cells
+ label nCells = 0;
+ labelList cellOffset(Pstream::nProcs(), 0);
+ forAll(allTgtCellIDs, procI)
+ {
+ cellOffset[procI] = nCells;
+ nCells += allTgtCellIDs[procI].size();
+ }
+
+ // Count any coupled faces
+ typedef FixedList<label, 3> label3;
+ typedef HashTable<label, label3, label3::Hash<> > procCoupleInfo;
+ procCoupleInfo procFaceToGlobalCell;
+
+ forAll(allNbrProcIDs, procI)
+ {
+ const labelList& nbrProcI = allNbrProcIDs[procI];
+ const labelList& localFaceI = allProcLocalFaceIDs[procI];
+
+ forAll(nbrProcI, i)
+ {
+ if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+ {
+ label3 key;
+ key[0] = min(procI, nbrProcI[i]);
+ key[1] = max(procI, nbrProcI[i]);
+ key[2] = localFaceI[i];
+
+ procCoupleInfo::const_iterator fnd =
+ procFaceToGlobalCell.find(key);
+
+ if (fnd == procFaceToGlobalCell.end())
+ {
+ procFaceToGlobalCell.insert(key, -1);
+ }
+ else
+ {
+ if (debug)
+ {
+ Pout<< "Additional internal face between procs:"
+ << key[0] << " and " << key[1]
+ << " across local face " << key[2] << endl;
+ }
+
+ allNIntCoupledFaces[procI]++;
+ }
+ }
+ }
+ }
+
+
+ // Starting offset for internal faces
+ label nIntFaces = 0;
+ label nFacesTotal = 0;
+ labelList internalFaceOffset(Pstream::nProcs(), 0);
+ forAll(allNIntCoupledFaces, procI)
+ {
+ label nCoupledFaces =
+ allNIntCoupledFaces[procI] - allNInternalFaces[procI];
+
+ internalFaceOffset[procI] = nIntFaces;
+ nIntFaces += allNIntCoupledFaces[procI];
+ nFacesTotal += allFaceOwners[procI].size() - nCoupledFaces;
+ }
+
+ tgtPoints.setSize(nPoints);
+ tgtFaces.setSize(nFacesTotal);
+ tgtFaceOwners.setSize(nFacesTotal);
+ tgtFaceNeighbours.setSize(nFacesTotal);
+ tgtCellIDs.setSize(nCells);
+
+ // Insert points
+ forAll(allPoints, procI)
+ {
+ const pointField& pts = allPoints[procI];
+ SubList<point>(tgtPoints, pts.size(), pointOffset[procI]).assign(pts);
+ }
+
+ // Insert cellIDs
+ forAll(allTgtCellIDs, procI)
+ {
+ const labelList& cellIDs = allTgtCellIDs[procI];
+ SubList<label>(tgtCellIDs, cellIDs.size(), cellOffset[procI]).assign
+ (
+ cellIDs
+ );
+ }
+
+
+ // Insert internal faces (from internal faces)
+ forAll(allFaces, procI)
+ {
+ const faceList& fcs = allFaces[procI];
+ const labelList& faceOs = allFaceOwners[procI];
+ const labelList& faceNs = allFaceNeighbours[procI];
+
+ SubList<face> slice
+ (
+ tgtFaces,
+ allNInternalFaces[procI],
+ internalFaceOffset[procI]
+ );
+ slice.assign(SubList<face>(fcs, allNInternalFaces[procI]));
+ forAll(slice, i)
+ {
+ add(slice[i], pointOffset[procI]);
+ }
+
+ SubField<label> ownSlice
+ (
+ tgtFaceOwners,
+ allNInternalFaces[procI],
+ internalFaceOffset[procI]
+ );
+ ownSlice.assign(SubField<label>(faceOs, allNInternalFaces[procI]));
+ add(ownSlice, cellOffset[procI]);
+
+ SubField<label> nbrSlice
+ (
+ tgtFaceNeighbours,
+ allNInternalFaces[procI],
+ internalFaceOffset[procI]
+ );
+ nbrSlice.assign(SubField<label>(faceNs, allNInternalFaces[procI]));
+ add(nbrSlice, cellOffset[procI]);
+
+ internalFaceOffset[procI] += allNInternalFaces[procI];
+ }
+
+
+ // Insert internal faces (from coupled face-pairs)
+ forAll(allNbrProcIDs, procI)
+ {
+ const labelList& nbrProcI = allNbrProcIDs[procI];
+ const labelList& localFaceI = allProcLocalFaceIDs[procI];
+ const labelList& faceOs = allFaceOwners[procI];
+ const faceList& fcs = allFaces[procI];
+
+ forAll(nbrProcI, i)
+ {
+ if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+ {
+ label3 key;
+ key[0] = min(procI, nbrProcI[i]);
+ key[1] = max(procI, nbrProcI[i]);
+ key[2] = localFaceI[i];
+
+ procCoupleInfo::iterator fnd = procFaceToGlobalCell.find(key);
+
+ if (fnd != procFaceToGlobalCell.end())
+ {
+ label tgtFaceI = fnd();
+ if (tgtFaceI == -1)
+ {
+ // on first visit store the new cell on this side
+ fnd() = cellOffset[procI] + faceOs[i];
+ }
+ else
+ {
+ // get owner and neighbour in new cell numbering
+ label newOwn = cellOffset[procI] + faceOs[i];
+ label newNbr = fnd();
+ label tgtFaceI = internalFaceOffset[procI]++;
+
+ if (debug)
+ {
+ Pout<< " proc " << procI
+ << "\tinserting face:" << tgtFaceI
+ << " connection between owner " << newOwn
+ << " and neighbour " << newNbr
+ << endl;
+ }
+
+ if (newOwn < newNbr)
+ {
+ // we have correct orientation
+ tgtFaces[tgtFaceI] = fcs[i];
+ tgtFaceOwners[tgtFaceI] = newOwn;
+ tgtFaceNeighbours[tgtFaceI] = newNbr;
+ }
+ else
+ {
+ // reverse orientation
+ tgtFaces[tgtFaceI] = fcs[i].reverseFace();
+ tgtFaceOwners[tgtFaceI] = newNbr;
+ tgtFaceNeighbours[tgtFaceI] = newOwn;
+ }
+
+ add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+ // mark with unique value
+ fnd() = -2;
+ }
+ }
+ }
+ }
+ }
+
+
+ forAll(allNbrProcIDs, procI)
+ {
+ const labelList& nbrProcI = allNbrProcIDs[procI];
+ const labelList& localFaceI = allProcLocalFaceIDs[procI];
+ const labelList& faceOs = allFaceOwners[procI];
+ const labelList& faceNs = allFaceNeighbours[procI];
+ const faceList& fcs = allFaces[procI];
+
+ forAll(nbrProcI, i)
+ {
+ // coupled boundary face
+ if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+ {
+ label3 key;
+ key[0] = min(procI, nbrProcI[i]);
+ key[1] = max(procI, nbrProcI[i]);
+ key[2] = localFaceI[i];
+
+ label tgtFaceI = procFaceToGlobalCell[key];
+
+ if (tgtFaceI == -1)
+ {
+ FatalErrorIn
+ (
+ "void Foam::meshToMeshNew::"
+ "distributeAndMergeCells"
+ "("
+ "const mapDistribute&, "
+ "const polyMesh&, "
+ "const globalIndex&, "
+ "pointField&, "
+ "faceList&, "
+ "labelList&, "
+ "labelList&, "
+ "labelList&"
+ ") const"
+ )
+ << "Unvisited " << key
+ << abort(FatalError);
+ }
+ else if (tgtFaceI != -2)
+ {
+ label newOwn = cellOffset[procI] + faceOs[i];
+ label tgtFaceI = nIntFaces++;
+
+ if (debug)
+ {
+ Pout<< " proc " << procI
+ << "\tinserting boundary face:" << tgtFaceI
+ << " from coupled face " << key
+ << endl;
+ }
+
+ tgtFaces[tgtFaceI] = fcs[i];
+ add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+ tgtFaceOwners[tgtFaceI] = newOwn;
+ tgtFaceNeighbours[tgtFaceI] = -1;
+ }
+ }
+ // normal boundary face
+ else
+ {
+ label own = faceOs[i];
+ label nbr = faceNs[i];
+ if ((own != -1) && (nbr == -1))
+ {
+ label newOwn = cellOffset[procI] + faceOs[i];
+ label tgtFaceI = nIntFaces++;
+
+ tgtFaces[tgtFaceI] = fcs[i];
+ add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+ tgtFaceOwners[tgtFaceI] = newOwn;
+ tgtFaceNeighbours[tgtFaceI] = -1;
+ }
+ }
+ }
+ }
+
+
+ if (debug)
+ {
+ // only merging points in debug mode
+
+ labelList oldToNew;
+ pointField newTgtPoints;
+ bool hasMerged = mergePoints
+ (
+ tgtPoints,
+ SMALL,
+ false,
+ oldToNew,
+ newTgtPoints
+ );
+
+ if (hasMerged)
+ {
+ if (debug)
+ {
+ Pout<< "Merged from " << tgtPoints.size()
+ << " down to " << newTgtPoints.size() << " points" << endl;
+ }
+
+ tgtPoints.transfer(newTgtPoints);
+ forAll(tgtFaces, i)
+ {
+ inplaceRenumber(oldToNew, tgtFaces[i]);
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
new file mode 100644
index 00000000000..4587b8ea60a
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
@@ -0,0 +1,537 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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 "fvMesh.H"
+#include "volFields.H"
+//#include "ops.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ //- Helper class for list
+ template<class Type>
+ class ListPlusEqOp
+ {
+ public:
+ void operator()(List<Type>& x, const List<Type> y) const
+ {
+ if (y.size())
+ {
+ if (x.size())
+ {
+ label sz = x.size();
+ x.setSize(sz + y.size());
+ forAll(y, i)
+ {
+ x[sz++] = y[i];
+ }
+ }
+ else
+ {
+ x = y;
+ }
+ }
+ }
+ };
+
+ //- Combine operator for maps/interpolations
+ template<class Type, class CombineOp>
+ class combineBinaryOp
+ {
+ const CombineOp& cop_;
+
+ public:
+
+ combineBinaryOp(const CombineOp& cop)
+ :
+ cop_(cop)
+ {}
+
+ void operator()
+ (
+ Type& x,
+ const label faceI,
+ const Type& y,
+ const scalar weight
+ ) const
+ {
+ cop_(x, weight*y);
+ }
+ };
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::meshToMeshNew::add
+(
+ UList<Type>& fld,
+ const label offset
+) const
+{
+ forAll(fld, i)
+ {
+ fld[i] += offset;
+ }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::mapSrcToTgt
+(
+ const UList<Type>& srcField,
+ const CombineOp& cop,
+ List<Type>& result
+) const
+{
+ if (result.size() != tgtToSrcCellAddr_.size())
+ {
+ FatalErrorIn
+ (
+ "void Foam::meshToMeshNew::mapSrcToTgt"
+ "("
+ "const UList<Type>&, "
+ "const CombineOp&, "
+ "List<Type>&"
+ ") const"
+ ) << "Supplied field size is not equal to target mesh size" << nl
+ << " source mesh = " << srcToTgtCellAddr_.size() << nl
+ << " target mesh = " << tgtToSrcCellAddr_.size() << nl
+ << " supplied field = " << result.size()
+ << abort(FatalError);
+ }
+
+ combineBinaryOp<Type, CombineOp> cbop(cop);
+
+ if (singleMeshProc_ == -1)
+ {
+ const mapDistribute& map = srcMapPtr_();
+
+ List<Type> work(srcField);
+ map.distribute(work);
+
+ forAll(result, cellI)
+ {
+ const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
+ const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
+
+ if (srcAddress.size())
+ {
+// result[cellI] = pTraits<Type>::zero;
+ result[cellI] *= (1.0 - sum(srcWeight));
+ forAll(srcAddress, i)
+ {
+ label srcI = srcAddress[i];
+ scalar w = srcWeight[i];
+ cbop(result[cellI], cellI, work[srcI], w);
+ }
+ }
+ }
+ }
+ else
+ {
+ forAll(result, cellI)
+ {
+ const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
+ const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
+
+ if (srcAddress.size())
+ {
+// result[cellI] = pTraits<Type>::zero;
+ result[cellI] *= (1.0 - sum(srcWeight));
+ forAll(srcAddress, i)
+ {
+ label srcI = srcAddress[i];
+ scalar w = srcWeight[i];
+ cbop(result[cellI], cellI, srcField[srcI], w);
+ }
+ }
+ }
+ }
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+ const Field<Type>& srcField,
+ const CombineOp& cop
+) const
+{
+ tmp<Field<Type> > tresult
+ (
+ new Field<Type>
+ (
+ tgtToSrcCellAddr_.size(),
+ pTraits<Type>::zero
+ )
+ );
+
+ mapSrcToTgt(srcField, cop, tresult());
+
+ return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+ const tmp<Field<Type> >& tsrcField,
+ const CombineOp& cop
+) const
+{
+ return mapSrcToTgt(tsrcField(), cop);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+ const Field<Type>& srcField
+) const
+{
+ return mapSrcToTgt(srcField, plusEqOp<Type>());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+ const tmp<Field<Type> >& tsrcField
+) const
+{
+ return mapSrcToTgt(tsrcField());
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::mapTgtToSrc
+(
+ const UList<Type>& tgtField,
+ const CombineOp& cop,
+ List<Type>& result
+) const
+{
+ if (result.size() != srcToTgtCellAddr_.size())
+ {
+ FatalErrorIn
+ (
+ "void Foam::meshToMeshNew::mapTgtToSrc"
+ "("
+ "const UList<Type>&, "
+ "const CombineOp&, "
+ "List<Type>&"
+ ") const"
+ ) << "Supplied field size is not equal to source mesh size" << nl
+ << " source mesh = " << srcToTgtCellAddr_.size() << nl
+ << " target mesh = " << tgtToSrcCellAddr_.size() << nl
+ << " supplied field = " << result.size()
+ << abort(FatalError);
+ }
+
+ combineBinaryOp<Type, CombineOp> cbop(cop);
+
+ if (singleMeshProc_ == -1)
+ {
+ const mapDistribute& map = tgtMapPtr_();
+
+ List<Type> work(tgtField);
+ map.distribute(work);
+
+ forAll(result, cellI)
+ {
+ const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
+ const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
+
+ if (tgtAddress.size())
+ {
+// result[cellI] = pTraits<Type>::zero;
+ result[cellI] *= (1.0 - sum(tgtWeight));
+ forAll(tgtAddress, i)
+ {
+ label tgtI = tgtAddress[i];
+ scalar w = tgtWeight[i];
+ cbop(result[cellI], cellI, work[tgtI], w);
+ }
+ }
+ }
+ }
+ else
+ {
+ forAll(result, cellI)
+ {
+ const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
+ const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
+
+ if (tgtAddress.size())
+ {
+// result[cellI] = pTraits<Type>::zero;
+ result[cellI] *= (1.0 - sum(tgtWeight));
+ forAll(tgtAddress, i)
+ {
+ label tgtI = tgtAddress[i];
+ scalar w = tgtWeight[i];
+ cbop(result[cellI], cellI, tgtField[tgtI], w);
+ }
+ }
+ }
+ }
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+ const Field<Type>& tgtField,
+ const CombineOp& cop
+) const
+{
+ tmp<Field<Type> > tresult
+ (
+ new Field<Type>
+ (
+ srcToTgtCellAddr_.size(),
+ pTraits<Type>::zero
+ )
+ );
+
+ mapTgtToSrc(tgtField, cop, tresult());
+
+ return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+ const tmp<Field<Type> >& ttgtField,
+ const CombineOp& cop
+) const
+{
+ return mapTgtToSrc(ttgtField(), cop);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+ const Field<Type>& tgtField
+) const
+{
+ return mapTgtToSrc(tgtField, plusEqOp<Type>());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+ const tmp<Field<Type> >& ttgtField
+) const
+{
+ return mapTgtToSrc(ttgtField(), plusEqOp<Type>());
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::interpolate
+(
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const CombineOp& cop,
+ GeometricField<Type, fvPatchField, volMesh>& result,
+ const bool interpPatches
+) const
+{
+ const fvMesh& mesh = field.mesh();
+
+ if (mesh.name() == srcRegionName_)
+ {
+ mapSrcToTgt(field, cop, result.internalField());
+ }
+ else if (mesh.name() == tgtRegionName_)
+ {
+ mapTgtToSrc(field, cop, result.internalField());
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "void Foam::meshToMeshNew::interpolate"
+ "("
+ "const GeometricField<Type, fvPatchField, volMesh>&, "
+ "const CombineOp&, "
+ "GeometricField<Type, fvPatchField, volMesh>&, "
+ "const bool"
+ ") const"
+ )
+ << "Supplied field " << field.name() << " did not originate from "
+ << "either the source or target meshes used to create this "
+ << "interpolation object"
+ << abort(FatalError);
+ }
+
+ if (interpPatches)
+ {
+ if (directMapping_)
+ {
+ result.boundaryField() == field.boundaryField();
+ }
+ else
+ {
+ notImplemented
+ (
+ "void Foam::meshToMeshNew::interpolate"
+ "("
+ "const GeometricField<Type, fvPatchField, volMesh>&, "
+ "const CombineOp&, "
+ "GeometricField<Type, fvPatchField, volMesh>&, "
+ "const bool"
+ ") const - non-conformal patches"
+ )
+
+ // do something...
+ }
+ }
+}
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const CombineOp& cop,
+ const bool interpPatches
+) const
+{
+ typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+ const fvMesh& mesh = field.mesh();
+
+ tmp<fieldType> tresult;
+
+ if (mesh.name() == srcRegionName_)
+ {
+ const fvMesh& tgtMesh =
+ mesh.time().lookupObject<fvMesh>(tgtRegionName_);
+
+ tresult =
+ new fieldType
+ (
+ IOobject
+ (
+ type() + "::interpolate(" + field.name() + ")",
+ tgtMesh.time().timeName(),
+ tgtMesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ tgtMesh,
+ dimensioned<Type>
+ (
+ "zero",
+ field.dimensions(),
+ pTraits<Type>::zero
+ )
+ );
+
+ interpolate(field, cop, tresult(), interpPatches);
+ }
+ else if (mesh.name() == tgtRegionName_)
+ {
+ const fvMesh& srcMesh =
+ mesh.time().lookupObject<fvMesh>(srcRegionName_);
+
+ tresult =
+ new fieldType
+ (
+ IOobject
+ (
+ type() + "::interpolate(" + field.name() + ")",
+ srcMesh.time().timeName(),
+ srcMesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ srcMesh,
+ dimensioned<Type>
+ (
+ "zero",
+ field.dimensions(),
+ pTraits<Type>::zero
+ )
+ );
+
+ interpolate(field, cop, tresult(), interpPatches);
+ }
+
+ return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+ const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
+ const CombineOp& cop,
+ const bool interpPatches
+) const
+{
+ return
+ interpolate
+ (
+ tfield(),
+ combineBinaryOp<Type, CombineOp>(cop),
+ interpPatches
+ );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+ const GeometricField<Type, fvPatchField, volMesh>& field,
+ const bool interpPatches
+) const
+{
+ return interpolate(field, plusEqOp<Type>(), interpPatches);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+ const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
+ const bool interpPatches
+) const
+{
+ return interpolate(tfield(), plusEqOp<Type>(), interpPatches);
+}
+
+
+// ************************************************************************* //
--
GitLab