src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.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 "isoSurfaceCell.H"
  27 #include "polyMesh.H"
  28 #include "tetMatcher.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 template<class Type>
  33 Type Foam::isoSurfaceCell::generatePoint
  34 (
  35     const DynamicList<Type>& snappedPoints,
  36
  37     const scalar s0,
  38     const Type& p0,
  39     const label p0Index,
  40
  41     const scalar s1,
  42     const Type& p1,
  43     const label p1Index
  44 ) const
  45 {
  46     scalar d = s1-s0;
  47
  48     if (mag(d) > VSMALL)
  49     {
  50         scalar s = (iso_-s0)/d;
  51
  52         if (s >= 0.5 && s <= 1 && p1Index != -1)
  53         {
  54             return snappedPoints[p1Index];
  55         }
  56         else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
  57         {
  58             return snappedPoints[p0Index];
  59         }
  60         else
  61         {
  62             return s*p1 + (1.0-s)*p0;
  63         }
  64     }
  65     else
  66     {
  67         scalar s = 0.4999;
  68
  69         return s*p1 + (1.0-s)*p0;
  70     }
  71 }
  72
  73
  74 template<class Type>
  75 void Foam::isoSurfaceCell::generateTriPoints
  76 (
  77     const DynamicList<Type>& snapped,
  78
  79     const scalar s0,
  80     const Type& p0,
  81     const label p0Index,
  82
  83     const scalar s1,
  84     const Type& p1,
  85     const label p1Index,
  86
  87     const scalar s2,
  88     const Type& p2,
  89     const label p2Index,
  90
  91     const scalar s3,
  92     const Type& p3,
  93     const label p3Index,
  94
  95     DynamicList<Type>& points
  96 ) const
  97 {
  98     int triIndex = 0;
  99     if (s0 < iso_)
 100     {
 101         triIndex |= 1;
 102     }
 103     if (s1 < iso_)
 104     {
 105         triIndex |= 2;
 106     }
 107     if (s2 < iso_)
 108     {
 109         triIndex |= 4;
 110     }
 111     if (s3 < iso_)
 112     {
 113         triIndex |= 8;
 114     }
 115
 116     /* Form the vertices of the triangles for each case */
 117     switch (triIndex)
 118     {
 119         case 0x00:
 120         case 0x0F:
 121         break;
 122
 123         case 0x0E:
 124         case 0x01:
 125             points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
 126             points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
 127             points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
 128         break;
 129
 130         case 0x0D:
 131         case 0x02:
 132             points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
 133             points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
 134             points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
 135         break;
 136
 137         case 0x0C:
 138         case 0x03:
 139         {
 140             Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
 141             Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 142
 143             points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
 144             points.append(tp1);
 145             points.append(tp2);
 146             points.append(tp2);
 147             points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
 148             points.append(tp1);
 149         }
 150         break;
 151
 152         case 0x0B:
 153         case 0x04:
 154         {
 155             points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
 156             points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
 157             points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
 158         }
 159         break;
 160
 161         case 0x0A:
 162         case 0x05:
 163         {
 164             Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
 165             Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 166
 167             points.append(tp0);
 168             points.append(tp1);
 169             points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
 170             points.append(tp0);
 171             points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
 172             points.append(tp1);
 173         }
 174         break;
 175
 176         case 0x09:
 177         case 0x06:
 178         {
 179             Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
 180             Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 181
 182             points.append(tp0);
 183             points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
 184             points.append(tp1);
 185             points.append(tp0);
 186             points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
 187             points.append(tp1);
 188         }
 189         break;
 190
 191         case 0x07:
 192         case 0x08:
 193             points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
 194             points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
 195             points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
 196         break;
 197     }
 198 }
 199
 200
 201 template<class Type>
 202 void Foam::isoSurfaceCell::generateTriPoints
 203 (
 204     const scalarField& cVals,
 205     const scalarField& pVals,
 206
 207     const Field<Type>& cCoords,
 208     const Field<Type>& pCoords,
 209
 210     const DynamicList<Type>& snappedPoints,
 211     const labelList& snappedCc,
 212     const labelList& snappedPoint,
 213
 214     DynamicList<Type>& triPoints,
 215     DynamicList<label>& triMeshCells
 216 ) const
 217 {
 218     tetMatcher tet;
 219
 220     forAll(mesh_.cells(), cellI)
 221     {
 222         if (cellCutType_[cellI] != NOTCUT)
 223         {
 224             label oldNPoints = triPoints.size();
 225
 226             const cell& cFaces = mesh_.cells()[cellI];
 227
 228             if (tet.isA(mesh_, cellI))
 229             {
 230                 // For tets don't do cell-centre decomposition, just use the
 231                 // tet points and values
 232
 233                 const face& f0 = mesh_.faces()[cFaces[0]];
 234
 235                 // Get the other point
 236                 const face& f1 = mesh_.faces()[cFaces[1]];
 237                 label oppositeI = -1;
 238                 forAll(f1, fp)
 239                 {
 240                     oppositeI = f1[fp];
 241
 242                     if (findIndex(f0, oppositeI) == -1)
 243                     {
 244                         break;
 245                     }
 246                 }
 247
 248                 generateTriPoints
 249                 (
 250                     snappedPoints,
 251                     pVals[f0[0]],
 252                     pCoords[f0[0]],
 253                     snappedPoint[f0[0]],
 254
 255                     pVals[f0[1]],
 256                     pCoords[f0[1]],
 257                     snappedPoint[f0[1]],
 258
 259                     pVals[f0[2]],
 260                     pCoords[f0[2]],
 261                     snappedPoint[f0[2]],
 262
 263                     pVals[oppositeI],
 264                     pCoords[oppositeI],
 265                     snappedPoint[oppositeI],
 266
 267                     triPoints
 268                 );
 269             }
 270             else
 271             {
 272                 const cell& cFaces = mesh_.cells()[cellI];
 273
 274                 forAll(cFaces, cFaceI)
 275                 {
 276                     label faceI = cFaces[cFaceI];
 277                     const face& f = mesh_.faces()[faceI];
 278
 279                     for (label fp = 1; fp < f.size() - 1; fp++)
 280                     {
 281                         triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
 282                     //List<triFace> tris(triangulate(f));
 283                     //forAll(tris, i)
 284                     //{
 285                     //    const triFace& tri = tris[i];
 286
 287
 288                         generateTriPoints
 289                         (
 290                             snappedPoints,
 291
 292                             pVals[tri[0]],
 293                             pCoords[tri[0]],
 294                             snappedPoint[tri[0]],
 295
 296                             pVals[tri[1]],
 297                             pCoords[tri[1]],
 298                             snappedPoint[tri[1]],
 299
 300                             pVals[tri[2]],
 301                             pCoords[tri[2]],
 302                             snappedPoint[tri[2]],
 303
 304                             cVals[cellI],
 305                             cCoords[cellI],
 306                             snappedCc[cellI],
 307
 308                             triPoints
 309                         );
 310                     }
 311                 }
 312             }
 313
 314
 315             // Every three triPoints is a cell
 316             label nCells = (triPoints.size()-oldNPoints)/3;
 317             for (label i = 0; i < nCells; i++)
 318             {
 319                 triMeshCells.append(cellI);
 320             }
 321         }
 322     }
 323
 324     triPoints.shrink();
 325     triMeshCells.shrink();
 326 }
 327
 328
 329 template <class Type>
 330 Foam::tmp<Foam::Field<Type> >
 331 Foam::isoSurfaceCell::interpolate
 332 (
 333     const scalarField& cVals,
 334     const scalarField& pVals,
 335     const Field<Type>& cCoords,
 336     const Field<Type>& pCoords
 337 ) const
 338 {
 339     DynamicList<Type> triPoints(nCutCells_);
 340     DynamicList<label> triMeshCells(nCutCells_);
 341
 342     // Dummy snap data
 343     DynamicList<Type> snappedPoints;
 344     labelList snappedCc(mesh_.nCells(), -1);
 345     labelList snappedPoint(mesh_.nPoints(), -1);
 346
 347
 348     generateTriPoints
 349     (
 350         cVals,
 351         pVals,
 352
 353         cCoords,
 354         pCoords,
 355
 356         snappedPoints,
 357         snappedCc,
 358         snappedPoint,
 359
 360         triPoints,
 361         triMeshCells
 362     );
 363
 364
 365     // One value per point
 366     tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
 367     Field<Type>& values = tvalues();
 368
 369     forAll(triPoints, i)
 370     {
 371         label mergedPointI = triPointMergeMap_[i];
 372
 373         if (mergedPointI >= 0)
 374         {
 375             values[mergedPointI] = triPoints[i];
 376         }
 377     }
 378
 379     return tvalues;
 380 }
 381
 382
 383 template <class Type>
 384 Foam::tmp<Foam::Field<Type> >
 385 Foam::isoSurfaceCell::interpolate
 386 (
 387     const Field<Type>& cCoords,
 388     const Field<Type>& pCoords
 389 ) const
 390 {
 391     DynamicList<Type> triPoints(nCutCells_);
 392     DynamicList<label> triMeshCells(nCutCells_);
 393
 394     // Dummy snap data
 395     DynamicList<Type> snappedPoints;
 396     labelList snappedCc(mesh_.nCells(), -1);
 397     labelList snappedPoint(mesh_.nPoints(), -1);
 398
 399
 400     generateTriPoints
 401     (
 402         cVals_,
 403         pVals_,
 404
 405         cCoords,
 406         pCoords,
 407
 408         snappedPoints,
 409         snappedCc,
 410         snappedPoint,
 411
 412         triPoints,
 413         triMeshCells
 414     );
 415
 416
 417     // One value per point
 418     tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
 419     Field<Type>& values = tvalues();
 420
 421     forAll(triPoints, i)
 422     {
 423         label mergedPointI = triPointMergeMap_[i];
 424
 425         if (mergedPointI >= 0)
 426         {
 427             values[mergedPointI] = triPoints[i];
 428         }
 429     }
 430
 431     return tvalues;
 432 }
 433
 434
 435 // ************************************************************************* //