src/sampling/sampledSurface/thresholdCellFaces/thresholdCellFaces.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 "thresholdCellFaces.H"
  27
  28 #include "polyMesh.H"
  29 #include "DynamicList.H"
  30
  31 #include "emptyPolyPatch.H"
  32 #include "processorPolyPatch.H"
  33
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::thresholdCellFaces, 0);
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 void Foam::thresholdCellFaces::calculate
  42 (
  43     const scalarField& field,
  44     const scalar lowerThreshold,
  45     const scalar upperThreshold,
  46     const bool triangulate
  47 )
  48 {
  49     const labelList& own = mesh_.faceOwner();
  50     const labelList& nei = mesh_.faceNeighbour();
  51
  52     const faceList& origFaces = mesh_.faces();
  53     const pointField& origPoints = mesh_.points();
  54
  55     const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
  56
  57
  58     surfZoneList surfZones(bMesh.size()+1);
  59
  60     surfZones[0] = surfZone
  61     (
  62         "internalMesh",
  63         0,  // size
  64         0,  // start
  65         0   // index
  66     );
  67
  68     forAll(bMesh, patchI)
  69     {
  70         surfZones[patchI+1] = surfZone
  71         (
  72             bMesh[patchI].name(),
  73             0,        // size
  74             0,        // start
  75             patchI+1  // index
  76         );
  77     }
  78
  79
  80     label nFaces = 0;
  81     label nPoints = 0;
  82
  83
  84     meshCells_.clear();
  85
  86     DynamicList<face>  surfFaces(0.5 * mesh_.nFaces());
  87     DynamicList<label> surfCells(surfFaces.size());
  88
  89     labelList oldToNewPoints(origPoints.size(), -1);
  90
  91
  92     // internal faces only
  93     for (label faceI = 0; faceI < mesh_.nInternalFaces(); ++faceI)
  94     {
  95         int side = 0;
  96
  97         // check lowerThreshold
  98         if (field[own[faceI]] > lowerThreshold)
  99         {
 100             if (field[nei[faceI]] < lowerThreshold)
 101             {
 102                 side = +1;
 103             }
 104         }
 105         else if (field[nei[faceI]] > lowerThreshold)
 106         {
 107             side = -1;
 108         }
 109
 110         // check upperThreshold
 111         if (field[own[faceI]] < upperThreshold)
 112         {
 113             if (field[nei[faceI]] > upperThreshold)
 114             {
 115                 side = +1;
 116             }
 117         }
 118         else if (field[nei[faceI]] < upperThreshold)
 119         {
 120             side = -1;
 121         }
 122
 123
 124         if (side)
 125         {
 126             const face& f = origFaces[faceI];
 127
 128             forAll(f, fp)
 129             {
 130                 if (oldToNewPoints[f[fp]] == -1)
 131                 {
 132                     oldToNewPoints[f[fp]] = nPoints++;
 133                 }
 134             }
 135
 136
 137             label cellId;
 138             face  surfFace;
 139
 140             if (side > 0)
 141             {
 142                 surfFace = f;
 143                 cellId = own[faceI];
 144             }
 145             else
 146             {
 147                 surfFace = f.reverseFace();
 148                 cellId = nei[faceI];
 149             }
 150
 151
 152             if (triangulate)
 153             {
 154                 label count = surfFace.triangles(origPoints, surfFaces);
 155                 while (count-- > 0)
 156                 {
 157                     surfCells.append(cellId);
 158                 }
 159             }
 160             else
 161             {
 162                 surfFaces.append(surfFace);
 163                 surfCells.append(cellId);
 164             }
 165         }
 166     }
 167
 168     surfZones[0].size() = surfFaces.size();
 169
 170
 171     // nothing special for processor patches?
 172     forAll(bMesh, patchI)
 173     {
 174         const polyPatch& p = bMesh[patchI];
 175         surfZone& zone = surfZones[patchI+1];
 176
 177         zone.start() = nFaces;
 178
 179         if
 180         (
 181             isA<emptyPolyPatch>(p)
 182          || (Pstream::parRun() && isA<processorPolyPatch>(p))
 183         )
 184         {
 185             continue;
 186         }
 187
 188         label faceI = p.start();
 189
 190         // patch faces
 191         forAll(p, localFaceI)
 192         {
 193             if
 194             (
 195                 field[own[faceI]] > lowerThreshold
 196              && field[own[faceI]] < upperThreshold
 197             )
 198             {
 199                 const face& f = origFaces[faceI];
 200                 forAll(f, fp)
 201                 {
 202                     if (oldToNewPoints[f[fp]] == -1)
 203                     {
 204                         oldToNewPoints[f[fp]] = nPoints++;
 205                     }
 206                 }
 207
 208                 label cellId = own[faceI];
 209
 210                 if (triangulate)
 211                 {
 212                     label count = f.triangles(origPoints, surfFaces);
 213                     while (count-- > 0)
 214                     {
 215                         surfCells.append(cellId);
 216                     }
 217                 }
 218                 else
 219                 {
 220                     surfFaces.append(f);
 221                     surfCells.append(cellId);
 222                 }
 223             }
 224
 225             ++faceI;
 226         }
 227
 228         zone.size() = surfFaces.size() - zone.start();
 229     }
 230
 231
 232     surfFaces.shrink();
 233     surfCells.shrink();
 234
 235     // renumber
 236     forAll(surfFaces, faceI)
 237     {
 238         inplaceRenumber(oldToNewPoints, surfFaces[faceI]);
 239     }
 240
 241
 242     pointField surfPoints(nPoints);
 243     nPoints = 0;
 244     forAll(oldToNewPoints, pointI)
 245     {
 246         if (oldToNewPoints[pointI] >= 0)
 247         {
 248             surfPoints[oldToNewPoints[pointI]] = origPoints[pointI];
 249             nPoints++;
 250         }
 251     }
 252     surfPoints.setSize(nPoints);
 253
 254     this->storedPoints().transfer(surfPoints);
 255     this->storedFaces().transfer(surfFaces);
 256     this->storedZones().transfer(surfZones);
 257
 258     meshCells_.transfer(surfCells);
 259 }
 260
 261
 262 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 263
 264 Foam::thresholdCellFaces::thresholdCellFaces
 265 (
 266     const polyMesh& mesh,
 267     const scalarField& field,
 268     const scalar lowerThreshold,
 269     const scalar upperThreshold,
 270     const bool triangulate
 271 )
 272 :
 273     mesh_(mesh)
 274 {
 275
 276     if (lowerThreshold > upperThreshold)
 277     {
 278         WarningIn("thresholdCellFaces::thresholdCellFaces(...)")
 279             << "lower > upper limit!  "
 280             << lowerThreshold << " > " << upperThreshold << endl;
 281     }
 282
 283     calculate(field, lowerThreshold, upperThreshold, triangulate);
 284 }
 285
 286
 287 // ************************************************************************* //