src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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 "error.H"
  27 #include "blockDescriptor.H"
  28
  29 #include "lineEdge.H"
  30 #include "lineDivide.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     //! \cond fileScope
  37     //  Calculate the geometric expension factor from the expansion ratio
  38     inline scalar calcGexp(const scalar expRatio, const label dim)
  39     {
  40         return dim > 1 ? pow(expRatio, 1.0/(dim - 1)) : 0.0;
  41     }
  42     //! \endcond
  43
  44 } // End namespace Foam
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 void Foam::blockDescriptor::makeBlockEdges()
  50 {
  51     const label ni = meshDensity_.x();
  52     const label nj = meshDensity_.y();
  53     const label nk = meshDensity_.z();
  54
  55     // these edges correspond to the "hex" cellModel
  56
  57     // x-direction
  58     setEdge(0,  0, 1, ni);
  59     setEdge(1,  3, 2, ni);
  60     setEdge(2,  7, 6, ni);
  61     setEdge(3,  4, 5, ni);
  62
  63     // y-direction
  64     setEdge(4,  0, 3, nj);
  65     setEdge(5,  1, 2, nj);
  66     setEdge(6,  5, 6, nj);
  67     setEdge(7,  4, 7, nj);
  68
  69     // z-direction
  70     setEdge(8,  0, 4, nk);
  71     setEdge(9,  1, 5, nk);
  72     setEdge(10, 2, 6, nk);
  73     setEdge(11, 3, 7, nk);
  74 }
  75
  76
  77 void Foam::blockDescriptor::setEdge
  78 (
  79     label edgeI,
  80     label start,
  81     label end,
  82     label dim
  83 )
  84 {
  85     // set reference to the list of labels defining the block
  86     const labelList& blockLabels = blockShape_;
  87
  88     // set reference to global list of points
  89     const pointField blockPoints = blockShape_.points(blockPointField_);
  90
  91     // Set the edge points/weights
  92     // The edge is a straight-line if it is not in the list of curvedEdges
  93
  94     // calc geometric expension factor from the expansion ratio
  95     const scalar gExp = calcGexp(expand_[edgeI], dim);
  96
  97     forAll(curvedEdges_, cedgeI)
  98     {
  99         const curvedEdge& cedge = curvedEdges_[cedgeI];
 100
 101         int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
 102
 103         if (cmp)
 104         {
 105             if (cmp > 0)
 106             {
 107                 // curve has the same orientation
 108
 109                 // divide the line
 110                 lineDivide divEdge(cedge, dim, gExp);
 111
 112                 edgePoints_[edgeI]  = divEdge.points();
 113                 edgeWeights_[edgeI] = divEdge.lambdaDivisions();
 114             }
 115             else
 116             {
 117                 // curve has the opposite orientation
 118
 119                 // divide the line
 120                 lineDivide divEdge(cedge, dim, 1.0/(gExp+SMALL));
 121
 122                 const pointField& p = divEdge.points();
 123                 const scalarList& d = divEdge.lambdaDivisions();
 124
 125                 edgePoints_[edgeI].setSize(p.size());
 126                 edgeWeights_[edgeI].setSize(d.size());
 127
 128                 label pMax = p.size() - 1;
 129                 forAll(p, pI)
 130                 {
 131                     edgePoints_[edgeI][pI]  = p[pMax - pI];
 132                     edgeWeights_[edgeI][pI] = 1.0 - d[pMax - pI];
 133                 }
 134
 135             }
 136
 137             // found curved-edge: done
 138             return;
 139         }
 140     }
 141
 142
 143     // not found: divide the edge as a straight line
 144
 145     lineDivide divEdge
 146     (
 147         lineEdge(blockPoints, start, end),
 148         dim,
 149         gExp
 150     );
 151
 152     edgePoints_[edgeI]  = divEdge.points();
 153     edgeWeights_[edgeI] = divEdge.lambdaDivisions();
 154 }
 155
 156
 157 // ************************************************************************* //