applications/utilities/mesh/generation/blockMesh/blockPoints.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 Description
  26     private member of block. Creates vertices for cells filling the block.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "error.H"
  31 #include "block.H"
  32
  33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  34
  35 void Foam::block::blockPoints()
  36 {
  37     // set local variables for mesh specification
  38     const label ni = blockDef_.n().x();
  39     const label nj = blockDef_.n().y();
  40     const label nk = blockDef_.n().z();
  41
  42     const point p000 = blockDef_.points()[blockDef_.blockShape()[0]];
  43     const point p100 = blockDef_.points()[blockDef_.blockShape()[1]];
  44     const point p110 = blockDef_.points()[blockDef_.blockShape()[2]];
  45     const point p010 = blockDef_.points()[blockDef_.blockShape()[3]];
  46
  47     const point p001 = blockDef_.points()[blockDef_.blockShape()[4]];
  48     const point p101 = blockDef_.points()[blockDef_.blockShape()[5]];
  49     const point p111 = blockDef_.points()[blockDef_.blockShape()[6]];
  50     const point p011 = blockDef_.points()[blockDef_.blockShape()[7]];
  51
  52     // list of edge point and weighting factors
  53     const List<List<point> >& p = blockDef_.blockEdgePoints();
  54     const scalarListList& w = blockDef_.blockEdgeWeights();
  55
  56     // generate vertices
  57
  58     for (label k = 0; k <= nk; k++)
  59     {
  60         for (label j = 0; j <= nj; j++)
  61         {
  62             for (label i = 0; i <= ni; i++)
  63             {
  64                 label vertexNo = vtxLabel(i, j, k);
  65
  66                 // points on edges
  67                 vector edgex1 = p000 + (p100 - p000)*w[0][i];
  68                 vector edgex2 = p010 + (p110 - p010)*w[1][i];
  69                 vector edgex3 = p011 + (p111 - p011)*w[2][i];
  70                 vector edgex4 = p001 + (p101 - p001)*w[3][i];
  71
  72                 vector edgey1 = p000 + (p010 - p000)*w[4][j];
  73                 vector edgey2 = p100 + (p110 - p100)*w[5][j];
  74                 vector edgey3 = p101 + (p111 - p101)*w[6][j];
  75                 vector edgey4 = p001 + (p011 - p001)*w[7][j];
  76
  77                 vector edgez1 = p000 + (p001 - p000)*w[8][k];
  78                 vector edgez2 = p100 + (p101 - p100)*w[9][k];
  79                 vector edgez3 = p110 + (p111 - p110)*w[10][k];
  80                 vector edgez4 = p010 + (p011 - p010)*w[11][k];
  81
  82                 // calculate the importance factors for all edges
  83                 // x - direction
  84                 scalar impx1 =
  85                 (
  86                     (1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
  87                   + w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
  88                 );
  89
  90                 scalar impx2 =
  91                 (
  92                     (1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
  93                   + w[1][i]*w[5][j]*(1.0 - w[10][k])
  94                 );
  95
  96                 scalar impx3 =
  97                 (
  98                      (1.0 - w[2][i])*w[7][j]*w[11][k]
  99                    + w[2][i]*w[6][j]*w[10][k]
 100                 );
 101
 102
 103                 scalar impx4 =
 104                 (
 105                     (1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
 106                   + w[3][i]*(1.0 - w[6][j])*w[9][k]
 107                 );
 108
 109                 scalar magImpx = impx1 + impx2 + impx3 + impx4;
 110                 impx1 /= magImpx;
 111                 impx2 /= magImpx;
 112                 impx3 /= magImpx;
 113                 impx4 /= magImpx;
 114
 115
 116                 // y - direction
 117                 scalar impy1 =
 118                 (
 119                     (1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
 120                   + w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
 121                 );
 122
 123                 scalar impy2 =
 124                 (
 125                     (1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
 126                   + w[5][j]*w[1][i]*(1.0 - w[10][k])
 127                 );
 128
 129                 scalar impy3 =
 130                 (
 131                     (1.0 - w[6][j])*w[3][i]*w[9][k]
 132                   + w[6][j]*w[2][i]*w[10][k]
 133                 );
 134
 135                 scalar impy4 =
 136                 (
 137                     (1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
 138                   + w[7][j]*(1.0 - w[2][i])*w[11][k]
 139                 );
 140
 141                 scalar magImpy = impy1 + impy2 + impy3 + impy4;
 142                 impy1 /= magImpy;
 143                 impy2 /= magImpy;
 144                 impy3 /= magImpy;
 145                 impy4 /= magImpy;
 146
 147
 148                 // z - direction
 149                 scalar impz1 =
 150                 (
 151                     (1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
 152                   + w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
 153                 );
 154
 155                 scalar impz2 =
 156                 (
 157                     (1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
 158                   + w[9][k]*w[3][i]*(1.0 - w[6][j])
 159                 );
 160
 161                 scalar impz3 =
 162                 (
 163                     (1.0 - w[10][k])*w[1][i]*w[5][j]
 164                   + w[10][k]*w[2][i]*w[6][j]
 165                 );
 166
 167                 scalar impz4 =
 168                 (
 169                     (1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
 170                   + w[11][k]*(1.0 - w[2][i])*w[7][j]
 171                 );
 172
 173                 scalar magImpz = impz1 + impz2 + impz3 + impz4;
 174                 impz1 /= magImpz;
 175                 impz2 /= magImpz;
 176                 impz3 /= magImpz;
 177                 impz4 /= magImpz;
 178
 179                 // calculate the correction vectors
 180                 vector corx1 = impx1*(p[0][i] - edgex1);
 181                 vector corx2 = impx2*(p[1][i] - edgex2);
 182                 vector corx3 = impx3*(p[2][i] - edgex3);
 183                 vector corx4 = impx4*(p[3][i] - edgex4);
 184
 185                 vector cory1 = impy1*(p[4][j] - edgey1);
 186                 vector cory2 = impy2*(p[5][j] - edgey2);
 187                 vector cory3 = impy3*(p[6][j] - edgey3);
 188                 vector cory4 = impy4*(p[7][j] - edgey4);
 189
 190                 vector corz1 = impz1*(p[8][k] - edgez1);
 191                 vector corz2 = impz2*(p[9][k] - edgez2);
 192                 vector corz3 = impz3*(p[10][k] - edgez3);
 193                 vector corz4 = impz4*(p[11][k] - edgez4);
 194
 195
 196                 // multiply by the importance factor
 197
 198                 // x - direction
 199                 edgex1 *= impx1;
 200                 edgex2 *= impx2;
 201                 edgex3 *= impx3;
 202                 edgex4 *= impx4;
 203
 204                 // y - direction
 205                 edgey1 *= impy1;
 206                 edgey2 *= impy2;
 207                 edgey3 *= impy3;
 208                 edgey4 *= impy4;
 209
 210                 // z - direction
 211                 edgez1 *= impz1;
 212                 edgez2 *= impz2;
 213                 edgez3 *= impz3;
 214                 edgez4 *= impz4;
 215
 216
 217                 // add the contributions
 218                 vertices_[vertexNo] = edgex1 + edgex2 + edgex3 + edgex4;
 219                 vertices_[vertexNo] += edgey1 + edgey2 + edgey3 + edgey4;
 220                 vertices_[vertexNo] += edgez1 + edgez2 + edgez3 + edgez4;
 221
 222                 vertices_[vertexNo] /= 3.0;
 223
 224                 vertices_[vertexNo] += corx1 + corx2 + corx3 + corx4;
 225                 vertices_[vertexNo] += cory1 + cory2 + cory3 + cory4;
 226                 vertices_[vertexNo] += corz1 + corz2 + corz3 + corz4;
 227             }
 228         }
 229     }
 230 }
 231
 232 // ************************************************************************* //