applications/utilities/preProcessing/mapFields/mapLagrangian.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "MapLagrangianFields.H"
  27 #include "CloudTemplate.H"
  28 #include "passiveParticle.H"
  29 #include "meshSearch.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 static const scalar perturbFactor = 1E-6;
  37
  38
  39 // Special version of findCell that generates a cell guaranteed to be
  40 // compatible with tracking.
  41 static label findCell(const meshSearch& meshSearcher, const point& pt)
  42 {
  43     const polyMesh& mesh = meshSearcher.mesh();
  44
  45     // Use tracking to find cell containing pt
  46     label cellI = meshSearcher.findCell(pt);
  47
  48     if (cellI >= 0)
  49     {
  50         return cellI;
  51     }
  52     else
  53     {
  54         // See if particle on face by finding nearest face and shifting
  55         // particle.
  56
  57         label faceI = meshSearcher.findNearestBoundaryFace(pt);
  58
  59         if (faceI >= 0)
  60         {
  61             const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
  62             const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
  63
  64             return meshSearcher.findCell(perturbPt);
  65         }
  66     }
  67     return -1;
  68 }
  69
  70
  71 void mapLagrangian(const meshToMesh& meshToMeshInterp)
  72 {
  73     // Determine which particles are in meshTarget
  74     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  75
  76     // target to source cell map
  77     const labelList& cellAddressing = meshToMeshInterp.cellAddressing();
  78
  79     // Invert celladdressing to get source to target(s).
  80     // Note: could use sparse addressing but that is too storage inefficient
  81     // (Map<labelList>)
  82     labelListList sourceToTargets
  83     (
  84         invertOneToMany(meshToMeshInterp.fromMesh().nCells(), cellAddressing)
  85     );
  86
  87     const fvMesh& meshSource = meshToMeshInterp.fromMesh();
  88     const fvMesh& meshTarget = meshToMeshInterp.toMesh();
  89     const pointField& targetCc = meshTarget.cellCentres();
  90
  91
  92     fileNameList cloudDirs
  93     (
  94         readDir
  95         (
  96             meshSource.time().timePath()/cloud::prefix,
  97             fileName::DIRECTORY
  98         )
  99     );
 100
 101     forAll(cloudDirs, cloudI)
 102     {
 103         // Search for list of lagrangian objects for this time
 104         IOobjectList objects
 105         (
 106             meshSource,
 107             meshSource.time().timeName(),
 108             cloud::prefix/cloudDirs[cloudI]
 109         );
 110
 111         IOobject* positionsPtr = objects.lookup("positions");
 112
 113         if (positionsPtr)
 114         {
 115             Info<< nl << "    processing cloud " << cloudDirs[cloudI] << endl;
 116
 117             // Read positions & cell
 118             Cloud<passiveParticle> sourceParcels
 119             (
 120                 meshSource,
 121                 cloudDirs[cloudI],
 122                 false
 123             );
 124             Info<< "    read " << sourceParcels.size()
 125                 << " parcels from source mesh." << endl;
 126
 127             // Construct empty target cloud
 128             Cloud<passiveParticle> targetParcels
 129             (
 130                 meshTarget,
 131                 cloudDirs[cloudI],
 132                 IDLList<passiveParticle>()
 133             );
 134
 135             label sourceParticleI = 0;
 136
 137             // Indices of source particles that get added to targetParcels
 138             DynamicList<label> addParticles(sourceParcels.size());
 139
 140             // Unmapped particles
 141             labelHashSet unmappedSource(sourceParcels.size());
 142
 143
 144             // Initial: track from fine-mesh cell centre to particle position
 145             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 146             // This requires there to be no boundary in the way.
 147
 148
 149             forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
 150             {
 151                 bool foundCell = false;
 152
 153                 // Assume that cell from read parcel is the correct one...
 154                 if (iter().cell() >= 0)
 155                 {
 156                     const labelList& targetCells =
 157                         sourceToTargets[iter().cell()];
 158
 159                     // Particle probably in one of the targetcells. Try
 160                     // all by tracking from their cell centre to the parcel
 161                     // position.
 162
 163                     forAll(targetCells, i)
 164                     {
 165                         // Track from its cellcentre to position to make sure.
 166                         autoPtr<passiveParticle> newPtr
 167                         (
 168                             new passiveParticle
 169                             (
 170                                 targetParcels,
 171                                 targetCc[targetCells[i]],
 172                                 targetCells[i]
 173                             )
 174                         );
 175                         passiveParticle& newP = newPtr();
 176
 177                         scalar fraction = 0;
 178                         label faceI = newP.track(iter().position(), fraction);
 179
 180                         if (faceI < 0 && newP.cell() >= 0)
 181                         {
 182                             // Hit position.
 183                             foundCell = true;
 184                             addParticles.append(sourceParticleI);
 185                             targetParcels.addParticle(newPtr.ptr());
 186                             break;
 187                         }
 188                     }
 189                 }
 190
 191                 if (!foundCell)
 192                 {
 193                     // Store for closer analysis
 194                     unmappedSource.insert(sourceParticleI);
 195                 }
 196
 197                 sourceParticleI++;
 198             }
 199
 200             Info<< "    after meshToMesh addressing found "
 201                 << targetParcels.size()
 202                 << " parcels in target mesh." << endl;
 203
 204
 205             // Do closer inspection for unmapped particles
 206             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 207
 208             if (unmappedSource.size())
 209             {
 210                 meshSearch targetSearcher(meshTarget, false);
 211
 212                 sourceParticleI = 0;
 213
 214                 forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
 215                 {
 216                     if (unmappedSource.found(sourceParticleI))
 217                     {
 218                         label targetCell =
 219                             findCell(targetSearcher, iter().position());
 220
 221                         if (targetCell >= 0)
 222                         {
 223                             unmappedSource.erase(sourceParticleI);
 224                             addParticles.append(sourceParticleI);
 225                             iter().cell()=targetCell;
 226                             targetParcels.addParticle
 227                             (
 228                                 sourceParcels.remove(&iter())
 229                             );
 230                         }
 231                     }
 232                     sourceParticleI++;
 233                 }
 234             }
 235             addParticles.shrink();
 236
 237             Info<< "    after additional mesh searching found "
 238                 << targetParcels.size() << " parcels in target mesh." << endl;
 239
 240             if (addParticles.size())
 241             {
 242                 IOPosition<passiveParticle>(targetParcels).write();
 243
 244                 // addParticles now contains the indices of the sourceMesh
 245                 // particles that were appended to the target mesh.
 246
 247                 // Map lagrangian fields
 248                 // ~~~~~~~~~~~~~~~~~~~~~
 249
 250                 MapLagrangianFields<label>
 251                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 252                 MapLagrangianFields<scalar>
 253                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 254                 MapLagrangianFields<vector>
 255                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 256                 MapLagrangianFields<sphericalTensor>
 257                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 258                 MapLagrangianFields<symmTensor>
 259                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 260                 MapLagrangianFields<tensor>
 261                 (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
 262             }
 263         }
 264     }
 265 }
 266
 267
 268 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 269
 270 } // End namespace Foam
 271
 272 // ************************************************************************* //