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