src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "faceDecompIsoSurfaceCuts.H"
  30 #include "volPointInterpolation.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 void Foam::faceDecompIsoSurfaceCuts::constructEdgeCuts
  35 (
  36     const volScalarField& volField,
  37     const pointScalarField& ptField,
  38     const scalar isoVal,
  39     const scalar tol
  40 )
  41 {
  42     const primitiveMesh& mesh = volField.mesh();
  43
  44     // Intermediate storage for labels of cut edges and cut points
  45     label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
  46
  47     labelHashSet cutCells(nBFaces);
  48     labelHashSet cutVerts(nBFaces);
  49     labelHashSet cutFaceCentres(nBFaces);
  50     labelHashSet cutCellCentres(nBFaces);
  51
  52     DynamicList<label> meshCutEdges(nBFaces);
  53     DynamicList<scalar> meshEdgeWeights(nBFaces);
  54
  55     DynamicList<pyramidEdge> cutPyrEdges(nBFaces);
  56     DynamicList<scalar> pyrEdgeWeights(nBFaces);
  57
  58     DynamicList<faceEdge> cutFaceEdges(nBFaces);
  59     DynamicList<scalar> faceEdgeWeights(nBFaces);
  60
  61     DynamicList<centreEdge> cutCentreEdges(nBFaces);
  62     DynamicList<scalar> centreEdgeWeights(nBFaces);
  63
  64
  65     // Walk over faces
  66     forAll(mesh.faces(), faceI)
  67     {
  68         const face& f = mesh.faces()[faceI];
  69
  70         // Face centre value
  71         scalar fcVal = 0;
  72
  73         forAll(f, fp)
  74         {
  75             fcVal += ptField[f[fp]];
  76         }
  77         fcVal /= f.size();
  78
  79
  80         label ownerI = mesh.faceOwner()[faceI];
  81         label neighbourI = -1;
  82         if (mesh.isInternalFace(faceI))
  83         {
  84             neighbourI = mesh.faceNeighbour()[faceI];
  85         }
  86
  87         scalar weight;
  88
  89         //
  90         // Check faceCentre-cellCentre crossing for owner and neigbour
  91         //
  92
  93         if (crosses(isoVal, fcVal, volField[ownerI], weight))
  94         {
  95             mark(ownerI, cutCells);
  96
  97             if (weight < tol)
  98             {
  99                 mark(faceI, cutFaceCentres);
 100             }
 101             else if (weight > 1-tol)
 102             {
 103                 mark(ownerI, cutCellCentres);
 104             }
 105             else
 106             {
 107                 cutCentreEdges.append(centreEdge(faceI, true));
 108                 centreEdgeWeights.append(weight);
 109             }
 110         }
 111
 112         if (neighbourI != -1)
 113         {
 114             if (crosses(isoVal, fcVal, volField[neighbourI], weight))
 115             {
 116                 mark(neighbourI, cutCells);
 117
 118                 if (weight < tol)
 119                 {
 120                     mark(faceI, cutFaceCentres);
 121                 }
 122                 else if (weight > 1-tol)
 123                 {
 124                     mark(neighbourI, cutCellCentres);
 125                 }
 126                 else
 127                 {
 128                     cutCentreEdges.append(centreEdge(faceI, false));
 129                     centreEdgeWeights.append(weight);
 130                 }
 131             }
 132         }
 133
 134
 135         forAll(f, fp)
 136         {
 137             //
 138             // Check face internal (faceEdge) edge crossing
 139             //
 140
 141             if (crosses(isoVal, ptField[f[fp]], fcVal, weight))
 142             {
 143                 mark(ownerI, cutCells);
 144                 if (neighbourI != -1)
 145                 {
 146                     mark(neighbourI, cutCells);
 147                 }
 148
 149                 if (weight < tol)
 150                 {
 151                     mark(f[fp], cutVerts);
 152                 }
 153                 else if (weight > 1-tol)
 154                 {
 155                     mark(faceI, cutFaceCentres);
 156                 }
 157                 else
 158                 {
 159                     cutFaceEdges.append(faceEdge(faceI, fp));
 160                     faceEdgeWeights.append(weight);
 161                 }
 162             }
 163         }
 164     }
 165
 166
 167     // Walk over all mesh points
 168     forAll(mesh.points(), pointI)
 169     {
 170         const labelList& myCells = mesh.pointCells()[pointI];
 171
 172         forAll(myCells, myCellI)
 173         {
 174             label cellI = myCells[myCellI];
 175
 176             //
 177             // Check pyramid edge crossing
 178             //
 179
 180             scalar weight;
 181
 182             if (crosses(isoVal, ptField[pointI], volField[cellI], weight))
 183             {
 184                 mark(cellI, cutCells);
 185
 186                 if (weight < tol)
 187                 {
 188                     mark(pointI, cutVerts);
 189                 }
 190                 else if (weight > 1-tol)
 191                 {
 192                     mark(cellI, cutCellCentres);
 193                 }
 194                 else
 195                 {
 196                     cutPyrEdges.append(pyramidEdge(pointI, cellI));
 197                     pyrEdgeWeights.append(weight);
 198                 }
 199             }
 200         }
 201     }
 202
 203
 204     // Walk over all mesh edges
 205     forAll(mesh.edges(), edgeI)
 206     {
 207         const edge& e = mesh.edges()[edgeI];
 208
 209         scalar weight;
 210
 211         if (crosses(isoVal, ptField[e.start()], ptField[e.end()], weight))
 212         {
 213             mark(mesh.edgeCells()[edgeI], cutCells);
 214
 215             if (weight < tol)
 216             {
 217                 mark(e.start(), cutVerts);
 218             }
 219             else if (weight > 1-tol)
 220             {
 221                 mark(e.end(), cutVerts);
 222             }
 223             else
 224             {
 225                 meshCutEdges.append(edgeI);
 226                 meshEdgeWeights.append(weight);
 227             }
 228         }
 229     }
 230
 231     meshCutEdges.shrink();
 232     meshEdgeWeights.shrink();
 233
 234     cutPyrEdges.shrink();
 235     pyrEdgeWeights.shrink();
 236
 237     cutFaceEdges.shrink();
 238     faceEdgeWeights.shrink();
 239
 240     cutCentreEdges.shrink();
 241     centreEdgeWeights.shrink();
 242
 243
 244
 245     // Tranfer lists to faceDecompCuts
 246
 247     cells_ = cutCells.toc();
 248     meshVerts_ = cutVerts.toc();
 249     meshFaceCentres_ = cutFaceCentres.toc();
 250     meshCellCentres_ = cutCellCentres.toc();
 251
 252     meshEdges_.transfer(meshCutEdges);
 253     meshEdgeWeights_.transfer(meshEdgeWeights);
 254
 255     pyrEdges_.transfer(cutPyrEdges);
 256     pyrEdgeWeights_.transfer(pyrEdgeWeights);
 257
 258     faceEdges_.transfer(cutFaceEdges);
 259     faceEdgeWeights_.transfer(faceEdgeWeights);
 260
 261     centreEdges_.transfer(cutCentreEdges);
 262     centreEdgeWeights_.transfer(centreEdgeWeights);
 263 }
 264
 265
 266 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 267
 268 Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts
 269 (
 270     const volScalarField& volField,
 271     const pointScalarField& ptField,
 272     const scalar isoVal,
 273     const scalar tol
 274 )
 275 :
 276     faceDecompCuts
 277     (
 278         volField.mesh(),
 279         labelList(),
 280
 281         labelList(),
 282         labelList(),
 283         labelList(),
 284
 285         labelList(),            // mesh edges
 286         scalarField(),
 287
 288         List<pyramidEdge>(),    // pyramid edges
 289         scalarField(),
 290
 291         List<faceEdge>(),       // face edges
 292         scalarField(),
 293
 294         List<centreEdge>(),     // cell centre edges
 295         scalarField()
 296     )
 297 {
 298     constructEdgeCuts(volField, ptField, isoVal, tol);
 299 }
 300
 301
 302 // Construct from interpolator and field (does interpolation)
 303 Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts
 304 (
 305     const volScalarField& volField,
 306     const volPointInterpolation& pInterp,
 307     const scalar isoVal,
 308     const scalar tol
 309 )
 310 :
 311     faceDecompCuts
 312     (
 313         volField.mesh(),
 314         labelList(),
 315
 316         labelList(),
 317         labelList(),
 318         labelList(),
 319
 320         labelList(),            // mesh edges
 321         scalarField(),
 322
 323         List<pyramidEdge>(),    // pyramid edges
 324         scalarField(),
 325
 326         List<faceEdge>(),       // face edges
 327         scalarField(),
 328
 329         List<centreEdge>(),     // cell centre edges
 330         scalarField()
 331     )
 332 {
 333     // Get field on vertices
 334     tmp<pointScalarField> tptField = pInterp.interpolate(volField);
 335     const pointScalarField& ptField = tptField();
 336
 337
 338     constructEdgeCuts(volField, ptField, isoVal, tol);
 339 }
 340
 341
 342 // ************************************************************************* //