From 03b0619ee1cf07f1f263de748fab795016967a23 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 11 Oct 2022 08:09:25 +0100 Subject: [PATCH] lagrangian: Support meshToMesh mapping Lagrangian is now compatible with the meshToMesh topology changer. If a cloud is being simulated and this topology changer is active, then the cloud data will be automatically mapped between the specified sequence of meshes in the same way as the finite volume data. This works both for serial and parallel simulations. In addition, mapFieldsPar now also supports mapping of Lagrangian data when run in parallel. --- .../preProcessing/mapFields/Make/files | 1 + .../preProcessing/mapFields/Make/options | 4 +- .../preProcessing/mapFields/meshToMesh0.C | 762 +++++++++++++++ .../preProcessing/mapFields}/meshToMesh0.H | 0 .../mapFields}/meshToMesh0Templates.C | 0 .../preProcessing/mapFieldsPar/Make/files | 4 +- .../preProcessing/mapFieldsPar/Make/options | 4 +- .../mapFieldsPar/MapLagrangianFields.H | 184 ---- .../preProcessing/mapFieldsPar/mapClouds.C | 287 ++++++ .../mapFieldsPar/{UnMapped.H => mapClouds.H} | 34 +- .../preProcessing/mapFieldsPar/mapFieldsPar.C | 94 +- .../{MapVolFields.H => mapGeometricFields.C} | 91 +- .../{mapLagrangian.H => mapGeometricFields.H} | 22 +- .../mapFieldsPar/mapLagrangian.C | 319 ------- .../preProcessing/mapFieldsPar/mapMeshes.C | 169 ---- .../meshes/polyMesh/polyMeshMap/polyMeshMap.C | 9 +- .../meshes/polyMesh/polyMeshMap/polyMeshMap.H | 13 +- src/finiteVolume/Make/files | 3 + .../fvMeshToFvMesh/fvMeshToFvMesh.C | 50 +- .../fvMeshToFvMesh/fvMeshToFvMesh.H | 109 +++ .../fvMeshToFvMesh/fvMeshToFvMeshTemplates.C | 171 ++++ .../patchToPatchFvPatchFieldMapper.C | 0 .../patchToPatchFvPatchFieldMapper.H | 0 .../patchToPatchFvPatchFieldMapperTemplates.C | 0 .../meshToMesh/Make/options | 2 - .../meshToMesh/MeshToMeshMapGeometricFields.H | 23 +- .../meshToMesh/fvMeshTopoChangersMeshToMesh.C | 18 +- src/lagrangian/basic/Cloud/Cloud.C | 105 ++- src/meshTools/Make/files | 10 + .../cellVolumeWeight/cellVolumeWeightMethod.C | 0 .../cellVolumeWeight/cellVolumeWeightMethod.H | 0 .../calcMethod/direct/directMethod.C | 0 .../calcMethod/direct/directMethod.H | 0 .../calcMethod/mapNearest/mapNearestMethod.C | 0 .../calcMethod/mapNearest/mapNearestMethod.H | 0 .../meshToMeshMethod/meshToMeshMethod.C | 0 .../meshToMeshMethod/meshToMeshMethod.H | 0 .../meshToMeshMethod/meshToMeshMethodNew.C | 0 src/meshTools/meshToMesh/meshToMesh.C | 414 +++++++++ .../meshToMesh/meshToMesh.H | 188 ++-- .../meshToMesh/meshToMeshI.H | 77 +- .../meshToMesh/meshToMeshParallelOps.C | 875 ++++++++++++++++++ .../meshToMesh/meshToMeshTemplates.C | 200 ++++ .../intersection/intersectionPatchToPatch.C | 11 +- .../intersection/intersectionPatchToPatch.H | 4 +- .../inverseDistancePatchToPatch.C | 11 +- .../inverseDistancePatchToPatch.H | 4 +- .../nearest/nearestPatchToPatch.C | 11 +- .../nearest/nearestPatchToPatch.H | 4 +- .../patchToPatch/patchToPatch/patchToPatch.C | 212 +---- .../patchToPatch/patchToPatch/patchToPatch.H | 83 +- .../patchToPatch/patchToPatch/patchToPatchI.H | 74 +- .../patchToPatch/patchToPatchParallelOps.C | 272 +----- .../patchToPatch/rays/raysPatchToPatch.C | 10 +- .../patchToPatch/rays/raysPatchToPatch.H | 2 +- .../patchToPatchTools/patchToPatchTools.C | 402 ++++++++ .../patchToPatchTools/patchToPatchTools.H | 153 +++ .../patchToPatchToolsTemplates.C | 85 ++ src/sampling/Make/files | 14 - src/sampling/meshToMesh/meshToMesh.C | 462 --------- .../meshToMesh/meshToMeshParallelOps.C | 833 ----------------- src/sampling/meshToMesh/meshToMeshTemplates.C | 344 ------- .../calculateMeshToMesh0Addressing.C | 344 ------- .../meshToMesh0/calculateMeshToMesh0Weights.C | 271 ------ src/sampling/meshToMesh0/meshToMesh0.C | 202 ---- .../lagrangian/particleFoam/hopper/Allclean | 13 +- .../lagrangian/particleFoam/hopper/Allrun | 29 +- .../hopper/hopperEmptying/{0.orig => 0}/U | 0 .../hopper/hopperEmptying/{0.orig => 0}/mu | 0 .../hopper/hopperEmptying/{0.orig => 0}/rho | 0 70 files changed, 4079 insertions(+), 4013 deletions(-) create mode 100644 applications/utilities/preProcessing/mapFields/meshToMesh0.C rename {src/sampling/meshToMesh0 => applications/utilities/preProcessing/mapFields}/meshToMesh0.H (100%) rename {src/sampling/meshToMesh0 => applications/utilities/preProcessing/mapFields}/meshToMesh0Templates.C (100%) delete mode 100644 applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H create mode 100644 applications/utilities/preProcessing/mapFieldsPar/mapClouds.C rename applications/utilities/preProcessing/mapFieldsPar/{UnMapped.H => mapClouds.H} (73%) rename applications/utilities/preProcessing/mapFieldsPar/{MapVolFields.H => mapGeometricFields.C} (64%) rename applications/utilities/preProcessing/mapFieldsPar/{mapLagrangian.H => mapGeometricFields.H} (80%) delete mode 100644 applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C delete mode 100644 applications/utilities/preProcessing/mapFieldsPar/mapMeshes.C rename applications/utilities/preProcessing/mapFieldsPar/mapMeshes.H => src/finiteVolume/fvMeshToFvMesh/fvMeshToFvMesh.C (55%) create mode 100644 src/finiteVolume/fvMeshToFvMesh/fvMeshToFvMesh.H create mode 100644 src/finiteVolume/fvMeshToFvMesh/fvMeshToFvMeshTemplates.C rename src/{sampling/meshToMesh => finiteVolume/fvMeshToFvMesh}/patchToPatchFvPatchFieldMapper.C (100%) rename src/{sampling/meshToMesh => finiteVolume/fvMeshToFvMesh}/patchToPatchFvPatchFieldMapper.H (100%) rename src/{sampling/meshToMesh => finiteVolume/fvMeshToFvMesh}/patchToPatchFvPatchFieldMapperTemplates.C (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/direct/directMethod.C (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/direct/directMethod.H (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H (100%) rename src/{sampling => meshTools}/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C (100%) create mode 100644 src/meshTools/meshToMesh/meshToMesh.C rename src/{sampling => meshTools}/meshToMesh/meshToMesh.H (56%) rename src/{sampling => meshTools}/meshToMesh/meshToMeshI.H (66%) create mode 100644 src/meshTools/meshToMesh/meshToMeshParallelOps.C create mode 100644 src/meshTools/meshToMesh/meshToMeshTemplates.C create mode 100644 src/meshTools/patchToPatchTools/patchToPatchTools.C create mode 100644 src/meshTools/patchToPatchTools/patchToPatchTools.H create mode 100644 src/meshTools/patchToPatchTools/patchToPatchToolsTemplates.C delete mode 100644 src/sampling/meshToMesh/meshToMesh.C delete mode 100644 src/sampling/meshToMesh/meshToMeshParallelOps.C delete mode 100644 src/sampling/meshToMesh/meshToMeshTemplates.C delete mode 100644 src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C delete mode 100644 src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C delete mode 100644 src/sampling/meshToMesh0/meshToMesh0.C rename tutorials/lagrangian/particleFoam/hopper/hopperEmptying/{0.orig => 0}/U (100%) rename tutorials/lagrangian/particleFoam/hopper/hopperEmptying/{0.orig => 0}/mu (100%) rename tutorials/lagrangian/particleFoam/hopper/hopperEmptying/{0.orig => 0}/rho (100%) 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