src/sampling/sampledSet/patchCloud/patchCloudSet.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011-2011 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 "patchCloudSet.H"
  27 #include "polyMesh.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "pointIndexHit.H"
  30 #include "Tuple2.H"
  31 #include "treeBoundBox.H"
  32 #include "indexedOctree.H"
  33 #include "treeDataFace.H"
  34 #include "Time.H"
  35 #include "meshTools.H"
  36 #include "wordReList.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(patchCloudSet, 0);
  43     addToRunTimeSelectionTable(sampledSet, patchCloudSet, word);
  44
  45
  46     //- Helper class for finding nearest
  47     // Nearest:
  48     //  - point+local index
  49     //  - sqr(distance)
  50     //  - processor
  51     typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
  52
  53     class nearestEqOp
  54     {
  55
  56     public:
  57
  58         void operator()(nearInfo& x, const nearInfo& y) const
  59         {
  60             if (y.first().hit())
  61             {
  62                 if (!x.first().hit())
  63                 {
  64                     x = y;
  65                 }
  66                 else if (y.second().first() < x.second().first())
  67                 {
  68                     x = y;
  69                 }
  70             }
  71         }
  72     };
  73
  74 }
  75
  76
  77 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  78
  79 void Foam::patchCloudSet::calcSamples
  80 (
  81     DynamicList<point>& samplingPts,
  82     DynamicList<label>& samplingCells,
  83     DynamicList<label>& samplingFaces,
  84     DynamicList<label>& samplingSegments,
  85     DynamicList<scalar>& samplingCurveDist
  86 ) const
  87 {
  88     if (debug)
  89     {
  90         Info<< "patchCloudSet : sampling on patches :" << endl;
  91     }
  92
  93     // Construct search tree for all patch faces.
  94     label sz = 0;
  95     forAllConstIter(labelHashSet, patchSet_, iter)
  96     {
  97         const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
  98
  99         sz += pp.size();
 100
 101         if (debug)
 102         {
 103             Info<< "    " << pp.name() << " size " << pp.size() << endl;
 104         }
 105     }
 106
 107     labelList patchFaces(sz);
 108     sz = 0;
 109     treeBoundBox bb(point::max, point::min);
 110     forAllConstIter(labelHashSet, patchSet_, iter)
 111     {
 112         const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
 113
 114         forAll(pp, i)
 115         {
 116             patchFaces[sz++] = pp.start()+i;
 117         }
 118
 119         const boundBox patchBb(pp.localPoints());
 120
 121         bb.min() = min(bb.min(), patchBb.min());
 122         bb.max() = max(bb.max(), patchBb.max());
 123     }
 124
 125     // Not very random
 126     Random rndGen(123456);
 127     // Make bb asymetric just to avoid problems on symmetric meshes
 128     bb = bb.extend(rndGen, 1E-4);
 129
 130     bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 131     bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 132
 133
 134     indexedOctree<treeDataFace> patchTree
 135     (
 136         treeDataFace    // all information needed to search faces
 137         (
 138             false,      // do not cache bb
 139             mesh(),
 140             patchFaces  // boundary faces only
 141         ),
 142         bb,             // overall search domain
 143         8,              // maxLevel
 144         10,             // leafsize
 145         3.0             // duplicity
 146     );
 147
 148
 149
 150     // All the info for nearest. Construct to miss
 151     List<nearInfo> nearest(sampleCoords_.size());
 152
 153     forAll(sampleCoords_, sampleI)
 154     {
 155         const point& sample = sampleCoords_[sampleI];
 156
 157         pointIndexHit& nearInfo = nearest[sampleI].first();
 158
 159         // Find the nearest locally
 160         nearInfo = patchTree.findNearest(sample, magSqr(bb.span()));
 161
 162         // Fill in the distance field and the processor field
 163         if (!nearInfo.hit())
 164         {
 165             nearest[sampleI].second().first() = Foam::sqr(GREAT);
 166             nearest[sampleI].second().second() = Pstream::myProcNo();
 167         }
 168         else
 169         {
 170             // Set nearest to mesh face label
 171             nearInfo.setIndex(patchFaces[nearInfo.index()]);
 172
 173             nearest[sampleI].second().first() = magSqr
 174             (
 175                 nearInfo.hitPoint()
 176               - sample
 177             );
 178             nearest[sampleI].second().second() = Pstream::myProcNo();
 179         }
 180     }
 181
 182
 183     // Find nearest.
 184     Pstream::listCombineGather(nearest, nearestEqOp());
 185     Pstream::listCombineScatter(nearest);
 186
 187
 188     if (debug && Pstream::master())
 189     {
 190         OFstream str
 191         (
 192             mesh().time().path()
 193           / name()
 194           + "_nearest.obj"
 195         );
 196         Info<< "Dumping mapping as lines from supplied points to"
 197             << " nearest patch face to file " << str.name() << endl;
 198
 199         label vertI = 0;
 200
 201         forAll(nearest, i)
 202         {
 203             meshTools::writeOBJ(str, sampleCoords_[i]);
 204             vertI++;
 205             meshTools::writeOBJ(str, nearest[i].first().hitPoint());
 206             vertI++;
 207             str << "l " << vertI-1 << ' ' << vertI << nl;
 208         }
 209     }
 210
 211
 212     // Store the sampling locations on the nearest processor
 213     forAll(nearest, sampleI)
 214     {
 215         if (nearest[sampleI].second().second() == Pstream::myProcNo())
 216         {
 217             const pointIndexHit& nearInfo = nearest[sampleI].first();
 218
 219             label faceI = nearInfo.index();
 220
 221             samplingPts.append(nearInfo.hitPoint());
 222             samplingCells.append(mesh().faceOwner()[faceI]);
 223             samplingFaces.append(faceI);
 224             samplingSegments.append(0);
 225             samplingCurveDist.append(1.0 * sampleI);
 226         }
 227     }
 228 }
 229
 230
 231 void Foam::patchCloudSet::genSamples()
 232 {
 233     // Storage for sample points
 234     DynamicList<point> samplingPts;
 235     DynamicList<label> samplingCells;
 236     DynamicList<label> samplingFaces;
 237     DynamicList<label> samplingSegments;
 238     DynamicList<scalar> samplingCurveDist;
 239
 240     calcSamples
 241     (
 242         samplingPts,
 243         samplingCells,
 244         samplingFaces,
 245         samplingSegments,
 246         samplingCurveDist
 247     );
 248
 249     samplingPts.shrink();
 250     samplingCells.shrink();
 251     samplingFaces.shrink();
 252     samplingSegments.shrink();
 253     samplingCurveDist.shrink();
 254
 255     setSamples
 256     (
 257         samplingPts,
 258         samplingCells,
 259         samplingFaces,
 260         samplingSegments,
 261         samplingCurveDist
 262     );
 263 }
 264
 265
 266 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 267
 268 Foam::patchCloudSet::patchCloudSet
 269 (
 270     const word& name,
 271     const polyMesh& mesh,
 272     meshSearch& searchEngine,
 273     const word& axis,
 274     const List<point>& sampleCoords,
 275     const labelHashSet& patchSet
 276 )
 277 :
 278     sampledSet(name, mesh, searchEngine, axis),
 279     sampleCoords_(sampleCoords),
 280     patchSet_(patchSet)
 281 {
 282     genSamples();
 283
 284     if (debug)
 285     {
 286         write(Info);
 287     }
 288 }
 289
 290
 291 Foam::patchCloudSet::patchCloudSet
 292 (
 293     const word& name,
 294     const polyMesh& mesh,
 295     meshSearch& searchEngine,
 296     const dictionary& dict
 297 )
 298 :
 299     sampledSet(name, mesh, searchEngine, dict),
 300     sampleCoords_(dict.lookup("points")),
 301     patchSet_
 302     (
 303         mesh.boundaryMesh().patchSet
 304         (
 305             wordList(dict.lookup("patches"))
 306         )
 307     )
 308 {
 309     genSamples();
 310
 311     if (debug)
 312     {
 313         write(Info);
 314     }
 315 }
 316
 317
 318 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 319
 320 Foam::patchCloudSet::~patchCloudSet()
 321 {}
 322
 323
 324 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 325
 326 Foam::point Foam::patchCloudSet::getRefPoint(const List<point>& pts) const
 327 {
 328     if (pts.size())
 329     {
 330         // Use first samplePt as starting point
 331         return pts[0];
 332     }
 333     else
 334     {
 335         return vector::zero;
 336     }
 337 }
 338
 339
 340 // ************************************************************************* //