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