src/sampling/probes/patchProbes.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 "patchProbes.H"
  27 #include "volFields.H"
  28 #include "IOmanip.H"
  29 // For 'nearInfo' helper class only
  30 #include "directMappedPatchBase.H"
  31 //#include "meshSearch.H"
  32 #include "treeBoundBox.H"
  33 #include "treeDataFace.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(patchProbes, 0);
  40 }
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 void Foam::patchProbes::findElements(const fvMesh& mesh)
  45 {
  46
  47     const polyBoundaryMesh& bm = mesh.boundaryMesh();
  48
  49     label patchI = bm.findPatchID(patchName_);
  50
  51     if (patchI == -1)
  52     {
  53         FatalErrorIn
  54         (
  55             " Foam::patchProbes::findElements(const fvMesh&)"
  56         )   << " Unknown patch name "
  57             << patchName_ << endl
  58             << exit(FatalError);
  59     }
  60
  61      // All the info for nearest. Construct to miss
  62     List<directMappedPatchBase::nearInfo> nearest(this->size());
  63
  64     const polyPatch& pp = bm[patchI];
  65
  66     if (pp.size() > 0)
  67     {
  68         labelList bndFaces(pp.size());
  69         forAll(bndFaces, i)
  70         {
  71             bndFaces[i] =  pp.start() + i;
  72         }
  73
  74         treeBoundBox overallBb(pp.points());
  75         Random rndGen(123456);
  76         overallBb = overallBb.extend(rndGen, 1E-4);
  77         overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
  78         overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
  79
  80         const indexedOctree<treeDataFace> boundaryTree
  81         (
  82             treeDataFace    // all information needed to search faces
  83             (
  84                 false,                      // do not cache bb
  85                 mesh,
  86                 bndFaces                    // patch faces only
  87             ),
  88             overallBb,                      // overall search domain
  89             8,                              // maxLevel
  90             10,                             // leafsize
  91             3.0                             // duplicity
  92         );
  93
  94
  95         forAll(probeLocations(), probeI)
  96         {
  97             const point sample = probeLocations()[probeI];
  98
  99             scalar span = boundaryTree.bb().mag();
 100
 101             pointIndexHit info = boundaryTree.findNearest
 102             (
 103                 sample,
 104                 Foam::sqr(span)
 105             );
 106
 107             if (!info.hit())
 108             {
 109                 info = boundaryTree.findNearest
 110                 (
 111                     sample,
 112                     Foam::sqr(GREAT)
 113                 );
 114             }
 115
 116             label faceI = boundaryTree.shapes().faceLabels()[info.index()];
 117
 118             const label patchi = bm.whichPatch(faceI);
 119
 120             if (isA<emptyPolyPatch>(bm[patchi]))
 121             {
 122                 WarningIn
 123                 (
 124                     " Foam::patchProbes::findElements(const fvMesh&)"
 125                 )
 126                 << " The sample point: " << sample
 127                 << " belongs to " << patchi
 128                 << " which is an empty patch. This is not permitted. "
 129                 << " This sample will not be included "
 130                 << endl;
 131             }
 132             else
 133             {
 134                 const point& fc = mesh.faceCentres()[faceI];
 135
 136                 directMappedPatchBase::nearInfo sampleInfo;
 137
 138                 sampleInfo.first() = pointIndexHit
 139                 (
 140                     true,
 141                     fc,
 142                     faceI
 143                 );
 144
 145                 sampleInfo.second().first() = magSqr(fc-sample);
 146                 sampleInfo.second().second() = Pstream::myProcNo();
 147
 148                 nearest[probeI]= sampleInfo;
 149             }
 150         }
 151     }
 152
 153
 154     // Find nearest.
 155     Pstream::listCombineGather(nearest, directMappedPatchBase::nearestEqOp());
 156     Pstream::listCombineScatter(nearest);
 157
 158     if (debug)
 159     {
 160         Info<< "patchProbes::findElements" << " : " << endl;
 161         forAll(nearest, sampleI)
 162         {
 163             label procI = nearest[sampleI].second().second();
 164             label localI = nearest[sampleI].first().index();
 165
 166             Info<< "    " << sampleI << " coord:"<< operator[](sampleI)
 167                 << " found on processor:" << procI
 168                 << " in local cell/face:" << localI
 169                 << " with fc:" << nearest[sampleI].first().rawPoint() << endl;
 170         }
 171     }
 172
 173
 174     // Extract any local faces to sample
 175     elementList_.setSize(nearest.size(), -1);
 176
 177     forAll(nearest, sampleI)
 178     {
 179         if (nearest[sampleI].second().second() == Pstream::myProcNo())
 180         {
 181             // Store the face to sample
 182             elementList_[sampleI] = nearest[sampleI].first().index();
 183         }
 184     }
 185 }
 186
 187
 188 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 189
 190 Foam::patchProbes::patchProbes
 191 (
 192     const word& name,
 193     const objectRegistry& obr,
 194     const dictionary& dict,
 195     const bool loadFromFiles
 196 )
 197 :
 198     probes(name, obr, dict, loadFromFiles)
 199 {
 200     // When constructing probes above it will have called the
 201     // probes::findElements (since the virtual mechanism not yet operating).
 202     // Not easy to workaround (apart from feeding through flag into constructor)
 203     // so clear out any cells found for now.
 204     elementList_.clear();
 205
 206     read(dict);
 207 }
 208
 209
 210 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 211
 212 Foam::patchProbes::~patchProbes()
 213 {}
 214
 215
 216 void Foam::patchProbes::write()
 217 {
 218     if (this->size() && prepare())
 219     {
 220         sampleAndWrite(scalarFields_);
 221         sampleAndWrite(vectorFields_);
 222         sampleAndWrite(sphericalTensorFields_);
 223         sampleAndWrite(symmTensorFields_);
 224         sampleAndWrite(tensorFields_);
 225     }
 226 }
 227
 228 void Foam::patchProbes::read(const dictionary& dict)
 229 {
 230     dict.lookup("patchName") >> patchName_;
 231     probes::read(dict);
 232 }
 233
 234
 235 // ************************************************************************* //