diff --git a/applications/utilities/preProcessing/mapFields/Make/files b/applications/utilities/preProcessing/mapFields/Make/files index 4503a5fa98..6b11f80616 100644 --- a/applications/utilities/preProcessing/mapFields/Make/files +++ b/applications/utilities/preProcessing/mapFields/Make/files @@ -1,5 +1,6 @@ mapLagrangian.C mapMeshes.C mapFields.C +meshToMesh0.C EXE = $(FOAM_APPBIN)/mapFields diff --git a/applications/utilities/preProcessing/mapFields/Make/options b/applications/utilities/preProcessing/mapFields/Make/options index 2f93ae2d61..518b236496 100644 --- a/applications/utilities/preProcessing/mapFields/Make/options +++ b/applications/utilities/preProcessing/mapFields/Make/options @@ -2,12 +2,10 @@ EXE_INC = \ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -ldecompositionMethods \ - -lsampling \ -lmeshTools \ -llagrangian \ -lfiniteVolume \ diff --git a/applications/utilities/preProcessing/mapFields/meshToMesh0.C b/applications/utilities/preProcessing/mapFields/meshToMesh0.C new file mode 100644 index 0000000000..2569e11ea0 --- /dev/null +++ b/applications/utilities/preProcessing/mapFields/meshToMesh0.C @@ -0,0 +1,762 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2022 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 . + +\*---------------------------------------------------------------------------*/ + +#include "meshToMesh0.H" +#include "processorFvPatch.H" +#include "demandDrivenData.H" +#include "treeDataCell.H" +#include "treeDataFace.H" +#include "tetOverlapVolume.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +defineTypeNameAndDebug(meshToMesh0, 0); +} + +const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-5; + + + + + + +void Foam::meshToMesh0::calcAddressing() +{ + if (debug) + { + InfoInFunction + << "Calculating mesh-to-mesh cell addressing" << endl; + } + + // Set reference to cells + const cellList& fromCells = fromMesh_.cells(); + const pointField& fromPoints = fromMesh_.points(); + + // In an attempt to preserve the efficiency of linear search and prevent + // failure, a RESCUE mechanism will be set up. Here, we shall mark all + // cells next to the solid boundaries. If such a cell is found as the + // closest, the relationship between the origin and cell will be examined. + // If the origin is outside the cell, a global n-squared search is + // triggered. + + // SETTING UP RESCUE + + // Visit all boundaries and mark the cell next to the boundary. + + if (debug) + { + InfoInFunction << "Setting up rescue" << endl; + } + + List boundaryCell(fromCells.size(), false); + + // Set reference to boundary + const polyPatchList& patchesFrom = fromMesh_.boundaryMesh(); + + forAll(patchesFrom, patchi) + { + // Get reference to cells next to the boundary + const labelUList& bCells = patchesFrom[patchi].faceCells(); + + forAll(bCells, facei) + { + boundaryCell[bCells[facei]] = true; + } + } + + treeBoundBox meshBb(fromPoints); + + scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size()))); + + treeBoundBox shiftedBb + ( + meshBb.min(), + meshBb.max() + vector(typDim, typDim, typDim) + ); + + if (debug) + { + Info<< "\nMesh\n" + << " bounding box : " << meshBb << nl + << " bounding box (shifted) : " << shiftedBb << nl + << " typical dimension :" << shiftedBb.typDim() << endl; + } + + indexedOctree oc + ( + treeDataCell(false, fromMesh_, polyMesh::CELL_TETS), + shiftedBb, // overall bounding box + 8, // maxLevel + 10, // leafsize + 6.0 // duplicity + ); + + if (debug) + { + oc.print(Pout, false, 0); + } + + cellAddresses + ( + cellAddressing_, + toMesh_.cellCentres(), + fromMesh_, + boundaryCell, + oc + ); + + forAll(toMesh_.boundaryMesh(), patchi) + { + const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi]; + + if (cuttingPatches_.found(toPatch.name())) + { + boundaryAddressing_[patchi].setSize(toPatch.size()); + + cellAddresses + ( + boundaryAddressing_[patchi], + toPatch.faceCentres(), + fromMesh_, + boundaryCell, + oc + ); + } + else if + ( + patchMap_.found(toPatch.name()) + && fromMeshPatches_.found(patchMap_.find(toPatch.name())()) + ) + { + const polyPatch& fromPatch = fromMesh_.boundaryMesh() + [ + fromMeshPatches_.find(patchMap_.find(toPatch.name())())() + ]; + + if (fromPatch.empty()) + { + WarningInFunction + << "Source patch " << fromPatch.name() + << " has no faces. Not performing mapping for it." + << endl; + boundaryAddressing_[patchi].setSize(toPatch.size()); + boundaryAddressing_[patchi] = -1; + } + else + { + treeBoundBox wallBb(fromPatch.localPoints()); + scalar typDim = + wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size()))); + + treeBoundBox shiftedBb + ( + wallBb.min(), + wallBb.max() + vector(typDim, typDim, typDim) + ); + + // Note: allow more levels than in meshSearch. Assume patch + // is not as big as all boundary faces + indexedOctree oc + ( + treeDataFace(false, fromPatch), + shiftedBb, // overall search domain + 12, // maxLevel + 10, // leafsize + 6.0 // duplicity + ); + + const vectorField::subField centresToBoundary = + toPatch.faceCentres(); + + boundaryAddressing_[patchi].setSize(toPatch.size()); + + const scalar distSqr = sqr(wallBb.mag()); + + forAll(toPatch, toi) + { + boundaryAddressing_[patchi][toi] = oc.findNearest + ( + centresToBoundary[toi], + distSqr + ).index(); + } + } + } + } + + if (debug) + { + InfoInFunction + << "Finished calculating mesh-to-mesh cell addressing" << endl; + } +} + + +void Foam::meshToMesh0::cellAddresses +( + labelList& cellAddressing_, + const pointField& points, + const fvMesh& fromMesh, + const List& boundaryCell, + const indexedOctree& oc +) const +{ + // The implemented search method is a simple neighbour array search. + // It starts from a cell zero, searches its neighbours and finds one + // which is nearer to the target point than the current position. + // The location of the "current position" is reset to that cell and + // search through the neighbours continues. The search is finished + // when all the neighbours of the cell are farther from the target + // point than the current cell + + // Set curCell label to zero (start) + label curCell = 0; + + // Set reference to cell to cell addressing + const vectorField& centresFrom = fromMesh.cellCentres(); + const labelListList& cc = fromMesh.cellCells(); + + forAll(points, toi) + { + // Pick up target position + const vector& p = points[toi]; + + // Set the sqr-distance + scalar distSqr = magSqr(p - centresFrom[curCell]); + + bool closer; + + do + { + closer = false; + + // Set the current list of neighbouring cells + const labelList& neighbours = cc[curCell]; + + forAll(neighbours, ni) + { + const scalar curDistSqr = + magSqr(p - centresFrom[neighbours[ni]]); + + // Search through all the neighbours. + // If the cell is closer, reset current cell and distance + if (curDistSqr < (1 - small)*distSqr) + { + curCell = neighbours[ni]; + distSqr = curDistSqr; + closer = true; // A closer neighbour has been found + } + } + } while (closer); + + cellAddressing_[toi] = -1; + + // Check point is actually in the nearest cell + if (fromMesh.pointInCell(p, curCell)) + { + cellAddressing_[toi] = curCell; + } + else + { + // If curCell is a boundary cell then the point maybe either outside + // the domain or in an other region of the doamin, either way use + // the octree search to find it. + if (boundaryCell[curCell]) + { + cellAddressing_[toi] = oc.findInside(p); + if (cellAddressing_[toi] != -1) + { + curCell = cellAddressing_[toi]; + } + } + else + { + // If not on the boundary search the neighbours + bool found = false; + + // Set the current list of neighbouring cells + const labelList& neighbours = cc[curCell]; + + forAll(neighbours, ni) + { + // Search through all the neighbours. + // If point is in neighbour reset current cell + if (fromMesh.pointInCell(p, neighbours[ni])) + { + cellAddressing_[toi] = neighbours[ni]; + found = true; + break; + } + } + + if (!found) + { + // If still not found search the neighbour-neighbours + + // Set the current list of neighbouring cells + const labelList& neighbours = cc[curCell]; + + forAll(neighbours, ni) + { + // Set the current list of neighbour-neighbouring cells + const labelList& nn = cc[neighbours[ni]]; + + forAll(nn, ni) + { + // Search through all the neighbours. + // If point is in neighbour reset current cell + if (fromMesh.pointInCell(p, nn[ni])) + { + cellAddressing_[toi] = nn[ni]; + found = true; + break; + } + } + if (found) break; + } + } + + if (!found) + { + // Still not found so use the octree + cellAddressing_[toi] = oc.findInside(p); + + if (cellAddressing_[toi] != -1) + { + curCell = cellAddressing_[toi]; + } + } + } + } + } +} + + +void Foam::meshToMesh0::calculateInverseDistanceWeights() const +{ + if (debug) + { + InfoInFunction + << "Calculating inverse distance weighting factors" << endl; + } + + if (inverseDistanceWeightsPtr_) + { + FatalErrorInFunction + << "weighting factors already calculated" + << exit(FatalError); + } + + //- Initialise overlap volume to zero + V_ = 0.0; + + inverseDistanceWeightsPtr_ = new scalarListList(toMesh_.nCells()); + scalarListList& invDistCoeffs = *inverseDistanceWeightsPtr_; + + // get reference to source mesh data + const labelListList& cc = fromMesh_.cellCells(); + const vectorField& centreFrom = fromMesh_.C(); + const vectorField& centreTo = toMesh_.C(); + + forAll(cellAddressing_, celli) + { + if (cellAddressing_[celli] != -1) + { + const vector& target = centreTo[celli]; + scalar m = mag(target - centreFrom[cellAddressing_[celli]]); + + const labelList& neighbours = cc[cellAddressing_[celli]]; + + // if the nearest cell is a boundary cell or there is a direct hit, + // pick up the value + label directCelli = -1; + if (m < directHitTol || neighbours.empty()) + { + directCelli = celli; + } + else + { + forAll(neighbours, ni) + { + scalar nm = mag(target - centreFrom[neighbours[ni]]); + if (nm < directHitTol) + { + directCelli = neighbours[ni]; + break; + } + } + } + + + if (directCelli != -1) + { + // Direct hit + invDistCoeffs[directCelli].setSize(1); + invDistCoeffs[directCelli][0] = 1.0; + V_ += fromMesh_.V()[cellAddressing_[directCelli]]; + } + else + { + invDistCoeffs[celli].setSize(neighbours.size() + 1); + + // The first coefficient corresponds to the centre cell. + // The rest is ordered in the same way as the cellCells list. + scalar invDist = 1.0/m; + invDistCoeffs[celli][0] = invDist; + scalar sumInvDist = invDist; + + // now add the neighbours + forAll(neighbours, ni) + { + invDist = 1.0/mag(target - centreFrom[neighbours[ni]]); + invDistCoeffs[celli][ni + 1] = invDist; + sumInvDist += invDist; + } + + // divide by the total inverse-distance + forAll(invDistCoeffs[celli], i) + { + invDistCoeffs[celli][i] /= sumInvDist; + } + + + V_ += + invDistCoeffs[celli][0] + *fromMesh_.V()[cellAddressing_[celli]]; + for (label i = 1; i < invDistCoeffs[celli].size(); i++) + { + V_ += + invDistCoeffs[celli][i]*fromMesh_.V()[neighbours[i-1]]; + } + } + } + } +} + + +void Foam::meshToMesh0::calculateInverseVolumeWeights() const +{ + if (debug) + { + InfoInFunction + << "Calculating inverse volume weighting factors" << endl; + } + + if (inverseVolumeWeightsPtr_) + { + FatalErrorInFunction + << "weighting factors already calculated" + << exit(FatalError); + } + + //- Initialise overlap volume to zero + V_ = 0.0; + + inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells()); + scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_; + + const labelListList& cellToCell = cellToCellAddressing(); + + tetOverlapVolume overlapEngine; + + forAll(cellToCell, celli) + { + const labelList& overlapCells = cellToCell[celli]; + + if (overlapCells.size() > 0) + { + invVolCoeffs[celli].setSize(overlapCells.size()); + + forAll(overlapCells, j) + { + label cellFrom = overlapCells[j]; + treeBoundBox bbFromMesh + ( + pointField + ( + fromMesh_.points(), + fromMesh_.cellPoints()[cellFrom] + ) + ); + + scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp + ( + toMesh_, + celli, + + fromMesh_, + cellFrom, + bbFromMesh + ); + invVolCoeffs[celli][j] = v/toMesh_.V()[celli]; + + V_ += v; + } + } + } +} + + +void Foam::meshToMesh0::calculateCellToCellAddressing() const +{ + if (debug) + { + InfoInFunction + << "Calculating cell to cell addressing" << endl; + } + + if (cellToCellAddressingPtr_) + { + FatalErrorInFunction + << "addressing already calculated" + << exit(FatalError); + } + + //- Initialise overlap volume to zero + V_ = 0.0; + + tetOverlapVolume overlapEngine; + + cellToCellAddressingPtr_ = new labelListList(toMesh_.nCells()); + labelListList& cellToCell = *cellToCellAddressingPtr_; + + + forAll(cellToCell, iTo) + { + const labelList overLapCells = + overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo); + if (overLapCells.size() > 0) + { + cellToCell[iTo].setSize(overLapCells.size()); + forAll(overLapCells, j) + { + cellToCell[iTo][j] = overLapCells[j]; + V_ += fromMesh_.V()[overLapCells[j]]; + } + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::meshToMesh0::meshToMesh0 +( + const fvMesh& meshFrom, + const fvMesh& meshTo, + const HashTable& patchMap, + const wordList& cuttingPatchNames +) +: + fromMesh_(meshFrom), + toMesh_(meshTo), + patchMap_(patchMap), + cellAddressing_(toMesh_.nCells()), + boundaryAddressing_(toMesh_.boundaryMesh().size()), + inverseDistanceWeightsPtr_(nullptr), + inverseVolumeWeightsPtr_(nullptr), + cellToCellAddressingPtr_(nullptr), + V_(0.0) +{ + forAll(fromMesh_.boundaryMesh(), patchi) + { + fromMeshPatches_.insert + ( + fromMesh_.boundaryMesh()[patchi].name(), + patchi + ); + } + + forAll(toMesh_.boundaryMesh(), patchi) + { + toMeshPatches_.insert + ( + toMesh_.boundaryMesh()[patchi].name(), + patchi + ); + } + + forAll(cuttingPatchNames, i) + { + if (toMeshPatches_.found(cuttingPatchNames[i])) + { + cuttingPatches_.insert + ( + cuttingPatchNames[i], + toMeshPatches_.find(cuttingPatchNames[i])() + ); + } + else + { + WarningInFunction + << "Cannot find cutting-patch " << cuttingPatchNames[i] + << " in destination mesh" << endl; + } + } + + forAll(toMesh_.boundaryMesh(), patchi) + { + // Add the processor patches in the toMesh to the cuttingPatches list + if (isA(toMesh_.boundaryMesh()[patchi])) + { + cuttingPatches_.insert + ( + toMesh_.boundaryMesh()[patchi].name(), + patchi + ); + } + } + + calcAddressing(); +} + + +Foam::meshToMesh0::meshToMesh0 +( + const fvMesh& meshFrom, + const fvMesh& meshTo +) +: + fromMesh_(meshFrom), + toMesh_(meshTo), + cellAddressing_(toMesh_.nCells()), + boundaryAddressing_(toMesh_.boundaryMesh().size()), + inverseDistanceWeightsPtr_(nullptr), + inverseVolumeWeightsPtr_(nullptr), + cellToCellAddressingPtr_(nullptr), + V_(0.0) +{ + // check whether both meshes have got the same number + // of boundary patches + if (fromMesh_.boundary().size() != toMesh_.boundary().size()) + { + FatalErrorInFunction + << "Incompatible meshes: different number of patches, " + << "fromMesh = " << fromMesh_.boundary().size() + << ", toMesh = " << toMesh_.boundary().size() + << exit(FatalError); + } + + forAll(fromMesh_.boundaryMesh(), patchi) + { + if + ( + fromMesh_.boundaryMesh()[patchi].name() + != toMesh_.boundaryMesh()[patchi].name() + ) + { + FatalErrorInFunction + << "Incompatible meshes: different patch names for patch " + << patchi + << ", fromMesh = " << fromMesh_.boundary()[patchi].name() + << ", toMesh = " << toMesh_.boundary()[patchi].name() + << exit(FatalError); + } + + if + ( + fromMesh_.boundaryMesh()[patchi].type() + != toMesh_.boundaryMesh()[patchi].type() + ) + { + FatalErrorInFunction + << "Incompatible meshes: different patch types for patch " + << patchi + << ", fromMesh = " << fromMesh_.boundary()[patchi].type() + << ", toMesh = " << toMesh_.boundary()[patchi].type() + << exit(FatalError); + } + + fromMeshPatches_.insert + ( + fromMesh_.boundaryMesh()[patchi].name(), + patchi + ); + + toMeshPatches_.insert + ( + toMesh_.boundaryMesh()[patchi].name(), + patchi + ); + + patchMap_.insert + ( + toMesh_.boundaryMesh()[patchi].name(), + fromMesh_.boundaryMesh()[patchi].name() + ); + } + + calcAddressing(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::meshToMesh0::~meshToMesh0() +{ + deleteDemandDrivenData(inverseDistanceWeightsPtr_); + deleteDemandDrivenData(inverseVolumeWeightsPtr_); + deleteDemandDrivenData(cellToCellAddressingPtr_); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::scalarListList& Foam::meshToMesh0::inverseDistanceWeights() const +{ + if (!inverseDistanceWeightsPtr_) + { + calculateInverseDistanceWeights(); + } + + return *inverseDistanceWeightsPtr_; +} + + +const Foam::scalarListList& Foam::meshToMesh0::inverseVolumeWeights() const +{ + if (!inverseVolumeWeightsPtr_) + { + calculateInverseVolumeWeights(); + } + + return *inverseVolumeWeightsPtr_; +} + + +const Foam::labelListList& Foam::meshToMesh0::cellToCellAddressing() const +{ + if (!cellToCellAddressingPtr_) + { + calculateCellToCellAddressing(); + } + + return *cellToCellAddressingPtr_; +} + + +// ************************************************************************* // diff --git a/src/sampling/meshToMesh0/meshToMesh0.H b/applications/utilities/preProcessing/mapFields/meshToMesh0.H similarity index 100% rename from src/sampling/meshToMesh0/meshToMesh0.H rename to applications/utilities/preProcessing/mapFields/meshToMesh0.H diff --git a/src/sampling/meshToMesh0/meshToMesh0Templates.C b/applications/utilities/preProcessing/mapFields/meshToMesh0Templates.C similarity index 100% rename from src/sampling/meshToMesh0/meshToMesh0Templates.C rename to applications/utilities/preProcessing/mapFields/meshToMesh0Templates.C diff --git a/applications/utilities/preProcessing/mapFieldsPar/Make/files b/applications/utilities/preProcessing/mapFieldsPar/Make/files index c2519449c4..83b9c4c454 100644 --- a/applications/utilities/preProcessing/mapFieldsPar/Make/files +++ b/applications/utilities/preProcessing/mapFieldsPar/Make/files @@ -1,5 +1,5 @@ -mapMeshes.C -mapLagrangian.C +mapGeometricFields.C +mapClouds.C mapFieldsPar.C EXE = $(FOAM_APPBIN)/mapFieldsPar diff --git a/applications/utilities/preProcessing/mapFieldsPar/Make/options b/applications/utilities/preProcessing/mapFieldsPar/Make/options index 7bd964094e..f597efdd05 100644 --- a/applications/utilities/preProcessing/mapFieldsPar/Make/options +++ b/applications/utilities/preProcessing/mapFieldsPar/Make/options @@ -1,11 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ - -lsampling \ -lmeshTools \ -llagrangian \ -lfiniteVolume \ diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H b/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H deleted file mode 100644 index 21327e69b2..0000000000 --- a/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H +++ /dev/null @@ -1,184 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 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 . - -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 -void MapLagrangianFields -( - const string& cloudName, - const IOobjectList& objects, - const polyMesh& meshTarget, - const labelList& addParticles -) -{ - { - IOobjectList fields = objects.lookupClass(IOField::typeName); - - forAllIter(IOobjectList, fields, fieldIter) - { - const word& fieldName = fieldIter()->name(); - - Info<< " mapping lagrangian field " << fieldName << endl; - - // Read field (does not need mesh) - IOField fieldSource(*fieldIter()); - - // Map - IOField 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>::typeName); - - forAllIter(IOobjectList, fieldFields, fieldIter) - { - const word& fieldName = fieldIter()->name(); - - Info<< " mapping lagrangian fieldField " << fieldName << endl; - - // Read field (does not need mesh) - IOField> fieldSource(*fieldIter()); - - // Map - use CompactIOField to automatically write in - // compact form for binary format. - CompactIOField> 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>::typeName); - - forAllIter(IOobjectList, fieldFields, fieldIter) - { - Info<< " mapping lagrangian fieldField " - << fieldIter()->name() << endl; - - // Read field (does not need mesh) - CompactIOField> fieldSource(*fieldIter()); - - // Map - CompactIOField> 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 - -// ************************************************************************* // diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapClouds.C b/applications/utilities/preProcessing/mapFieldsPar/mapClouds.C new file mode 100644 index 0000000000..d1bed95210 --- /dev/null +++ b/applications/utilities/preProcessing/mapFieldsPar/mapClouds.C @@ -0,0 +1,287 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2022 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 . + +\*---------------------------------------------------------------------------*/ + +#include "mapClouds.H" +#include "fvMeshToFvMesh.H" +#include "IOobjectList.H" +#include "OSspecific.H" +#include "passiveParticleCloud.H" +#include "patchToPatchTools.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template +List gatherAndFlatten(const List& l) +{ + List> procLs(Pstream::nProcs()); + procLs[Pstream::myProcNo()] = l; + + Pstream::gatherList(procLs); + Pstream::scatterList(procLs); + + return + HashSet + ( + ListListOps::combine> + ( + procLs, + accessOp>() + ) + ).sortedToc(); +} + + +template +void mapCloudTypeFields +( + const fileName& cloudDir, + const IOobjectList& objects, + const polyMesh& srcMesh, + const polyMesh& tgtMesh, + const distributionMap& map +) +{ + const wordList fieldNames = + gatherAndFlatten(objects.lookupClass(ReadIOField::typeName).names()); + + forAll(fieldNames, fieldi) + { + const word& fieldName = fieldNames[fieldi]; + + Info<< " mapping lagrangian field " << fieldName << endl; + + const IOobject* fieldIOPtr = objects.lookup(fieldName); + + // Read the field + ReadIOField field + ( + fieldIOPtr != nullptr + ? *fieldIOPtr + : IOobject + ( + fieldName, + srcMesh.time().timeName(), + cloud::prefix/cloudDir, + srcMesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ) + ); + + // Distribute the field + map.distribute(field); + + // Write it out into the target case + if (field.size()) + { + WriteIOField + ( + IOobject + ( + fieldName, + tgtMesh.time().timeName(), + cloud::prefix/cloudDir, + tgtMesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + move(field) + ).write(); + } + } +} + + +template +void mapCloudTypeFieldsAndFieldFields +( + const fileName& cloudDir, + const IOobjectList& objects, + const polyMesh& srcMesh, + const polyMesh& tgtMesh, + const distributionMap& map +) +{ + mapCloudTypeFields + < + IOField, + IOField + > + ( + cloudDir, + objects, + srcMesh, + tgtMesh, + map + ); + + mapCloudTypeFields + < + IOField>, + CompactIOField> + > + ( + cloudDir, + objects, + srcMesh, + tgtMesh, + map + ); + + mapCloudTypeFields + < + CompactIOField>, + CompactIOField> + > + ( + cloudDir, + objects, + srcMesh, + tgtMesh, + map + ); +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void Foam::mapClouds(const fvMeshToFvMesh& interp) +{ + const polyMesh& srcMesh = interp.srcMesh(); + const polyMesh& tgtMesh = interp.tgtMesh(); + + // Determine the clouds present in this mesh + const fileNameList cloudDirs = + gatherAndFlatten + ( + readDir + ( + srcMesh.time().timePath()/cloud::prefix, + fileType::directory + ) + ); + + forAll(cloudDirs, cloudi) + { + Info<< nl << "Mapping cloud " << cloudDirs[cloudi] << endl; + + // Read the source cloud positions + passiveParticleCloud srcCloud(srcMesh, cloudDirs[cloudi], false); + + Info<< " read " + << returnReduce(srcCloud.size(), sumOp