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
00027
00028
00029 #include "polyMesh.H"
00030 #include <OpenFOAM/Time.H>
00031 #include <OpenFOAM/primitiveMesh.H>
00032 #include <OpenFOAM/DynamicList.H>
00033
00034
00035
00036 Foam::labelListList Foam::polyMesh::cellShapePointCells
00037 (
00038 const cellShapeList& c
00039 ) const
00040 {
00041 List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
00042 pc(points().size());
00043
00044
00045 forAll(c, i)
00046 {
00047
00048 const labelList& labels = c[i];
00049
00050 forAll(labels, j)
00051 {
00052
00053 label curPoint = labels[j];
00054 DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
00055 pc[curPoint];
00056
00057
00058 curPointCells.append(i);
00059 }
00060 }
00061
00062 labelListList pointCellAddr(pc.size());
00063
00064 forAll (pc, pointI)
00065 {
00066 pointCellAddr[pointI].transfer(pc[pointI]);
00067 }
00068
00069 return pointCellAddr;
00070 }
00071
00072
00073 Foam::labelList Foam::polyMesh::facePatchFaceCells
00074 (
00075 const faceList& patchFaces,
00076 const labelListList& pointCells,
00077 const faceListList& cellsFaceShapes,
00078 const label patchID
00079 ) const
00080 {
00081 register bool found;
00082
00083 labelList FaceCells(patchFaces.size());
00084
00085 forAll(patchFaces, fI)
00086 {
00087 found = false;
00088
00089 const face& curFace = patchFaces[fI];
00090 const labelList& facePoints = patchFaces[fI];
00091
00092 forAll(facePoints, pointI)
00093 {
00094 const labelList& facePointCells = pointCells[facePoints[pointI]];
00095
00096 forAll(facePointCells, cellI)
00097 {
00098 faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
00099
00100 forAll(cellFaces, cellFace)
00101 {
00102 if (cellFaces[cellFace] == curFace)
00103 {
00104
00105 FaceCells[fI] = facePointCells[cellI];
00106
00107 found = true;
00108 }
00109 if (found) break;
00110 }
00111 if (found) break;
00112 }
00113 if (found) break;
00114 }
00115
00116 if (!found)
00117 {
00118 FatalErrorIn
00119 (
00120 "polyMesh::facePatchFaceCells(const faceList& patchFaces,"
00121 "const labelListList& pointCells,"
00122 "const faceListList& cellsFaceShapes,"
00123 "const label patchID)"
00124 ) << "face " << fI << " in patch " << patchID
00125 << " does not have neighbour cell"
00126 << " face: " << patchFaces[fI]
00127 << abort(FatalError);
00128 }
00129 }
00130
00131 return FaceCells;
00132 }
00133
00134
00135 Foam::polyMesh::polyMesh
00136 (
00137 const IOobject& io,
00138 const Xfer<pointField>& points,
00139 const cellShapeList& cellsAsShapes,
00140 const faceListList& boundaryFaces,
00141 const wordList& boundaryPatchNames,
00142 const wordList& boundaryPatchTypes,
00143 const word& defaultBoundaryPatchName,
00144 const word& defaultBoundaryPatchType,
00145 const wordList& boundaryPatchPhysicalTypes,
00146 const bool syncPar
00147 )
00148 :
00149 objectRegistry(io),
00150 primitiveMesh(),
00151 points_
00152 (
00153 IOobject
00154 (
00155 "points",
00156 instance(),
00157 meshSubDir,
00158 *this,
00159 IOobject::NO_READ,
00160 IOobject::AUTO_WRITE
00161 ),
00162 points
00163 ),
00164 faces_
00165 (
00166 IOobject
00167 (
00168 "faces",
00169 instance(),
00170 meshSubDir,
00171 *this,
00172 IOobject::NO_READ,
00173 IOobject::AUTO_WRITE
00174 ),
00175 0
00176 ),
00177 owner_
00178 (
00179 IOobject
00180 (
00181 "owner",
00182 instance(),
00183 meshSubDir,
00184 *this,
00185 IOobject::NO_READ,
00186 IOobject::AUTO_WRITE
00187 ),
00188 0
00189 ),
00190 neighbour_
00191 (
00192 IOobject
00193 (
00194 "neighbour",
00195 instance(),
00196 meshSubDir,
00197 *this,
00198 IOobject::NO_READ,
00199 IOobject::AUTO_WRITE
00200 ),
00201 0
00202 ),
00203 clearedPrimitives_(false),
00204 boundary_
00205 (
00206 IOobject
00207 (
00208 "boundary",
00209 instance(),
00210 meshSubDir,
00211 *this,
00212 IOobject::NO_READ,
00213 IOobject::AUTO_WRITE
00214 ),
00215 *this,
00216 boundaryFaces.size() + 1
00217 ),
00218 bounds_(points_, syncPar),
00219 geometricD_(Vector<label>::zero),
00220 solutionD_(Vector<label>::zero),
00221 pointZones_
00222 (
00223 IOobject
00224 (
00225 "pointZones",
00226 instance(),
00227 meshSubDir,
00228 *this,
00229 IOobject::NO_READ,
00230 IOobject::NO_WRITE
00231 ),
00232 *this,
00233 0
00234 ),
00235 faceZones_
00236 (
00237 IOobject
00238 (
00239 "faceZones",
00240 instance(),
00241 meshSubDir,
00242 *this,
00243 IOobject::NO_READ,
00244 IOobject::NO_WRITE
00245 ),
00246 *this,
00247 0
00248 ),
00249 cellZones_
00250 (
00251 IOobject
00252 (
00253 "cellZones",
00254 instance(),
00255 meshSubDir,
00256 *this,
00257 IOobject::NO_READ,
00258 IOobject::NO_WRITE
00259 ),
00260 *this,
00261 0
00262 ),
00263 globalMeshDataPtr_(NULL),
00264 moving_(false),
00265 curMotionTimeIndex_(time().timeIndex()),
00266 oldPointsPtr_(NULL)
00267 {
00268 if (debug)
00269 {
00270 Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
00271 }
00272
00273
00274 removeFiles(instance());
00275
00276
00277
00278 label maxFaces = 0;
00279
00280
00281 faceListList cellsFaceShapes(cellsAsShapes.size());
00282 cellList cells(cellsAsShapes.size());
00283
00284 forAll(cellsFaceShapes, cellI)
00285 {
00286 cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
00287
00288 cells[cellI].setSize(cellsFaceShapes[cellI].size());
00289
00290
00291 static_cast<labelList&>(cells[cellI]) = -1;
00292
00293
00294 maxFaces += cellsFaceShapes[cellI].size();
00295 }
00296
00297
00298 faces_.setSize(maxFaces);
00299
00300
00301 label nFaces = 0;
00302
00303
00304 labelListList PointCells = cellShapePointCells(cellsAsShapes);
00305
00306 bool found = false;
00307
00308 forAll(cells, cellI)
00309 {
00310
00311
00312
00313
00314
00315
00316 const faceList& curFaces = cellsFaceShapes[cellI];
00317
00318
00319 labelList neiCells(curFaces.size(), -1);
00320
00321
00322 labelList faceOfNeiCell(curFaces.size(), -1);
00323
00324 label nNeighbours = 0;
00325
00326
00327 forAll(curFaces, faceI)
00328 {
00329
00330 if (cells[cellI][faceI] >= 0) continue;
00331
00332 found = false;
00333
00334 const face& curFace = curFaces[faceI];
00335
00336
00337 const labelList& curPoints = curFace;
00338
00339
00340 forAll(curPoints, pointI)
00341 {
00342
00343 const labelList& curNeighbours =
00344 PointCells[curPoints[pointI]];
00345
00346
00347 forAll(curNeighbours, neiI)
00348 {
00349 label curNei = curNeighbours[neiI];
00350
00351
00352 if (curNei > cellI)
00353 {
00354
00355 const faceList& searchFaces = cellsFaceShapes[curNei];
00356
00357 forAll(searchFaces, neiFaceI)
00358 {
00359 if (searchFaces[neiFaceI] == curFace)
00360 {
00361
00362 found = true;
00363
00364
00365 neiCells[faceI] = curNei;
00366 faceOfNeiCell[faceI] = neiFaceI;
00367 nNeighbours++;
00368
00369 break;
00370 }
00371 }
00372 if (found) break;
00373 }
00374 if (found) break;
00375 }
00376 if (found) break;
00377 }
00378 }
00379
00380
00381 for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
00382 {
00383
00384 label nextNei = -1;
00385 label minNei = cells.size();
00386
00387 forAll (neiCells, ncI)
00388 {
00389 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
00390 {
00391 nextNei = ncI;
00392 minNei = neiCells[ncI];
00393 }
00394 }
00395
00396 if (nextNei > -1)
00397 {
00398
00399 faces_[nFaces] = curFaces[nextNei];
00400
00401
00402 cells[cellI][nextNei] = nFaces;
00403 cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
00404
00405
00406 neiCells[nextNei] = -1;
00407
00408
00409 nFaces++;
00410 }
00411 else
00412 {
00413 FatalErrorIn
00414 (
00415 "polyMesh::polyMesh\n"
00416 "(\n"
00417 " const IOobject&,\n"
00418 " const Xfer<pointField>&,\n"
00419 " const cellShapeList& cellsAsShapes,\n"
00420 " const faceListList& boundaryFaces,\n"
00421 " const wordList& boundaryPatchTypes,\n"
00422 " const wordList& boundaryPatchNames,\n"
00423 " const word& defaultBoundaryPatchType\n"
00424 ")"
00425 ) << "Error in internal face insertion"
00426 << abort(FatalError);
00427 }
00428 }
00429 }
00430
00431
00432
00433 labelList patchSizes(boundaryFaces.size(), -1);
00434 labelList patchStarts(boundaryFaces.size(), -1);
00435
00436 forAll (boundaryFaces, patchI)
00437 {
00438 const faceList& patchFaces = boundaryFaces[patchI];
00439
00440 labelList curPatchFaceCells =
00441 facePatchFaceCells
00442 (
00443 patchFaces,
00444 PointCells,
00445 cellsFaceShapes,
00446 patchI
00447 );
00448
00449
00450 label curPatchStart = nFaces;
00451
00452 forAll (patchFaces, faceI)
00453 {
00454 const face& curFace = patchFaces[faceI];
00455
00456 const label cellInside = curPatchFaceCells[faceI];
00457
00458 faces_[nFaces] = curFace;
00459
00460
00461 const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
00462
00463 bool found = false;
00464
00465 forAll (facesOfCellInside, cellFaceI)
00466 {
00467 if (facesOfCellInside[cellFaceI] == curFace)
00468 {
00469 if (cells[cellInside][cellFaceI] >= 0)
00470 {
00471 FatalErrorIn
00472 (
00473 "polyMesh::polyMesh\n"
00474 "(\n"
00475 " const IOobject&,\n"
00476 " const Xfer<pointField>&,\n"
00477 " const cellShapeList& cellsAsShapes,\n"
00478 " const faceListList& boundaryFaces,\n"
00479 " const wordList& boundaryPatchTypes,\n"
00480 " const wordList& boundaryPatchNames,\n"
00481 " const word& defaultBoundaryPatchType\n"
00482 ")"
00483 ) << "Trying to specify a boundary face " << curFace
00484 << " on the face on cell " << cellInside
00485 << " which is either an internal face or already "
00486 << "belongs to some other patch. This is face "
00487 << faceI << " of patch "
00488 << patchI << " named "
00489 << boundaryPatchNames[patchI] << "."
00490 << abort(FatalError);
00491 }
00492
00493 found = true;
00494
00495 cells[cellInside][cellFaceI] = nFaces;
00496
00497 break;
00498 }
00499 }
00500
00501 if (!found)
00502 {
00503 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
00504 << "face " << faceI << " of patch " << patchI
00505 << " does not seem to belong to cell " << cellInside
00506 << " which, according to the addressing, "
00507 << "should be next to it."
00508 << abort(FatalError);
00509 }
00510
00511
00512 nFaces++;
00513 }
00514
00515 patchSizes[patchI] = nFaces - curPatchStart;
00516 patchStarts[patchI] = curPatchStart;
00517 }
00518
00519
00520
00521 label defaultPatchStart = nFaces;
00522
00523 forAll(cells, cellI)
00524 {
00525 labelList& curCellFaces = cells[cellI];
00526
00527 forAll(curCellFaces, faceI)
00528 {
00529 if (curCellFaces[faceI] == -1)
00530 {
00531 curCellFaces[faceI] = nFaces;
00532 faces_[nFaces] = cellsFaceShapes[cellI][faceI];
00533
00534 nFaces++;
00535 }
00536 }
00537 }
00538
00539
00540 faces_.setSize(nFaces);
00541
00542
00543
00544 forAll (boundaryFaces, patchI)
00545 {
00546
00547 boundary_.set
00548 (
00549 patchI,
00550 polyPatch::New
00551 (
00552 boundaryPatchTypes[patchI],
00553 boundaryPatchNames[patchI],
00554 patchSizes[patchI],
00555 patchStarts[patchI],
00556 patchI,
00557 boundary_
00558 )
00559 );
00560
00561 if
00562 (
00563 boundaryPatchPhysicalTypes.size()
00564 && boundaryPatchPhysicalTypes[patchI].size()
00565 )
00566 {
00567 boundary_[patchI].physicalType() =
00568 boundaryPatchPhysicalTypes[patchI];
00569 }
00570 }
00571
00572 label nAllPatches = boundaryFaces.size();
00573
00574 if (nFaces > defaultPatchStart)
00575 {
00576 WarningIn("polyMesh::polyMesh(... construct from shapes...)")
00577 << "Found " << nFaces - defaultPatchStart
00578 << " undefined faces in mesh; adding to default patch." << endl;
00579
00580 boundary_.set
00581 (
00582 nAllPatches,
00583 polyPatch::New
00584 (
00585 defaultBoundaryPatchType,
00586 defaultBoundaryPatchName,
00587 nFaces - defaultPatchStart,
00588 defaultPatchStart,
00589 boundary_.size() - 1,
00590 boundary_
00591 )
00592 );
00593
00594 nAllPatches++;
00595 }
00596
00597
00598 boundary_.setSize(nAllPatches);
00599
00600
00601 initMesh(cells);
00602
00603 if (syncPar)
00604 {
00605
00606 boundary_.updateMesh();
00607
00608
00609 boundary_.calcGeometry();
00610 }
00611
00612 if (debug)
00613 {
00614 if (checkMesh())
00615 {
00616 Info << "Mesh OK" << endl;
00617 }
00618 }
00619 }
00620
00621
00622