src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "sampledTriSurfaceMesh.H"
  27 #include "treeDataPoint.H"
  28 #include "meshSearch.H"
  29 #include "Tuple2.H"
  30 #include "globalIndex.H"
  31
  32 #include "addToRunTimeSelectionTable.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(sampledTriSurfaceMesh, 0);
  39     addToRunTimeSelectionTable
  40     (
  41         sampledSurface,
  42         sampledTriSurfaceMesh,
  43         word
  44     );
  45
  46     //- Private class for finding nearest
  47     //  - global index
  48     //  - sqr(distance)
  49     typedef Tuple2<scalar, label> nearInfo;
  50
  51     class nearestEqOp
  52     {
  53
  54     public:
  55
  56         void operator()(nearInfo& x, const nearInfo& y) const
  57         {
  58             if (y.first() < x.first())
  59             {
  60                 x = y;
  61             }
  62         }
  63     };
  64 }
  65
  66
  67 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  68
  69 Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
  70 (
  71     const word& name,
  72     const polyMesh& mesh,
  73     const word& surfaceName
  74 )
  75 :
  76     sampledSurface(name, mesh),
  77     surface_
  78     (
  79         IOobject
  80         (
  81             name,
  82             mesh.time().constant(), // instance
  83             "triSurface",           // local
  84             mesh,                   // registry
  85             IOobject::MUST_READ,
  86             IOobject::NO_WRITE,
  87             false
  88         )
  89     ),
  90     needsUpdate_(true),
  91     cellLabels_(0),
  92     pointToFace_(0)
  93 {}
  94
  95
  96 Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
  97 (
  98     const word& name,
  99     const polyMesh& mesh,
 100     const dictionary& dict
 101 )
 102 :
 103     sampledSurface(name, mesh, dict),
 104     surface_
 105     (
 106         IOobject
 107         (
 108             dict.lookup("surface"),
 109             mesh.time().constant(), // instance
 110             "triSurface",           // local
 111             mesh,                   // registry
 112             IOobject::MUST_READ,
 113             IOobject::NO_WRITE,
 114             false
 115         )
 116     ),
 117     needsUpdate_(true),
 118     cellLabels_(0),
 119     pointToFace_(0)
 120 {}
 121
 122
 123 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 124
 125 Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
 126 {}
 127
 128
 129 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 130
 131 bool Foam::sampledTriSurfaceMesh::needsUpdate() const
 132 {
 133     return needsUpdate_;
 134 }
 135
 136
 137 bool Foam::sampledTriSurfaceMesh::expire()
 138 {
 139     // already marked as expired
 140     if (needsUpdate_)
 141     {
 142         return false;
 143     }
 144
 145     sampledSurface::clearGeom();
 146     MeshStorage::clear();
 147     cellLabels_.clear();
 148     pointToFace_.clear();
 149
 150     needsUpdate_ = true;
 151     return true;
 152 }
 153
 154
 155 bool Foam::sampledTriSurfaceMesh::update()
 156 {
 157     if (!needsUpdate_)
 158     {
 159         return false;
 160     }
 161
 162
 163     // Find the cells the triangles of the surface are in.
 164     // Does approximation by looking at the face centres only
 165     const pointField& fc = surface_.faceCentres();
 166
 167     meshSearch meshSearcher(mesh(), false);
 168
 169     const indexedOctree<treeDataPoint>& cellCentreTree =
 170         meshSearcher.cellCentreTree();
 171
 172
 173     // Global numbering for cells - only used to uniquely identify local cells.
 174     globalIndex globalCells(mesh().nCells());
 175     List<nearInfo> nearest(fc.size());
 176     forAll(nearest, i)
 177     {
 178         nearest[i].first() = GREAT;
 179         nearest[i].second() = labelMax;
 180     }
 181
 182     // Search triangles using nearest on local mesh
 183     forAll(fc, triI)
 184     {
 185         pointIndexHit nearInfo = cellCentreTree.findNearest
 186         (
 187             fc[triI],
 188             sqr(GREAT)
 189         );
 190         if (nearInfo.hit())
 191         {
 192             nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
 193             nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
 194         }
 195     }
 196
 197     // See which processor has the nearest.
 198     Pstream::listCombineGather(nearest, nearestEqOp());
 199     Pstream::listCombineScatter(nearest);
 200
 201     boolList include(surface_.size(), false);
 202
 203     cellLabels_.setSize(fc.size());
 204     cellLabels_ = -1;
 205
 206     label nFound = 0;
 207     forAll(nearest, triI)
 208     {
 209         if (nearest[triI].second() == labelMax)
 210         {
 211             // Not found on any processor. How to map?
 212         }
 213         else if (globalCells.isLocal(nearest[triI].second()))
 214         {
 215             cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
 216
 217             include[triI] = true;
 218             nFound++;
 219         }
 220     }
 221
 222
 223     if (debug)
 224     {
 225         Pout<< "Local out of faces:" << cellLabels_.size()
 226             << " keeping:" << nFound << endl;
 227     }
 228
 229     // Now subset the surface. Do not use triSurface::subsetMesh since requires
 230     // original surface to be in compact numbering.
 231
 232     const triSurface& s = surface_;
 233
 234     // Compact to original triangle
 235     labelList faceMap(s.size());
 236     // Compact to original points
 237     labelList pointMap(s.points().size());
 238     // From original point to compact points
 239     labelList reversePointMap(s.points().size(), -1);
 240
 241     {
 242         label newPointI = 0;
 243         label newTriI = 0;
 244
 245         forAll(s, triI)
 246         {
 247             if (include[triI])
 248             {
 249                 faceMap[newTriI++] = triI;
 250
 251                 const labelledTri& f = s[triI];
 252                 forAll(f, fp)
 253                 {
 254                     if (reversePointMap[f[fp]] == -1)
 255                     {
 256                         pointMap[newPointI] = f[fp];
 257                         reversePointMap[f[fp]] = newPointI++;
 258                     }
 259                 }
 260             }
 261         }
 262         faceMap.setSize(newTriI);
 263         pointMap.setSize(newPointI);
 264     }
 265
 266     // Subset cellLabels
 267     cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
 268
 269     // Store any face per point
 270     pointToFace_.setSize(pointMap.size());
 271
 272     // Create faces and points for subsetted surface
 273     faceList& faces = this->storedFaces();
 274     faces.setSize(faceMap.size());
 275     forAll(faceMap, i)
 276     {
 277         const triFace& f = s[faceMap[i]];
 278         triFace newF
 279         (
 280             reversePointMap[f[0]],
 281             reversePointMap[f[1]],
 282             reversePointMap[f[2]]
 283         );
 284         faces[i] = newF.triFaceFace();
 285
 286         forAll(newF, fp)
 287         {
 288             pointToFace_[newF[fp]] = i;
 289         }
 290     }
 291
 292     this->storedPoints() = pointField(s.points(), pointMap);
 293
 294     if (debug)
 295     {
 296         print(Pout);
 297         Pout<< endl;
 298     }
 299
 300     needsUpdate_ = false;
 301     return true;
 302 }
 303
 304
 305 Foam::tmp<Foam::scalarField>
 306 Foam::sampledTriSurfaceMesh::sample
 307 (
 308     const volScalarField& vField
 309 ) const
 310 {
 311     return sampleField(vField);
 312 }
 313
 314
 315 Foam::tmp<Foam::vectorField>
 316 Foam::sampledTriSurfaceMesh::sample
 317 (
 318     const volVectorField& vField
 319 ) const
 320 {
 321     return sampleField(vField);
 322 }
 323
 324 Foam::tmp<Foam::sphericalTensorField>
 325 Foam::sampledTriSurfaceMesh::sample
 326 (
 327     const volSphericalTensorField& vField
 328 ) const
 329 {
 330     return sampleField(vField);
 331 }
 332
 333
 334 Foam::tmp<Foam::symmTensorField>
 335 Foam::sampledTriSurfaceMesh::sample
 336 (
 337     const volSymmTensorField& vField
 338 ) const
 339 {
 340     return sampleField(vField);
 341 }
 342
 343
 344 Foam::tmp<Foam::tensorField>
 345 Foam::sampledTriSurfaceMesh::sample
 346 (
 347     const volTensorField& vField
 348 ) const
 349 {
 350     return sampleField(vField);
 351 }
 352
 353
 354 Foam::tmp<Foam::scalarField>
 355 Foam::sampledTriSurfaceMesh::interpolate
 356 (
 357     const interpolation<scalar>& interpolator
 358 ) const
 359 {
 360     return interpolateField(interpolator);
 361 }
 362
 363
 364 Foam::tmp<Foam::vectorField>
 365 Foam::sampledTriSurfaceMesh::interpolate
 366 (
 367     const interpolation<vector>& interpolator
 368 ) const
 369 {
 370     return interpolateField(interpolator);
 371 }
 372
 373 Foam::tmp<Foam::sphericalTensorField>
 374 Foam::sampledTriSurfaceMesh::interpolate
 375 (
 376     const interpolation<sphericalTensor>& interpolator
 377 ) const
 378 {
 379     return interpolateField(interpolator);
 380 }
 381
 382
 383 Foam::tmp<Foam::symmTensorField>
 384 Foam::sampledTriSurfaceMesh::interpolate
 385 (
 386     const interpolation<symmTensor>& interpolator
 387 ) const
 388 {
 389     return interpolateField(interpolator);
 390 }
 391
 392
 393 Foam::tmp<Foam::tensorField>
 394 Foam::sampledTriSurfaceMesh::interpolate
 395 (
 396     const interpolation<tensor>& interpolator
 397 ) const
 398 {
 399     return interpolateField(interpolator);
 400 }
 401
 402
 403 void Foam::sampledTriSurfaceMesh::print(Ostream& os) const
 404 {
 405     os  << "sampledTriSurfaceMesh: " << name() << " :"
 406         << "  surface:" << surface_.objectRegistry::name()
 407         << "  faces:" << faces().size()
 408         << "  points:" << points().size();
 409 }
 410
 411
 412 // ************************************************************************* //