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 "localPointRegion.H"
00027 #include <OpenFOAM/syncTools.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mapPolyMesh.H>
00030 #include <OpenFOAM/globalIndex.H>
00031 #include <OpenFOAM/indirectPrimitivePatch.H>
00032
00033
00034
00035 namespace Foam
00036 {
00037
00038 defineTypeNameAndDebug(localPointRegion, 0);
00039
00040
00041 class minEqOpFace
00042 {
00043 public:
00044
00045 void operator()(face& x, const face& y) const
00046 {
00047 if (x.size())
00048 {
00049 label j = 0;
00050 forAll(x, i)
00051 {
00052 x[i] = min(x[i], y[j]);
00053
00054 j = y.rcIndex(j);
00055 }
00056 }
00057 };
00058 };
00059
00060
00061
00062 void transformList
00063 (
00064 const tensorField& rotTensor,
00065 UList<face>& field
00066 )
00067 {};
00068
00069 }
00070
00071
00072
00073
00074 bool Foam::localPointRegion::isDuplicate
00075 (
00076 const face& f0,
00077 const face& f1,
00078 const bool forward
00079 )
00080 {
00081 label fp1 = findIndex(f1, f0[0]);
00082
00083 if (fp1 == -1)
00084 {
00085 return false;
00086 }
00087
00088 forAll(f0, fp0)
00089 {
00090 if (f0[fp0] != f1[fp1])
00091 {
00092 return false;
00093 }
00094
00095 if (forward)
00096 {
00097 fp1 = f1.fcIndex(fp1);
00098 }
00099 else
00100 {
00101 fp1 = f1.rcIndex(fp1);
00102 }
00103 }
00104 return true;
00105 }
00106
00107
00108
00109 void Foam::localPointRegion::countPointRegions
00110 (
00111 const polyMesh& mesh,
00112 const boolList& candidatePoint,
00113 const Map<label>& candidateFace,
00114 faceList& minRegion
00115 )
00116 {
00117
00118
00119 labelList minPointRegion(mesh.nPoints(), -1);
00120
00121 meshPointMap_.resize(candidateFace.size()/100);
00122
00123 DynamicList<labelList> pointRegions(meshPointMap_.size());
00124
00125
00126 meshFaceMap_.resize(meshPointMap_.size());
00127
00128 forAllConstIter(Map<label>, candidateFace, iter)
00129 {
00130 label faceI = iter.key();
00131
00132 if (!mesh.isInternalFace(faceI))
00133 {
00134 const face& f = mesh.faces()[faceI];
00135
00136 if (minRegion[faceI].empty())
00137 {
00138 FatalErrorIn("localPointRegion::countPointRegions(..)")
00139 << "Face from candidateFace without minRegion set." << endl
00140 << "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
00141 << " verts:" << f << abort(FatalError);
00142 }
00143
00144 forAll(f, fp)
00145 {
00146 label pointI = f[fp];
00147
00148
00149
00150
00151 if (candidatePoint[pointI])
00152 {
00153 label region = minRegion[faceI][fp];
00154
00155 if (minPointRegion[pointI] == -1)
00156 {
00157 minPointRegion[pointI] = region;
00158 }
00159 else if (minPointRegion[pointI] != region)
00160 {
00161
00162 Map<label>::iterator iter = meshPointMap_.find(pointI);
00163 if (iter != meshPointMap_.end())
00164 {
00165 labelList& regions = pointRegions[iter()];
00166 if (findIndex(regions, region) == -1)
00167 {
00168 label sz = regions.size();
00169 regions.setSize(sz+1);
00170 regions[sz] = region;
00171 }
00172 }
00173 else
00174 {
00175 label localPointI = meshPointMap_.size();
00176 meshPointMap_.insert(pointI, localPointI);
00177 labelList regions(2);
00178 regions[0] = minPointRegion[pointI];
00179 regions[1] = region;
00180 pointRegions.append(regions);
00181 }
00182
00183 label meshFaceMapI = meshFaceMap_.size();
00184 meshFaceMap_.insert(faceI, meshFaceMapI);
00185 }
00186 }
00187 }
00188 }
00189 }
00190 minPointRegion.clear();
00191
00192
00193
00194 forAllConstIter(Map<label>, candidateFace, iter)
00195 {
00196 label faceI = iter.key();
00197
00198 if (mesh.isInternalFace(faceI))
00199 {
00200 const face& f = mesh.faces()[faceI];
00201
00202 forAll(f, fp)
00203 {
00204
00205
00206 if (candidatePoint[f[fp]] && meshPointMap_.found(f[fp]))
00207 {
00208 label meshFaceMapI = meshFaceMap_.size();
00209 meshFaceMap_.insert(faceI, meshFaceMapI);
00210 }
00211 }
00212 }
00213 }
00214
00215
00216
00217 pointRegions.shrink();
00218 pointRegions_.setSize(pointRegions.size());
00219 forAll(pointRegions, i)
00220 {
00221 pointRegions_[i].transfer(pointRegions[i]);
00222 }
00223
00224
00225 faceRegions_.setSize(meshFaceMap_.size());
00226 forAllConstIter(Map<label>, meshFaceMap_, iter)
00227 {
00228 faceRegions_[iter()].labelList::transfer(minRegion[iter.key()]);
00229
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 }
00244
00245
00246
00247 }
00248
00249
00250 void Foam::localPointRegion::calcPointRegions
00251 (
00252 const polyMesh& mesh,
00253 boolList& candidatePoint
00254 )
00255 {
00256 label nBnd = mesh.nFaces()-mesh.nInternalFaces();
00257 const labelList& faceOwner = mesh.faceOwner();
00258 const labelList& faceNeighbour = mesh.faceNeighbour();
00259
00260
00261 syncTools::syncPointList
00262 (
00263 mesh,
00264 candidatePoint,
00265 orEqOp<bool>(),
00266 false,
00267 false
00268 );
00269
00270
00271
00272
00273
00274 Map<label> candidateFace(2*nBnd);
00275 label candidateFaceI = 0;
00276
00277 Map<label> candidateCell(nBnd);
00278 label candidateCellI = 0;
00279
00280 forAll(mesh.faces(), faceI)
00281 {
00282 const face& f = mesh.faces()[faceI];
00283
00284 forAll(f, fp)
00285 {
00286 if (candidatePoint[f[fp]])
00287 {
00288
00289 if (candidateFace.insert(faceI, candidateFaceI))
00290 {
00291 candidateFaceI++;
00292 }
00293
00294
00295 if (candidateCell.insert(faceOwner[faceI], candidateCellI))
00296 {
00297 candidateCellI++;
00298 }
00299
00300 if (mesh.isInternalFace(faceI))
00301 {
00302 label nei = faceNeighbour[faceI];
00303 if (candidateCell.insert(nei, candidateCellI))
00304 {
00305 candidateCellI++;
00306 }
00307 }
00308
00309 break;
00310 }
00311 }
00312 }
00313
00314
00315
00316
00317 globalIndex globalCells(mesh.nCells());
00318
00319
00320
00321
00322
00323
00324 faceList minRegion(mesh.nFaces());
00325 forAllConstIter(Map<label>, candidateFace, iter)
00326 {
00327 label faceI = iter.key();
00328 const face& f = mesh.faces()[faceI];
00329
00330 if (mesh.isInternalFace(faceI))
00331 {
00332 label globOwn = globalCells.toGlobal(faceOwner[faceI]);
00333 label globNei = globalCells.toGlobal(faceNeighbour[faceI]);
00334 minRegion[faceI].setSize(f.size(), min(globOwn, globNei));
00335 }
00336 else
00337 {
00338 label globOwn = globalCells.toGlobal(faceOwner[faceI]);
00339 minRegion[faceI].setSize(f.size(), globOwn);
00340 }
00341 }
00342
00343
00344
00345
00346
00347 while (true)
00348 {
00349
00350
00351
00352 Map<label> minPointValue(100);
00353 label nChanged = 0;
00354 forAllConstIter(Map<label>, candidateCell, iter)
00355 {
00356 minPointValue.clear();
00357
00358 label cellI = iter.key();
00359 const cell& cFaces = mesh.cells()[cellI];
00360
00361
00362 forAll(cFaces, cFaceI)
00363 {
00364 label faceI = cFaces[cFaceI];
00365
00366 if (minRegion[faceI].size())
00367 {
00368 const face& f = mesh.faces()[faceI];
00369
00370 forAll(f, fp)
00371 {
00372 label pointI = f[fp];
00373 Map<label>::iterator iter = minPointValue.find(pointI);
00374
00375 if (iter == minPointValue.end())
00376 {
00377 minPointValue.insert(pointI, minRegion[faceI][fp]);
00378 }
00379 else
00380 {
00381 label currentMin = iter();
00382 iter() = min(currentMin, minRegion[faceI][fp]);
00383 }
00384 }
00385 }
00386 }
00387
00388
00389 forAll(cFaces, cFaceI)
00390 {
00391 label faceI = cFaces[cFaceI];
00392
00393 if (minRegion[faceI].size())
00394 {
00395 const face& f = mesh.faces()[faceI];
00396
00397 forAll(f, fp)
00398 {
00399 label minVal = minPointValue[f[fp]];
00400
00401 if (minVal != minRegion[faceI][fp])
00402 {
00403 minRegion[faceI][fp] = minVal;
00404 nChanged++;
00405 }
00406 }
00407 }
00408 }
00409 }
00410
00411
00412
00413 if (returnReduce(nChanged, sumOp<label>()) == 0)
00414 {
00415 break;
00416 }
00417
00418
00419
00420
00421
00422 syncTools::syncFaceList
00423 (
00424 mesh,
00425 minRegion,
00426 minEqOpFace(),
00427 false
00428 );
00429 }
00430
00431
00432
00433 countPointRegions(mesh, candidatePoint, candidateFace, minRegion);
00434 minRegion.clear();
00435
00436
00438
00439
00440
00441
00442
00443
00444 }
00445
00446
00447
00448
00449 Foam::localPointRegion::localPointRegion(const polyMesh& mesh)
00450 :
00451
00452 meshPointMap_(0),
00453 pointRegions_(0),
00454 meshFaceMap_(0),
00455 faceRegions_(0)
00456 {
00457 const polyBoundaryMesh& patches = mesh.boundaryMesh();
00458
00459
00460 boolList candidatePoint(mesh.nPoints(), false);
00461
00462 forAll(patches, patchI)
00463 {
00464 if (!patches[patchI].coupled())
00465 {
00466 const polyPatch& pp = patches[patchI];
00467
00468 forAll(pp.meshPoints(), i)
00469 {
00470 candidatePoint[pp.meshPoints()[i]] = true;
00471 }
00472 }
00473 }
00474
00475 calcPointRegions(mesh, candidatePoint);
00476 }
00477
00478
00479 Foam::localPointRegion::localPointRegion
00480 (
00481 const polyMesh& mesh,
00482 const labelList& candidatePoints
00483 )
00484 :
00485
00486 meshPointMap_(0),
00487 pointRegions_(0),
00488 meshFaceMap_(0),
00489 faceRegions_(0)
00490 {
00491 boolList candidatePoint(mesh.nPoints(), false);
00492
00493 forAll(candidatePoints, i)
00494 {
00495 candidatePoint[candidatePoints[i]] = true;
00496 }
00497
00498 calcPointRegions(mesh, candidatePoint);
00499 }
00500
00501
00502
00503
00504
00505
00506 Foam::labelList Foam::localPointRegion::findDuplicateFaces
00507 (
00508 const primitiveMesh& mesh,
00509 const labelList& boundaryFaces
00510 )
00511 {
00512
00513 indirectPrimitivePatch allPatch
00514 (
00515 IndirectList<face>(mesh.faces(), boundaryFaces),
00516 mesh.points()
00517 );
00518
00519 labelList duplicateFace(allPatch.size(), -1);
00520 label nDuplicateFaces = 0;
00521
00522
00523 forAll(allPatch, bFaceI)
00524 {
00525 const face& f = allPatch.localFaces()[bFaceI];
00526
00527
00528
00529 const labelList& pFaces = allPatch.pointFaces()[f[0]];
00530
00531 forAll(pFaces, i)
00532 {
00533 label otherFaceI = pFaces[i];
00534
00535 if (otherFaceI > bFaceI)
00536 {
00537 const face& otherF = allPatch.localFaces()[otherFaceI];
00538
00539 if (isDuplicate(f, otherF, true))
00540 {
00541 FatalErrorIn
00542 (
00543 "findDuplicateFaces(const primitiveMesh&"
00544 ", const labelList&)"
00545 ) << "Face:" << bFaceI + mesh.nInternalFaces()
00546 << " has local points:" << f
00547 << " which are in same order as face:"
00548 << otherFaceI + mesh.nInternalFaces()
00549 << " with local points:" << otherF
00550 << abort(FatalError);
00551 }
00552 else if (isDuplicate(f, otherF, false))
00553 {
00554 label meshFace0 = bFaceI + mesh.nInternalFaces();
00555 label meshFace1 = otherFaceI + mesh.nInternalFaces();
00556
00557 if
00558 (
00559 duplicateFace[bFaceI] != -1
00560 || duplicateFace[otherFaceI] != -1
00561 )
00562 {
00563 FatalErrorIn
00564 (
00565 "findDuplicateFaces(const primitiveMesh&"
00566 ", const labelList&)"
00567 ) << "One of two duplicate faces already marked"
00568 << " as duplicate." << nl
00569 << "This means that three or more faces share"
00570 << " the same points and this is illegal." << nl
00571 << "Face:" << meshFace0
00572 << " with local points:" << f
00573 << " which are in same order as face:"
00574 << meshFace1
00575 << " with local points:" << otherF
00576 << abort(FatalError);
00577 }
00578
00579 duplicateFace[bFaceI] = otherFaceI;
00580 duplicateFace[otherFaceI] = bFaceI;
00581 nDuplicateFaces++;
00582 }
00583 }
00584 }
00585 }
00586
00587 return duplicateFace;
00588 }
00589
00590
00591 void Foam::localPointRegion::updateMesh(const mapPolyMesh& map)
00592 {
00593 {
00594 Map<label> newMap(meshFaceMap_.size());
00595
00596 forAllConstIter(Map<label>, meshFaceMap_, iter)
00597 {
00598 label newFaceI = map.reverseFaceMap()[iter.key()];
00599
00600 if (newFaceI >= 0)
00601 {
00602 newMap.insert(newFaceI, iter());
00603 }
00604 }
00605 meshFaceMap_.transfer(newMap);
00606 }
00607 {
00608 Map<label> newMap(meshPointMap_.size());
00609
00610 forAllConstIter(Map<label>, meshPointMap_, iter)
00611 {
00612 label newPointI = map.reversePointMap()[iter.key()];
00613
00614 if (newPointI >= 0)
00615 {
00616 newMap.insert(newPointI, iter());
00617 }
00618 }
00619
00620 meshPointMap_.transfer(newMap);
00621 }
00622 }
00623
00624
00625