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