00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "hexBlock.H"
00027
00028
00029
00030 namespace Foam
00031 {
00032
00033
00034
00035 label hexBlock::vtxLabel(label a, label b, label c) const
00036 {
00037 return (a + b*(xDim_ + 1) + c*(xDim_ + 1)*(yDim_ + 1));
00038 }
00039
00040
00041
00042
00043
00044 void hexBlock::setHandedness()
00045 {
00046 const pointField& p = points_;
00047
00048 for (label k = 0; k <= zDim_ - 1; k++)
00049 {
00050 for (label j = 0; j <= yDim_ - 1; j++)
00051 {
00052 for (label i = 0; i <= xDim_ - 1; i++)
00053 {
00054 vector x = p[vtxLabel(i+1, j, k)] - p[vtxLabel(i, j, k)];
00055 vector y = p[vtxLabel(i, j+1, k)] - p[vtxLabel(i, j, k)];
00056 vector z = p[vtxLabel(i, j, k+1)] - p[vtxLabel(i, j, k)];
00057
00058 if (mag(x) > SMALL && mag(y) > SMALL && mag(z) > SMALL)
00059 {
00060 Info<< "Looking at cell "
00061 << i << ' ' << j << ' ' << k
00062 << " to determine orientation."
00063 << endl;
00064
00065 if (((x ^ y) & z) > 0)
00066 {
00067 Info<< "Right-handed block." << endl;
00068 blockHandedness_ = right;
00069 }
00070 else
00071 {
00072 Info << "Left-handed block." << endl;
00073 blockHandedness_ = left;
00074 }
00075 return;
00076 }
00077 else
00078 {
00079 Info<< "Cannot determine orientation of cell "
00080 << i << ' ' << j << ' ' << k
00081 << " since has base vectors " << x << y << z << endl;
00082 }
00083 }
00084 }
00085 }
00086
00087 if (blockHandedness_ == noPoints)
00088 {
00089 WarningIn("hexBlock::hexBlock::setHandedness()")
00090 << "Cannot determine orientation of block."
00091 << " Continuing as if right handed." << endl;
00092 blockHandedness_ = right;
00093 }
00094 }
00095
00096
00097
00098
00099
00100 hexBlock::hexBlock(const label nx, const label ny, const label nz)
00101 :
00102 xDim_(nx - 1),
00103 yDim_(ny - 1),
00104 zDim_(nz - 1),
00105 blockHandedness_(noPoints),
00106 points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
00107 {}
00108
00109
00110
00111
00112 void hexBlock::readPoints
00113 (
00114 const bool readBlank,
00115 const scalar twoDThicknes,
00116 Istream& is
00117 )
00118 {
00119 scalar iBlank;
00120
00121 label nPoints = points_.size();
00122
00123 if (twoDThicknes > 0)
00124 {
00125 nPoints /= 2;
00126 }
00127
00128 Info<< "Reading " << nPoints << " x coordinates..." << endl;
00129 for (label i=0; i < nPoints; i++)
00130 {
00131 is >> points_[i].x();
00132 }
00133
00134 Info<< "Reading " << nPoints << " y coordinates..." << endl;
00135 for (label i=0; i < nPoints; i++)
00136 {
00137 is >> points_[i].y();
00138 }
00139
00140 if (twoDThicknes > 0)
00141 {
00142 Info<< "Extruding " << nPoints << " points in z direction..." << endl;
00143
00144 for (label i=0; i < nPoints; i++)
00145 {
00146 points_[i+nPoints] = points_[i];
00147 }
00148 for (label i=0; i < nPoints; i++)
00149 {
00150 points_[i].z() = 0;
00151 points_[i+nPoints].z() = twoDThicknes;
00152 }
00153 }
00154 else
00155 {
00156 Info<< "Reading " << nPoints << " z coordinates..." << endl;
00157 for (label i=0; i < nPoints; i++)
00158 {
00159 is >> points_[i].z();
00160 }
00161 }
00162
00163
00164 if (readBlank)
00165 {
00166 Info<< "Reading " << nPoints << " blanks..." << endl;
00167 for (label i=0; i < nPoints; i++)
00168 {
00169 is >> iBlank;
00170 }
00171 }
00172
00173
00174 setHandedness();
00175 }
00176
00177
00178 labelListList hexBlock::blockCells() const
00179 {
00180 labelListList result(xDim_*yDim_*zDim_);
00181
00182 label cellNo = 0;
00183
00184 if (blockHandedness_ == right)
00185 {
00186 for (label k = 0; k <= zDim_ - 1; k++)
00187 {
00188 for (label j = 0; j <= yDim_ - 1; j++)
00189 {
00190 for (label i = 0; i <= xDim_ - 1; i++)
00191 {
00192 labelList& hexLabels = result[cellNo];
00193 hexLabels.setSize(8);
00194
00195 hexLabels[0] = vtxLabel(i, j, k);
00196 hexLabels[1] = vtxLabel(i+1, j, k);
00197 hexLabels[2] = vtxLabel(i+1, j+1, k);
00198 hexLabels[3] = vtxLabel(i, j+1, k);
00199 hexLabels[4] = vtxLabel(i, j, k+1);
00200 hexLabels[5] = vtxLabel(i+1, j, k+1);
00201 hexLabels[6] = vtxLabel(i+1, j+1, k+1);
00202 hexLabels[7] = vtxLabel(i, j+1, k+1);
00203
00204 cellNo++;
00205 }
00206 }
00207 }
00208 }
00209 else if (blockHandedness_ == left)
00210 {
00211 for (label k = 0; k <= zDim_ - 1; k++)
00212 {
00213 for (label j = 0; j <= yDim_ - 1; j++)
00214 {
00215 for (label i = 0; i <= xDim_ - 1; i++)
00216 {
00217 labelList& hexLabels = result[cellNo];
00218 hexLabels.setSize(8);
00219
00220 hexLabels[0] = vtxLabel(i, j, k+1);
00221 hexLabels[1] = vtxLabel(i+1, j, k+1);
00222 hexLabels[2] = vtxLabel(i+1, j+1, k+1);
00223 hexLabels[3] = vtxLabel(i, j+1, k+1);
00224 hexLabels[4] = vtxLabel(i, j, k);
00225 hexLabels[5] = vtxLabel(i+1, j, k);
00226 hexLabels[6] = vtxLabel(i+1, j+1, k);
00227 hexLabels[7] = vtxLabel(i, j+1, k);
00228
00229 cellNo++;
00230 }
00231 }
00232 }
00233 }
00234 else
00235 {
00236 FatalErrorIn("hexBlock::cellShapes()")
00237 << "Unable to determine block handedness as points "
00238 << "have not been read in yet"
00239 << abort(FatalError);
00240 }
00241
00242 return result;
00243 }
00244
00245
00246
00247
00248
00249
00250
00251 faceList hexBlock::patchFaces(const label direc, const labelList& range) const
00252 {
00253 if (range.size() != 6)
00254 {
00255 FatalErrorIn
00256 (
00257 "patchFaces(const label direc, const labelList& range) const"
00258 ) << "Invalid size of the range array: " << range.size()
00259 << ". Should be 6 (xMin, xMax, yMin, yMax, zMin, zMax"
00260 << abort(FatalError);
00261 }
00262
00263 label xMinRange = range[0];
00264 label xMaxRange = range[1];
00265 label yMinRange = range[2];
00266 label yMaxRange = range[3];
00267 label zMinRange = range[4];
00268 label zMaxRange = range[5];
00269
00270 faceList result(0);
00271
00272 switch (direc)
00273 {
00274 case 1:
00275 {
00276
00277
00278 result.setSize
00279 (
00280 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
00281 );
00282
00283 label p = 0;
00284 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00285 {
00286 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00287 {
00288 result[p].setSize(4);
00289
00290
00291 result[p][0] = vtxLabel(xDim_, j, k);
00292 result[p][1] = vtxLabel(xDim_, j+1, k);
00293 result[p][2] = vtxLabel(xDim_, j+1, k+1);
00294 result[p][3] = vtxLabel(xDim_, j, k+1);
00295
00296 p++;
00297 }
00298 }
00299
00300 result.setSize(p);
00301 break;
00302 }
00303
00304 case 2:
00305 {
00306
00307 result.setSize
00308 (
00309 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
00310 );
00311
00312 label p = 0;
00313 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00314 {
00315 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00316 {
00317 result[p].setSize(4);
00318
00319
00320 result[p][0] = vtxLabel(i, yDim_, k);
00321 result[p][1] = vtxLabel(i, yDim_, k + 1);
00322 result[p][2] = vtxLabel(i + 1, yDim_, k + 1);
00323 result[p][3] = vtxLabel(i + 1, yDim_, k);
00324
00325 p++;
00326 }
00327 }
00328
00329 result.setSize(p);
00330 break;
00331 }
00332
00333 case 3:
00334 {
00335
00336 result.setSize
00337 (
00338 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
00339 );
00340
00341 label p = 0;
00342 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00343 {
00344 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00345 {
00346 result[p].setSize(4);
00347
00348
00349 result[p][0] = vtxLabel(i, j, zDim_);
00350 result[p][1] = vtxLabel(i + 1, j, zDim_);
00351 result[p][2] = vtxLabel(i + 1, j + 1, zDim_);
00352 result[p][3] = vtxLabel(i, j + 1, zDim_);
00353
00354 p++;
00355 }
00356 }
00357
00358 result.setSize(p);
00359 break;
00360 }
00361
00362 case 4:
00363 {
00364
00365 result.setSize
00366 (
00367 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
00368 );
00369
00370 label p = 0;
00371 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00372 {
00373 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00374 {
00375 result[p].setSize(4);
00376
00377
00378 result[p][0] = vtxLabel(0, j, k);
00379 result[p][1] = vtxLabel(0, j, k + 1);
00380 result[p][2] = vtxLabel(0, j + 1, k + 1);
00381 result[p][3] = vtxLabel(0, j + 1, k);
00382
00383 p++;
00384 }
00385 }
00386
00387 result.setSize(p);
00388 break;
00389 }
00390
00391 case 5:
00392 {
00393
00394 result.setSize
00395 (
00396 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
00397 );
00398
00399 label p = 0;
00400 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00401 {
00402 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00403 {
00404 result[p].setSize(4);
00405
00406
00407 result[p][0] = vtxLabel(i, 0, k);
00408 result[p][1] = vtxLabel(i + 1, 0, k);
00409 result[p][2] = vtxLabel(i + 1, 0, k + 1);
00410 result[p][3] = vtxLabel(i, 0, k + 1);
00411
00412 p++;
00413 }
00414 }
00415
00416 result.setSize(p);
00417 break;
00418 }
00419
00420 case 6:
00421 {
00422
00423 result.setSize
00424 (
00425 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
00426 );
00427
00428 label p = 0;
00429 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00430 {
00431 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00432 {
00433 result[p].setSize(4);
00434
00435
00436 result[p][0] = vtxLabel(i, j, 0);
00437 result[p][1] = vtxLabel(i, j + 1, 0);
00438 result[p][2] = vtxLabel(i + 1, j + 1, 0);
00439 result[p][3] = vtxLabel(i + 1, j, 0);
00440
00441 p++;
00442 }
00443 }
00444
00445 result.setSize(p);
00446 break;
00447 }
00448
00449 default:
00450 {
00451 FatalErrorIn
00452 (
00453 "patchFaces(const label direc, const labelList& range) const"
00454 ) << "direction out of range (1 to 6): " << direc
00455 << abort(FatalError);
00456 }
00457 }
00458
00459
00460
00461 if (blockHandedness_ == noPoints)
00462 {
00463 FatalErrorIn
00464 (
00465 "patchFaces(const label direc, const labelList& range) const"
00466 ) << "Unable to determine block handedness as points "
00467 << "have not been read in yet"
00468 << abort(FatalError);
00469 }
00470 else if (blockHandedness_ == left)
00471 {
00472
00473 forAll (result, faceI)
00474 {
00475 result[faceI] = result[faceI].reverseFace();
00476 }
00477 }
00478
00479 return result;
00480 }
00481
00482
00483
00484
00485 }
00486
00487