src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "cellModel.H"
  27 #include "pyramid.H"
  28
  29 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  30
  31 Foam::vector Foam::cellModel::centre
  32 (
  33     const labelList& pointLabels,
  34     const pointField& points
  35 ) const
  36 {
  37     // Estimate centre of cell
  38     vector cEst = vector::zero;
  39
  40     // Sum the points idicated by the label list
  41     forAll(pointLabels, i)
  42     {
  43         cEst += points[pointLabels[i]];
  44     }
  45
  46     // Average by dividing by the number summed over.
  47     cEst /= scalar(pointLabels.size());
  48
  49
  50     // Calculate the centre by breaking the cell into pyramids and
  51     // volume-weighted averaging their centres
  52     scalar sumV = 0.0;
  53     vector sumVc = vector::zero;
  54
  55     const faceList cellFaces = faces(pointLabels);
  56
  57     forAll(cellFaces, i)
  58     {
  59         const face& curFace = cellFaces[i];
  60
  61         scalar pyrVol =
  62             pyramid<point, const point&, const face&>
  63             (
  64                 curFace,
  65                 cEst
  66             ).mag(points);
  67
  68         if (pyrVol > SMALL)
  69         {
  70             WarningIn("cellModel::centre(const labelList&, const pointField&)")
  71                 << "zero or negative pyramid volume: " << -pyrVol
  72                 << " for face " << i
  73                 << endl;
  74         }
  75
  76         sumVc -=
  77             pyrVol
  78            *pyramid<point, const point&, const face&>(curFace, cEst)
  79            .centre(points);
  80
  81         sumV -= pyrVol;
  82     }
  83
  84     return sumVc/(sumV + VSMALL);
  85 }
  86
  87
  88 Foam::scalar Foam::cellModel::mag
  89 (
  90     const labelList& pointLabels,
  91     const pointField& points
  92 ) const
  93 {
  94     // Estimate centre of cell
  95     vector cEst = vector::zero;
  96
  97     // Sum the points idicated by the label list
  98     forAll(pointLabels, i)
  99     {
 100         cEst += points[pointLabels[i]];
 101     }
 102
 103     // Average by dividing by the number summed over.
 104     cEst /= scalar(pointLabels.size());
 105
 106
 107     // Calculate the magnitude by summing the -mags of the pyramids
 108     // The sign change is because the faces point outwards
 109     // and a pyramid is constructed from an inward pointing face
 110     // and the base centre-apex vector
 111     scalar v = 0;
 112
 113     const faceList cellFaces = faces(pointLabels);
 114
 115     forAll(cellFaces, i)
 116     {
 117         const face& curFace =cellFaces[i];
 118
 119         scalar pyrVol =
 120             pyramid<point, const point&, const face&>
 121             (
 122                 curFace,
 123                 cEst
 124             ).mag(points);
 125
 126         if (pyrVol > SMALL)
 127         {
 128             WarningIn("cellModel::mag(const labelList&, const pointField&)")
 129                 << "zero or negative pyramid volume: " << -pyrVol
 130                 << " for face " << i
 131                 << endl;
 132         }
 133
 134         v -= pyrVol;
 135     }
 136
 137     return v;
 138 }
 139
 140 // ************************************************************************* //