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
00030
00031
00032
00033
00034
00035
00036
00037 #include "SundanceHNMesh2D.hpp"
00038
00039 #include "SundanceMeshType.hpp"
00040 #include "SundanceCellJacobianBatch.hpp"
00041 #include "SundanceMaximalCofacetBatch.hpp"
00042 #include "SundanceMeshSource.hpp"
00043 #include "SundanceDebug.hpp"
00044 #include "SundanceOut.hpp"
00045 #include "PlayaMPIContainerComm.hpp"
00046 #include "Teuchos_Time.hpp"
00047 #include "Teuchos_TimeMonitor.hpp"
00048 #include "SundanceObjectWithVerbosity.hpp"
00049 #include "SundanceCollectiveExceptionCheck.hpp"
00050
00051 using namespace Sundance;
00052 using namespace Teuchos;
00053 using namespace std;
00054 using Playa::MPIComm;
00055 using Playa::MPIContainerComm;
00056
00057 int HNMesh2D::offs_Points_x_[4] = {0, 1, 0, 1};
00058
00059 int HNMesh2D::offs_Points_y_[4] = {0, 0, 1, 1};
00060
00061 int HNMesh2D::edge_Points_localIndex[4][2] = { {0,1} , {0,2} , {1,3} , {2,3} };
00062
00063
00064 HNMesh2D::HNMesh2D(int dim, const MPIComm& comm ,
00065 const MeshEntityOrder& order)
00066 : MeshBase(dim, comm , order),_dimension(dim), _comm(comm)
00067 {
00068 setVerb(0);
00069
00070 nrProc_ = MPIComm::world().getNProc();
00071 myRank_ = MPIComm::world().getRank();
00072
00073 points_.resize(0);
00074 nrElem_.resize(3,0);
00075 nrElemOwned_.resize(3,0);
00076
00077 cellsPoints_.resize(0);
00078 cellsEdges_.resize(0);
00079 isCellOut_.resize(0);
00080 edgePoints_.resize(0);
00081 edgeVertex_.resize(0);
00082
00083 edgeMaxCoF_.resize(0);
00084 pointMaxCoF_.resize(0);
00085
00086 elementOwner_.resize(3); elementOwner_[0].resize(0); elementOwner_[1].resize(0); elementOwner_[2].resize(0);
00087
00088 edgeIsProcBonudary_.resize(0);
00089 edgeIsMeshDomainBonudary_.resize(0);
00090
00091 indexInParent_.resize(0);
00092 parentCellLID_.resize(0);
00093 cellLevel_.resize(0);
00094 isCellLeaf_.resize(0);
00095
00096 isPointHanging_.resize(0);
00097 isEdgeHanging_.resize(0);
00098
00099 hangElmStore_.resize(2);
00100 hangElmStore_[0] = Hashtable< int, Array<int> >();
00101 hangElmStore_[1] = Hashtable< int, Array<int> >();
00102 refineCell_.resize(0);
00103
00104 nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0; nrVertexLeafGID_ = 0;
00105 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrEdgeLeafLID_ = 0;
00106 }
00107
00108 int HNMesh2D::numCells(int dim) const {
00109 SUNDANCE_MSG3(verb(),"HNMesh2D::numCells(int dim): dim:" << dim );
00110 switch (dim){
00111 case 0: return nrVertexLeafLID_;
00112 case 1: return nrEdgeLeafLID_;
00113 case 2: return nrCellLeafLID_;
00114 }
00115 return 0;
00116 }
00117
00118 Point HNMesh2D::nodePosition(int i) const {
00119 SUNDANCE_MSG3(verb(),"HNMesh2D::nodePosition(int i) i:"<< i);
00120
00121 return points_[vertexLeafToLIDMapping_[i]];
00122 }
00123
00124 const double* HNMesh2D::nodePositionView(int i) const {
00125 SUNDANCE_MSG3(verb(),"HNMesh2D::nodePositionView(int i) i:" << i);
00126
00127 return &(points_[vertexLeafToLIDMapping_[i]][0]);;
00128 }
00129
00130 void HNMesh2D::getJacobians(int cellDim, const Array<int>& cellLID,
00131 CellJacobianBatch& jBatch) const
00132 {
00133 SUNDANCE_MSG3(verb(),"HNMesh2D::getJacobians cellDim:"<<cellDim<<" _x:"<<_ofs_x<<" _y:"<<_ofs_y);
00134 SUNDANCE_VERB_HIGH("getJacobians()");
00135 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00136 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00137 int nCells = cellLID.size();
00138 int LID;
00139 Point pnt(0.0,0.0);
00140 jBatch.resize(cellLID.size(), spatialDim(), cellDim);
00141 if (cellDim < spatialDim())
00142 {
00143 for (int i=0; i<nCells; i++)
00144 {
00145 double* detJ = jBatch.detJ(i);
00146 switch(cellDim)
00147 {
00148 case 0: *detJ = 1.0;
00149 break;
00150 case 1:
00151 LID = edgeLeafToLIDMapping_[cellLID[i]];
00152 pnt = (points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]]);
00153 *detJ = sqrt(pnt * pnt);
00154 break;
00155 default:
00156 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh2D::getJacobians()");
00157 }
00158 }
00159 }else{
00160
00161 SUNDANCE_VERB_HIGH("cellDim == spatialDim()");
00162 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00163 {
00164 double* J = jBatch.jVals(i);
00165 switch(cellDim)
00166 {
00167 case 2:
00168 LID = cellLeafToLIDMapping_[cellLID[i]];
00169
00170 J[0] = points_[cellsPoints_[LID][1]][0] - points_[cellsPoints_[LID][0]][0];
00171 J[1] = 0.0;
00172 J[2] = 0.0;
00173 J[3] = points_[cellsPoints_[LID][2]][1] - points_[cellsPoints_[LID][0]][1];
00174 SUNDANCE_MSG3(verb() , "HNMesh2D::getJacobians X:" << J[0] << " Y:" << J[3] );
00175 break;
00176 default:
00177 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh2D::getJacobians()");
00178 }
00179 }
00180 }
00181 }
00182
00183 void HNMesh2D::getCellDiameters(int cellDim, const Array<int>& cellLID,
00184 Array<double>& cellDiameters) const {
00185 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00186 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00187 SUNDANCE_VERB_HIGH("getCellDiameters()");
00188 cellDiameters.resize(cellLID.size());
00189 Point pnt(0.0,0.0);
00190 int LID;
00191 if (cellDim < spatialDim())
00192 {
00193 SUNDANCE_MSG3(verb(),"HNMesh2D::getCellDiameters(), cellDim < spatialDim() ");
00194 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00195 {
00196 switch(cellDim)
00197 {
00198 case 0:
00199 cellDiameters[i] = 1.0;
00200 break;
00201 case 1:
00202 LID = edgeLeafToLIDMapping_[cellLID[i]];
00203 pnt = (points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]]);
00204 cellDiameters[i] = sqrt(pnt * pnt);
00205 break;
00206 default:
00207 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh2D::getCellDiameters()");
00208 }
00209 }
00210 }
00211 else
00212 {
00213 SUNDANCE_MSG3(verb(),"HNMesh2D::getCellDiameters(), cellDim == spatialDim() ");
00214 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00215 {
00216 switch(cellDim)
00217 {
00218 case 2:
00219 LID = cellLeafToLIDMapping_[cellLID[i]];
00220 pnt = points_[cellsPoints_[LID][3]] - points_[cellsPoints_[LID][0]];
00221 cellDiameters[i] = sqrt(pnt * pnt);
00222 break;
00223 default:
00224 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00225 "cellDim=" << cellDim << " in HNMesh2D::getCellDiameters()");
00226 }
00227 }
00228 }
00229 }
00230
00231 void HNMesh2D::pushForward(int cellDim, const Array<int>& cellLID,
00232 const Array<Point>& refQuadPts,
00233 Array<Point>& physQuadPts) const {
00234
00235 SUNDANCE_MSG3(verb(),"HNMesh2D::pushForward cellDim:"<<cellDim);
00236 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00237 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00238
00239 int nQuad = refQuadPts.size();
00240 Point pnt( 0.0 , 0.0 );
00241 Point pnt1( 0.0 , 0.0 );
00242
00243 if (physQuadPts.size() > 0) physQuadPts.resize(0);
00244 physQuadPts.reserve(cellLID.size() * refQuadPts.size());
00245 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00246 {
00247 switch(cellDim)
00248 {
00249 case 0:
00250 physQuadPts.append(pnt);
00251 break;
00252 case 1:{
00253 int LID = edgeLeafToLIDMapping_[cellLID[i]];
00254 pnt = points_[edgePoints_[LID][0]];
00255 pnt1 = points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]];
00256 for (int q=0; q<nQuad; q++) {
00257 physQuadPts.append(pnt + (pnt1)*refQuadPts[q][0]);
00258 }
00259 break;}
00260 case 2:{
00261 int LID = cellLeafToLIDMapping_[cellLID[i]];
00262 pnt = points_[cellsPoints_[LID][0]];
00263
00264 pnt1 = points_[cellsPoints_[LID][3]] - points_[cellsPoints_[LID][0]];
00265 for (int q=0; q<nQuad; q++) {
00266 physQuadPts.append( pnt + Point( refQuadPts[q][0] * pnt1[0] , refQuadPts[q][1] * pnt1[1] ) );
00267 }
00268 break;}
00269 default:
00270 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "in HNMesh2D::getJacobians()");
00271 }
00272 }
00273 }
00274
00275 int HNMesh2D::ownerProcID(int cellDim, int cellLID) const {
00276 int ID = -1;
00277 if (cellDim == 0) ID = vertexLeafToLIDMapping_[cellLID];
00278 if (cellDim == 1) ID = edgeLeafToLIDMapping_[cellLID];
00279 if (cellDim == 2) ID = cellLeafToLIDMapping_[cellLID];
00280 SUNDANCE_MSG3(verb() , " HNMesh2D::ownerProcID ,cellDim:" << cellDim << ", cellLID:"
00281 << cellLID <<" ,ID:" << ID << ", ret:"<< elementOwner_[cellDim][ID] );
00282 return elementOwner_[cellDim][ID];
00283 }
00284
00285
00286 int HNMesh2D::numFacets(int cellDim, int cellLID,
00287 int facetDim) const {
00288
00289 if (cellDim==1) {
00290 return 2;
00291 }
00292 else if (cellDim==2) {
00293 return 4;
00294 }
00295 return -1;
00296 }
00297
00298 int HNMesh2D::facetLID(int cellDim, int cellLID,
00299 int facetDim, int facetIndex,
00300 int& facetOrientation) const {
00301
00302 facetOrientation = 1;
00303 SUNDANCE_MSG3(verb(),"HNMesh2D::facetLID cellDim:"<<cellDim<<", cellLID:"<<cellLID<<", facetDim:"<<facetDim<< ", facetIndex:" << facetIndex);
00304 int rnt = -1 , LID=-1 , tmp=-1;
00305 if (facetDim == 0){
00306 if (cellDim == 2 ){
00307 LID = cellLeafToLIDMapping_[cellLID];
00308 rnt = cellsPoints_[LID][facetIndex]; tmp = rnt;
00309 rnt = vertexLIDToLeafMapping_[rnt];
00310 }
00311 else if ((cellDim==1)){
00312 LID = edgeLeafToLIDMapping_[cellLID];
00313 rnt = edgePoints_[LID][facetIndex]; tmp = rnt;
00314 rnt = vertexLIDToLeafMapping_[rnt];
00315 }
00316 } else if (facetDim == 1){
00317 LID = cellLeafToLIDMapping_[cellLID];
00318 rnt = cellsEdges_[LID][facetIndex]; tmp = rnt;
00319 rnt = edgeLIDToLeafMapping_[rnt];
00320 }
00321 SUNDANCE_MSG3(verb()," RET = " << rnt << ", LID:" << LID << ", tmp:" << tmp);
00322 return rnt;
00323 }
00324
00325
00326 void HNMesh2D::getFacetLIDs(int cellDim,
00327 const Array<int>& cellLID,
00328 int facetDim,
00329 Array<int>& facetLID,
00330 Array<int>& facetSign) const {
00331
00332 SUNDANCE_MSG3(verb(),"HNMesh2D::getFacetLIDs() cellDim:"<<cellDim<<" cellLID.size():"<<cellLID.size()<<" facetDim:" <<facetDim);
00333 int LID = 0 , cLID , facetOrientation ;
00334 int ptr = 0;
00335
00336 int nf = numFacets(cellDim, cellLID[0], facetDim);
00337 facetLID.resize(cellLID.size() * nf);
00338 facetSign.resize(cellLID.size() * nf);
00339
00340 for (unsigned int i = 0 ; i < (unsigned int)cellLID.size() ; i++){
00341 cLID = cellLID[i];
00342 for (int f=0; f<nf; f++, ptr++) {
00343
00344 LID = this->facetLID( cellDim, cLID, facetDim, f , facetOrientation);
00345 facetLID[ptr] = LID;
00346 facetSign[ptr] = facetOrientation;
00347 }
00348 }
00349 SUNDANCE_MSG3(verb() ,"HNMesh2D::getFacetLIDs() DONE. ");
00350 }
00351
00352
00353 const int* HNMesh2D::elemZeroFacetView(int cellLID) const {
00354 int LID = cellLeafToLIDMapping_[cellLID];
00355 SUNDANCE_MSG3(verb() , "HNMesh2D::elemZeroFacetView ");
00356 return (const int*)(&cellsPoints_[LID]);
00357 }
00358
00359
00360 int HNMesh2D::numMaxCofacets(int cellDim, int cellLID) const {
00361
00362 SUNDANCE_MSG3(verb() , "HNMesh2D::numMaxCofacets(): cellDim:"<<cellDim<<" cellLID:"<<cellLID );
00363 int rnt = -1;
00364 if (cellDim==0) {
00365 int LID = vertexLeafToLIDMapping_[cellLID];
00366 int sum = 0;
00367 for (int i = 0 ; i < 4 ; i++)
00368 if ( (pointMaxCoF_[LID][i] >= 0) && ( cellLIDToLeafMapping_[pointMaxCoF_[LID][i]] >= 0) )
00369 sum++;
00370
00371 rnt = sum;
00372 }
00373 else if (cellDim==1) {
00374 int LID = edgeLeafToLIDMapping_[cellLID];
00375 if ((edgeMaxCoF_[LID][0] >= 0) && ( cellLIDToLeafMapping_[edgeMaxCoF_[LID][0]] >= 0) &&
00376 (edgeMaxCoF_[LID][1] >= 0) && ( cellLIDToLeafMapping_[edgeMaxCoF_[LID][1]] >= 0))
00377 rnt = 2;
00378 else
00379 rnt = 1;
00380 }
00381 SUNDANCE_MSG3(verb() ," RET = " << rnt );
00382 return rnt;
00383 }
00384
00385
00386 int HNMesh2D::maxCofacetLID(int cellDim, int cellLID,
00387 int cofacetIndex,
00388 int& facetIndex) const {
00389
00390 SUNDANCE_MSG3(verb() ,"HNMesh2D::maxCofacetLID() cellDim:"<<cellDim<<" cellLID:"<<cellLID<<" cofacetIndex:"<<cofacetIndex<< " facetIndex:"
00391 << facetIndex);
00392 int rnt =-1;
00393 if (cellDim==0) {
00394
00395 int actCoFacetIndex = 0;
00396 int LID = vertexLeafToLIDMapping_[cellLID];
00397 for (int ii = 0 ; ii < 4 ; ii++){
00398
00399 if ( (pointMaxCoF_[LID][ii] >= 0) && (cellLIDToLeafMapping_[pointMaxCoF_[LID][ii]] >= 0) ){
00400 if ( actCoFacetIndex < cofacetIndex ){
00401 actCoFacetIndex++;
00402 }else{
00403 facetIndex = ii;
00404 rnt = pointMaxCoF_[LID][ii];
00405 break;
00406 }
00407 }
00408 }
00409 }
00410 else if (cellDim==1) {
00411 int orientation = 0;
00412 int addFakt = 0;
00413 int maxCoFacet = 0;
00414 int LID = edgeLeafToLIDMapping_[cellLID];
00415
00416 if ( edgeVertex_[LID] == 0 ){
00417 orientation = 0; addFakt = 2;
00418 }else{
00419 orientation = 1; addFakt = 1;
00420 }
00421 SUNDANCE_MSG3(verb() ," HNMesh2D::maxCofacetLID() , orientation:" <<orientation<< " addFakt:" << addFakt );
00422
00423 int actCoFacetIndex = 0;
00424 for (int ii = 0 ; ii < 2 ; ii++){
00425
00426 if ( (edgeMaxCoF_[LID][ii] >= 0) && (cellLIDToLeafMapping_[edgeMaxCoF_[LID][ii]] >= 0) ){
00427 if ( actCoFacetIndex < cofacetIndex ){
00428 actCoFacetIndex++;
00429 }else{
00430 facetIndex = 1-ii;
00431 maxCoFacet = edgeMaxCoF_[LID][ii];
00432 break;
00433 }
00434 }
00435 }
00436 SUNDANCE_MSG3(verb() ,"HNMesh2D::maxCofacetLID() , facetIndex:" << facetIndex <<", _edgeMaxCoFacet[0]:" << edgeMaxCoF_[LID][0] <<
00437 "_edgeMaxCoFacet[1]:" << edgeMaxCoF_[LID][1] );
00438
00439 if ( orientation == 0 ){
00440 facetIndex = facetIndex + facetIndex*addFakt;
00441 } else {
00442 facetIndex = facetIndex + addFakt;
00443 }
00444 SUNDANCE_MSG3(verb() ," maxCoFacet = " << maxCoFacet );
00445 rnt = ( maxCoFacet );
00446 }
00447
00448 rnt = cellLIDToLeafMapping_[ rnt ];
00449 SUNDANCE_MSG3(verb() ," RET = " << rnt << ", facetIndex:" << facetIndex);
00450 return rnt;
00451 }
00452
00453 void HNMesh2D::getCofacets(int cellDim, int cellLID,
00454 int cofacetDim, Array<int>& cofacetLIDs) const {
00455
00456 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::getCofacets() not implemented");
00457 }
00458
00459
00460 void HNMesh2D::getMaxCofacetLIDs(const Array<int>& cellLIDs,
00461 MaximalCofacetBatch& cofacets) const {
00462
00463 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::getMaxCofacetLIDs() not implemented");
00464 }
00465
00466
00467 int HNMesh2D::mapGIDToLID(int cellDim, int globalIndex) const {
00468
00469 switch (cellDim){
00470 case 0:{
00471 int ID = vertexLeafToGIDMapping_[globalIndex];
00472 SUNDANCE_MSG3(verb() , " HNMesh2D::mapGIDToLID 0 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< vertexLIDToLeafMapping_[ID]);
00473 return vertexLIDToLeafMapping_[ID];
00474 break;}
00475 case 1:{
00476 int ID = edgeLeafToGIDMapping_[globalIndex];
00477 SUNDANCE_MSG3(verb() , " HNMesh2D::mapGIDToLID 1 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< edgeLIDToLeafMapping_[ID]);
00478 return edgeLIDToLeafMapping_[ID];
00479 break;}
00480 case 2:{
00481 int ID = cellLeafToGIDMapping_[globalIndex];
00482 SUNDANCE_MSG3(verb() , " HNMesh2D::mapGIDToLID 2 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< cellLIDToLeafMapping_[ID]);
00483 return cellLIDToLeafMapping_[ID];
00484 break;}
00485 }
00486 return -1;
00487 }
00488
00489
00490 bool HNMesh2D::hasGID(int cellDim, int globalIndex) const {
00491
00492
00493 return true;
00494 }
00495
00496
00497 int HNMesh2D::mapLIDToGID(int cellDim, int localIndex) const {
00498
00499 switch (cellDim){
00500 case 0:{
00501 int ID = vertexLeafToLIDMapping_[localIndex];
00502 SUNDANCE_MSG3(verb() , " HNMesh2D::mapLIDToGID 0 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< vertexGIDToLeafMapping_[ID]);
00503 return vertexGIDToLeafMapping_[ID];
00504 break;}
00505 case 1:{
00506 int ID = edgeLeafToLIDMapping_[localIndex];
00507 SUNDANCE_MSG3(verb() , " HNMesh2D::mapLIDToGID 1 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< edgeGIDToLeafMapping_[ID]);
00508 return edgeGIDToLeafMapping_[ID];
00509 break;}
00510 case 2:{
00511 int ID = cellLeafToLIDMapping_[localIndex];
00512 SUNDANCE_MSG3(verb() , " HNMesh2D::mapLIDToGID 2 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< cellGIDToLeafMapping_[ID]);
00513 return cellGIDToLeafMapping_[ID];
00514 break;}
00515 }
00516 return -1;
00517 }
00518
00519
00520 CellType HNMesh2D::cellType(int cellDim) const {
00521 switch(cellDim)
00522 {
00523 case 0: return PointCell;
00524 case 1: return LineCell;
00525 case 2: return QuadCell;
00526 default:
00527 return NullCell;
00528 }
00529 }
00530
00531
00532 int HNMesh2D::label(int cellDim, int cellLID) const {
00533
00534 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::label() not implemented yet");
00535 return 0;
00536 }
00537
00538
00539 void HNMesh2D::getLabels(int cellDim, const Array<int>& cellLID,
00540 Array<int>& labels) const {
00541
00542 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::getLabels() not implemented yet");
00543 }
00544
00545 Set<int> HNMesh2D::getAllLabelsForDimension(int cellDim) const {
00546 Set<int> rtn;
00547
00548 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::getAllLabelsForDimension() not implemented yet");
00549 return rtn;
00550 }
00551
00552 void HNMesh2D::getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const {
00553
00554 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::getLIDsForLabel() not implemented yet");
00555 }
00556
00557 void HNMesh2D::setLabel(int cellDim, int cellLID, int label) {
00558
00559 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh2D::setLabel() not implemented yet");
00560 }
00561
00562
00563 void HNMesh2D::assignIntermediateCellGIDs(int cellDim) {
00564
00565
00566 }
00567
00568
00569 bool HNMesh2D::hasIntermediateGIDs(int dim) const {
00570
00571
00572 return true;
00573 }
00574
00575
00576
00577 bool HNMesh2D::isElementHangingNode(int cellDim , int cellLID) const {
00578 SUNDANCE_MSG3(verb() ,"HNMesh2D::isElementHangingNode cellDim:"<<cellDim<<" LID:"<< cellLID);
00579 if (cellDim==0) {
00580 int LID = vertexLeafToLIDMapping_[cellLID];
00581 return (isPointHanging_[LID]);
00582 }
00583 else if (cellDim==1) {
00584 int LID = edgeLeafToLIDMapping_[cellLID];
00585 return (isEdgeHanging_[LID]);
00586 } else {
00587 return false;
00588 }
00589 return false;
00590 }
00591
00592 int HNMesh2D::indexInParent(int maxCellLID) const {
00593
00594 int mapMyChilderIndexToStandardDoFMAP[9] = {0,5,6,1,4,7,2,3,8};
00595 int LID = cellLeafToLIDMapping_[maxCellLID];
00596 int indexInPar = indexInParent_[LID];
00597
00598 return mapMyChilderIndexToStandardDoFMAP[indexInPar];
00599 }
00600
00601 void HNMesh2D::returnParentFacets( int childCellLID , int dimFacets ,
00602 Array<int> &facetsLIDs , int &parentCellLIDs ) const {
00603 int LID = cellLeafToLIDMapping_[childCellLID];
00604 parentCellLIDs = parentCellLID_[LID];
00605 facetsLIDs.resize(4);
00606 SUNDANCE_MSG3(verb() , "HNMesh2D::returnParentFacets childCellLID:"<<childCellLID<<" dimFacets:"<<dimFacets<<
00607 " parentCellLIDs:"<< parentCellLIDs);
00608
00609 facetsLIDs[0] = facetLID_tree( 2 , parentCellLIDs , dimFacets , 0 );
00610 facetsLIDs[1] = facetLID_tree( 2 , parentCellLIDs , dimFacets , 1 );
00611 facetsLIDs[2] = facetLID_tree( 2 , parentCellLIDs , dimFacets , 2 );
00612 facetsLIDs[3] = facetLID_tree( 2 , parentCellLIDs , dimFacets , 3 );
00613
00614 parentCellLIDs = cellLIDToLeafMapping_[parentCellLIDs];
00615 }
00616
00617
00618 int HNMesh2D::facetLID_tree(int cellDim, int cellLID,
00619 int facetDim, int facetIndex) const{
00620 int rnt = -1;
00621 if (facetDim == 0){
00622 rnt = cellsPoints_[cellLID][facetIndex];
00623 rnt = vertexLIDToLeafMapping_[rnt];
00624
00625 } else if (facetDim == 1){
00626 rnt = cellsEdges_[cellLID][facetIndex];
00627 rnt = edgeLIDToLeafMapping_[rnt];
00628
00629 }
00630 SUNDANCE_MSG3(verb() , "HNMesh2D::facetLID_tree cellDim:"<<cellDim<<", cellLID:"<<cellLID<<", facetDim:"<<facetDim<<
00631 ", facetIndex:"<<facetIndex<<" RET = " << rnt );
00632 return rnt;
00633 }
00634
00635
00636
00637
00638 void HNMesh2D::addVertex(int vertexLID , int ownerProc , bool isHanging ,
00639 double coordx , double coordy , const Array<int> &maxCoF){
00640
00641
00642 if (points_.size() <= vertexLID){
00643 TEUCHOS_TEST_FOR_EXCEPTION(vertexLID != nrElem_[0] , std::logic_error ,"HNMesh2D::addVertex " <<
00644 " vertexLID:" << vertexLID << " nrElem_[0]:" << nrElem_[0] );
00645 Point pt(coordx,coordy);
00646 points_.append( pt );
00647 pointMaxCoF_.append( maxCoF );
00648 isPointHanging_.append( isHanging );
00649 elementOwner_[0].append( ownerProc );
00650 SUNDANCE_MSG3(verb() , "HNMesh2D::addVertex: " << nrElem_[0] << " P:" << pt << " , maxCoF:" << maxCoF);
00651 SUNDANCE_MSG3(verb() , " ownerProc:" << ownerProc << " , isHanging:" << isHanging);
00652 nrElem_[0]++;
00653 }
00654 }
00655
00656
00657 void HNMesh2D::addEdge(int edgeLID , int ownerProc , bool isHanging , int edgeVertex ,
00658 bool isProcBoundary , bool isMeshBoundary ,
00659 const Array<int> &vertexLIDs , const Array<int> &maxCoF){
00660
00661 if (edgePoints_.size() <= edgeLID ){
00662 TEUCHOS_TEST_FOR_EXCEPTION(edgeLID != nrElem_[1], std::logic_error, "HNMesh2D::addEdge edgeLID != nrElem_[1]");
00663 edgePoints_.append( vertexLIDs );
00664 edgeVertex_.append( edgeVertex );
00665 edgeMaxCoF_.append( maxCoF );
00666 edgeIsProcBonudary_.append( isProcBoundary );
00667 edgeIsMeshDomainBonudary_.append( isMeshBoundary );
00668 isEdgeHanging_.append(isHanging);
00669 elementOwner_[1].append( ownerProc );
00670 SUNDANCE_MSG3(verb() , "HNMesh2D::addEdge: " << nrElem_[1] << " vertexLIDs:" << vertexLIDs << " , maxCoF:" << maxCoF );
00671 SUNDANCE_MSG3(verb() , " ownerProc:" << ownerProc << ", isHanging:" << isHanging << ", edgeVertex:" << edgeVertex <<
00672 ", isProcBoundary:" << isProcBoundary << ", isMeshBoundary:" << isMeshBoundary);
00673 nrElem_[1]++;
00674 }
00675 }
00676
00677
00678 void HNMesh2D::addCell(int cellLID , int ownerProc ,
00679 int indexInParent , int parentCellLID , int level ,
00680 const Array<int> &edgeLIDs , const Array<int> &vertexLIDs){
00681
00682 if (cellsPoints_.size() <= cellLID ) {
00683 TEUCHOS_TEST_FOR_EXCEPTION(cellLID != nrElem_[2], std::logic_error, "HNMesh2D::cellLID cellLID != nrElem_[2]");
00684 cellsPoints_.append( vertexLIDs );
00685 cellsEdges_.append( edgeLIDs );
00686 indexInParent_.append( indexInParent );
00687 parentCellLID_.append( parentCellLID );
00688 cellLevel_.append( level );
00689 isCellLeaf_.append( true );
00690 refineCell_.append( 0 );
00691 cellsChildren_.append( tuple(1) );
00692 elementOwner_[2].append( ownerProc );
00693 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell: " << nrElem_[2] << " vertexLIDs:" << vertexLIDs << " edgeLIDs:" << edgeLIDs);
00694 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell p0:" << points_[vertexLIDs[0]] );
00695 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell p1:" << points_[vertexLIDs[1]] );
00696 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell p2:" << points_[vertexLIDs[2]] );
00697 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell p3:" << points_[vertexLIDs[3]] );
00698
00699 isCellOut_.append( !( meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[0]]) ||
00700 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[1]]) ||
00701 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[2]]) ||
00702 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[3]]) ) );
00703 SUNDANCE_MSG3(verb() , "HNMesh2D::addCell IN DOMAIN:" << isCellOut_[nrElem_[2]] );
00704 nrElem_[2]++;
00705 }
00706 }
00707
00708
00709
00710
00711
00712 void HNMesh2D::createMesh(
00713 double position_x,
00714 double position_y,
00715 double offset_x,
00716 double offset_y,
00717 int resolution_x,
00718 int resolution_y,
00719 const RefinementClass& refineClass ,
00720 const MeshDomainDef& meshDomain
00721 ){
00722
00723 setVerb(0);
00724
00725
00726 _pos_x = position_x; _pos_y = position_y;
00727 _ofs_x = offset_x; _ofs_y = offset_y;
00728 _res_x = resolution_x; _res_y = resolution_y;
00729 refineClass_ = refineClass;
00730 meshDomain_ = meshDomain;
00731
00732
00733 createCoarseMesh();
00734
00735
00736 bool doRefinement = true;
00737 while (doRefinement){
00738 doRefinement = oneRefinementIteration();
00739 }
00740
00741
00742
00743 createLeafNumbering_sophisticated();
00744
00745 }
00746
00747 void HNMesh2D::createCoarseMesh(){
00748
00749
00750
00751
00752 int nrCoarseCell = _res_x*_res_y;
00753 int nrCoarsePoints = (_res_x+1)*(_res_y+1);
00754 int nrCoarseEdge = (_res_x+1)*_res_y + _res_x*(_res_y+1);
00755 Array<int> coarseProcessCell( nrCoarseCell , -1 );
00756 Array<int> coarseCellOwner( nrCoarseCell , -1 );
00757 Array<int> coarsePointOwner( nrCoarsePoints , -1 );
00758 Array<int> coarseEdgeOwner( nrCoarseEdge , -1 );
00759 Array<int> coarseCellLID( nrCoarseCell , -1 );
00760 Array<int> coarsePointLID( nrCoarsePoints , -1 );
00761 Array<int> coarseEdgeLID( nrCoarseEdge , -1 );
00762 Array<int> coarseCellsLoad( nrCoarseCell , 0 );
00763 int totalLoad = 0;
00764
00765 SUNDANCE_MSG3(verb() , "HNMesh2D::createMesh nrCoarseCell:" << nrCoarseCell << " nrCoarsePoints:" << nrCoarsePoints
00766 << " nrCoarseEdge:" << nrCoarseEdge << " nrProc_:" << nrProc_ << " myRank_:" << myRank_);
00767 TEUCHOS_TEST_FOR_EXCEPTION( nrCoarseCell < nrProc_ , std::logic_error," HNMesh2D::createMesh nrCoarseCell < nrProc_ ");
00768
00769
00770
00771 double h[2];
00772 h[0] = _ofs_x/(double)_res_x; h[1] = _ofs_y/(double)_res_y;
00773 Point pos(h[0],h[1]);
00774 Point res(h[0],h[1]);
00775
00776 for (int i=0; i < nrCoarseCell; i++){
00777
00778 pos[0] = _pos_x + (double)(i / _res_y)*h[0] + 0.5*h[0];
00779 pos[1] = _pos_y + (double)(i % _res_y)*h[1] + 0.5*h[1];
00780
00781 coarseCellsLoad[i] = refineClass_.estimateRefinementLevel( pos , res );
00782 totalLoad += coarseCellsLoad[i];
00783 }
00784
00785 double loadPerProc = (double)totalLoad / (double)nrProc_;
00786 int actualProc=0;
00787 totalLoad = 0;
00788 Array<int> ind(2);
00789
00790 int vertexOffs[4] = {0 , _res_y + 1 , 1 , _res_y + 2};
00791 int edgeOffs[4] = {_res_y , 0 , 2*_res_y + 1 , _res_y + 1};
00792
00793
00794 SUNDANCE_MSG3(verb() , "Processor asign, loadPerProc:" << loadPerProc );
00795 double diff_load = 0.0;
00796 for (int i=0; i < nrCoarseCell; i++){
00797 ind[0] = (i / _res_y);
00798 ind[1] = (i % _res_y);
00799 int vertexInd = (_res_y+1)*ind[0] + ind[1];
00800 int edgeInd = (2*_res_y+1)*ind[0] + ind[1];
00801
00802
00803
00804 for (int jj = 0 ; jj < 4 ; jj++){
00805 if (coarsePointOwner[vertexInd+vertexOffs[jj]] < 0){
00806 coarsePointOwner[vertexInd+vertexOffs[jj]] = actualProc;
00807 SUNDANCE_MSG3(verb() , "Vertex CPU assign " << vertexInd+vertexOffs[jj] << " ->" << actualProc );
00808 }
00809 }
00810
00811 for (int jj = 0 ; jj < 4 ; jj++){
00812 if (coarseEdgeOwner[edgeInd+edgeOffs[jj]] < 0){
00813 coarseEdgeOwner[edgeInd+edgeOffs[jj]] = actualProc;
00814
00815 }
00816 }
00817
00818 coarseCellOwner[i] = actualProc;
00819 totalLoad += coarseCellsLoad[i];
00820 SUNDANCE_MSG3(verb() , "Cell CPU assign " << i << " ->" << actualProc <<
00821 ", totalLoad:" << totalLoad << " loadPerProc:" << loadPerProc);
00822
00823 if (((double)totalLoad >= (loadPerProc - 1e-8 - diff_load)) && ( actualProc < nrProc_-1 )){
00824 SUNDANCE_MSG3(verb() , "Increase CPU , totalLoad:" << totalLoad << " loadPerProc:" << loadPerProc );
00825
00826 diff_load = totalLoad - loadPerProc;
00827 actualProc = actualProc + 1;
00828 totalLoad = 0;
00829 }
00830 }
00831
00832
00833 for (int i=0; i < nrCoarseCell; i++){
00834 ind[0] = (i / _res_y);
00835 ind[1] = (i % _res_y);
00836 coarseProcessCell[i] = 1;
00837 }
00838
00839
00840 SUNDANCE_MSG3(verb() , " Process Cells:" << nrCoarseCell );
00841 for (int i=0; i < nrCoarseCell; i++){
00842 ind[0] = (i / _res_y);
00843 ind[1] = (i % _res_y);
00844 pos[0] = _pos_x + (double)(ind[0])*h[0];
00845 pos[1] = _pos_y + (double)(ind[1])*h[1];
00846 SUNDANCE_MSG3(verb() , "PCell ID:" << i << " coarseProcessCell[i]:" <<coarseProcessCell[i] << " pos:" << pos );
00847 SUNDANCE_MSG3(verb() , "PCell, actual index" << ind );
00848
00849 if (coarseProcessCell[i] > 0) {
00850 int vertexInd = (_res_y+1)*ind[0] + ind[1];
00851 int edgeInd = (2*_res_y+1)*ind[0] + ind[1];
00852 int cellLID = coarseCellLID[i];
00853 Array<int> vertexLID(4,-1) , vertexMaxCoF(4,-1) , vLID(4,-1);;
00854 Array<int> edgeLID(4,-1) , edgeVert(2,-1) , edgeMaxCoef(2,-1);
00855
00856
00857
00858 if (coarseCellLID[i] < 0 ){
00859 coarseCellLID[i] = nrElem_[2];
00860 cellLID = coarseCellLID[i];
00861 }
00862
00863 for (int jj = 0 ; jj < 4 ; jj++){
00864 vLID[jj] = coarsePointLID[vertexInd+vertexOffs[jj]];
00865
00866
00867
00868 if (coarsePointLID[vertexInd+vertexOffs[jj]] < 0){
00869 coarsePointLID[vertexInd+vertexOffs[jj]] = nrElem_[0];
00870 }
00871 vLID[jj] = coarsePointLID[vertexInd+vertexOffs[jj]];
00872
00873
00874
00875 addVertex( vLID[jj] , coarsePointOwner[vertexInd+vertexOffs[jj]] , false ,
00876 pos[0] + ((double)offs_Points_x_[jj])*h[0] , pos[1] + ((double)offs_Points_y_[jj])*h[1] ,
00877 vertexMaxCoF );
00878 }
00879
00880 for (int jj = 0 ; jj < 4 ; jj++){
00881 edgeLID[jj] = coarseEdgeLID[edgeInd+edgeOffs[jj]];
00882 if (coarseEdgeLID[edgeInd+edgeOffs[jj]] < 0){
00883 coarseEdgeLID[edgeInd+edgeOffs[jj]] = nrElem_[1];
00884 }
00885 edgeLID[jj] = coarseEdgeLID[edgeInd+edgeOffs[jj]];
00886 edgeVert[0] = vLID[edge_Points_localIndex[jj][0]];
00887 edgeVert[1] = vLID[edge_Points_localIndex[jj][1]];
00888
00889 addEdge( edgeLID[jj] , coarseEdgeOwner[edgeInd+edgeOffs[jj]] , false , (jj==0 || jj==3) ? 0 : 1 ,
00890 false , false ,
00891 edgeVert , edgeMaxCoef );
00892 }
00893
00894 addCell( cellLID , coarseCellOwner[i] , 0 , cellLID , 0 , edgeLID , vLID);
00895 }
00896 }
00897
00898
00899 SUNDANCE_MSG3(verb() , " Process maxCofacets:" << nrCoarseCell );
00900 for (int i=0; i < nrCoarseCell; i++){
00901
00902 if (coarseProcessCell[i] > 0){
00903 Array<int> vLID(4,-1) , eLID(4,-1) , maxVertexCoF(4,-1);
00904 ind[0] = (i / _res_y);
00905 ind[1] = (i % _res_y);
00906 int vertexInd = (_res_y+1)*ind[0] + ind[1];
00907 int edgeInd = (2*_res_y+1)*ind[0] + ind[1];
00908 int cellLID = coarseCellLID[i];
00909 SUNDANCE_MSG3(verb() , "MCell cellLID:" << cellLID << " coarseProcessCell[i]:" <<coarseProcessCell[i] );
00910
00911
00912
00913 for (int jj = 0 ; jj < 4 ; jj++){
00914 vLID[jj] = coarsePointLID[vertexInd+vertexOffs[jj]];
00915
00916 pointMaxCoF_[vLID[jj]][jj] = cellLID;
00917 SUNDANCE_MSG3(verb() , "Vertex MaxCoFacet vLID[jj]:" << vLID[jj] << " jj:" << jj << " cellLID:" << cellLID );
00918 }
00919
00920 for (int jj = 0 ; jj < 4 ; jj++){
00921 eLID[jj] = coarseEdgeLID[edgeInd+edgeOffs[jj]];
00922
00923 edgeMaxCoF_[eLID[jj]][(3-jj)/2] = cellLID;
00924 SUNDANCE_MSG3(verb() , "Edge MaxCoFacet eLID[jj]:" << eLID[jj] << " (3-jj)/2:" << (3-jj)/2 << " cellLID:" << cellLID );
00925
00926 int orientation = ( (jj==0) || (jj==3)) ? 1 : 0;
00927
00928 if (orientation == 0){
00929
00930 if ( (ind[0] == 0) && (jj==1)) edgeIsMeshDomainBonudary_[eLID[jj]] = true;
00931 if ( (ind[0] == _res_x - 1 ) && (jj==2)) edgeIsMeshDomainBonudary_[eLID[jj]] = true;
00932
00933 if ( edgeIsMeshDomainBonudary_[eLID[jj]] != true){
00934
00935 edgeIsProcBonudary_[eLID[jj]] = false;
00936 }
00937 }else{
00938
00939 if ( (ind[1] == 0) && (jj==0)) { edgeIsMeshDomainBonudary_[eLID[jj]] = true; }
00940 if ( (ind[1] == _res_y - 1 ) && (jj==3)) { edgeIsMeshDomainBonudary_[eLID[jj]] = true; }
00941
00942 if ( edgeIsMeshDomainBonudary_[eLID[jj]] != true){
00943
00944 edgeIsProcBonudary_[eLID[jj]] = false;
00945 }
00946 SUNDANCE_MSG3(verb() , "Neighb. H Cell DONE ");
00947 }
00948 }
00949 }
00950 }
00951
00952 }
00953
00954
00955 bool HNMesh2D::oneRefinementIteration(){
00956
00957 Array<int> ind(2);
00958 int nrActualCell = nrElem_[2];
00959 bool rtn = false;
00960 SUNDANCE_MSG3(verb() , " HNMesh2D::oneRefinementIteration, start one refinement iteration cycle ");
00961
00962 for (int i=0 ; i < nrActualCell ; i++){
00963
00964 ind[0] = (i / _res_y);
00965 ind[1] = (i % _res_y);
00966
00967
00968 SUNDANCE_MSG3(verb() , " Test cell " << i << ", elementOwner_[2][i]:" << elementOwner_[2][i] <<
00969 ", isCellLeaf_[i]:" << isCellLeaf_[i] << ", out:" << (!isCellOut_[i]));
00970
00971 if ( (isCellLeaf_[i] == true) )
00972
00973
00974
00975 {
00976
00977 Array<int>& cellsEdges = cellsEdges_[i];
00978 bool doRefined = true;
00979 bool refFunction = false;
00980 for (int jj = 0 ; jj < cellsEdges.size() ; jj++){
00981 SUNDANCE_MSG3(verb() , " eLID: " << cellsEdges[jj] << ", isHanging:" << isEdgeHanging_[cellsEdges[jj]]);
00982 doRefined = ( doRefined && ((isEdgeHanging_[cellsEdges[jj]]) == false) );
00983 }
00984
00985 SUNDANCE_MSG3(verb() , " is possible to refine cell nr: " << i << ", doRefined:" << doRefined);
00986
00987 Array<int>& cellsVertex = cellsPoints_[i];
00988 Point h = points_[cellsVertex[3]] - points_[cellsVertex[0]];
00989 Point p2 = points_[cellsVertex[0]] + 0.5*h;
00990 refFunction = ((refineCell_[i] == 1) || refineClass_.refine( cellLevel_[i] , p2 , h ));
00991
00992
00993 SUNDANCE_MSG3( verb() , " execute refinement on cell nr: " << i << ", refFunction:" << refFunction);
00994 if (doRefined && refFunction)
00995 {
00996
00997 refineCell(i);
00998 rtn = true;
00999 refineCell_[i] = 0;
01000 }
01001 else
01002 {
01003
01004
01005 if (refFunction){
01006
01007 for (int jj = 0 ; jj < cellsEdges.size() ; jj++){
01008 if (isEdgeHanging_[cellsEdges[jj]]){
01009
01010
01011 int cLevel = cellLevel_[i];
01012 int parentID = parentCellLID_[i];
01013 int parentCEdge = cellsEdges_[parentID][jj];
01014 int refCell = -1;
01015
01016 if ( (jj == 0) || ( jj == 1) )
01017 refCell = edgeMaxCoF_[parentCEdge][0];
01018 else
01019 refCell = edgeMaxCoF_[parentCEdge][1];
01020
01021
01022
01023
01024 if ( (refCell >= 0) && (cellLevel_[refCell] < cLevel) && (isCellLeaf_[refCell])){
01025
01026 refineCell_[refCell] = 1;
01027 rtn = true;
01028
01029
01030 }
01031 }
01032 }
01033 }
01034 }
01035 }
01036 }
01037
01038
01039
01040
01041 SUNDANCE_MSG3(verb() , " HNMesh2D::oneRefinementIteration DONE with one refinement iteration");
01042
01043 return rtn;
01044 }
01045
01046
01047 void HNMesh2D::refineCell(int cellLID){
01048
01049 Array<int>& cellsEdges = cellsEdges_[cellLID];
01050 Array<int>& cellsVertex = cellsPoints_[cellLID];
01051 int cellOwner = elementOwner_[2][cellLID];
01052 Point p = points_[cellsVertex[0]];
01053 Point h = points_[cellsVertex[3]] - points_[cellsVertex[0]];
01054
01055 double vertex_X[16] = { 0.0 , 0.0 , 0.0 , 0.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
01056 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 1.0 , 1.0 , 1.0 , 1.0 };
01057 double vertex_Y[16] = { 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
01058 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 };
01059
01060 int hanginEdgeI[4][3] = { {3,10,17} , {0,1,2} , {21,22,23} , {6,13,20} };
01061
01062 int hanginVertexI[4][2] = { {4,8} , {1,2} , {13,14} , {7,11} };
01063
01064 int edgeVertex[24] = { 1,1,1 , 0,0,0,0 , 1,1,1 , 0,0,0,0 , 1,1,1 , 0,0,0,0 , 1,1,1};
01065
01066 int VertexIndex[4] = { 0 , 1 , 1 , 0 };
01067
01068 int VertexI[4] = { 0 , 0 , 1 , 2 };
01069
01070
01071 int vertexOwner[16] = { elementOwner_[0][cellsVertex[0]] , elementOwner_[1][cellsEdges[1]] ,
01072 elementOwner_[1][cellsEdges[1]] , elementOwner_[0][cellsVertex[2]] ,
01073 elementOwner_[1][cellsEdges[0]] , cellOwner , cellOwner , elementOwner_[1][cellsEdges[3]] ,
01074 elementOwner_[1][cellsEdges[0]] , cellOwner , cellOwner , elementOwner_[1][cellsEdges[3]] ,
01075 elementOwner_[0][cellsVertex[1]] , elementOwner_[1][cellsEdges[2]] ,
01076 elementOwner_[1][cellsEdges[2]] , elementOwner_[0][cellsVertex[3]] };
01077
01078 bool vertexIsHanging[16] = { false , !edgeIsMeshDomainBonudary_[cellsEdges[1]] , !edgeIsMeshDomainBonudary_[cellsEdges[1]], false ,
01079 !edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , !edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01080 !edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , !edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01081 false , !edgeIsMeshDomainBonudary_[cellsEdges[2]] , !edgeIsMeshDomainBonudary_[cellsEdges[2]], false };
01082
01083 int edgeOwner[24] = { elementOwner_[1][cellsEdges[1]] , elementOwner_[1][cellsEdges[1]] , elementOwner_[1][cellsEdges[1]] ,
01084 elementOwner_[1][cellsEdges[0]] , cellOwner , cellOwner , elementOwner_[1][cellsEdges[3]] ,
01085 cellOwner , cellOwner , cellOwner ,
01086 elementOwner_[1][cellsEdges[0]] , cellOwner , cellOwner , elementOwner_[1][cellsEdges[3]] ,
01087 cellOwner , cellOwner , cellOwner ,
01088 elementOwner_[1][cellsEdges[0]] , cellOwner , cellOwner , elementOwner_[1][cellsEdges[3]] ,
01089 elementOwner_[1][cellsEdges[2]] , elementOwner_[1][cellsEdges[2]] , elementOwner_[1][cellsEdges[2]] };
01090
01091 bool edgeIsHanging[24] = { !edgeIsMeshDomainBonudary_[cellsEdges[1]] , !edgeIsMeshDomainBonudary_[cellsEdges[1]] , !edgeIsMeshDomainBonudary_[cellsEdges[1]] ,
01092 !edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , !edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01093 false , false , false ,
01094 !edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , !edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01095 false , false , false ,
01096 !edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , !edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01097 !edgeIsMeshDomainBonudary_[cellsEdges[2]] , !edgeIsMeshDomainBonudary_[cellsEdges[2]] , !edgeIsMeshDomainBonudary_[cellsEdges[2]] };
01098
01099 int edgeVertexes[24][2] = { {0,1} , {1,2} , {2,3} , {0,4} , {1,5} , {2,6} , {3,7} , {4,5} , {5,6} , {6,7} ,
01100 {4,8} , {5,9} , {6,10} , {7,11} , {8,9} , {9,10} , {10,11} ,
01101 {8,12} , {9,13} , {10,14} , {11,15} , {12,13} , {13,14} , {14,15}};
01102
01103 bool edgePBnd[24] = { edgeIsProcBonudary_[cellsEdges[1]] , edgeIsProcBonudary_[cellsEdges[1]] , edgeIsProcBonudary_[cellsEdges[1]] ,
01104 edgeIsProcBonudary_[cellsEdges[0]] , false , false , edgeIsProcBonudary_[cellsEdges[3]] ,
01105 false , false , false ,
01106 edgeIsProcBonudary_[cellsEdges[0]] , false , false , edgeIsProcBonudary_[cellsEdges[3]] ,
01107 false , false , false ,
01108 edgeIsProcBonudary_[cellsEdges[0]] , false , false , edgeIsProcBonudary_[cellsEdges[3]] ,
01109 edgeIsProcBonudary_[cellsEdges[2]] , edgeIsProcBonudary_[cellsEdges[2]] , edgeIsProcBonudary_[cellsEdges[2]] };
01110
01111 bool edgeMBnd[24] = { edgeIsMeshDomainBonudary_[cellsEdges[1]] , edgeIsMeshDomainBonudary_[cellsEdges[1]] , edgeIsMeshDomainBonudary_[cellsEdges[1]] ,
01112 edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01113 false , false , false ,
01114 edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01115 false , false , false ,
01116 edgeIsMeshDomainBonudary_[cellsEdges[0]] , false , false , edgeIsMeshDomainBonudary_[cellsEdges[3]] ,
01117 edgeIsMeshDomainBonudary_[cellsEdges[2]] , edgeIsMeshDomainBonudary_[cellsEdges[2]] , edgeIsMeshDomainBonudary_[cellsEdges[2]] };
01118
01119 int refinedCellsVertexes[9][4] = { {0,4,1,5} , {1,5,2,6} , {2,6,3,7} ,
01120 {4,8,5,9} , {5,9,6,10} , {6,10,7,11} ,
01121 {8,12,9,13} , {9,13,10,14} , {10,14,11,15} };
01122
01123 int refinedCellsEdges[9][4] = { {3,0,7,4} , {4,1,8,5} , {5,2,9,6} ,
01124 {10,7,14,11} , {11,8,15,12} , {12,9,16,13} ,
01125 {17,14,21,18} , {18,15,22,19} , {19,16,23,20} };
01126
01127 int hVertexIndex[16] = { -1 , 0 , 1 , -1 , 0 , -1 , -1 , 0 , 1 , -1 , -1 , 1 , -1 , 0 , 1 , -1 };
01128 int hEdgeIndex[24] = { 0 , 1 , 2 , 0 , -1 , -1 , 0 , -1 , -1 , -1 , 1 , -1 , -1 , 1 , -1 , -1 , -1 , 2 , -1 , -1 , 2 , 0 , 1 , 2 };
01129 int hVertexVertexIndex[16] = { -1 , 1 , 1 , -1 , 0 , -1 , -1 , 3 , 0 , -1 , -1 , 3 , -1 , 2 , 2 , -1 };
01130 int hEdgeVertexIndex[24] = { 1 , 1 , 1 , 0 , -1 , -1 , 3 , -1 , -1 , -1 , 0 , -1 , -1 , 3 , -1 , -1 , -1 , 0 , -1 , -1 , 3 , 2 , 2 , 2 };
01131
01132
01133 int vertexMaxCoFacet[16][4] = { {0,-1,-1,-1} , {1,-1,0,-1} , {2 ,-1,1,-1} , {-1,-1,2,-1} ,
01134 {3,0,-1,-1} , {4,1,3,0} , {5 ,2,4,1} , {-1,-1,5,2} ,
01135 {6,3,-1,-1} , {7,4,6,3} , {8 ,5,7,4} , {-1,-1,8,5} ,
01136 {-1,6,-1,-1} , {-1,7,-1,6} , {-1,8,-1,7} , {-1,-1,-1,8} };
01137
01138 int edgeMaxCoFacet[24][2] = { {-1,0} , {-1,1} , {-1,2} ,
01139 {-1,0} , {0,1} , {1,2} , {2,-1} ,
01140 {0,3} , {1,4} , {2,5} ,
01141 {-1,3} , {3,4} , {4,5} , {5,-1} ,
01142 {3,6} , {4,7} , {5,8} ,
01143 {-1,6} , {6,7} , {7,8} , {8,-1} ,
01144 {6,-1} , {7,-1} , {8,-1} };
01145
01146 int indexEdgeParentMaxCoFacet[24] = { 1 , 1, 1 , 0 , -1 , -1 , 3 , -1 , -1 , - 1, 0 , -1 , -1 , 3 , -1 , -1 , -1 , 0 ,-1,-1,3, 2,2,2 };
01147
01148 int offsEdgeParentMaxCoFacet[24] = { 0 , 0, 0 , 0 , -1 , -1 , 1 , -1 , -1 , - 1, 0 , -1 , -1 , 1 , -1 , -1 , -1 , 0 ,-1,-1,1, 1,1,1 };
01149
01150
01151 Array< Array<int> > hangingInfo(4);
01152 hangingInfo[0].resize(0); hangingInfo[1].resize(0); hangingInfo[2].resize(0); hangingInfo[3].resize(0);
01153
01154
01155 Array<int> vertexLIDs(16,-1);
01156 Array<int> edgeLIDs(24,-1);
01157 Array<int> cellLIDs(9,-1);
01158
01159 SUNDANCE_MSG3(verb() , " ========== HNMesh2D::refineCell, cellLID:" << cellLID << " ==========");
01160
01161
01162 vertexLIDs[0] = cellsVertex[0]; vertexLIDs[3] = cellsVertex[2];
01163 vertexLIDs[12] = cellsVertex[1]; vertexLIDs[15] = cellsVertex[3];
01164
01165
01166 for (int v = 0 ; v < 4 ; v++){
01167
01168
01169
01170 if (hangElmStore_[VertexIndex[v]].containsKey(cellsVertex[VertexI[v]])){
01171 const Array<int>& hnInfo = hangElmStore_[VertexIndex[v]].get(cellsVertex[VertexI[v]]);
01172
01173
01174
01175 vertexLIDs[hanginVertexI[v][0]] = hnInfo[0]; isPointHanging_[hnInfo[0]] = false;
01176 vertexLIDs[hanginVertexI[v][1]] = hnInfo[1]; isPointHanging_[hnInfo[1]] = false;
01177 edgeLIDs[hanginEdgeI[v][0]] = hnInfo[2];
01178 isEdgeHanging_[hnInfo[2]] = false;
01179 edgeLIDs[hanginEdgeI[v][1]] = hnInfo[3];
01180 isEdgeHanging_[hnInfo[3]] = false;
01181 edgeLIDs[hanginEdgeI[v][2]] = hnInfo[4]; isEdgeHanging_[hnInfo[4]] = false;
01182
01183 hangElmStore_[VertexIndex[v]].remove(cellsVertex[VertexI[v]]);
01184 }
01185 }
01186 SUNDANCE_MSG3(verb() , " refineCell bef. vertexLIDs:" << vertexLIDs );
01187 SUNDANCE_MSG3(verb() , " refineCell bef. edgeLIDs:" << edgeLIDs );
01188
01189
01190 for (int v = 0 ; v < 16 ; v++){
01191
01192 if (vertexLIDs[v] < 0){
01193
01194 Array<int> maxCoFacets(4,-1);
01195
01196 int vLID = nrElem_[0];
01197 addVertex( nrElem_[0] , vertexOwner[v] , vertexIsHanging[v] ,
01198 p[0] + vertex_X[v]*h[0] , p[1] + vertex_Y[v]*h[1] , maxCoFacets );
01199 vertexLIDs[v] = vLID;
01200
01201 if (vertexIsHanging[v]){
01202 Array<int>& hinf = hangingInfo[hVertexVertexIndex[v]];
01203 if (hinf.size() < 1) hinf.resize(5,-1);
01204 hinf[hVertexIndex[v]] = vLID;
01205
01206 if ( (v == 1) || ( v==2 )) { pointMaxCoF_[vertexLIDs[v]][1] = edgeMaxCoF_[cellsEdges[1]][0];
01207 pointMaxCoF_[vertexLIDs[v]][3] = edgeMaxCoF_[cellsEdges[1]][0]; }
01208 if ( (v == 4) || ( v==8 )) { pointMaxCoF_[vertexLIDs[v]][2] = edgeMaxCoF_[cellsEdges[0]][0];
01209 pointMaxCoF_[vertexLIDs[v]][3] = edgeMaxCoF_[cellsEdges[0]][0]; }
01210 if ( (v == 7) || ( v==11 )) { pointMaxCoF_[vertexLIDs[v]][0] = edgeMaxCoF_[cellsEdges[3]][1];
01211 pointMaxCoF_[vertexLIDs[v]][1] = edgeMaxCoF_[cellsEdges[3]][1]; }
01212 if ( (v == 13) || ( v==14 )) { pointMaxCoF_[vertexLIDs[v]][0] = edgeMaxCoF_[cellsEdges[2]][1];
01213 pointMaxCoF_[vertexLIDs[v]][2] = edgeMaxCoF_[cellsEdges[2]][1]; }
01214
01215
01216
01217
01218
01219
01220 }
01221 }
01222
01223 for (int hh = 0 ; hh < 4 ; hh++){
01224 SUNDANCE_MSG3(verb() , " vertexMaxCoFacet[v][hh]:" << vertexMaxCoFacet[v][hh] << " , vertexLIDs[v]:"<< vertexLIDs[v] );
01225 if (vertexMaxCoFacet[v][hh] >= 0)
01226 pointMaxCoF_[vertexLIDs[v]][hh] = nrElem_[2] + vertexMaxCoFacet[v][hh];
01227 }
01228 SUNDANCE_MSG3(verb() , " MaxCoFac point:" << vertexLIDs[v] << " , MaxCoFac:"<< pointMaxCoF_[vertexLIDs[v]] );
01229 }
01230 SUNDANCE_MSG3(verb() , " refineCell after Vertex creation vertexLIDs:" << vertexLIDs );
01231
01232
01233 for (int e = 0 ; e < 24 ; e++){
01234
01235 if (edgeLIDs[e] < 0){
01236 Array<int> maxCoFacets(2,-1);
01237 Array<int> edgeVertexLIDs(2,-1);
01238 edgeVertexLIDs[0] = vertexLIDs[edgeVertexes[e][0]];
01239 edgeVertexLIDs[1] = vertexLIDs[edgeVertexes[e][1]];
01240 int eLID = nrElem_[1];
01241 addEdge( nrElem_[1] , edgeOwner[e] , edgeIsHanging[e] , edgeVertex[e] ,
01242 edgePBnd[e] , edgeMBnd[e] , edgeVertexLIDs , maxCoFacets );
01243 edgeLIDs[e] = eLID;
01244
01245 if (edgeIsHanging[e]){
01246 Array<int>& hinf = hangingInfo[hEdgeVertexIndex[e]];
01247 if (hinf.size() < 1) hinf.resize(5,-1);
01248 SUNDANCE_MSG3(verb() , " HNI e:" << e << ", hEdgeIndex[e]:" << hEdgeIndex[e] );
01249 hinf[2+hEdgeIndex[e]] = eLID;
01250
01251 if ( edgeMBnd[e] == false){
01252
01253 edgeMaxCoF_[edgeLIDs[e]][offsEdgeParentMaxCoFacet[e]] =
01254 edgeMaxCoF_[cellsEdges[indexEdgeParentMaxCoFacet[e]]][offsEdgeParentMaxCoFacet[e]];
01255 }
01256 }
01257 }
01258
01259 for (int hh = 0 ; hh < 2 ; hh++){
01260
01261 if (edgeMaxCoFacet[e][hh] >= 0)
01262 edgeMaxCoF_[edgeLIDs[e]][hh] = nrElem_[2] + edgeMaxCoFacet[e][hh];
01263 }
01264
01265 }
01266 SUNDANCE_MSG3(verb() , " refineCell after Edges creation edgeLIDs:" << edgeLIDs );
01267
01268
01269 for (int c = 0 ; c < 9 ; c++){
01270 Array<int> eLIDs(4,-1);
01271 Array<int> vLIDs(4,-1);
01272 vLIDs[0] = vertexLIDs[refinedCellsVertexes[c][0]]; vLIDs[1] = vertexLIDs[refinedCellsVertexes[c][1]];
01273 vLIDs[2] = vertexLIDs[refinedCellsVertexes[c][2]]; vLIDs[3] = vertexLIDs[refinedCellsVertexes[c][3]];
01274 eLIDs[0] = edgeLIDs[refinedCellsEdges[c][0]]; eLIDs[1] = edgeLIDs[refinedCellsEdges[c][1]];
01275 eLIDs[2] = edgeLIDs[refinedCellsEdges[c][2]]; eLIDs[3] = edgeLIDs[refinedCellsEdges[c][3]];
01276
01277 cellLIDs[c] = nrElem_[2];
01278 addCell( nrElem_[2] , cellOwner , c , cellLID , cellLevel_[cellLID] + 1 , eLIDs , vLIDs );
01279 }
01280
01281 SUNDANCE_MSG3(verb() , " ---- Adding hanging node information ----- " );
01282
01283 if ( (hangElmStore_[0].containsKey(cellsPoints_[cellLID][0]) == false) && ( hangingInfo[0].size() > 0 ))
01284 {
01285 hangElmStore_[0].put(cellsPoints_[cellLID][0],hangingInfo[0]);
01286
01287 }
01288 if ( (hangElmStore_[1].containsKey(cellsPoints_[cellLID][0]) == false) && ( hangingInfo[1].size() > 0 ))
01289 {
01290 hangElmStore_[1].put(cellsPoints_[cellLID][0],hangingInfo[1]);
01291
01292 }
01293 if ( (hangElmStore_[1].containsKey(cellsPoints_[cellLID][1]) == false) && ( hangingInfo[2].size() > 0 ))
01294 {
01295 hangElmStore_[1].put(cellsPoints_[cellLID][1],hangingInfo[2]);
01296
01297 }
01298 if ( (hangElmStore_[0].containsKey(cellsPoints_[cellLID][2]) == false) && ( hangingInfo[3].size() > 0 ))
01299 {
01300 hangElmStore_[0].put(cellsPoints_[cellLID][2],hangingInfo[3]);
01301
01302 }
01303
01304 SUNDANCE_MSG3(verb() , " ---- Refinement cell DONE , mark cell as not leaf and store children LID ----");
01305 isCellLeaf_[cellLID] = false;
01306
01307 cellsChildren_[cellLID] = cellLIDs;
01308 }
01309
01310
01311 void HNMesh2D::createLeafNumbering(){
01312
01313
01314
01315
01316
01317
01318
01319 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering nrPoint:" << nrElem_[0] << " , nrEdge:" << nrElem_[1] << ", nrCell:" << nrElem_[2]);
01320
01321 vertexGIDToLeafMapping_.resize(nrElem_[0],-1);
01322 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexGIDToLeafMapping_[dd] = -1;
01323 vertexLeafToGIDMapping_.resize(nrElem_[0],-1);
01324 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToGIDMapping_[dd] = -1;
01325
01326 edgeGIDToLeafMapping_.resize(nrElem_[1],-1);
01327 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeGIDToLeafMapping_[dd] = -1;
01328 edgeLeafToGIDMapping_.resize(nrElem_[1],-1);
01329 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToGIDMapping_[dd] = -1;
01330
01331 cellGIDToLeafMapping_.resize(nrElem_[2],-1);
01332 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellGIDToLeafMapping_[dd] = -1;
01333 cellLeafToGIDMapping_.resize(nrElem_[2],-1);
01334 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLeafToGIDMapping_[dd] = -1;
01335
01336 nrVertexLeafGID_ = 0; nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0;
01337
01338 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrEdgeLeafLID_ = 0;
01339 vertexLIDToLeafMapping_.resize(nrElem_[0],-1);
01340 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLIDToLeafMapping_[dd] = -1;
01341 vertexLeafToLIDMapping_.resize(nrElem_[0],-1);
01342 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToLIDMapping_[dd] = -1;
01343
01344 edgeLIDToLeafMapping_.resize(nrElem_[1],-1);
01345 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLIDToLeafMapping_[dd] = -1;
01346 edgeLeafToLIDMapping_.resize(nrElem_[1],-1);
01347 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToLIDMapping_[dd] = -1;
01348
01349 cellLIDToLeafMapping_.resize(nrElem_[2],-1);
01350 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLIDToLeafMapping_[dd] = -1;
01351 cellLeafToLIDMapping_.resize(nrElem_[2],-1);
01352 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLeafToLIDMapping_[dd] = -1;
01353
01354 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering , start assigning leaf numbers");
01355
01356 for (int ind = 0 ; ind < nrElem_[0] ; ind++){
01357
01358 }
01359
01360
01361
01362 Array<bool> hasCellLID(nrElem_[2],false);
01363
01364 for (int ind = 0 ; ind < nrElem_[2] ; ind++){
01365 Array<int>& vertexIDs = cellsPoints_[ind];
01366 hasCellLID[ind] = false;
01367 for (int v = 0 ; v < 4 ; v++){
01368 Array<int>& maxCoFacet = pointMaxCoF_[vertexIDs[v]];
01369 hasCellLID[ind] = ( hasCellLID[ind]
01370 || ( (maxCoFacet[0] >= 0) && (elementOwner_[2][maxCoFacet[0]] == myRank_) )
01371 || ( (maxCoFacet[1] >= 0) && (elementOwner_[2][maxCoFacet[1]] == myRank_) )
01372 || ( (maxCoFacet[2] >= 0) && (elementOwner_[2][maxCoFacet[2]] == myRank_) )
01373 || ( (maxCoFacet[3] >= 0) && (elementOwner_[2][maxCoFacet[3]] == myRank_) ) ) ;
01374
01375
01376
01377
01378
01379 if ( (hasCellLID[ind] == false) && (isPointHanging_[vertexIDs[v]] == true)){
01380 int parentID = parentCellLID_[ind];
01381 Array<int>& parentVertexIDs = cellsPoints_[parentID];
01382 hasCellLID[ind] = hasCellLID[ind] || (elementOwner_[0][parentVertexIDs[v]] == myRank_);
01383 }
01384 }
01385 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering Cell ID :" << ind << " should be LID: " << hasCellLID[ind] <<
01386 " ,isCellLeaf_[ind]:" << isCellLeaf_[ind]);
01387 }
01388
01389
01390
01391
01392
01393 bool check_Ghost_cells = true;
01394 while (check_Ghost_cells){
01395 check_Ghost_cells = false;
01396 for (int ind = 0 ; ind < nrElem_[2] ; ind++){
01397 if ( (hasCellLID[ind] == true) && (elementOwner_[2][ind] != myRank_ ) ){
01398
01399 Array<int>& edgeIDs = cellsEdges_[ind];
01400
01401 for (int ii = 0 ; ii < 4 ; ii++ ){
01402
01403 if (isEdgeHanging_[edgeIDs[ii]] && ( elementOwner_[1][edgeIDs[ii]] != myRank_)){
01404
01405 int parentCell = parentCellLID_[ind];
01406 Array<int>& parentEdgesIDs = cellsEdges_[parentCell];
01407 for (int f = 0 ; f < 2 ; f++)
01408 if ( ( edgeMaxCoF_[parentEdgesIDs[ii]][f] >= 0 ) &&
01409 ( elementOwner_[2][ edgeMaxCoF_[parentEdgesIDs[ii]][f] ] != myRank_ ) &&
01410 ( hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] == false) &&
01411 ( isCellLeaf_[edgeMaxCoF_[parentEdgesIDs[ii]][f]] )
01412 ){
01413 hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] = true;
01414 check_Ghost_cells = true;
01415 }
01416 }
01417 }
01418 }
01419 }
01420 }
01421
01422
01423 for (int ind = 0 ; ind < nrElem_[2] ; ind++)
01424 {
01425
01426
01427 if ( (isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
01428 {
01429 Array<int>& vertexIDs = cellsPoints_[ind];
01430 for (int v = 0; v < 4 ; v++)
01431 {
01432 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering vertexIDs[v]:" << vertexIDs[v] );
01433 if (vertexGIDToLeafMapping_[vertexIDs[v]] < 0)
01434 {
01435 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafGID_:" << nrVertexLeafGID_ );
01436 vertexLeafToGIDMapping_[nrVertexLeafGID_] = vertexIDs[v];
01437 vertexGIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafGID_;
01438 nrVertexLeafGID_++;
01439 }
01440 }
01441 Array<int>& edgeIDs = cellsEdges_[ind];
01442
01443 for (int e = 0; e < 4 ; e++)
01444 {
01445
01446 if (edgeGIDToLeafMapping_[edgeIDs[e]] < 0)
01447 {
01448 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafGID_:" << nrEdgeLeafGID_ );
01449
01450 edgeLeafToGIDMapping_[nrEdgeLeafGID_] = edgeIDs[e];
01451 edgeGIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafGID_;
01452 nrEdgeLeafGID_++;
01453 }
01454 }
01455
01456 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering CELL cellID:" << ind << " , nrCellLeafGID_:" << nrCellLeafGID_ );
01457 cellLeafToGIDMapping_[nrCellLeafGID_] = ind;
01458 cellGIDToLeafMapping_[ind] = nrCellLeafGID_;
01459 nrCellLeafGID_++;
01460
01461
01462
01463 if (hasCellLID[ind]){
01464
01465 for (int v = 0; v < 4 ; v++)
01466 {
01467 if (vertexLIDToLeafMapping_[vertexIDs[v]] < 0)
01468 {
01469 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafLID_:" << nrVertexLeafLID_ );
01470 vertexLeafToLIDMapping_[nrVertexLeafLID_] = vertexIDs[v];
01471 vertexLIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafLID_;
01472 nrVertexLeafLID_++;
01473 }
01474 }
01475
01476 for (int e = 0; e < 4 ; e++)
01477 {
01478 if (edgeLIDToLeafMapping_[edgeIDs[e]] < 0)
01479 {
01480 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafLID_:" << nrEdgeLeafLID_ );
01481 edgeLeafToLIDMapping_[nrEdgeLeafLID_] = edgeIDs[e];
01482 edgeLIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafLID_;
01483 nrEdgeLeafLID_++;
01484 }
01485 }
01486
01487 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering CELL cellID:" << ind << " , nrCellLeafLID_:" << nrCellLeafLID_ );
01488 cellLeafToLIDMapping_[nrCellLeafLID_] = ind;
01489 cellLIDToLeafMapping_[ind] = nrCellLeafLID_;
01490 nrCellLeafLID_++;
01491 }
01492 }
01493 }
01494 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering , DONE");
01495 }
01496
01497
01498
01499
01500 int HNMesh2D::estimateCellLoad(int ID){
01501 int rtn = 0;
01502 if (isCellLeaf_[ID]){
01503 if (!isCellOut_[ID]) rtn = 1;
01504 } else {
01505
01506 for (int r = 0 ; r < (int)cellsChildren_[ID].size() ; r++){
01507 rtn = rtn + estimateCellLoad(cellsChildren_[ID][r]);
01508 }
01509 }
01510 return rtn;
01511 }
01512
01513
01514 void HNMesh2D::markCellsAndFacets(int cellID , int procID){
01515
01516 if (elementOwner_[2][cellID] < 0) { elementOwner_[2][cellID] = procID; }
01517
01518 if (elementOwner_[1][cellsEdges_[cellID][0]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][0]] = procID;}
01519 if (elementOwner_[1][cellsEdges_[cellID][1]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][1]] = procID;}
01520 if (elementOwner_[1][cellsEdges_[cellID][2]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][2]] = procID;}
01521 if (elementOwner_[1][cellsEdges_[cellID][3]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][3]] = procID;}
01522
01523
01524 if (elementOwner_[0][cellsPoints_[cellID][0]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][0]] = procID;}
01525 if (elementOwner_[0][cellsPoints_[cellID][1]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][1]] = procID;}
01526 if (elementOwner_[0][cellsPoints_[cellID][2]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][2]] = procID;}
01527 if (elementOwner_[0][cellsPoints_[cellID][3]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][3]] = procID;}
01528
01529
01530 if (!isCellLeaf_[cellID]){
01531
01532 for (int r = 0 ; r < (int)cellsChildren_[cellID].size() ; r++){
01533 markCellsAndFacets(cellsChildren_[cellID][r] , procID);
01534 }
01535 }
01536 }
01537
01538 void HNMesh2D::createLeafNumbering_sophisticated(){
01539
01540
01541 Array<bool> hasCellLID(nrElem_[2],false);
01542 double total_load = 0.0;
01543 int nrCoarseCell = _res_x * _res_y;
01544 Array<int> coarseCellLoad( _res_x * _res_y , 1 );
01545
01546
01547
01548
01549
01550
01551
01552
01553 for (int ind = 0 ; ind < nrElem_[2] ; ind++){
01554 if (ind < nrCoarseCell) {
01555
01556 coarseCellLoad[ind] = estimateCellLoad(ind);
01557 }
01558 if ((isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
01559 { total_load = total_load + 1 ; }
01560 }
01561
01562 SUNDANCE_MSG3(verb() , "total_load = " << total_load << " , nrCell = " << nrElem_[2]);
01563
01564
01565
01566 int levelM = ::ceil( ::fmax( ::log2(_res_x) , ::log2(_res_y ) ) );
01567
01568 Array<int> vectX1(4), vectY1(4), vectX2(4), vectY2(4);
01569 vectX1[0] = 0; vectX1[1] = (int)::pow(2,levelM-1); vectX1[2] = 0; vectX1[3] = (int)::pow(2,levelM-1);
01570 vectY1[0] = 0; vectY1[1] = 0; vectY1[2] = (int)::pow(2,levelM-1); vectY1[3] = (int)::pow(2,levelM-1);
01571 vectX2[0] = 0; vectX2[1] = (int)::pow(2,levelM-1); vectX2[2] = 0; vectX2[3] = (int)::pow(2,levelM-1);
01572 vectY2[0] = 0; vectY2[1] = 0; vectY2[2] = (int)::pow(2,levelM-1); vectY2[3] = (int)::pow(2,levelM-1);
01573 int addX[4] = { 0 , 1 , 0 , 1};
01574 int addY[4] = { 0 , 0 , 1 , 1};
01575 Array<int> *inX = &vectX1 , *inY = &vectY1 , *outX = &vectX2 , *outY = &vectY2 , *tmpVectP;
01576 int levelActual = levelM - 2;
01577
01578 while (levelActual >= 0){
01579 outX->resize( 4 * inX->size() );
01580 outY->resize( 4 * inY->size() );
01581 int cI = 0 , addO = (int)::pow(2,levelActual);
01582 SUNDANCE_MSG3(verb() , " outX->size():" << outX->size() << ", levelActual:" << levelActual << " , addO:" << addO);
01583
01584 for (int ce = 0 ; ce < inX->size() ; ce++){
01585 (*outX)[cI+0] = (*inX)[ce] + addO*addX[0];
01586 (*outX)[cI+1] = (*inX)[ce] + addO*addX[1];
01587 (*outX)[cI+2] = (*inX)[ce] + addO*addX[2];
01588 (*outX)[cI+3] = (*inX)[ce] + addO*addX[3];
01589 (*outY)[cI+0] = (*inY)[ce] + addO*addY[0];
01590 (*outY)[cI+1] = (*inY)[ce] + addO*addY[1];
01591 (*outY)[cI+2] = (*inY)[ce] + addO*addY[2];
01592 (*outY)[cI+3] = (*inY)[ce] + addO*addY[3];
01593 cI = cI + 4;
01594 }
01595 SUNDANCE_MSG3(verb() , " EX: " << (*outX)[0] << " , " << (*outX)[1] << " , " << (*outX)[2]);
01596 SUNDANCE_MSG3(verb() , " EY: " << (*outY)[0] << " , " << (*outY)[1] << " , " << (*outY)[2]);
01597
01598 levelActual = levelActual - 1;
01599 tmpVectP = inX; inX = outX; outX = tmpVectP;
01600 tmpVectP = inY; inY = outY; outY = tmpVectP;
01601 }
01602
01603 tmpVectP = inX; inX = outX; outX = tmpVectP;
01604 tmpVectP = inY; inY = outY; outY = tmpVectP;
01605
01606
01607 for (int tmp = 0 ; tmp < nrElem_[0] ; tmp++ ){ elementOwner_[0][tmp] = -1; }
01608 for (int tmp = 0 ; tmp < nrElem_[1] ; tmp++ ){ elementOwner_[1][tmp] = -1; }
01609 for (int tmp = 0 ; tmp < nrElem_[2] ; tmp++ ){ elementOwner_[2][tmp] = -1; }
01610
01611
01612 int coarseCellID , actProcID = 0 , actualLoad = 0;
01613 double loadPerProc = (double)total_load / (double)nrProc_ , diff_load = 0.0;
01614 for (int ind = 0 ; ind < outX->size() ; ind++){
01615
01616 if ( ((*outX)[ind] < _res_x) && ((*outY)[ind] < _res_y) ){
01617
01618 coarseCellID = ((*outX)[ind])*_res_y + ((*outY)[ind]);
01619 SUNDANCE_MSG3(verb(),"Z-curve trav. ind:" << ind << " , coarseCellID:" << coarseCellID << " , indX:" << (*outX)[ind] << " , indY:" << (*outY)[ind]);
01620
01621 TEUCHOS_TEST_FOR_EXCEPTION( cellLevel_[coarseCellID] > 0 , std::logic_error, " coarseCellID:" << coarseCellID << " has level:" << cellLevel_[coarseCellID] );
01622 markCellsAndFacets( coarseCellID , actProcID);
01623 actualLoad = actualLoad + coarseCellLoad[coarseCellID];
01624
01625 if (((double)actualLoad >= (loadPerProc - 1e-8 - diff_load)) && ( actProcID < nrProc_-1 )){
01626 SUNDANCE_MSG3(verb() , "Increase CPU , actualLoad:" << actualLoad << " loadPerProc:" << loadPerProc );
01627
01628 diff_load = actualLoad - loadPerProc;
01629 actProcID = actProcID + 1;
01630 actualLoad = 0;
01631 }
01632 }
01633 }
01634
01635
01636 SUNDANCE_MSG3(verb()," nrElem_[0]:" << nrElem_[0] << " , nrElem_[1]:" << nrElem_[1] << " , nrElem_[2]" << nrElem_[2]);
01637 for (int tmp = 0 ; tmp < nrElem_[0] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[0][tmp] < 0 , std::logic_error, " 0 tmp:" << tmp); }
01638 for (int tmp = 0 ; tmp < nrElem_[1] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[1][tmp] < 0 , std::logic_error, " 1 tmp:" << tmp); }
01639 for (int tmp = 0 ; tmp < nrElem_[2] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[2][tmp] < 0 , std::logic_error, " 2 tmp:" << tmp); }
01640
01641
01642
01643
01644
01645
01646
01647 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering nrPoint:" << nrElem_[0] << " , nrEdge:" << nrElem_[1] << ", nrCell:" << nrElem_[2]);
01648
01649 vertexGIDToLeafMapping_.resize(nrElem_[0],-1);
01650 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexGIDToLeafMapping_[dd] = -1;
01651 vertexLeafToGIDMapping_.resize(nrElem_[0],-1);
01652 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToGIDMapping_[dd] = -1;
01653
01654 edgeGIDToLeafMapping_.resize(nrElem_[1],-1);
01655 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeGIDToLeafMapping_[dd] = -1;
01656 edgeLeafToGIDMapping_.resize(nrElem_[1],-1);
01657 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToGIDMapping_[dd] = -1;
01658
01659 cellGIDToLeafMapping_.resize(nrElem_[2],-1);
01660 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellGIDToLeafMapping_[dd] = -1;
01661 cellLeafToGIDMapping_.resize(nrElem_[2],-1);
01662 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLeafToGIDMapping_[dd] = -1;
01663
01664 nrVertexLeafGID_ = 0; nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0;
01665
01666 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrEdgeLeafLID_ = 0;
01667 vertexLIDToLeafMapping_.resize(nrElem_[0],-1);
01668 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLIDToLeafMapping_[dd] = -1;
01669 vertexLeafToLIDMapping_.resize(nrElem_[0],-1);
01670 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToLIDMapping_[dd] = -1;
01671
01672 edgeLIDToLeafMapping_.resize(nrElem_[1],-1);
01673 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLIDToLeafMapping_[dd] = -1;
01674 edgeLeafToLIDMapping_.resize(nrElem_[1],-1);
01675 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToLIDMapping_[dd] = -1;
01676
01677 cellLIDToLeafMapping_.resize(nrElem_[2],-1);
01678 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLIDToLeafMapping_[dd] = -1;
01679 cellLeafToLIDMapping_.resize(nrElem_[2],-1);
01680 for (int dd = 0 ; dd < nrElem_[2] ; dd++) cellLeafToLIDMapping_[dd] = -1;
01681
01682 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering , start assigning leaf numbers");
01683
01684 for (int ind = 0 ; ind < nrElem_[2] ; ind++){
01685 Array<int>& vertexIDs = cellsPoints_[ind];
01686 hasCellLID[ind] = false;
01687 for (int v = 0 ; v < 4 ; v++){
01688 Array<int>& maxCoFacet = pointMaxCoF_[vertexIDs[v]];
01689 hasCellLID[ind] = ( hasCellLID[ind]
01690 || ( (maxCoFacet[0] >= 0) && (elementOwner_[2][maxCoFacet[0]] == myRank_) )
01691 || ( (maxCoFacet[1] >= 0) && (elementOwner_[2][maxCoFacet[1]] == myRank_) )
01692 || ( (maxCoFacet[2] >= 0) && (elementOwner_[2][maxCoFacet[2]] == myRank_) )
01693 || ( (maxCoFacet[3] >= 0) && (elementOwner_[2][maxCoFacet[3]] == myRank_) ) ) ;
01694
01695
01696
01697
01698
01699 if ( (hasCellLID[ind] == false) && (isPointHanging_[vertexIDs[v]] == true)){
01700 int parentID = parentCellLID_[ind];
01701 Array<int>& parentVertexIDs = cellsPoints_[parentID];
01702 hasCellLID[ind] = hasCellLID[ind] || (elementOwner_[0][parentVertexIDs[v]] == myRank_);
01703 }
01704 }
01705 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering Cell ID :" << ind << " should be LID: " << hasCellLID[ind] <<
01706 " ,isCellLeaf_[ind]:" << isCellLeaf_[ind]);
01707 }
01708
01709
01710
01711
01712
01713 bool check_Ghost_cells = true;
01714 while (check_Ghost_cells){
01715 check_Ghost_cells = false;
01716 for (int ind = 0 ; ind < nrElem_[2] ; ind++){
01717 if ( (hasCellLID[ind] == true) && (elementOwner_[2][ind] != myRank_ ) ){
01718
01719 Array<int>& edgeIDs = cellsEdges_[ind];
01720
01721 for (int ii = 0 ; ii < 4 ; ii++ ){
01722
01723 if (isEdgeHanging_[edgeIDs[ii]] && ( elementOwner_[1][edgeIDs[ii]] != myRank_)){
01724
01725 int parentCell = parentCellLID_[ind];
01726 Array<int>& parentEdgesIDs = cellsEdges_[parentCell];
01727 for (int f = 0 ; f < 2 ; f++)
01728 if ( ( edgeMaxCoF_[parentEdgesIDs[ii]][f] >= 0 ) &&
01729 ( elementOwner_[2][ edgeMaxCoF_[parentEdgesIDs[ii]][f] ] != myRank_ ) &&
01730 ( hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] == false) &&
01731 ( isCellLeaf_[edgeMaxCoF_[parentEdgesIDs[ii]][f]] )
01732 ){
01733 hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] = true;
01734 check_Ghost_cells = true;
01735 }
01736 }
01737 }
01738 }
01739 }
01740 }
01741
01742
01743 for (int ind = 0 ; ind < nrElem_[2] ; ind++)
01744 {
01745
01746
01747 if ( (isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
01748 {
01749 Array<int>& vertexIDs = cellsPoints_[ind];
01750 for (int v = 0; v < 4 ; v++)
01751 {
01752 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering vertexIDs[v]:" << vertexIDs[v] );
01753 if (vertexGIDToLeafMapping_[vertexIDs[v]] < 0)
01754 {
01755 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafGID_:" << nrVertexLeafGID_ );
01756 vertexLeafToGIDMapping_[nrVertexLeafGID_] = vertexIDs[v];
01757 vertexGIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafGID_;
01758 nrVertexLeafGID_++;
01759 }
01760 }
01761 Array<int>& edgeIDs = cellsEdges_[ind];
01762
01763 for (int e = 0; e < 4 ; e++)
01764 {
01765
01766 if (edgeGIDToLeafMapping_[edgeIDs[e]] < 0)
01767 {
01768 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafGID_:" << nrEdgeLeafGID_ );
01769
01770 edgeLeafToGIDMapping_[nrEdgeLeafGID_] = edgeIDs[e];
01771 edgeGIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafGID_;
01772 nrEdgeLeafGID_++;
01773 }
01774 }
01775
01776 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering CELL cellID:" << ind << " , nrCellLeafGID_:" << nrCellLeafGID_ );
01777 cellLeafToGIDMapping_[nrCellLeafGID_] = ind;
01778 cellGIDToLeafMapping_[ind] = nrCellLeafGID_;
01779 nrCellLeafGID_++;
01780
01781
01782
01783 if (hasCellLID[ind]){
01784
01785 for (int v = 0; v < 4 ; v++)
01786 {
01787 if (vertexLIDToLeafMapping_[vertexIDs[v]] < 0)
01788 {
01789 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafLID_:" << nrVertexLeafLID_ );
01790 vertexLeafToLIDMapping_[nrVertexLeafLID_] = vertexIDs[v];
01791 vertexLIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafLID_;
01792 nrVertexLeafLID_++;
01793 }
01794 }
01795
01796 for (int e = 0; e < 4 ; e++)
01797 {
01798 if (edgeLIDToLeafMapping_[edgeIDs[e]] < 0)
01799 {
01800 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafLID_:" << nrEdgeLeafLID_ );
01801 edgeLeafToLIDMapping_[nrEdgeLeafLID_] = edgeIDs[e];
01802 edgeLIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafLID_;
01803 nrEdgeLeafLID_++;
01804 }
01805 }
01806
01807 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering CELL cellID:" << ind << " , nrCellLeafLID_:" << nrCellLeafLID_ );
01808 cellLeafToLIDMapping_[nrCellLeafLID_] = ind;
01809 cellLIDToLeafMapping_[ind] = nrCellLeafLID_;
01810 nrCellLeafLID_++;
01811 }
01812 }
01813 }
01814 SUNDANCE_MSG3(verb() , "HNMesh2D::createLeafNumbering , DONE");
01815 }