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