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