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 "SundanceHNMesh3D.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 HNMesh3D::offs_Points_x_[8] = {0, 1, 0, 1 , 0 , 1 , 0 , 1};
00058
00059 int HNMesh3D::offs_Points_y_[8] = {0, 0, 1, 1 , 0 , 0 , 1 , 1};
00060
00061 int HNMesh3D::offs_Points_z_[8] = {0, 0, 0, 0 , 1 , 1 , 1 , 1 };
00062
00063 int HNMesh3D::edge_Points_localIndex[12][2] = { {0,1} , {0,2} , {0,4} , {1,3} , {1,5} , {2,3} , {2,6} , {3,7} ,
00064 {4,5} , {4,6} , {5,7} , {6,7} };
00065
00066 int HNMesh3D::edge_Orientation[12] = { 0, 1, 2, 1, 2, 0, 2, 2, 0, 1, 1, 0 };
00067 int HNMesh3D::edge_MaxCofacetIndex[3][4] = { {0,5,8,11},{1,3,9,10},{2,4,6,7} };
00068 int HNMesh3D::edge_MaxCof[12] = { 0,0,0, 1, 1, 1, 2, 3, 2, 2, 3, 3 };
00069
00070 int HNMesh3D::face_Points_localIndex[6][4] = { {0,1,2,3} , {0,1,4,5} , {0,2,4,6} , {1,3,5,7} , {2,3,6,7} , {4,5,6,7} };
00071
00072 int HNMesh3D::face_Edges_localIndex[6][4]= { {0,1,3,5} , {0,2,4,8} , {1,2,6,9} , {3,4,7,10}, {5,6,7,11}, {8,9,10,11}};
00073
00074 int HNMesh3D::face_Orientation[6] = { 0,1,2,2,1,0 };
00075 int HNMesh3D::face_MaxCofacetIndex[3][2] = { {0,5},{1,4},{2,3}};
00076 int HNMesh3D::face_MaxCof[6] = { 0,0,0,1,1,1 };
00077
00078
00079 int HNMesh3D::vInd[8];
00080 int HNMesh3D::eInd[12];
00081 int HNMesh3D::fInd[6];
00082
00083
00084 double HNMesh3D::vertex_X[64] =
00085 { 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00086 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00087 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00088 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00089 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00090 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00091 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 ,
00092 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 , 0.0 , 1.0/3.0 , 2.0/3.0 , 1.0 };
00093
00094 double HNMesh3D::vertex_Y[64] =
00095 { 0.0 , 0.0 , 0.0 , 0.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00096 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 1.0 , 1.0 , 1.0 , 1.0 ,
00097 0.0 , 0.0 , 0.0 , 0.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00098 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 1.0 , 1.0 , 1.0 , 1.0 ,
00099 0.0 , 0.0 , 0.0 , 0.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00100 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 1.0 , 1.0 , 1.0 , 1.0 ,
00101 0.0 , 0.0 , 0.0 , 0.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00102 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 1.0 , 1.0 , 1.0 , 1.0 };
00103 double HNMesh3D::vertex_Z[64] =
00104 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
00105 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
00106 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00107 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 , 1.0/3.0 ,
00108 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 ,
00109 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 , 2.0/3.0 ,
00110 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 ,
00111 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 , 1.0 };
00112
00113 int HNMesh3D::vertexToParentEdge[64] =
00114 { -1, 0, 0, -1, 1, 20, 20, 3, 1, 20, 20, 3, -1, 5, 5, -1, 2, 21, 21, 4,
00115 22, -1, -1, 23, 22, -1, -1, 23, 6, 24, 24, 7, 2, 21, 21, 4, 22, -1, -1, 23,
00116 22, -1, -1, 23, 6, 24, 24, 7, -1, 8, 8, -1, 9, 25, 25, 10, 9, 25, 25, 10,
00117 -1, 11, 11, -1 };
00118
00119 int HNMesh3D::vertexInParentIndex[64] =
00120 { -1, 0, 1, -1, 0, 20, 21, 0, 1, 22, 23, 1, -1, 0, 1, -1, 0, 20, 21, 0,
00121 20, -1, -1, 20, 21, -1, -1, 21, 0, 20, 21, 0, 1, 22, 23, 1, 22, -1, -1, 22,
00122 23, -1, -1, 23, 1, 22, 23, 1, -1, 0, 1, -1, 0, 20, 21, 0, 1, 22, 23, 1,
00123 -1, 0, 1, -1 };
00124
00125 int HNMesh3D::edgeToParentEdge[144] =
00126 { 0, 0, 0, 1, 20, 20, 3, 20, 20, 20, 1, 20, 20, 3, 20, 20, 20, 1, 20, 20,
00127 3, 5, 5, 5, 2, 21, 21, 4, 22, -1, -1, 23, 22, -1, -1, 23, 6, 24, 24, 7,
00128 21, 21, 21, 22, -1, -1, 23, -1, -1, -1, 22, -1, -1, 23, -1, -1, -1, 22, -1, -1,
00129 23, 24, 24, 24, 2, 21, 21, 4, 22, -1, -1, 23, 22, -1, -1, 23, 6, 24, 24, 7,
00130 21, 21, 21, 22, -1, -1, 23, -1, -1, -1, 22, -1, -1, 23, -1, -1, -1, 22, -1, -1,
00131 23, 24, 24, 24, 2, 21, 21, 4, 22, -1, -1, 23, 22, -1, -1, 23, 6, 24, 24, 7,
00132 8, 8, 8, 9, 25, 25, 10, 25, 25, 25, 9, 25, 25, 10, 25, 25, 25, 9, 25, 25,
00133 10, 11, 11, 11 };
00134
00135 int HNMesh3D::edgeInParentIndex[144] =
00136 { 0, 1, 2, 0, 20, 21, 0, 22, 23, 24, 1, 25, 26, 1, 27, 28, 29, 2, 30, 31,
00137 2, 0, 1, 2, 0, 20, 21, 0, 20, -1, -1, 20, 21, -1, -1, 21, 0, 20, 21, 0,
00138 22, 23, 24, 22, -1, -1, 22, -1, -1, -1, 23, -1, -1, 23, -1, -1, -1, 24, -1, -1,
00139 24, 22, 23, 24, 1, 25, 26, 1, 25, -1, -1, 25, 26, -1, -1, 26, 1, 25, 26, 1,
00140 27, 28, 29, 27, -1, -1, 27, -1, -1, -1, 28, -1, -1, 28, -1, -1, -1, 29, -1, -1,
00141 29, 27, 28, 29, 2, 30, 31, 2, 30, -1, -1, 30, 31, -1, -1, 31, 2, 30, 31, 2,
00142 0, 1, 2, 0, 20, 21, 0, 22, 23, 24, 1, 25, 26, 1, 27, 28, 29, 2, 30, 31,
00143 2, 0, 1, 2 };
00144
00145 int HNMesh3D::faceToParentFace[108] =
00146 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, -1, -1, 3, -1, -1, -1, 2,
00147 -1, -1, 3, -1, -1, -1, 2, -1, -1, 3, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1,
00148 -1, -1, 1, 1, 1, 2, -1, -1, 3, -1, -1, -1, 2, -1, -1, 3, -1, -1, -1, 2,
00149 -1, -1, 3, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 2, -1,
00150 -1, 3, -1, -1, -1, 2, -1, -1, 3, -1, -1, -1, 2, -1, -1, 3, 4, 4, 4, 5,
00151 5, 5, 5, 5, 5, 5, 5, 5 };
00152
00153 int HNMesh3D::faceInParentIndex[108] =
00154 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 0, -1, -1, 0, -1, -1, -1, 1,
00155 -1, -1, 1, -1, -1, -1, 2, -1, -1, 2, 0, 1, 2, -1, -1, -1, -1, -1, -1, -1,
00156 -1, -1, 3, 4, 5, 3, -1, -1, 3, -1, -1, -1, 4, -1, -1, 4, -1, -1, -1, 5,
00157 -1, -1, 5, 3, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, 7, 8, 6, -1,
00158 -1, 6, -1, -1, -1, 7, -1, -1, 7, -1, -1, -1, 8, -1, -1, 8, 6, 7, 8, 0,
00159 1, 2, 3, 4, 5, 6, 7, 8 };
00160
00161
00162 HNMesh3D::HNMesh3D(int dim, const MPIComm& comm ,
00163 const MeshEntityOrder& order)
00164 : MeshBase(dim, comm , order), _comm(comm)
00165 {
00166 setVerb(0);
00167
00168
00169 nrProc_ = MPIComm::world().getNProc();
00170 myRank_ = MPIComm::world().getRank();
00171
00172 points_.resize(0);
00173 nrElem_.resize(4,0);
00174 nrElemOwned_.resize(4,0);
00175
00176 cellsPoints_.resize(0);
00177 cellsEdges_.resize(0);
00178 cellsFaces_.resize(0);
00179 isCellOut_.resize(0);
00180 faceEdges_.resize(0);
00181 facePoints_.resize(0);
00182 edgePoints_.resize(0);
00183 edgeOrientation_.resize(0);
00184 faceOrientation_.resize(0);
00185
00186 faceMaxCoF_.resize(0);
00187 edgeMaxCoF_.resize(0);
00188 pointMaxCoF_.resize(0);
00189
00190 elementOwner_.resize(4); elementOwner_[0].resize(0); elementOwner_[1].resize(0); elementOwner_[2].resize(0); elementOwner_[3].resize(0);
00191
00192 indexInParent_.resize(0);
00193 parentCellLID_.resize(0);
00194 cellLevel_.resize(0);
00195 isCellLeaf_.resize(0);
00196
00197 isPointHanging_.resize(0);
00198 isEdgeHanging_.resize(0);
00199
00200 edgeHangingElmStore_ = Hashtable< int, Array<int> >();
00201 hangingAccessCount_.resize(0);
00202 faceHangingElmStore_ = Hashtable< int, Array<int> >();
00203 refineCell_.resize(0);
00204
00205 nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0; nrFaceLeafGID_ = 0; nrVertexLeafGID_ = 0;
00206 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrFaceLeafLID_ = 0; nrEdgeLeafLID_ = 0;
00207 }
00208
00209 int HNMesh3D::numCells(int dim) const {
00210 SUNDANCE_MSG3(verb(),"HNMesh3D::numCells(int dim): dim:" << dim );
00211 switch (dim){
00212 case 0: return nrVertexLeafLID_;
00213 case 1: return nrEdgeLeafLID_;
00214 case 2: return nrFaceLeafLID_;
00215 case 3: return nrCellLeafLID_;
00216 }
00217 return 0;
00218 }
00219
00220 Point HNMesh3D::nodePosition(int i) const {
00221 SUNDANCE_MSG3(verb(),"HNMesh3D::nodePosition(int i) i:"<< i);
00222
00223 return points_[vertexLeafToLIDMapping_[i]];
00224 }
00225
00226 const double* HNMesh3D::nodePositionView(int i) const {
00227 SUNDANCE_MSG3(verb(),"HNMesh3D::nodePositionView(int i) i:" << i);
00228
00229 return &(points_[vertexLeafToLIDMapping_[i]][0]);;
00230 }
00231
00232 void HNMesh3D::getJacobians(int cellDim, const Array<int>& cellLID,
00233 CellJacobianBatch& jBatch) const
00234 {
00235
00236 SUNDANCE_MSG3(verb(),"HNMesh3D::getJacobians cellDim:"<<cellDim<<" _x:"<<_ofs_x<<" _y:"<<_ofs_y<<" _z:"<<_ofs_z);
00237 SUNDANCE_VERB_HIGH("getJacobians()");
00238 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00239 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00240 int nCells = cellLID.size();
00241 int LID;
00242 Point pnt(0.0,0.0,0.0);
00243 jBatch.resize(cellLID.size(), spatialDim(), cellDim);
00244 if (cellDim < spatialDim())
00245 {
00246 for (int i=0; i<nCells; i++)
00247 {
00248 double* detJ = jBatch.detJ(i);
00249 switch(cellDim)
00250 {
00251 case 0:{ *detJ = 1.0;
00252 break;}
00253 case 1:{
00254 LID = edgeLeafToLIDMapping_[cellLID[i]];
00255 pnt = (points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]]);
00256 *detJ = sqrt(pnt * pnt);
00257 break;}
00258 case 2:{
00259 LID = faceLeafToLIDMapping_[cellLID[i]];
00260 int a = facePoints_[LID][0];
00261 int b = facePoints_[LID][1];
00262 int c = facePoints_[LID][2];
00263 const Point& pa = points_[a];
00264 const Point& pb = points_[b];
00265 const Point& pc = points_[c];
00266 Point directedArea = cross( pb - pa , pc - pa );
00267 *detJ = sqrt(directedArea * directedArea);
00268 break;}
00269 default:
00270 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh3D::getJacobians()");
00271 }
00272 }
00273 }else{
00274
00275 SUNDANCE_VERB_HIGH("cellDim == spatialDim()");
00276 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00277 {
00278 double* J = jBatch.jVals(i);
00279 switch(cellDim)
00280 {
00281 case 3:{
00282 LID = cellLeafToLIDMapping_[cellLID[i]];
00283
00284 J[0] = points_[cellsPoints_[LID][1]][0] - points_[cellsPoints_[LID][0]][0];
00285 J[1] = 0.0; J[2] = 0.0; J[3] = 0.0;
00286 J[4] = points_[cellsPoints_[LID][2]][1] - points_[cellsPoints_[LID][0]][1];
00287 J[5] = 0.0; J[6] = 0.0; J[7] = 0.0;
00288 J[8] = points_[cellsPoints_[LID][4]][2] - points_[cellsPoints_[LID][0]][2];
00289 SUNDANCE_MSG3(verb() , "HNMesh3D::getJacobians LID:" << LID << " X:" << J[0] << " Y:" << J[4] << " Z:" << J[8]);
00290
00291
00292
00293
00294 break;}
00295 default:
00296 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh3D::getJacobians()");
00297 }
00298 }
00299 }
00300 }
00301
00302 void HNMesh3D::getCellDiameters(int cellDim, const Array<int>& cellLID,
00303 Array<double>& cellDiameters) const {
00304
00305 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00306 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00307 SUNDANCE_VERB_HIGH("getCellDiameters()");
00308 cellDiameters.resize(cellLID.size());
00309 Point pnt(0.0 , 0.0 , 0.0 );
00310 int LID;
00311 if (cellDim < spatialDim())
00312 {
00313 SUNDANCE_MSG3(verb(),"HNMesh3D::getCellDiameters(), cellDim < spatialDim() ");
00314 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00315 {
00316 switch(cellDim)
00317 {
00318 case 0:
00319 cellDiameters[i] = 1.0;
00320 break;
00321 case 1:
00322 LID = edgeLeafToLIDMapping_[cellLID[i]];
00323 pnt = (points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]]);
00324 cellDiameters[i] = sqrt(pnt * pnt);
00325 break;
00326 case 2:
00327 LID = faceLeafToLIDMapping_[cellLID[i]];
00328 pnt = (points_[facePoints_[LID][3]] - points_[facePoints_[LID][0]]);
00329 cellDiameters[i] = sqrt(pnt * pnt);
00330 break;
00331 default:
00332 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "cellDim=" << cellDim << " in HNMesh3D::getCellDiameters()");
00333 }
00334 }
00335 }
00336 else
00337 {
00338 SUNDANCE_MSG3(verb(),"HNMesh3D::getCellDiameters(), cellDim == spatialDim() ");
00339 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00340 {
00341 switch(cellDim)
00342 {
00343 case 3:
00344 LID = cellLeafToLIDMapping_[cellLID[i]];
00345 pnt = points_[cellsPoints_[LID][7]] - points_[cellsPoints_[LID][0]];
00346 cellDiameters[i] = sqrt(pnt * pnt);
00347 break;
00348 default:
00349 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00350 "cellDim=" << cellDim << " in HNMesh3D::getCellDiameters()");
00351 }
00352 }
00353 }
00354 }
00355
00356 void HNMesh3D::pushForward(int cellDim, const Array<int>& cellLID,
00357 const Array<Point>& refQuadPts,
00358 Array<Point>& physQuadPts) const {
00359
00360 SUNDANCE_MSG3(verb(),"HNMesh3D::pushForward cellDim:"<<cellDim);
00361 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00362 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00363
00364 int nQuad = refQuadPts.size();
00365 Point pnt( 0.0 , 0.0 , 0.0);
00366 Point pnt1( 0.0 , 0.0 , 0.0);
00367
00368 if (physQuadPts.size() > 0) physQuadPts.resize(0);
00369 physQuadPts.reserve(cellLID.size() * refQuadPts.size());
00370 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00371 {
00372 switch(cellDim)
00373 {
00374 case 0:
00375 physQuadPts.append(pnt);
00376 break;
00377 case 1:{
00378 int LID = edgeLeafToLIDMapping_[cellLID[i]];
00379 pnt = points_[edgePoints_[LID][0]];
00380 pnt1 = points_[edgePoints_[LID][1]] - points_[edgePoints_[LID][0]];
00381 for (int q=0; q<nQuad; q++) {
00382 physQuadPts.append(pnt + (pnt1)*refQuadPts[q][0]);
00383 }
00384 break;}
00385 case 2:{
00386 int LID = faceLeafToLIDMapping_[cellLID[i]];
00387 pnt = points_[facePoints_[LID][0]];
00388
00389 pnt1 = points_[facePoints_[LID][3]] - points_[facePoints_[LID][0]];
00390 for (int q=0; q<nQuad; q++) {
00391 physQuadPts.append( pnt + Point( refQuadPts[q][0] * pnt1[0] , refQuadPts[q][1] * pnt1[1] , refQuadPts[q][2] * pnt1[2]) );
00392 }
00393 break;}
00394 case 3:{
00395 int LID = cellLeafToLIDMapping_[cellLID[i]];
00396 pnt = points_[cellsPoints_[LID][0]];
00397
00398 pnt1 = points_[cellsPoints_[LID][7]] - points_[cellsPoints_[LID][0]];
00399 for (int q=0; q<nQuad; q++) {
00400 physQuadPts.append( pnt + Point( refQuadPts[q][0] * pnt1[0] , refQuadPts[q][1] * pnt1[1] , refQuadPts[q][2] * pnt1[2] ) );
00401 }
00402 break;}
00403 default:
00404 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value " "in HNMesh3D::getJacobians()");
00405 }
00406 }
00407 }
00408
00409 int HNMesh3D::ownerProcID(int cellDim, int cellLID) const {
00410 int ID = -1;
00411 if (cellDim == 0) ID = vertexLeafToLIDMapping_[cellLID];
00412 if (cellDim == 1) ID = edgeLeafToLIDMapping_[cellLID];
00413 if (cellDim == 2) ID = faceLeafToLIDMapping_[cellLID];
00414 if (cellDim == 3) ID = cellLeafToLIDMapping_[cellLID];
00415 SUNDANCE_MSG3(verb() , " HNMesh3D::ownerProcID ,cellDim:" << cellDim << ", cellLID:"
00416 << cellLID <<" ,ID:" << ID << ", ret:"<< elementOwner_[cellDim][ID] );
00417 return elementOwner_[cellDim][ID];
00418 }
00419
00420
00421 int HNMesh3D::numFacets(int cellDim, int cellLID,
00422 int facetDim) const {
00423
00424 if (cellDim==1) {
00425 return 2;
00426 }
00427 else if (cellDim==2) {
00428 return 4;
00429 }
00430 else if (cellDim==3) {
00431 if (facetDim == 0) return 8;
00432 if (facetDim == 1) return 12;
00433 if (facetDim == 2) return 6;
00434 }
00435 return -1;
00436 }
00437
00438 int HNMesh3D::facetLID(int cellDim, int cellLID,
00439 int facetDim, int facetIndex,
00440 int& facetOrientation) const {
00441
00442
00443 facetOrientation = 1;
00444 SUNDANCE_MSG3(verb(),"HNMesh3D::facetLID cellDim:"<<cellDim<<", cellLID:"<<cellLID<<", facetDim:"<<facetDim<< ", facetIndex:" << facetIndex);
00445 int rnt = -1 , LID=-1 , tmp=-1;
00446 if (facetDim == 0){
00447 if (cellDim == 3 ){
00448 LID = cellLeafToLIDMapping_[cellLID];
00449 rnt = cellsPoints_[LID][facetIndex]; tmp = rnt;
00450 rnt = vertexLIDToLeafMapping_[rnt];
00451 }
00452 else if ((cellDim==2)){
00453 LID = faceLeafToLIDMapping_[cellLID];
00454 rnt = facePoints_[LID][facetIndex]; tmp = rnt;
00455 rnt = vertexLIDToLeafMapping_[rnt];
00456 }
00457 else if ((cellDim==1)){
00458 LID = edgeLeafToLIDMapping_[cellLID];
00459 rnt = edgePoints_[LID][facetIndex]; tmp = rnt;
00460 rnt = vertexLIDToLeafMapping_[rnt];
00461 }
00462 } else if (facetDim == 1){
00463 if (cellDim == 3 ){
00464 LID = cellLeafToLIDMapping_[cellLID];
00465 rnt = cellsEdges_[LID][facetIndex]; tmp = rnt;
00466 rnt = edgeLIDToLeafMapping_[rnt];
00467 } else if ((cellDim==2)){
00468 LID = faceLeafToLIDMapping_[cellLID];
00469 rnt = faceEdges_[LID][facetIndex]; tmp = rnt;
00470 rnt = edgeLIDToLeafMapping_[rnt];
00471 }
00472 } else if (facetDim == 2){
00473 if (cellDim == 3 ){
00474 LID = cellLeafToLIDMapping_[cellLID];
00475 rnt = cellsFaces_[LID][facetIndex]; tmp = rnt;
00476 rnt = faceLIDToLeafMapping_[rnt];
00477 }
00478 }
00479 SUNDANCE_MSG3(verb()," RET = " << rnt << ", LID:" << LID << ", tmp:" << tmp);
00480 return rnt;
00481 }
00482
00483
00484 void HNMesh3D::getFacetLIDs(int cellDim,
00485 const Array<int>& cellLID,
00486 int facetDim,
00487 Array<int>& facetLID,
00488 Array<int>& facetSign) const {
00489
00490 SUNDANCE_MSG3(verb(),"HNMesh3D::getFacetLIDs() cellDim:"<<cellDim<<" cellLID.size():"<<cellLID.size()<<" facetDim:" <<facetDim);
00491 int LID = 0 , cLID , facetOrientation ;
00492 int ptr = 0;
00493
00494 int nf = numFacets(cellDim, cellLID[0], facetDim);
00495 facetLID.resize(cellLID.size() * nf);
00496 facetSign.resize(cellLID.size() * nf);
00497
00498 for (unsigned int i = 0 ; i < (unsigned int)cellLID.size() ; i++){
00499 cLID = cellLID[i];
00500 for (int f=0; f<nf; f++, ptr++) {
00501
00502 LID = this->facetLID( cellDim, cLID, facetDim, f , facetOrientation);
00503 facetLID[ptr] = LID;
00504 facetSign[ptr] = facetOrientation;
00505 }
00506 }
00507 SUNDANCE_MSG3(verb() ,"HNMesh3D::getFacetLIDs() DONE. ");
00508 }
00509
00510
00511 const int* HNMesh3D::elemZeroFacetView(int cellLID) const {
00512 int LID = cellLeafToLIDMapping_[cellLID];
00513 SUNDANCE_MSG3(verb() , "HNMesh3D::elemZeroFacetView ");
00514 return (const int*)(&cellsPoints_[LID]);
00515 }
00516
00517
00518 int HNMesh3D::numMaxCofacets(int cellDim, int cellLID) const {
00519
00520 SUNDANCE_MSG3(verb() , "HNMesh3D::numMaxCofacets(): cellDim:"<<cellDim<<" cellLID:"<<cellLID );
00521 int rnt = -1;
00522
00523 if (cellDim==0) {
00524 int LID = vertexLeafToLIDMapping_[cellLID];
00525 int sum = 0;
00526 SUNDANCE_MSG3(verb() ," pointMaxCoF_[LID] = " << pointMaxCoF_[LID] );
00527 for (int i = 0 ; i < 8 ; i++)
00528 if ( (pointMaxCoF_[LID][i] >= 0) && ( cellLIDToLeafMapping_[pointMaxCoF_[LID][i]] >= 0) )
00529 sum++;
00530
00531 rnt = sum;
00532 }
00533 else if (cellDim==1) {
00534 int LID = edgeLeafToLIDMapping_[cellLID];
00535 int sum = 0;
00536 SUNDANCE_MSG3(verb() ," edgeMaxCoF_[LID] = " << edgeMaxCoF_[LID] );
00537 for (int i = 0 ; i < 4 ; i++)
00538 if ( (edgeMaxCoF_[LID][i] >= 0) && ( cellLIDToLeafMapping_[edgeMaxCoF_[LID][i]] >= 0) )
00539 sum++;
00540
00541 rnt = sum;
00542 }
00543 else if (cellDim==2) {
00544 int LID = faceLeafToLIDMapping_[cellLID];
00545 int sum = 0;
00546 SUNDANCE_MSG3(verb() ," faceMaxCoF_[LID] = " << faceMaxCoF_[LID] );
00547 for (int i = 0 ; i < 2 ; i++)
00548 if ( (faceMaxCoF_[LID][i] >= 0) && ( cellLIDToLeafMapping_[faceMaxCoF_[LID][i]] >= 0) )
00549 sum++;
00550
00551 rnt = sum;
00552 }
00553 SUNDANCE_MSG3(verb() ," RET = " << rnt );
00554 return rnt;
00555 }
00556
00557
00558 int HNMesh3D::maxCofacetLID(int cellDim, int cellLID,
00559 int cofacetIndex,
00560 int& facetIndex) const {
00561
00562 SUNDANCE_MSG3(verb() ,"HNMesh3D::maxCofacetLID() cellDim:"<<cellDim<<" cellLID:"<<cellLID<<" cofacetIndex:"<<cofacetIndex<< " facetIndex:"
00563 << facetIndex);
00564 int rnt =-1;
00565
00566 if (cellDim==0) {
00567
00568 int actCoFacetIndex = 0;
00569 int LID = vertexLeafToLIDMapping_[cellLID];
00570 for (int ii = 0 ; ii < 8 ; ii++){
00571
00572 if ( (pointMaxCoF_[LID][ii] >= 0) && (cellLIDToLeafMapping_[pointMaxCoF_[LID][ii]] >= 0) ){
00573 if ( actCoFacetIndex < cofacetIndex ){
00574 actCoFacetIndex++;
00575 }else{
00576 facetIndex = ii;
00577 rnt = pointMaxCoF_[LID][ii];
00578 break;
00579 }
00580 }
00581 }
00582 }
00583 else if (cellDim==1) {
00584 int maxCoFacet = 0;
00585 int LID = edgeLeafToLIDMapping_[cellLID];
00586 int orientation = edgeOrientation_[LID];
00587 SUNDANCE_MSG3(verb() ," HNMesh3D::maxCofacetLID() 1 , orientation:" << orientation );
00588
00589 int actCoFacetIndex = 0;
00590 for (int ii = 0 ; ii < 4 ; ii++){
00591
00592 if ( (edgeMaxCoF_[LID][ii] >= 0) && (cellLIDToLeafMapping_[edgeMaxCoF_[LID][ii]] >= 0) ){
00593 if ( actCoFacetIndex < cofacetIndex ){
00594 actCoFacetIndex++;
00595 }else{
00596 facetIndex = ii;
00597 maxCoFacet = edgeMaxCoF_[LID][ii];
00598 break;
00599 }
00600 }
00601 }
00602
00603 facetIndex = edge_MaxCofacetIndex[orientation][facetIndex];
00604 SUNDANCE_MSG3(verb() ,"HNMesh3D::maxCofacetLID() 1 , facetIndex:" << facetIndex << " maxCoFacet = " << maxCoFacet );
00605 rnt = ( maxCoFacet );
00606 }
00607 else if (cellDim==2) {
00608 int maxCoFacet = 0;
00609 int LID = faceLeafToLIDMapping_[cellLID];
00610 int orientation = faceOrientation_[LID];
00611 SUNDANCE_MSG3(verb() ," HNMesh3D::maxCofacetLID() 2 , orientation:" << orientation );
00612
00613 int actCoFacetIndex = 0;
00614 for (int ii = 0 ; ii < 2 ; ii++){
00615
00616 if ( (faceMaxCoF_[LID][ii] >= 0) && (cellLIDToLeafMapping_[faceMaxCoF_[LID][ii]] >= 0) ){
00617 if ( actCoFacetIndex < cofacetIndex ){
00618 actCoFacetIndex++;
00619 }else{
00620 facetIndex = ii;
00621 maxCoFacet = faceMaxCoF_[LID][ii];
00622 break;
00623 }
00624 }
00625 }
00626
00627 facetIndex = face_MaxCofacetIndex[orientation][facetIndex];
00628 SUNDANCE_MSG3(verb() ,"HNMesh3D::maxCofacetLID() 2 , facetIndex:" << facetIndex << " maxCoFacet = " << maxCoFacet );
00629 rnt = ( maxCoFacet );
00630 }
00631
00632 rnt = cellLIDToLeafMapping_[ rnt ];
00633
00634 SUNDANCE_MSG3(verb() ," RET = " << rnt << ", facetIndex:" << facetIndex);
00635 return rnt;
00636 }
00637
00638 void HNMesh3D::getCofacets(int cellDim, int cellLID,
00639 int cofacetDim, Array<int>& cofacetLIDs) const {
00640
00641 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::getCofacets() not implemented");
00642 }
00643
00644
00645 void HNMesh3D::getMaxCofacetLIDs(const Array<int>& cellLIDs,
00646 MaximalCofacetBatch& cofacets) const {
00647
00648 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::getMaxCofacetLIDs() not implemented");
00649 }
00650
00651
00652 int HNMesh3D::mapGIDToLID(int cellDim, int globalIndex) const {
00653
00654 switch (cellDim){
00655 case 0:{
00656 int ID = vertexLeafToGIDMapping_[globalIndex];
00657 SUNDANCE_MSG3(verb() , " HNMesh3D::mapGIDToLID 0 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< vertexLIDToLeafMapping_[ID]);
00658 return vertexLIDToLeafMapping_[ID];
00659 break;}
00660 case 1:{
00661 int ID = edgeLeafToGIDMapping_[globalIndex];
00662 SUNDANCE_MSG3(verb() , " HNMesh3D::mapGIDToLID 1 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< edgeLIDToLeafMapping_[ID]);
00663 return edgeLIDToLeafMapping_[ID];
00664 break;}
00665 case 2:{
00666 int ID = faceLeafToGIDMapping_[globalIndex];
00667 SUNDANCE_MSG3(verb() , " HNMesh3D::mapGIDToLID 2 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< faceLIDToLeafMapping_[ID]);
00668 return faceLIDToLeafMapping_[ID];
00669 break;}
00670 case 3:{
00671 int ID = cellLeafToGIDMapping_[globalIndex];
00672 SUNDANCE_MSG3(verb() , " HNMesh3D::mapGIDToLID 3 , globalIndex:" << globalIndex << " ,ID:" << ID << ", ret:"<< cellLIDToLeafMapping_[ID]);
00673 return cellLIDToLeafMapping_[ID];
00674 break;}
00675 }
00676 return -1;
00677 }
00678
00679
00680 bool HNMesh3D::hasGID(int cellDim, int globalIndex) const {
00681
00682
00683 return true;
00684 }
00685
00686
00687 int HNMesh3D::mapLIDToGID(int cellDim, int localIndex) const {
00688
00689 switch (cellDim){
00690 case 0:{
00691 int ID = vertexLeafToLIDMapping_[localIndex];
00692 SUNDANCE_MSG3(verb() , " HNMesh3D::mapLIDToGID 0 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< vertexGIDToLeafMapping_[ID]);
00693 return vertexGIDToLeafMapping_[ID];
00694 break;}
00695 case 1:{
00696 int ID = edgeLeafToLIDMapping_[localIndex];
00697 SUNDANCE_MSG3(verb() , " HNMesh3D::mapLIDToGID 1 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< edgeGIDToLeafMapping_[ID]);
00698 return edgeGIDToLeafMapping_[ID];
00699 break;}
00700 case 2:{
00701 int ID = faceLeafToLIDMapping_[localIndex];
00702 SUNDANCE_MSG3(verb() , " HNMesh3D::mapLIDToGID 2 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< faceGIDToLeafMapping_[ID]);
00703 return faceGIDToLeafMapping_[ID];
00704 break;}
00705 case 3:{
00706 int ID = cellLeafToLIDMapping_[localIndex];
00707 SUNDANCE_MSG3(verb() , " HNMesh3D::mapLIDToGID 3 , localIndex:" << localIndex << " ,ID:" << ID << ", ret:"<< cellGIDToLeafMapping_[ID]);
00708 return cellGIDToLeafMapping_[ID];
00709 break;}
00710 }
00711 return -1;
00712 }
00713
00714
00715 CellType HNMesh3D::cellType(int cellDim) const {
00716 switch(cellDim)
00717 {
00718 case 0: return PointCell;
00719 case 1: return LineCell;
00720 case 2: return QuadCell;
00721 case 3: return BrickCell;
00722 default:
00723 return NullCell;
00724 }
00725 }
00726
00727
00728 int HNMesh3D::label(int cellDim, int cellLID) const {
00729
00730 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::label() not implemented yet");
00731 return 0;
00732 }
00733
00734
00735 void HNMesh3D::getLabels(int cellDim, const Array<int>& cellLID,
00736 Array<int>& labels) const {
00737
00738 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::getLabels() not implemented yet");
00739 }
00740
00741 Set<int> HNMesh3D::getAllLabelsForDimension(int cellDim) const {
00742 Set<int> rtn;
00743
00744 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::getAllLabelsForDimension() not implemented yet");
00745 return rtn;
00746 }
00747
00748 void HNMesh3D::getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const {
00749
00750 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::getLIDsForLabel() not implemented yet");
00751 }
00752
00753 void HNMesh3D::setLabel(int cellDim, int cellLID, int label) {
00754
00755 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," HNMesh3D::setLabel() not implemented yet");
00756 }
00757
00758
00759 void HNMesh3D::assignIntermediateCellGIDs(int cellDim) {
00760
00761 }
00762
00763
00764 bool HNMesh3D::hasIntermediateGIDs(int dim) const {
00765
00766 return true;
00767 }
00768
00769
00770
00771 bool HNMesh3D::isElementHangingNode(int cellDim , int cellLID) const {
00772 SUNDANCE_MSG3(verb() ,"HNMesh3D::isElementHangingNode cellDim:"<<cellDim<<" LID:"<< cellLID);
00773 if (cellDim==0) {
00774 int LID = vertexLeafToLIDMapping_[cellLID];
00775 return (isPointHanging_[LID]);
00776 }
00777 else if (cellDim==1) {
00778 int LID = edgeLeafToLIDMapping_[cellLID];
00779 return (isEdgeHanging_[LID]);
00780 }
00781 else if (cellDim==2)
00782 {
00783 int LID = faceLeafToLIDMapping_[cellLID];
00784 return (isFaceHanging_[LID] );
00785
00786
00787
00788 }
00789 return false;
00790 }
00791
00792 int HNMesh3D::indexInParent(int maxCellLID) const {
00793 int ID = cellLeafToLIDMapping_[maxCellLID];
00794 int indexInPar = indexInParent_[ID];
00795 return indexInPar;
00796
00797 }
00798
00799 void HNMesh3D::returnParentFacets( int childCellLID , int dimFacets ,
00800 Array<int> &facetsLIDs , int &parentCellLIDs ) const {
00801 int LID = cellLeafToLIDMapping_[childCellLID];
00802 parentCellLIDs = parentCellLID_[LID];
00803
00804 SUNDANCE_MSG3( verb() , "HNMesh3D::returnParentFacets childCellLID:"<<childCellLID<<" dimFacets:"<<dimFacets<<
00805 " parentCellLIDs:"<< parentCellLIDs);
00806
00807
00808 if (dimFacets == 0){
00809 facetsLIDs.resize(8);
00810 for (int kuku = 0; kuku < 8 ; kuku++)
00811 facetsLIDs[kuku] = facetLID_tree( 3 , parentCellLIDs , dimFacets , kuku );
00812 }
00813 else if (dimFacets == 1){
00814 facetsLIDs.resize(12);
00815 for (int kuku = 0; kuku < 12 ; kuku++)
00816 facetsLIDs[kuku] = facetLID_tree( 3 , parentCellLIDs , dimFacets , kuku );
00817 }
00818 else if (dimFacets == 2){
00819 facetsLIDs.resize(6);
00820 for (int kuku = 0; kuku < 6 ; kuku++)
00821 facetsLIDs[kuku] = facetLID_tree( 3 , parentCellLIDs , dimFacets , kuku );
00822 }
00823
00824
00825
00826 }
00827
00828
00829 int HNMesh3D::facetLID_tree(int cellDim, int cellLID,
00830 int facetDim, int facetIndex) const{
00831 int rnt = -1;
00832 if (facetDim == 0){
00833 rnt = cellsPoints_[cellLID][facetIndex];
00834 rnt = vertexLIDToLeafMapping_[rnt];
00835
00836 } else if (facetDim == 1){
00837 rnt = cellsEdges_[cellLID][facetIndex];
00838 rnt = edgeLIDToLeafMapping_[rnt];
00839
00840 } else if (facetDim == 2){
00841 rnt = cellsFaces_[cellLID][facetIndex];
00842 rnt = faceLIDToLeafMapping_[rnt];
00843
00844 }
00845 SUNDANCE_MSG3(verb() , "HNMesh3D::facetLID_tree cellDim:"<<cellDim<<", cellLID:"<<cellLID<<", facetDim:"<<facetDim<<
00846 ", facetIndex:"<<facetIndex<<" RET = " << rnt );
00847 return rnt;
00848 }
00849
00850
00851
00852
00853 void HNMesh3D::addVertex(int vertexLID , int ownerProc , bool isHanging ,
00854 double coordx , double coordy , double coordz , const Array<int> &maxCoF){
00855
00856 if (points_.size() <= vertexLID){
00857 TEUCHOS_TEST_FOR_EXCEPTION(vertexLID != nrElem_[0] , std::logic_error ,"HNMesh3D::addVertex " <<
00858 " vertexLID:" << vertexLID << " nrElem_[0]:" << nrElem_[0] );
00859 Point pt(coordx, coordy, coordz );
00860 points_.append( pt );
00861 pointMaxCoF_.append( maxCoF );
00862 isPointHanging_.append( isHanging );
00863 elementOwner_[0].append( (short int)ownerProc );
00864 SUNDANCE_MSG3(verb() , "HNMesh3D::addVertex: " << nrElem_[0] << " P:" << pt << " , maxCoF:" << maxCoF);
00865 SUNDANCE_MSG3(verb() , " ownerProc:" << ownerProc << " , isHanging:" << isHanging);
00866 nrElem_[0]++;
00867 }
00868 }
00869
00870
00871 void HNMesh3D::addEdge(int edgeLID , int ownerProc , bool isHanging , int edgeOrientation ,
00872 const Array<int> &vertexLIDs , const Array<int> &maxCoF){
00873
00874 SUNDANCE_MSG3(verb() , "HNMesh3D -- addEdge: " << edgeLID << " nrElem_[1]: " << nrElem_[1] << " edgePoints_.size():" << edgePoints_.size() );
00875 if (edgePoints_.size() <= edgeLID ){
00876 TEUCHOS_TEST_FOR_EXCEPTION(edgeLID != nrElem_[1], std::logic_error, "HNMesh3D::addEdge edgeLID != nrElem_[1]");
00877 edgePoints_.append( vertexLIDs );
00878 edgeOrientation_.append( (short int)edgeOrientation );
00879 edgeMaxCoF_.append( maxCoF );
00880 isEdgeHanging_.append(isHanging);
00881 hangingAccessCount_.append( (short int)0);
00882 elementOwner_[1].append( (short int)ownerProc );
00883 SUNDANCE_MSG3(verb() , "HNMesh3D::addEdge: " << nrElem_[1] << " vertexLIDs:" << vertexLIDs << " , maxCoF:" << maxCoF );
00884 SUNDANCE_MSG3(verb() , " ownerProc:" << ownerProc << ", isHanging:" << isHanging << ", edgeOrientation:" << edgeOrientation );
00885 nrElem_[1]++;
00886 }
00887 }
00888
00889 void HNMesh3D::addFace(int faceLID , int ownerProc , bool isHanging , int faceOrientation ,
00890 const Array<int> &vertexLIDs , const Array<int> &edgeLIDs ,
00891 const Array<int> &maxCoF){
00892
00893
00894 if (facePoints_.size() <= faceLID ){
00895 TEUCHOS_TEST_FOR_EXCEPTION(faceLID != nrElem_[2], std::logic_error, "HNMesh3D::addFace faceLID != nrElem_[2]");
00896 facePoints_.append( vertexLIDs );
00897 faceEdges_.append( edgeLIDs );
00898 faceOrientation_.append( (short int)faceOrientation );
00899 faceMaxCoF_.append( maxCoF );
00900 isFaceHanging_.append(isHanging);
00901 SUNDANCE_MSG3(verb() , "HNMesh3D::addFace: " << nrElem_[2] << " vertexLIDs:" << vertexLIDs << " , maxCoF:" << maxCoF );
00902 SUNDANCE_MSG3(verb() , "HNMesh3D::addFace , edgeLIDs:" << edgeLIDs );
00903 SUNDANCE_MSG3(verb() , " ownerProc:" << ownerProc << ", isHanging:" << isHanging << ", faceOrientation:" << faceOrientation);
00904 elementOwner_[2].append( (short int)ownerProc );
00905 nrElem_[2]++;
00906 }
00907 }
00908
00909
00910 void HNMesh3D::addCell(int cellLID , int ownerProc ,
00911 int indexInParent , int parentCellLID , int level ,
00912 const Array<int> &faceLIDs , const Array<int> &edgeLIDs ,
00913 const Array<int> &vertexLIDs)
00914 {
00915
00916 if (cellsPoints_.size() <= cellLID ) {
00917 TEUCHOS_TEST_FOR_EXCEPTION(cellLID != nrElem_[3], std::logic_error, "HNMesh3D::cellLID cellLID != nrElem_[3]");
00918 cellsFaces_.append( faceLIDs );
00919 cellsEdges_.append( edgeLIDs );
00920 cellsPoints_.append( vertexLIDs );
00921 indexInParent_.append( indexInParent );
00922 parentCellLID_.append( parentCellLID );
00923 cellLevel_.append( level );
00924 isCellLeaf_.append( true );
00925 refineCell_.append( 0 );
00926 cellsChildren_.append( tuple(1) );
00927 elementOwner_[3].append( (short int)ownerProc );
00928
00929 isCellOut_.append( !( meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[0]]) ||
00930 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[1]]) ||
00931 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[2]]) ||
00932 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[3]]) ||
00933 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[4]]) ||
00934 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[5]]) ||
00935 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[6]]) ||
00936 meshDomain_.isInsideComputationalDomain(points_[vertexLIDs[7]]) ) );
00937 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell: " << nrElem_[3] <<
00938 " vertexLIDs:" << vertexLIDs << " edgeLIDs:" << edgeLIDs << " faceLIDs:" << faceLIDs);
00939 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p0:" << points_[vertexLIDs[0]] );
00940 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p1:" << points_[vertexLIDs[1]] );
00941 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p2:" << points_[vertexLIDs[2]] );
00942 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p3:" << points_[vertexLIDs[3]] );
00943 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p4:" << points_[vertexLIDs[4]] );
00944 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p5:" << points_[vertexLIDs[5]] );
00945 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p6:" << points_[vertexLIDs[6]] );
00946 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell p7:" << points_[vertexLIDs[7]] );
00947 SUNDANCE_MSG3(verb() , "HNMesh3D::addCell IN DOMAIN:" << isCellOut_[nrElem_[3]] );
00948 nrElem_[3]++;
00949 }
00950 }
00951
00952
00953
00954
00955
00956 void HNMesh3D::createMesh(
00957 double position_x,
00958 double position_y,
00959 double position_z,
00960 double offset_x,
00961 double offset_y,
00962 double offset_z,
00963 int resolution_x,
00964 int resolution_y,
00965 int resolution_z,
00966 const RefinementClass& refineClass ,
00967 const MeshDomainDef& meshDomain
00968 ){
00969
00970 setVerb(0);
00971
00972
00973 _pos_x = position_x; _pos_y = position_y; _pos_z = position_z;
00974 _ofs_x = offset_x; _ofs_y = offset_y; _ofs_z = offset_z;
00975 _res_x = resolution_x; _res_y = resolution_y; _res_z = resolution_z;
00976 refineClass_ = refineClass;
00977 meshDomain_ = meshDomain;
00978
00979
00980 createCoarseMesh();
00981
00982
00983 bool doRefinement = true;
00984 while (doRefinement){
00985 doRefinement = oneRefinementIteration();
00986 }
00987
00988
00989
00990 createLeafNumbering_sophisticated();
00991
00992 }
00993
00994 void HNMesh3D::updateLocalCoarseNumbering(int ix, int iy , int iz , int Nx , int Ny){
00995
00996 vInd[0] = (Nx+1)*(Ny+1)*iz + (Nx+1)*iy + ix;
00997 vInd[1] = vInd[0] + 1;
00998 vInd[2] = vInd[0] + Nx+1;
00999 vInd[3] = vInd[0] + Nx+2;
01000 vInd[4] = vInd[0] + (Nx+1)*(Ny+1);
01001 vInd[5] = vInd[4] + 1;
01002 vInd[6] = vInd[4] + Nx+1;
01003 vInd[7] = vInd[4] + Nx+2;
01004
01005 eInd[0] = (Nx*(Ny+1) + Ny*(Nx+1) + (Nx+1)*(Ny+1))*iz + (Nx+1+Nx)*iy + ix;
01006 eInd[1] = eInd[0] + Nx;
01007 eInd[3] = eInd[0] + Nx+1;
01008 eInd[5] = eInd[0] + 2*Nx+1;
01009 int ed5 = (2*Nx+1)*(iy+1)+ix;
01010 eInd[2] = eInd[5]+(Nx*(Ny+1)+Ny*(Nx+1) - ed5) + (Nx+1)*iy+ix;
01011 eInd[4] = eInd[2]+1;
01012 eInd[6] = eInd[2]+Nx+1;
01013 eInd[7] = eInd[2]+Nx+2;
01014 int ed7 = (Nx+1)*(iy+1)+ix+1;
01015 eInd[8] = eInd[7] + ( (Nx+1)*(Ny+1) - ed7) + (2*Nx+1)*iy+ix;
01016 eInd[9] = eInd[8] + Nx;
01017 eInd[10] = eInd[8] + Nx+1;
01018 eInd[11] = eInd[8] + 2*Nx+1;
01019
01020 fInd[0] = (Nx*Ny+Nx*(Ny+1)+ Ny*(Nx+1))*iz+Nx*iy+ix;
01021 fInd[1] = fInd[0] + (Nx*Ny - (Nx*iy+ix)) + (iy*(2*Nx+1) + ix);
01022 fInd[2] = fInd[1] + Nx;
01023 fInd[3] = fInd[2] + 1;
01024 fInd[4] = fInd[3] + Nx;
01025 fInd[5] = fInd[4] + ((Nx+1)*Ny + (Ny+1)*Nx - (2*Nx+1)*(iy+1) - (ix) + Nx*iy + ix);
01026 }
01027
01028 void HNMesh3D::createCoarseMesh(){
01029
01030
01031
01032
01033
01034 int nrCoarseCell = _res_x * _res_y * _res_z;
01035 int nrCoarsePoints = (_res_x+1)*(_res_y+1)*(_res_z+1);
01036 int nrCoarseEdge = (_res_x+1)*(_res_y+1)*_res_z + _res_x*(_res_y+1)*(_res_z+1) + (_res_x+1)*_res_y*(_res_z+1);
01037 int nrCoarseFace = (_res_x+1)*_res_y*_res_z + _res_x*(_res_y+1)*_res_z+ + _res_x*_res_y*(_res_z+1);
01038 Array<int> coarseCellOwner( nrCoarseCell , -1 );
01039 Array<int> coarsePointOwner( nrCoarsePoints , -1 );
01040 Array<int> coarseEdgeOwner( nrCoarseEdge , -1 );
01041 Array<int> coarseFaceOwner( nrCoarseFace , -1 );
01042 Array<int> coarseCellLID( nrCoarseCell , -1 );
01043 Array<int> coarsePointLID( nrCoarsePoints , -1 );
01044 Array<int> coarseEdgeLID( nrCoarseEdge , -1 );
01045 Array<int> coarseFaceLID( nrCoarseFace , -1 );
01046 Array<int> coarseCellsLoad( nrCoarseCell , 0 );
01047 int totalLoad = 0;
01048
01049 SUNDANCE_MSG3(verb() , "HNMesh3D::createMesh nrCoarseCell:" << nrCoarseCell << " nrCoarsePoints:" << nrCoarsePoints
01050 << " nrCoarseEdge:" << nrCoarseEdge << " nrCoarseFace:" << nrCoarseFace << " nrProc_:" << nrProc_ << " myRank_:" << myRank_);
01051 TEUCHOS_TEST_FOR_EXCEPTION( nrCoarseCell < nrProc_ , std::logic_error," HNMesh3D::createMesh nrCoarseCell < nrProc_ ");
01052
01053
01054
01055 double h[3];
01056 Array<int> ind(3);
01057 h[0] = _ofs_x/(double)_res_x; h[1] = _ofs_y/(double)_res_y; h[2] = _ofs_z/(double)_res_z;
01058 Point pos(h[0],h[1],h[2]);
01059 Point res(h[0],h[1],h[2]);
01060
01061 for (int i=0; i < nrCoarseCell; i++){
01062
01063 ind[0] = ((i % (_res_x*_res_y)) % _res_x);
01064 ind[1] = ((i % (_res_x*_res_y)) / _res_x);
01065 ind[2] = (i / (_res_x*_res_y));
01066 pos[0] = _pos_x + (double)(ind[0])*h[0] + 0.5*h[0];
01067 pos[1] = _pos_y + (double)(ind[1])*h[1] + 0.5*h[1];
01068 pos[2] = _pos_z + (double)(ind[2])*h[2] + 0.5*h[2];
01069
01070 coarseCellsLoad[i] = refineClass_.estimateRefinementLevel( pos , res );
01071 totalLoad += coarseCellsLoad[i];
01072 }
01073
01074
01075 double loadPerProc = (double)totalLoad / (double)nrProc_;
01076 int actualProc=0;
01077 totalLoad = 0;
01078
01079 SUNDANCE_MSG3(verb() , "Processor asign, loadPerProc:" << loadPerProc );
01080 double diff_load = 0.0;
01081 for (int i=0; i < nrCoarseCell; i++){
01082 ind[0] = ((i % (_res_x*_res_y)) % _res_x);
01083 ind[1] = ((i % (_res_x*_res_y)) / _res_x);
01084 ind[2] = (i / (_res_x*_res_y));
01085
01086 updateLocalCoarseNumbering( ind[0] , ind[1] , ind[2] , _res_x , _res_y);
01087
01088
01089
01090 for (int jj = 0 ; jj < 8 ; jj++){
01091 if (coarsePointOwner[vInd[jj]] < 0){
01092 coarsePointOwner[vInd[jj]] = actualProc;
01093 SUNDANCE_MSG3(verb() , "Vertex CPU assign " << vInd[jj] << " ->" << actualProc );
01094 }
01095 }
01096
01097 for (int jj = 0 ; jj < 12 ; jj++){
01098 if (coarseEdgeOwner[eInd[jj]] < 0){
01099 coarseEdgeOwner[eInd[jj]] = actualProc;
01100 SUNDANCE_MSG3(verb() , "Edge CPU assign " << eInd[jj] << " ->" << actualProc );
01101 }
01102 }
01103
01104 for (int jj = 0 ; jj < 6 ; jj++){
01105 if (coarseFaceOwner[fInd[jj]] < 0){
01106 coarseFaceOwner[fInd[jj]] = actualProc;
01107 SUNDANCE_MSG3(verb() , "Face CPU assign " << fInd[jj] << " ->" << actualProc );
01108 }
01109 }
01110
01111 coarseCellOwner[i] = actualProc;
01112 totalLoad += coarseCellsLoad[i];
01113 SUNDANCE_MSG3(verb() , "Cell CPU assign " << i << " ->" << actualProc <<
01114 ", totalLoad:" << totalLoad << " loadPerProc:" << loadPerProc);
01115
01116 if (((double)totalLoad >= (loadPerProc - 1e-8 - diff_load)) && ( actualProc < nrProc_-1 )){
01117 SUNDANCE_MSG3(verb() , "Increase CPU , totalLoad:" << totalLoad << " loadPerProc:" << loadPerProc );
01118
01119 diff_load = totalLoad - loadPerProc;
01120 actualProc = actualProc + 1;
01121 totalLoad = 0;
01122 }
01123 }
01124
01125
01126 SUNDANCE_MSG3(verb() , " Process Cells:" << nrCoarseCell );
01127 for (int i=0; i < nrCoarseCell; i++)
01128 {
01129 ind[0] = ((i % (_res_x*_res_y)) % _res_x);
01130 ind[1] = ((i % (_res_x*_res_y)) / _res_x);
01131 ind[2] = (i / (_res_x*_res_y));
01132
01133 updateLocalCoarseNumbering( ind[0] , ind[1] , ind[2] , _res_x , _res_y );
01134 pos[0] = _pos_x + (double)(ind[0])*h[0];
01135 pos[1] = _pos_y + (double)(ind[1])*h[1];
01136 pos[2] = _pos_z + (double)(ind[2])*h[2];
01137 SUNDANCE_MSG3(verb() , "PCell ID:" << i <<" pos:"<<pos<<" _res_x:"<<_res_x<<" _res_y:"<<_res_y<<" _res_z:"<<_res_z);
01138 SUNDANCE_MSG3(verb() , "PCell, actual index" << ind );
01139
01140 int cellLID = coarseCellLID[i];
01141 Array<int> vLID(8,-1) , vertexMaxCoF(8,-1) ;
01142 Array<int> edgeLID(12,-1) , edgeVert(2,-1) , edgeMaxCoef(4,-1);
01143 Array<int> faceLID(6,-1) , faceVert(4,-1) , faceEdge(4,-1) , faceMaxCoef(2,-1);
01144
01145
01146
01147 if (coarseCellLID[i] < 0 ){
01148 coarseCellLID[i] = nrElem_[3];
01149 cellLID = coarseCellLID[i];
01150 }
01151
01152 for (int jj = 0 ; jj < 8 ; jj++){
01153 if (coarsePointLID[vInd[jj]] < 0){
01154 coarsePointLID[vInd[jj]] = nrElem_[0];
01155 }
01156 vLID[jj] = coarsePointLID[vInd[jj]];
01157
01158
01159
01160 addVertex( vLID[jj] , coarsePointOwner[vInd[jj]] , false ,
01161 pos[0] + ((double)offs_Points_x_[jj])*h[0] , pos[1] + ((double)offs_Points_y_[jj])*h[1] ,
01162 pos[2] + ((double)offs_Points_z_[jj])*h[2] ,
01163 vertexMaxCoF );
01164 }
01165
01166 for (int jj = 0 ; jj < 12 ; jj++){
01167 if (coarseEdgeLID[eInd[jj]] < 0){
01168 coarseEdgeLID[eInd[jj]] = nrElem_[1];
01169 }
01170 edgeLID[jj] = coarseEdgeLID[eInd[jj]];
01171 edgeVert[0] = vLID[edge_Points_localIndex[jj][0]];
01172 edgeVert[1] = vLID[edge_Points_localIndex[jj][1]];
01173 SUNDANCE_MSG3(verb() , "Edge local Index:" << eInd[jj] );
01174
01175 addEdge( edgeLID[jj] , coarseEdgeOwner[eInd[jj]] , false , edge_Orientation[jj] , edgeVert , edgeMaxCoef );
01176 }
01177
01178 for (int jj = 0 ; jj < 6 ; jj++){
01179 if (coarseFaceLID[fInd[jj]] < 0){
01180 coarseFaceLID[fInd[jj]] = nrElem_[2];
01181 }
01182 faceLID[jj] = coarseFaceLID[fInd[jj]];
01183
01184 faceVert[0] = vLID[face_Points_localIndex[jj][0]];
01185 faceVert[1] = vLID[face_Points_localIndex[jj][1]];
01186 faceVert[2] = vLID[face_Points_localIndex[jj][2]];
01187 faceVert[3] = vLID[face_Points_localIndex[jj][3]];
01188 faceEdge[0] = edgeLID[face_Edges_localIndex[jj][0]];
01189 faceEdge[1] = edgeLID[face_Edges_localIndex[jj][1]];
01190 faceEdge[2] = edgeLID[face_Edges_localIndex[jj][2]];
01191 faceEdge[3] = edgeLID[face_Edges_localIndex[jj][3]];
01192
01193 addFace( faceLID[jj] , coarseFaceOwner[fInd[jj]] , false , face_Orientation[jj] ,
01194 faceVert , faceEdge , faceMaxCoef );
01195 }
01196
01197 addCell( cellLID , coarseCellOwner[i] , 0 , cellLID , 0 , faceLID , edgeLID , vLID);
01198 }
01199
01200
01201 SUNDANCE_MSG3(verb() , " Process maxCofacets:" << nrCoarseCell );
01202 for (int i=0; i < nrCoarseCell; i++){
01203 Array<int> vLID(8,-1) , eLID(12,-1) ,fLID(6,-1) ;
01204 ind[0] = ((i % (_res_x*_res_y)) % _res_x);
01205 ind[1] = ((i % (_res_x*_res_y)) / _res_x);
01206 ind[2] = (i / (_res_x*_res_y));
01207
01208 updateLocalCoarseNumbering( ind[0] , ind[1] , ind[2] , _res_x , _res_y );
01209 int cellLID = coarseCellLID[i];
01210 SUNDANCE_MSG3(verb() , " MaxCoFs in Cell cellLID:" << cellLID );
01211
01212
01213
01214 for (int jj = 0 ; jj < 8 ; jj++){
01215 vLID[jj] = coarsePointLID[vInd[jj]];
01216
01217 pointMaxCoF_[vLID[jj]][jj] = cellLID;
01218 SUNDANCE_MSG3(verb() , "Vertex MaxCoFacet vLID[jj]:" << vLID[jj] << " jj:" << jj << " cellLID:" << cellLID );
01219 }
01220
01221 for (int jj = 0 ; jj < 12 ; jj++){
01222 eLID[jj] = coarseEdgeLID[eInd[jj]];
01223
01224 edgeMaxCoF_[eLID[jj]][edge_MaxCof[jj]] = cellLID;
01225 SUNDANCE_MSG3(verb() , "Edge MaxCoFacet eLID[jj]:" << eLID[jj] << " edge_MaxCof[jj]:" << edge_MaxCof[jj] << " cellLID:" << cellLID );
01226 }
01227
01228 for (int jj = 0 ; jj < 6 ; jj++){
01229 fLID[jj] = coarseFaceLID[fInd[jj]];
01230
01231 faceMaxCoF_[fLID[jj]][face_MaxCof[jj]] = cellLID;
01232 SUNDANCE_MSG3(verb() , "Face MaxCoFacet fLID[jj]:" << fLID[jj] << " face_MaxCof[jj]:" << face_MaxCof[jj] << " cellLID:" << cellLID );
01233 }
01234 }
01235
01236 }
01237
01238
01239 bool HNMesh3D::oneRefinementIteration(){
01240
01241 int nrActualCell = nrElem_[3];
01242 bool rtn = false;
01243 SUNDANCE_MSG3(verb() , " HNMesh3D::oneRefinementIteration, start one refinement iteration cycle ");
01244
01245 for (int i=0 ; i < nrActualCell ; i++){
01246
01247 SUNDANCE_MSG3(verb() , " Test cell " << i << ", elementOwner_[3][i]:" << elementOwner_[3][i] <<
01248 ", isCellLeaf_[i]:" << isCellLeaf_[i] << ", out:" << (!isCellOut_[i]));
01249
01250 if ( (isCellLeaf_[i] == true) )
01251
01252 {
01253
01254 Array<int>& cellsEdges = cellsEdges_[i];
01255 bool doRefined = true;
01256 bool refFunction = false;
01257 for (int jj = 0 ; jj < cellsEdges.size() ; jj++){
01258 SUNDANCE_MSG3(verb() , " eLID: " << cellsEdges[jj] << ", isHanging:" << isEdgeHanging_[cellsEdges[jj]]);
01259 doRefined = ( doRefined && ((isEdgeHanging_[cellsEdges[jj]]) == false) );
01260 }
01261
01262 SUNDANCE_MSG3(verb() , " is possible to refine cell nr: " << i << ", doRefined:" << doRefined);
01263
01264 Array<int>& cellsVertex = cellsPoints_[i];
01265 Point h = points_[cellsVertex[7]] - points_[cellsVertex[0]];
01266 Point p2 = points_[cellsVertex[0]] + 0.5*h;
01267 SUNDANCE_MSG3(verb() , " points_[cellsVertex[0]]: " << points_[cellsVertex[0]]);
01268 SUNDANCE_MSG3(verb() , " points_[cellsVertex[7]]: " << points_[cellsVertex[7]]);
01269 refFunction = ( (refineCell_[i] == 1) || refineClass_.refine( cellLevel_[i] , p2 , h ) );
01270
01271
01272
01273
01274 SUNDANCE_MSG3(verb() , " execute refinement on cell nr: " << i << ", refFunction:" << refFunction);
01275 if (doRefined && refFunction)
01276 {
01277
01278 refineCell(i);
01279 rtn = true;
01280 refineCell_[i] = 0;
01281 }
01282 else
01283 {
01284
01285
01286 if (refFunction) {
01287
01288 for (int jj = 0 ; jj < cellsEdges.size() ; jj++)
01289 if (isEdgeHanging_[cellsEdges[jj]]){
01290
01291 int cLevel = cellLevel_[i];
01292 int parentID = parentCellLID_[i];
01293 int parentCEdge = cellsEdges_[parentID][jj];
01294 int refCell = -1;
01295
01296 for (int coF = 0 ; coF < 4 ; coF++){
01297 refCell = edgeMaxCoF_[parentCEdge][coF];
01298
01299 if ( (refCell >= 0) && (cellLevel_[refCell] < cLevel) && (isCellLeaf_[refCell])){
01300
01301
01302 refineCell_[refCell] = 1;
01303 rtn = true;
01304
01305 }
01306 }
01307 }
01308 }
01309 }
01310 }
01311 }
01312
01313 SUNDANCE_MSG3(verb() , " HNMesh3D::oneRefinementIteration DONE with one refinement iteration");
01314
01315 return rtn;
01316 }
01317
01318
01319 void HNMesh3D::refineCell(int cellLID){
01320
01321
01322 int cellOwner = elementOwner_[3][cellLID];
01323 Point p = points_[cellsPoints_[cellLID][0]];
01324 Point h = points_[cellsPoints_[cellLID][7]] - points_[cellsPoints_[cellLID][0]];
01325
01326
01327 Array<int> vertexLIDs(64,-1);
01328
01329 vertexLIDs[0] = cellsPoints_[cellLID][0]; vertexLIDs[3] = cellsPoints_[cellLID][1];
01330 vertexLIDs[12] = cellsPoints_[cellLID][2]; vertexLIDs[15] = cellsPoints_[cellLID][3];
01331 vertexLIDs[48] = cellsPoints_[cellLID][4]; vertexLIDs[51] = cellsPoints_[cellLID][5];
01332 vertexLIDs[60] = cellsPoints_[cellLID][6]; vertexLIDs[63] = cellsPoints_[cellLID][7];
01333 Array<int> edgeLIDs(144,-1);
01334 Array<int> faceLIDs(108,-1);
01335
01336 Array<int> cellLIDs(27,-1);
01337 for (int c = 0; c < 27 ; cellLIDs[c] = nrElem_[3] + c , c++);
01338
01339
01340 isCellLeaf_[cellLID] = false;
01341
01342 SUNDANCE_MSG3(verb() , " ========== HNMesh3D::refineCell, cellLID:" << cellLID << " ==========");
01343
01344 Array<int> ind(3,-1);
01345 for (int childCell = 0 ; childCell < 27 ; childCell++){
01346 Array<int> currCellV(8,-1);
01347 Array<int> currCellE(12,-1);
01348 Array<int> currCellF(6,-1);
01349
01350
01351 ind[0] = ( (childCell % 9) % 3);
01352 ind[1] = ( (childCell % 9) / 3);
01353 ind[2] = ( childCell / 9);
01354
01355 updateLocalCoarseNumbering( ind[0] , ind[1] , ind[2] , 3 , 3);
01356
01357
01358 for (int v = 0; v < 8 ; v++){
01359 int vertexOwner = cellOwner;
01360 bool vertexIsHanging = false;
01361 bool useEdge = true;
01362 int pIndex = 0 , pOffset = 0 , pID = 0;
01363
01364 if ( vertexToParentEdge[vInd[v]] >= 0){
01365 if ( vertexToParentEdge[vInd[v]] > 19){
01366 useEdge = false;
01367 pIndex = vertexToParentEdge[vInd[v]] - 20;
01368 pID = cellsFaces_[cellLID][ pIndex ];
01369 pOffset = vertexInParentIndex[vInd[v]] - 20;
01370 SUNDANCE_MSG3(verb() , "Vertex -> Face vInd[v]:" << vInd[v] << " pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01371 SUNDANCE_MSG3(verb() , " cellsFaces:" << cellsFaces_[cellLID] << " elementOwner_[2][pID]:" << elementOwner_[2][pID]);
01372 vertexOwner = elementOwner_[2][pID];
01373 vertexIsHanging = ( numMaxCofacets_ID( 2 , pID ) > 1);
01374 } else {
01375 pIndex = vertexToParentEdge[vInd[v]];
01376 pID = cellsEdges_[cellLID][ pIndex ];
01377 pOffset = vertexInParentIndex[vInd[v]];
01378 SUNDANCE_MSG3(verb() , "Vertex -> Edge vInd[v]:" << vInd[v] << " pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01379 SUNDANCE_MSG3(verb() , " cellsEdges:" << cellsEdges_[cellLID] << " elementOwner_[1][pID]:" << elementOwner_[1][pID]);
01380 vertexOwner = elementOwner_[1][pID];
01381 vertexIsHanging = ( numMaxCofacets_ID( 1 , pID ) > 1);
01382 }
01383 }
01384
01385 if ( (vertexLIDs[vInd[v]] < 0) && (vertexToParentEdge[vInd[v]] >= 0) ){
01386 vertexLIDs[vInd[v]] = getHangingElement(0, useEdge, pID , pOffset );
01387 }
01388 SUNDANCE_MSG3(verb() , "Vertex vInd[v]:" << vInd[v] << " pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01389
01390 if ( vertexLIDs[vInd[v]] < 0){
01391 Array<int> maxCoF(8,-1);
01392 addVertex( nrElem_[0] , vertexOwner , vertexIsHanging ,
01393 p[0] + h[0]*vertex_X[vInd[v]] , p[1] + h[1]*vertex_Y[vInd[v]] , p[2] + h[2]*vertex_Z[vInd[v]] , maxCoF);
01394 vertexLIDs[vInd[v]] = nrElem_[0] - 1;
01395
01396 if (vertexIsHanging)
01397 addHangingElement(0 , vertexLIDs[vInd[v]] , useEdge , pID , pOffset);
01398 }
01399 currCellV[v] = vertexLIDs[vInd[v]];
01400
01401 SUNDANCE_MSG3(verb() , " vertexLIDs[vInd[v]]: " << vertexLIDs[vInd[v]] );
01402 pointMaxCoF_[ vertexLIDs[vInd[v]] ][v] = cellLIDs[childCell];
01403
01404 }
01405
01406
01407 for (int e = 0; e < 12 ; e++){
01408 int edgeOwner = cellOwner;
01409 bool edgeIsHanging = false;
01410 bool useEdge = true;
01411 int pIndex = 0 , pOffset = 0 , pID = 0;
01412
01413 if ( edgeToParentEdge[eInd[e]] >= 0){
01414 if ( edgeToParentEdge[eInd[e]] > 19){
01415 useEdge = false;
01416 pIndex = edgeToParentEdge[eInd[e]] - 20;
01417 pID = cellsFaces_[cellLID][ pIndex ];
01418 pOffset = edgeInParentIndex[eInd[e]] - 20;
01419 SUNDANCE_MSG3(verb() , "Edge -> Face eInd[e]:" << eInd[e] <<" pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01420 SUNDANCE_MSG3(verb() , " cellsFaces:" << cellsFaces_[cellLID] );
01421 edgeOwner = elementOwner_[2][pID];
01422 edgeIsHanging = ( numMaxCofacets_ID( 2 , pID ) > 1);
01423 } else {
01424 pIndex = edgeToParentEdge[eInd[e]];
01425 pID = cellsEdges_[cellLID][ pIndex ];
01426 pOffset = edgeInParentIndex[eInd[e]];
01427 SUNDANCE_MSG3(verb() , "Edge -> Edge eInd[e]:" << eInd[e] <<" pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01428
01429 edgeOwner = elementOwner_[1][pID];
01430 edgeIsHanging = ( numMaxCofacets_ID( 1 , pID ) > 1);
01431 }
01432 }
01433
01434 if ( (edgeLIDs[eInd[e]] < 0) && ( edgeToParentEdge[eInd[e]] >= 0) ){
01435 edgeLIDs[eInd[e]] = getHangingElement( 1 , useEdge, pID , pOffset );
01436 }
01437 SUNDANCE_MSG3(verb() , "Edge eInd[e]:" << eInd[e] <<" pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01438 if ( edgeLIDs[eInd[e]] < 0){
01439 Array<int> maxCoF(4,-1);
01440 Array<int> edgeVertexs(2,-1);
01441 edgeVertexs[0] = currCellV[edge_Points_localIndex[e][0]]; edgeVertexs[1] = currCellV[edge_Points_localIndex[e][1]];
01442
01443 addEdge( nrElem_[1] , edgeOwner , edgeIsHanging , edge_Orientation[e] , edgeVertexs , maxCoF );
01444 edgeLIDs[eInd[e]] = nrElem_[1] - 1;
01445
01446 if (edgeIsHanging)
01447 addHangingElement(1 , edgeLIDs[eInd[e]] , useEdge , pID , pOffset);
01448 }
01449 currCellE[e] = edgeLIDs[eInd[e]];
01450
01451 edgeMaxCoF_[edgeLIDs[eInd[e]]][edge_MaxCof[e]] = cellLIDs[childCell];
01452 SUNDANCE_MSG3(verb() , " edgeLIDs[eInd[e]]: " << edgeLIDs[eInd[e]] );
01453
01454 }
01455
01456
01457 for (int f = 0; f < 6 ; f++){
01458 int faceOwner = cellOwner;
01459 bool faceIsHanging = false;
01460 bool useEdge = false;
01461 int pIndex = 0 , pOffset = 0 , pID = 0;
01462
01463 if ( faceToParentFace[fInd[f]] >= 0 ){
01464 pIndex = faceToParentFace[fInd[f]];
01465 pID = cellsFaces_[cellLID][ pIndex ];
01466 pOffset = faceInParentIndex[fInd[f]];
01467 SUNDANCE_MSG3(verb() , " cellsFaces:" << cellsFaces_[cellLID] );
01468 faceOwner = elementOwner_[2][pID];
01469 faceIsHanging = ( numMaxCofacets_ID( 2 , pID ) > 1);
01470 }
01471
01472 if ( (faceLIDs[fInd[f]] < 0) && ( faceToParentFace[fInd[f]] >= 0 )){
01473 faceLIDs[fInd[f]] = getHangingElement( 2 , useEdge, pID , pOffset );
01474 }
01475 SUNDANCE_MSG3(verb() , "Face -> Face pIndex: " << pIndex << " pID:" << pID << " pOffset:" << pOffset);
01476
01477 if ( faceLIDs[fInd[f]] < 0){
01478 Array<int> maxCoF(2,-1);
01479 Array<int> faceVertexs(4,-1);
01480 faceVertexs[0] = currCellV[face_Points_localIndex[f][0]]; faceVertexs[1] = currCellV[face_Points_localIndex[f][1]];
01481 faceVertexs[2] = currCellV[face_Points_localIndex[f][2]]; faceVertexs[3] = currCellV[face_Points_localIndex[f][3]];
01482 Array<int> faceEdges(4,-1);
01483 faceEdges[0] = currCellE[face_Edges_localIndex[f][0]]; faceEdges[1] = currCellE[face_Edges_localIndex[f][1]];
01484 faceEdges[2] = currCellE[face_Edges_localIndex[f][2]]; faceEdges[3] = currCellE[face_Edges_localIndex[f][3]];
01485
01486 addFace( nrElem_[2] , faceOwner , faceIsHanging , face_Orientation[f] , faceVertexs , faceEdges , maxCoF );
01487 faceLIDs[fInd[f]] = nrElem_[2] - 1;
01488
01489 if (faceIsHanging)
01490 addHangingElement(2 , faceLIDs[fInd[f]] , useEdge , pID , pOffset);
01491 }
01492
01493 currCellF[f] = faceLIDs[fInd[f]];
01494 faceMaxCoF_[ currCellF[f] ][face_MaxCof[f]] = cellLIDs[childCell];
01495 SUNDANCE_MSG3(verb() , " faceLIDs[fInd[f]]: " << faceLIDs[fInd[f]] );
01496
01497 }
01498
01499
01500 addCell( nrElem_[3] , cellOwner , childCell , cellLID ,
01501 cellLevel_[cellLID] + 1 , currCellF , currCellE , currCellV );
01502
01503 cellsChildren_[cellLID].append(cellLIDs[childCell]);
01504 }
01505 }
01506
01507 void HNMesh3D::addHangingElement(int cellDim, int cellID ,bool useEdge,
01508 int parentID , int parentOffset){
01509
01510 if (useEdge){
01511 if (!edgeHangingElmStore_.containsKey(parentID)){
01512 Array<int> hnInfo(6,-1);
01513 hnInfo[5] = numMaxCofacets_ID( 1 , parentID );
01514 hangingAccessCount_[parentID] = hnInfo[5];
01515 edgeHangingElmStore_.put( parentID , hnInfo );
01516 }
01517 const Array<int>& hnInfo = edgeHangingElmStore_.get(parentID);
01518 Array<int> tmp_vect(hnInfo.size());
01519 for (int ii = 0; ii < hnInfo.size() ; ii++) tmp_vect[ii] = hnInfo[ii];
01520 int realOffset = 0;
01521
01522 if (cellDim == 0){ realOffset = 0;}
01523
01524 else{ realOffset = 2; }
01525
01526 SUNDANCE_MSG3(verb() ,"HNMesh3D::addHangingElement EDGE realOffset:" << realOffset << " parentOffset:" << parentOffset);
01527 SUNDANCE_MSG3(verb() ,"HNMesh3D::addHangingElement EDGE i:" << realOffset+parentOffset <<" hnInfo:" << hnInfo);
01528 tmp_vect[realOffset+parentOffset] = cellID;
01529
01530 edgeHangingElmStore_.remove( parentID );
01531 edgeHangingElmStore_.put( parentID , tmp_vect );
01532 }
01533
01534 else{
01535 if (!faceHangingElmStore_.containsKey(parentID)){
01536 Array<int> hnInfo(25,-1);
01537 faceHangingElmStore_.put( parentID , hnInfo );
01538 }
01539 const Array<int>& hnInfo = faceHangingElmStore_.get(parentID);
01540 Array<int> tmp_vect(hnInfo.size());
01541 for (int ii = 0; ii < hnInfo.size() ; ii++) tmp_vect[ii] = hnInfo[ii];
01542 int realOffset = 0;
01543
01544 if (cellDim == 0) { realOffset = 0; }
01545
01546 else if (cellDim == 1) { realOffset = 4; }
01547
01548 else { realOffset = 16; }
01549
01550 SUNDANCE_MSG3(verb() ,"HNMesh3D::addHangingElement FACE realOffset:" << realOffset << " parentOffset:" << parentOffset);
01551 SUNDANCE_MSG3(verb() ,"HNMesh3D::addHangingElement FACE i:" << realOffset+parentOffset <<" hnInfo:" << hnInfo);
01552 tmp_vect[realOffset+parentOffset] = cellID;
01553
01554 faceHangingElmStore_.remove( parentID );
01555 faceHangingElmStore_.put( parentID , tmp_vect );
01556 }
01557 }
01558
01559 int HNMesh3D::getHangingElement(int cellDim, bool useEdge, int parentID , int parentOffset){
01560 int rtn = -1;
01561 int realOffset = 0;
01562 if (useEdge){
01563
01564 if (cellDim == 0){ realOffset = 0;}
01565
01566 else{ realOffset = 2; }
01567 if (edgeHangingElmStore_.containsKey(parentID)){
01568 const Array<int>& hnInfo = edgeHangingElmStore_.get(parentID);
01569 rtn = hnInfo[realOffset+parentOffset];
01570 SUNDANCE_MSG3(verb() ,"HNMesh3D::getHangingElement EDGE rtn = " << rtn << " realOffset:" << realOffset << " parentOffset:" << parentOffset);
01571 SUNDANCE_MSG3(verb() ,"HNMesh3D::getHangingElement EDGE i:" << realOffset+parentOffset <<" hnInfo:" << hnInfo);
01572
01573 if ( rtn >= 0){
01574 if ((cellDim == 0) && (parentOffset == 1)) {
01575 if (hangingAccessCount_[parentID] <= 2){
01576 isPointHanging_[hnInfo[0]] = false; isPointHanging_[hnInfo[1]] = false;
01577 } else {
01578
01579 }
01580 }
01581 if ((cellDim == 1) && (parentOffset == 2)) {
01582 if (hangingAccessCount_[parentID] <= 2){
01583 isEdgeHanging_[hnInfo[2]] = false; isEdgeHanging_[hnInfo[3]] = false;
01584 isEdgeHanging_[hnInfo[4]] = false;
01585 } else {
01586 hangingAccessCount_[parentID]--;
01587 }
01588 }
01589 }
01590 }
01591 }
01592 else
01593 {
01594
01595 if (cellDim == 0) { realOffset = 0; }
01596
01597 else if (cellDim == 1) { realOffset = 4; }
01598
01599 else { realOffset = 16; }
01600 if (faceHangingElmStore_.containsKey(parentID)){
01601 const Array<int>& hnInfo = faceHangingElmStore_.get(parentID);
01602 rtn = hnInfo[realOffset+parentOffset];
01603 SUNDANCE_MSG3(verb() ,"HNMesh3D::getHangingElement FACE rtn = " << rtn << " realOffset:" << realOffset << " parentOffset:" << parentOffset);
01604 SUNDANCE_MSG3(verb() ,"HNMesh3D::getHangingElement FACE i:" << realOffset+parentOffset <<" hnInfo:" << hnInfo);
01605
01606 if ( rtn >= 0){
01607 if (cellDim == 0) { isPointHanging_[rtn] = false; }
01608 else if (cellDim == 1) { isEdgeHanging_[rtn] = false; }
01609 else { isFaceHanging_[rtn] = false; }
01610 }
01611 }
01612 }
01613 return rtn;
01614 }
01615
01616 int HNMesh3D::numMaxCofacets_ID(int cellDim, int cellID)
01617 {
01618 int rtn = -1;
01619 if (cellDim==1) {
01620 int LID = cellID;
01621 int sum = 0;
01622 SUNDANCE_MSG3(verb() ,"HNMesh3D::numMaxCofacets_ID ID:" << LID << " edgeMaxCoF_[LID] = " << edgeMaxCoF_[LID] );
01623 for (int i = 0 ; i < 4 ; i++){
01624 if (edgeMaxCoF_[LID][i] >= 0)
01625 sum++;
01626 }
01627
01628 rtn = sum;
01629 }
01630 else if (cellDim==2) {
01631 int LID = cellID;
01632 int sum = 0;
01633 SUNDANCE_MSG3(verb() ,"HNMesh3D::numMaxCofacets_ID ID:" << LID << " faceMaxCoF_[LID] = " << faceMaxCoF_[LID] );
01634 for (int i = 0 ; i < 2 ; i++){
01635 if (faceMaxCoF_[LID][i] >= 0)
01636 sum++;
01637 }
01638
01639 rtn = sum;
01640 }
01641 return rtn;
01642 }
01643
01644
01645 void HNMesh3D::createLeafNumbering(){
01646
01647
01648
01649
01650
01651
01652
01653 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering nrPoint:" << nrElem_[0] << " , nrEdge:"
01654 << nrElem_[1] << ", nrFace:" << nrElem_[2] << ", nrCell:" << nrElem_[3]);
01655
01656 vertexGIDToLeafMapping_.resize(nrElem_[0],-1);
01657 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexGIDToLeafMapping_[dd] = -1;
01658 vertexLeafToGIDMapping_.resize(nrElem_[0],-1);
01659 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToGIDMapping_[dd] = -1;
01660
01661 edgeGIDToLeafMapping_.resize(nrElem_[1],-1);
01662 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeGIDToLeafMapping_[dd] = -1;
01663 edgeLeafToGIDMapping_.resize(nrElem_[1],-1);
01664 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToGIDMapping_[dd] = -1;
01665
01666 faceGIDToLeafMapping_.resize(nrElem_[2],-1);
01667 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceGIDToLeafMapping_[dd] = -1;
01668 faceLeafToGIDMapping_.resize(nrElem_[2],-1);
01669 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLeafToGIDMapping_[dd] = -1;
01670
01671 cellGIDToLeafMapping_.resize(nrElem_[3],-1);
01672 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellGIDToLeafMapping_[dd] = -1;
01673 cellLeafToGIDMapping_.resize(nrElem_[3],-1);
01674 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLeafToGIDMapping_[dd] = -1;
01675
01676 nrVertexLeafGID_ = 0; nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0; nrFaceLeafGID_ = 0;
01677
01678 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrEdgeLeafLID_ = 0;
01679 vertexLIDToLeafMapping_.resize(nrElem_[0],-1);
01680 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLIDToLeafMapping_[dd] = -1;
01681 vertexLeafToLIDMapping_.resize(nrElem_[0],-1);
01682 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToLIDMapping_[dd] = -1;
01683
01684 edgeLIDToLeafMapping_.resize(nrElem_[1],-1);
01685 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLIDToLeafMapping_[dd] = -1;
01686 edgeLeafToLIDMapping_.resize(nrElem_[1],-1);
01687 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToLIDMapping_[dd] = -1;
01688
01689 faceLIDToLeafMapping_.resize(nrElem_[2],-1);
01690 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLIDToLeafMapping_[dd] = -1;
01691 faceLeafToLIDMapping_.resize(nrElem_[2],-1);
01692 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLeafToLIDMapping_[dd] = -1;
01693
01694 cellLIDToLeafMapping_.resize(nrElem_[3],-1);
01695 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLIDToLeafMapping_[dd] = -1;
01696 cellLeafToLIDMapping_.resize(nrElem_[3],-1);
01697 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLeafToLIDMapping_[dd] = -1;
01698
01699 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering , start assigning leaf numbers");
01700
01701
01702
01703 Array<bool> hasCellLID(nrElem_[3],false);
01704
01705 for (int ind = 0 ; ind < nrElem_[3] ; ind++){
01706 Array<int>& vertexIDs = cellsPoints_[ind];
01707 hasCellLID[ind] = false;
01708 for (int v = 0 ; v < 8 ; v++){
01709 Array<int>& maxCoFacet = pointMaxCoF_[vertexIDs[v]];
01710 hasCellLID[ind] = ( hasCellLID[ind]
01711 || ( (maxCoFacet[0] >= 0) && (elementOwner_[3][maxCoFacet[0]] == myRank_) )
01712 || ( (maxCoFacet[1] >= 0) && (elementOwner_[3][maxCoFacet[1]] == myRank_) )
01713 || ( (maxCoFacet[2] >= 0) && (elementOwner_[3][maxCoFacet[2]] == myRank_) )
01714 || ( (maxCoFacet[3] >= 0) && (elementOwner_[3][maxCoFacet[3]] == myRank_) )
01715 || ( (maxCoFacet[4] >= 0) && (elementOwner_[3][maxCoFacet[4]] == myRank_) )
01716 || ( (maxCoFacet[5] >= 0) && (elementOwner_[3][maxCoFacet[5]] == myRank_) )
01717 || ( (maxCoFacet[6] >= 0) && (elementOwner_[3][maxCoFacet[6]] == myRank_) )
01718 || ( (maxCoFacet[7] >= 0) && (elementOwner_[3][maxCoFacet[7]] == myRank_) )) ;
01719
01720
01721
01722
01723 if ( (hasCellLID[ind] == false) && (isPointHanging_[vertexIDs[v]] == true)){
01724 int parentID = parentCellLID_[ind];
01725 Array<int>& parentVertexIDs = cellsPoints_[parentID];
01726 hasCellLID[ind] = hasCellLID[ind] || (elementOwner_[0][parentVertexIDs[v]] == myRank_);
01727 }
01728 }
01729 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering Cell ID :" << ind << " should be LID: " << hasCellLID[ind] <<
01730 " ,isCellLeaf_[ind]:" << isCellLeaf_[ind]);
01731 }
01732
01733
01734
01735
01736
01737 bool check_Ghost_cells = true;
01738 while (check_Ghost_cells){
01739 check_Ghost_cells = false;
01740 for (int ind = 0 ; ind < nrElem_[3] ; ind++){
01741 if ( (hasCellLID[ind] == true) && (elementOwner_[3][ind] != myRank_ ) ){
01742 bool lookforEdges = true;
01743
01744 Array<int>& faceIDs = cellsFaces_[ind];
01745 for (int ii = 0 ; ii < 6 ; ii++ ){
01746
01747 if (isFaceHanging_[faceIDs[ii]] && ( elementOwner_[2][faceIDs[ii]] != myRank_)){
01748
01749 int parentCell = parentCellLID_[ind];
01750 Array<int>& parentfaceIDs = cellsFaces_[parentCell];
01751 for (int f = 0 ; f < 2 ; f++)
01752 if ( ( faceMaxCoF_[parentfaceIDs[ii]][f] >= 0 ) &&
01753 ( elementOwner_[3][ faceMaxCoF_[parentfaceIDs[ii]][f] ] != myRank_ ) &&
01754 ( hasCellLID[faceMaxCoF_[parentfaceIDs[ii]][f]] == false) &&
01755 ( isCellLeaf_[faceMaxCoF_[parentfaceIDs[ii]][f]] ) ){
01756 hasCellLID[faceMaxCoF_[parentfaceIDs[ii]][f]] = true;
01757 check_Ghost_cells = true;
01758 lookforEdges = false;
01759 }
01760 }
01761 }
01762
01763 Array<int>& edgeIDs = cellsEdges_[ind];
01764
01765 if (lookforEdges){
01766 for (int ii = 0 ; ii < 12 ; ii++ ){
01767
01768 if (isEdgeHanging_[edgeIDs[ii]] && ( elementOwner_[1][edgeIDs[ii]] != myRank_)){
01769
01770 int parentCell = parentCellLID_[ind];
01771 Array<int>& parentEdgesIDs = cellsEdges_[parentCell];
01772 for (int f = 0 ; f < 4 ; f++)
01773 if ( ( edgeMaxCoF_[parentEdgesIDs[ii]][f] >= 0 ) &&
01774 ( elementOwner_[3][ edgeMaxCoF_[parentEdgesIDs[ii]][f] ] != myRank_ ) &&
01775 ( hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] == false) &&
01776 ( isCellLeaf_[edgeMaxCoF_[parentEdgesIDs[ii]][f]] )
01777 ){
01778 hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] = true;
01779 check_Ghost_cells = true;
01780 }
01781 }
01782 }
01783 }
01784 }
01785 }
01786 }
01787
01788
01789 for (int ind = 0 ; ind < nrElem_[3] ; ind++)
01790 {
01791
01792
01793 if ( (isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
01794 {
01795 Array<int>& vertexIDs = cellsPoints_[ind];
01796 for (int v = 0; v < 8 ; v++)
01797 {
01798 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering vertexIDs[v]:" << vertexIDs[v] );
01799 if (vertexGIDToLeafMapping_[vertexIDs[v]] < 0)
01800 {
01801 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafGID_:" << nrVertexLeafGID_ );
01802 vertexLeafToGIDMapping_[nrVertexLeafGID_] = vertexIDs[v];
01803 vertexGIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafGID_;
01804 nrVertexLeafGID_++;
01805 }
01806 }
01807 Array<int>& edgeIDs = cellsEdges_[ind];
01808
01809 for (int e = 0; e < 12 ; e++)
01810 {
01811
01812 if (edgeGIDToLeafMapping_[edgeIDs[e]] < 0)
01813 {
01814 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafGID_:" << nrEdgeLeafGID_ );
01815
01816 edgeLeafToGIDMapping_[nrEdgeLeafGID_] = edgeIDs[e];
01817 edgeGIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafGID_;
01818 nrEdgeLeafGID_++;
01819 }
01820 }
01821 Array<int>& faceIDs = cellsFaces_[ind];
01822
01823 for (int f = 0; f < 6 ; f++)
01824 {
01825
01826 if (faceGIDToLeafMapping_[faceIDs[f]] < 0)
01827 {
01828 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> faceID:" << faceIDs[f] << " , nrFaceLeafGID_:" << nrFaceLeafGID_ );
01829
01830 faceLeafToGIDMapping_[nrFaceLeafGID_] = faceIDs[f];
01831 faceGIDToLeafMapping_[faceIDs[f]] = nrFaceLeafGID_;
01832 nrFaceLeafGID_++;
01833 }
01834 }
01835
01836 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering CELL cellID:" << ind << " , nrCellLeafGID_:" << nrCellLeafGID_ );
01837 cellLeafToGIDMapping_[nrCellLeafGID_] = ind;
01838 cellGIDToLeafMapping_[ind] = nrCellLeafGID_;
01839 nrCellLeafGID_++;
01840
01841
01842
01843 if (hasCellLID[ind]){
01844
01845 for (int v = 0; v < 8 ; v++)
01846 {
01847 if (vertexLIDToLeafMapping_[vertexIDs[v]] < 0)
01848 {
01849 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafLID_:" << nrVertexLeafLID_ );
01850 vertexLeafToLIDMapping_[nrVertexLeafLID_] = vertexIDs[v];
01851 vertexLIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafLID_;
01852 nrVertexLeafLID_++;
01853 }
01854 }
01855
01856 for (int e = 0; e < 12 ; e++)
01857 {
01858 if (edgeLIDToLeafMapping_[edgeIDs[e]] < 0)
01859 {
01860 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafLID_:" << nrEdgeLeafLID_ );
01861 edgeLeafToLIDMapping_[nrEdgeLeafLID_] = edgeIDs[e];
01862 edgeLIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafLID_;
01863 nrEdgeLeafLID_++;
01864 }
01865 }
01866
01867 for (int f = 0; f < 6 ; f++)
01868 {
01869 if (faceLIDToLeafMapping_[faceIDs[f]] < 0)
01870 {
01871 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> faceID:" << faceIDs[f] << " , nrFaceLeafLID_:" << nrFaceLeafLID_ );
01872 faceLeafToLIDMapping_[nrFaceLeafLID_] = faceIDs[f];
01873 faceLIDToLeafMapping_[faceIDs[f]] = nrFaceLeafLID_;
01874 nrFaceLeafLID_++;
01875 }
01876 }
01877
01878 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering CELL cellID:" << ind << " , nrCellLeafLID_:" << nrCellLeafLID_ );
01879 cellLeafToLIDMapping_[nrCellLeafLID_] = ind;
01880 cellLIDToLeafMapping_[ind] = nrCellLeafLID_;
01881 nrCellLeafLID_++;
01882 }
01883 }
01884 }
01885 SUNDANCE_MSG3(verb() , " nrVertexLeafGID_:" << nrVertexLeafGID_ << " nrEdgeLeafGID_:" << nrEdgeLeafGID_
01886 << " nrFaceLeafGID_:" << nrFaceLeafGID_ << " nrCellLeafGID_:" << nrCellLeafGID_ );
01887 SUNDANCE_MSG3(verb() , " nrVertexLeafLID_:" << nrVertexLeafLID_ << " nrEdgeLeafLID_:" << nrEdgeLeafLID_
01888 << " nrFaceLeafLID_:" <<nrFaceLeafLID_ << " nrCellLeafLID_:" << nrCellLeafLID_);
01889 SUNDANCE_MSG3(verb() , " vertexLIDToLeafMapping_: " << vertexLIDToLeafMapping_);
01890 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering , DONE");
01891 }
01892
01893
01894
01895
01896 int HNMesh3D::estimateCellLoad(int ID){
01897 int rtn = 0;
01898
01899 if (isCellLeaf_[ID]){
01900 if (!isCellOut_[ID]) rtn = 1;
01901 } else {
01902
01903 for (int r = 0 ; r < (int)cellsChildren_[ID].size() ; r++){
01904 if (cellLevel_[ID] < cellLevel_[cellsChildren_[ID][r]]){
01905 rtn = rtn + estimateCellLoad(cellsChildren_[ID][r]);
01906 }
01907 }
01908 }
01909 return rtn;
01910 }
01911
01912
01913 void HNMesh3D::markCellsAndFacets(int cellID , int procID){
01914
01915 if (elementOwner_[3][cellID] < 0) { elementOwner_[3][cellID] = procID; }
01916
01917 if (elementOwner_[2][cellsFaces_[cellID][0]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][0]] = procID;}
01918 if (elementOwner_[2][cellsFaces_[cellID][1]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][1]] = procID;}
01919 if (elementOwner_[2][cellsFaces_[cellID][2]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][2]] = procID;}
01920 if (elementOwner_[2][cellsFaces_[cellID][3]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][3]] = procID;}
01921 if (elementOwner_[2][cellsFaces_[cellID][4]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][4]] = procID;}
01922 if (elementOwner_[2][cellsFaces_[cellID][5]] < 0 ) { elementOwner_[2][cellsFaces_[cellID][5]] = procID;}
01923
01924 if (elementOwner_[1][cellsEdges_[cellID][0]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][0]] = procID;}
01925 if (elementOwner_[1][cellsEdges_[cellID][1]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][1]] = procID;}
01926 if (elementOwner_[1][cellsEdges_[cellID][2]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][2]] = procID;}
01927 if (elementOwner_[1][cellsEdges_[cellID][3]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][3]] = procID;}
01928 if (elementOwner_[1][cellsEdges_[cellID][4]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][4]] = procID;}
01929 if (elementOwner_[1][cellsEdges_[cellID][5]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][5]] = procID;}
01930 if (elementOwner_[1][cellsEdges_[cellID][6]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][6]] = procID;}
01931 if (elementOwner_[1][cellsEdges_[cellID][7]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][7]] = procID;}
01932 if (elementOwner_[1][cellsEdges_[cellID][8]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][8]] = procID;}
01933 if (elementOwner_[1][cellsEdges_[cellID][9]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][9]] = procID;}
01934 if (elementOwner_[1][cellsEdges_[cellID][10]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][10]] = procID;}
01935 if (elementOwner_[1][cellsEdges_[cellID][11]] < 0 ) { elementOwner_[1][cellsEdges_[cellID][11]] = procID;}
01936
01937 if (elementOwner_[0][cellsPoints_[cellID][0]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][0]] = procID;}
01938 if (elementOwner_[0][cellsPoints_[cellID][1]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][1]] = procID;}
01939 if (elementOwner_[0][cellsPoints_[cellID][2]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][2]] = procID;}
01940 if (elementOwner_[0][cellsPoints_[cellID][3]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][3]] = procID;}
01941 if (elementOwner_[0][cellsPoints_[cellID][4]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][4]] = procID;}
01942 if (elementOwner_[0][cellsPoints_[cellID][5]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][5]] = procID;}
01943 if (elementOwner_[0][cellsPoints_[cellID][6]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][6]] = procID;}
01944 if (elementOwner_[0][cellsPoints_[cellID][7]] < 0 ) { elementOwner_[0][cellsPoints_[cellID][7]] = procID;}
01945
01946 if (!isCellLeaf_[cellID]){
01947
01948 for (int r = 0 ; r < (int)cellsChildren_[cellID].size() ; r++){
01949 if (cellLevel_[cellID] < cellLevel_[cellsChildren_[cellID][r]]){
01950 markCellsAndFacets(cellsChildren_[cellID][r] , procID);
01951 }
01952 }
01953 }
01954 }
01955
01956 void HNMesh3D::createLeafNumbering_sophisticated(){
01957
01958
01959 Array<bool> hasCellLID(nrElem_[3],false);
01960 double total_load = 0.0;
01961
01962 Array<int> coarseCellLoad( _res_x * _res_y * _res_z , 1 );
01963
01964
01965
01966
01967
01968
01969
01970
01971 for (int ind = 0 ; ind < nrElem_[3] ; ind++){
01972 if (cellLevel_[ind] < 1) {
01973
01974 coarseCellLoad[ind] = estimateCellLoad(ind);
01975 }
01976 if ((isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
01977 { total_load = total_load + 1 ; }
01978 }
01979
01980 SUNDANCE_MSG3(verb() , "total_load = " << total_load << " , nrCell = " << nrElem_[3]);
01981
01982
01983
01984 int levelM = ::ceil( ::fmax( ::fmax( ::log2(_res_x) , ::log2(_res_y ) ) , ::log2(_res_z ) ) );
01985
01986 Array<int> vectX1(8), vectY1(8), vectZ1(8), vectX2(8), vectY2(8), vectZ2(8);
01987 vectX1[0] = 0; vectX1[1] = (int)::pow(2,levelM-1); vectX1[2] = 0; vectX1[3] = (int)::pow(2,levelM-1);
01988 vectX1[4] = 0; vectX1[5] = (int)::pow(2,levelM-1); vectX1[6] = 0; vectX1[7] = (int)::pow(2,levelM-1);
01989 vectY1[0] = 0; vectY1[1] = 0; vectY1[2] = (int)::pow(2,levelM-1); vectY1[3] = (int)::pow(2,levelM-1);
01990 vectY1[4] = 0; vectY1[5] = 0; vectY1[6] = (int)::pow(2,levelM-1); vectY1[7] = (int)::pow(2,levelM-1);
01991 vectZ1[0] = 0; vectZ1[1] = 0; vectZ1[2] = 0; vectZ1[3] = 0;
01992 vectZ1[4] = (int)::pow(2,levelM-1); vectZ1[5] = (int)::pow(2,levelM-1); vectZ1[6] = (int)::pow(2,levelM-1); vectZ1[7] = (int)::pow(2,levelM-1);
01993
01994 vectX2[0] = 0; vectX2[1] = (int)::pow(2,levelM-1); vectX2[2] = 0; vectX2[3] = (int)::pow(2,levelM-1);
01995 vectX2[4] = 0; vectX2[5] = (int)::pow(2,levelM-1); vectX2[6] = 0; vectX2[7] = (int)::pow(2,levelM-1);
01996 vectY2[0] = 0; vectY2[1] = 0; vectY2[2] = (int)::pow(2,levelM-1); vectY2[3] = (int)::pow(2,levelM-1);
01997 vectY2[4] = 0; vectY2[5] = 0; vectY2[6] = (int)::pow(2,levelM-1); vectY2[7] = (int)::pow(2,levelM-1);
01998 vectZ2[0] = 0; vectZ2[1] = 0; vectZ2[2] = 0; vectZ2[3] = 0;
01999 vectZ2[4] = (int)::pow(2,levelM-1); vectZ2[5] = (int)::pow(2,levelM-1); vectZ2[6] = (int)::pow(2,levelM-1); vectZ2[7] = (int)::pow(2,levelM-1);
02000
02001 int addX[8] = { 0 , 1 , 0 , 1 , 0 , 1 , 0 , 1};
02002 int addY[8] = { 0 , 0 , 1 , 1 , 0 , 0 , 1 , 1};
02003 int addZ[8] = { 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1};
02004 Array<int> *inX = &vectX1 , *inY = &vectY1 , *inZ = &vectZ1 ,*outX = &vectX2 , *outY = &vectY2 , *outZ = &vectZ2 , *tmpVectP;
02005 int levelActual = levelM - 2;
02006
02007 while (levelActual >= 0){
02008 outX->resize( 8 * inX->size() );
02009 outY->resize( 8 * inY->size() );
02010 outZ->resize( 8 * inZ->size() );
02011 int cI = 0 , addO = (int)::pow(2,levelActual);
02012 SUNDANCE_MSG3(verb() , " outX->size():" << outX->size() << ", levelActual:" << levelActual << " , addO:" << addO);
02013
02014 for (int ce = 0 ; ce < inX->size() ; ce++){
02015 (*outX)[cI+0] = (*inX)[ce] + addO*addX[0]; (*outY)[cI+0] = (*inY)[ce] + addO*addY[0]; (*outZ)[cI+0] = (*inZ)[ce] + addO*addZ[0];
02016 (*outX)[cI+1] = (*inX)[ce] + addO*addX[1]; (*outY)[cI+1] = (*inY)[ce] + addO*addY[1]; (*outZ)[cI+1] = (*inZ)[ce] + addO*addZ[1];
02017 (*outX)[cI+2] = (*inX)[ce] + addO*addX[2]; (*outY)[cI+2] = (*inY)[ce] + addO*addY[2]; (*outZ)[cI+2] = (*inZ)[ce] + addO*addZ[2];
02018 (*outX)[cI+3] = (*inX)[ce] + addO*addX[3]; (*outY)[cI+3] = (*inY)[ce] + addO*addY[3]; (*outZ)[cI+3] = (*inZ)[ce] + addO*addZ[3];
02019 (*outX)[cI+4] = (*inX)[ce] + addO*addX[4]; (*outY)[cI+4] = (*inY)[ce] + addO*addY[4]; (*outZ)[cI+4] = (*inZ)[ce] + addO*addZ[4];
02020 (*outX)[cI+5] = (*inX)[ce] + addO*addX[5]; (*outY)[cI+5] = (*inY)[ce] + addO*addY[5]; (*outZ)[cI+5] = (*inZ)[ce] + addO*addZ[5];
02021 (*outX)[cI+6] = (*inX)[ce] + addO*addX[6]; (*outY)[cI+6] = (*inY)[ce] + addO*addY[6]; (*outZ)[cI+6] = (*inZ)[ce] + addO*addZ[6];
02022 (*outX)[cI+7] = (*inX)[ce] + addO*addX[7]; (*outY)[cI+7] = (*inY)[ce] + addO*addY[7]; (*outZ)[cI+7] = (*inZ)[ce] + addO*addZ[7];
02023 cI = cI + 8;
02024 }
02025 SUNDANCE_MSG3(verb() , " EX: " << (*outX)[0] << " , " << (*outX)[1] << " , " << (*outX)[2]);
02026 SUNDANCE_MSG3(verb() , " EY: " << (*outY)[0] << " , " << (*outY)[1] << " , " << (*outY)[2]);
02027 SUNDANCE_MSG3(verb() , " EZ: " << (*outZ)[0] << " , " << (*outZ)[1] << " , " << (*outZ)[2]);
02028
02029 levelActual = levelActual - 1;
02030 tmpVectP = inX; inX = outX; outX = tmpVectP;
02031 tmpVectP = inY; inY = outY; outY = tmpVectP;
02032 tmpVectP = inZ; inZ = outZ; outZ = tmpVectP;
02033 }
02034
02035 tmpVectP = inX; inX = outX; outX = tmpVectP;
02036 tmpVectP = inY; inY = outY; outY = tmpVectP;
02037 tmpVectP = inZ; inZ = outZ; outZ = tmpVectP;
02038
02039
02040 for (int tmp = 0 ; tmp < nrElem_[0] ; tmp++ ){ elementOwner_[0][tmp] = -1; }
02041 for (int tmp = 0 ; tmp < nrElem_[1] ; tmp++ ){ elementOwner_[1][tmp] = -1; }
02042 for (int tmp = 0 ; tmp < nrElem_[2] ; tmp++ ){ elementOwner_[2][tmp] = -1; }
02043 for (int tmp = 0 ; tmp < nrElem_[3] ; tmp++ ){ elementOwner_[3][tmp] = -1; }
02044
02045
02046 int coarseCellID , actProcID = 0 , actualLoad = 0;
02047 double loadPerProc = (double)total_load / (double)nrProc_ , diff_load = 0.0;
02048 for (int ind = 0 ; ind < outX->size() ; ind++){
02049
02050 if ( ((*outX)[ind] < _res_x) && ((*outY)[ind] < _res_y) && ((*outZ)[ind] < _res_z) ){
02051
02052 coarseCellID = ((*outZ)[ind])*_res_x*_res_y + ((*outY)[ind])*_res_x + ((*outX)[ind]) ;
02053 SUNDANCE_MSG3(verb(),"Z-curve trav. ind:" << ind << " , coarseCellID:" << coarseCellID
02054 << " , indX:" << (*outX)[ind] << " , indY:" << (*outY)[ind] << " , indZ:" << (*outZ)[ind]);
02055
02056 TEUCHOS_TEST_FOR_EXCEPTION( cellLevel_[coarseCellID] > 0 , std::logic_error, " coarseCellID:" << coarseCellID << " has level:" << cellLevel_[coarseCellID] );
02057 markCellsAndFacets( coarseCellID , actProcID);
02058 actualLoad = actualLoad + coarseCellLoad[coarseCellID];
02059
02060 if (((double)actualLoad >= (loadPerProc - 1e-8 - diff_load)) && ( actProcID < nrProc_-1 )){
02061 SUNDANCE_MSG3(verb() , "Increase CPU , actualLoad:" << actualLoad << " loadPerProc:" << loadPerProc );
02062
02063 diff_load = actualLoad - loadPerProc;
02064 actProcID = actProcID + 1;
02065 actualLoad = 0;
02066 }
02067 }
02068 }
02069
02070
02071 SUNDANCE_MSG3(verb()," nrElem_[0]:" << nrElem_[0] << " , nrElem_[1]:" << nrElem_[1] << " , nrElem_[2]" << nrElem_[2]);
02072 for (int tmp = 0 ; tmp < nrElem_[0] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[0][tmp] < 0 , std::logic_error, " 0 tmp:" << tmp); }
02073 for (int tmp = 0 ; tmp < nrElem_[1] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[1][tmp] < 0 , std::logic_error, " 1 tmp:" << tmp); }
02074 for (int tmp = 0 ; tmp < nrElem_[2] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[2][tmp] < 0 , std::logic_error, " 2 tmp:" << tmp); }
02075 for (int tmp = 0 ; tmp < nrElem_[3] ; tmp++ ){ TEUCHOS_TEST_FOR_EXCEPTION( elementOwner_[3][tmp] < 0 , std::logic_error, " 3 tmp:" << tmp); }
02076
02077
02078
02079 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering nrPoint:" << nrElem_[0] << " , nrEdge:"
02080 << nrElem_[1] << ", nrFace:" << nrElem_[2] << ", nrCell:" << nrElem_[3]);
02081
02082 vertexGIDToLeafMapping_.resize(nrElem_[0],-1);
02083 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexGIDToLeafMapping_[dd] = -1;
02084 vertexLeafToGIDMapping_.resize(nrElem_[0],-1);
02085 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToGIDMapping_[dd] = -1;
02086
02087 edgeGIDToLeafMapping_.resize(nrElem_[1],-1);
02088 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeGIDToLeafMapping_[dd] = -1;
02089 edgeLeafToGIDMapping_.resize(nrElem_[1],-1);
02090 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToGIDMapping_[dd] = -1;
02091
02092 faceGIDToLeafMapping_.resize(nrElem_[2],-1);
02093 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceGIDToLeafMapping_[dd] = -1;
02094 faceLeafToGIDMapping_.resize(nrElem_[2],-1);
02095 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLeafToGIDMapping_[dd] = -1;
02096
02097 cellGIDToLeafMapping_.resize(nrElem_[3],-1);
02098 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellGIDToLeafMapping_[dd] = -1;
02099 cellLeafToGIDMapping_.resize(nrElem_[3],-1);
02100 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLeafToGIDMapping_[dd] = -1;
02101
02102 nrVertexLeafGID_ = 0; nrCellLeafGID_ = 0; nrEdgeLeafGID_ = 0; nrFaceLeafGID_ = 0;
02103
02104 nrVertexLeafLID_ = 0; nrCellLeafLID_ = 0; nrEdgeLeafLID_ = 0;
02105 vertexLIDToLeafMapping_.resize(nrElem_[0],-1);
02106 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLIDToLeafMapping_[dd] = -1;
02107 vertexLeafToLIDMapping_.resize(nrElem_[0],-1);
02108 for (int dd = 0 ; dd < nrElem_[0] ; dd++) vertexLeafToLIDMapping_[dd] = -1;
02109
02110 edgeLIDToLeafMapping_.resize(nrElem_[1],-1);
02111 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLIDToLeafMapping_[dd] = -1;
02112 edgeLeafToLIDMapping_.resize(nrElem_[1],-1);
02113 for (int dd = 0 ; dd < nrElem_[1] ; dd++) edgeLeafToLIDMapping_[dd] = -1;
02114
02115 faceLIDToLeafMapping_.resize(nrElem_[2],-1);
02116 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLIDToLeafMapping_[dd] = -1;
02117 faceLeafToLIDMapping_.resize(nrElem_[2],-1);
02118 for (int dd = 0 ; dd < nrElem_[2] ; dd++) faceLeafToLIDMapping_[dd] = -1;
02119
02120 cellLIDToLeafMapping_.resize(nrElem_[3],-1);
02121 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLIDToLeafMapping_[dd] = -1;
02122 cellLeafToLIDMapping_.resize(nrElem_[3],-1);
02123 for (int dd = 0 ; dd < nrElem_[3] ; dd++) cellLeafToLIDMapping_[dd] = -1;
02124
02125 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering , start assigning leaf numbers");
02126
02127
02128
02129 for (int ind = 0 ; ind < nrElem_[3] ; ind++){
02130 Array<int>& vertexIDs = cellsPoints_[ind];
02131 hasCellLID[ind] = false;
02132 for (int v = 0 ; v < 8 ; v++){
02133 Array<int>& maxCoFacet = pointMaxCoF_[vertexIDs[v]];
02134 hasCellLID[ind] = ( hasCellLID[ind]
02135 || ( (maxCoFacet[0] >= 0) && (elementOwner_[3][maxCoFacet[0]] == myRank_) )
02136 || ( (maxCoFacet[1] >= 0) && (elementOwner_[3][maxCoFacet[1]] == myRank_) )
02137 || ( (maxCoFacet[2] >= 0) && (elementOwner_[3][maxCoFacet[2]] == myRank_) )
02138 || ( (maxCoFacet[3] >= 0) && (elementOwner_[3][maxCoFacet[3]] == myRank_) )
02139 || ( (maxCoFacet[4] >= 0) && (elementOwner_[3][maxCoFacet[4]] == myRank_) )
02140 || ( (maxCoFacet[5] >= 0) && (elementOwner_[3][maxCoFacet[5]] == myRank_) )
02141 || ( (maxCoFacet[6] >= 0) && (elementOwner_[3][maxCoFacet[6]] == myRank_) )
02142 || ( (maxCoFacet[7] >= 0) && (elementOwner_[3][maxCoFacet[7]] == myRank_) )) ;
02143
02144
02145
02146
02147 if ( (hasCellLID[ind] == false) && (isPointHanging_[vertexIDs[v]] == true)){
02148 int parentID = parentCellLID_[ind];
02149 Array<int>& parentVertexIDs = cellsPoints_[parentID];
02150 hasCellLID[ind] = hasCellLID[ind] || (elementOwner_[0][parentVertexIDs[v]] == myRank_);
02151 }
02152 }
02153 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering Cell ID :" << ind << " should be LID: " << hasCellLID[ind] <<
02154 " ,isCellLeaf_[ind]:" << isCellLeaf_[ind]);
02155 }
02156
02157
02158
02159
02160
02161 bool check_Ghost_cells = true;
02162 while (check_Ghost_cells){
02163 check_Ghost_cells = false;
02164 for (int ind = 0 ; ind < nrElem_[3] ; ind++){
02165 if ( (hasCellLID[ind] == true) && (elementOwner_[3][ind] != myRank_ ) ){
02166 bool lookforEdges = true;
02167
02168 Array<int>& faceIDs = cellsFaces_[ind];
02169 for (int ii = 0 ; ii < 6 ; ii++ ){
02170
02171 if (isFaceHanging_[faceIDs[ii]] && ( elementOwner_[2][faceIDs[ii]] != myRank_)){
02172
02173 int parentCell = parentCellLID_[ind];
02174 Array<int>& parentfaceIDs = cellsFaces_[parentCell];
02175 for (int f = 0 ; f < 2 ; f++)
02176 if ( ( faceMaxCoF_[parentfaceIDs[ii]][f] >= 0 ) &&
02177 ( elementOwner_[3][ faceMaxCoF_[parentfaceIDs[ii]][f] ] != myRank_ ) &&
02178 ( hasCellLID[faceMaxCoF_[parentfaceIDs[ii]][f]] == false) &&
02179 ( isCellLeaf_[faceMaxCoF_[parentfaceIDs[ii]][f]] ) ){
02180 hasCellLID[faceMaxCoF_[parentfaceIDs[ii]][f]] = true;
02181 check_Ghost_cells = true;
02182 lookforEdges = false;
02183 }
02184 }
02185 }
02186
02187 Array<int>& edgeIDs = cellsEdges_[ind];
02188
02189 if (lookforEdges){
02190 for (int ii = 0 ; ii < 12 ; ii++ ){
02191
02192 if (isEdgeHanging_[edgeIDs[ii]] && ( elementOwner_[1][edgeIDs[ii]] != myRank_)){
02193
02194 int parentCell = parentCellLID_[ind];
02195 Array<int>& parentEdgesIDs = cellsEdges_[parentCell];
02196 for (int f = 0 ; f < 4 ; f++)
02197 if ( ( edgeMaxCoF_[parentEdgesIDs[ii]][f] >= 0 ) &&
02198 ( elementOwner_[3][ edgeMaxCoF_[parentEdgesIDs[ii]][f] ] != myRank_ ) &&
02199 ( hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] == false) &&
02200 ( isCellLeaf_[edgeMaxCoF_[parentEdgesIDs[ii]][f]] )
02201 ){
02202 hasCellLID[edgeMaxCoF_[parentEdgesIDs[ii]][f]] = true;
02203 check_Ghost_cells = true;
02204 }
02205 }
02206 }
02207 }
02208 }
02209 }
02210 }
02211
02212
02213 for (int ind = 0 ; ind < nrElem_[3] ; ind++)
02214 {
02215
02216
02217 if ( (isCellLeaf_[ind] == true) && (!isCellOut_[ind]) )
02218 {
02219 Array<int>& vertexIDs = cellsPoints_[ind];
02220 for (int v = 0; v < 8 ; v++)
02221 {
02222 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering vertexIDs[v]:" << vertexIDs[v] );
02223 if (vertexGIDToLeafMapping_[vertexIDs[v]] < 0)
02224 {
02225 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafGID_:" << nrVertexLeafGID_ );
02226 vertexLeafToGIDMapping_[nrVertexLeafGID_] = vertexIDs[v];
02227 vertexGIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafGID_;
02228 nrVertexLeafGID_++;
02229 }
02230 }
02231 Array<int>& edgeIDs = cellsEdges_[ind];
02232
02233 for (int e = 0; e < 12 ; e++)
02234 {
02235
02236 if (edgeGIDToLeafMapping_[edgeIDs[e]] < 0)
02237 {
02238 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafGID_:" << nrEdgeLeafGID_ );
02239
02240 edgeLeafToGIDMapping_[nrEdgeLeafGID_] = edgeIDs[e];
02241 edgeGIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafGID_;
02242 nrEdgeLeafGID_++;
02243 }
02244 }
02245 Array<int>& faceIDs = cellsFaces_[ind];
02246
02247 for (int f = 0; f < 6 ; f++)
02248 {
02249
02250 if (faceGIDToLeafMapping_[faceIDs[f]] < 0)
02251 {
02252 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering -> faceID:" << faceIDs[f] << " , nrFaceLeafGID_:" << nrFaceLeafGID_ );
02253
02254 faceLeafToGIDMapping_[nrFaceLeafGID_] = faceIDs[f];
02255 faceGIDToLeafMapping_[faceIDs[f]] = nrFaceLeafGID_;
02256 nrFaceLeafGID_++;
02257 }
02258 }
02259
02260 SUNDANCE_MSG3(verb() , " createLeafGIDNumbering CELL cellID:" << ind << " , nrCellLeafGID_:" << nrCellLeafGID_ );
02261 cellLeafToGIDMapping_[nrCellLeafGID_] = ind;
02262 cellGIDToLeafMapping_[ind] = nrCellLeafGID_;
02263 nrCellLeafGID_++;
02264
02265
02266
02267 if (hasCellLID[ind]){
02268
02269 for (int v = 0; v < 8 ; v++)
02270 {
02271 if (vertexLIDToLeafMapping_[vertexIDs[v]] < 0)
02272 {
02273 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> VertexID:" << vertexIDs[v] << " , nrVertexLeafLID_:" << nrVertexLeafLID_ );
02274 vertexLeafToLIDMapping_[nrVertexLeafLID_] = vertexIDs[v];
02275 vertexLIDToLeafMapping_[vertexIDs[v]] = nrVertexLeafLID_;
02276 nrVertexLeafLID_++;
02277 }
02278 }
02279
02280 for (int e = 0; e < 12 ; e++)
02281 {
02282 if (edgeLIDToLeafMapping_[edgeIDs[e]] < 0)
02283 {
02284 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> edgeID:" << edgeIDs[e] << " , nrEdgeLeafLID_:" << nrEdgeLeafLID_ );
02285 edgeLeafToLIDMapping_[nrEdgeLeafLID_] = edgeIDs[e];
02286 edgeLIDToLeafMapping_[edgeIDs[e]] = nrEdgeLeafLID_;
02287 nrEdgeLeafLID_++;
02288 }
02289 }
02290
02291 for (int f = 0; f < 6 ; f++)
02292 {
02293 if (faceLIDToLeafMapping_[faceIDs[f]] < 0)
02294 {
02295 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering -> faceID:" << faceIDs[f] << " , nrFaceLeafLID_:" << nrFaceLeafLID_ );
02296 faceLeafToLIDMapping_[nrFaceLeafLID_] = faceIDs[f];
02297 faceLIDToLeafMapping_[faceIDs[f]] = nrFaceLeafLID_;
02298 nrFaceLeafLID_++;
02299 }
02300 }
02301
02302 SUNDANCE_MSG3(verb() , " createLeafLIDNumbering CELL cellID:" << ind << " , nrCellLeafLID_:" << nrCellLeafLID_ );
02303 cellLeafToLIDMapping_[nrCellLeafLID_] = ind;
02304 cellLIDToLeafMapping_[ind] = nrCellLeafLID_;
02305 nrCellLeafLID_++;
02306 }
02307 }
02308 }
02309 SUNDANCE_MSG3(verb() , " nrVertexLeafGID_:" << nrVertexLeafGID_ << " nrEdgeLeafGID_:" << nrEdgeLeafGID_
02310 << " nrFaceLeafGID_:" << nrFaceLeafGID_ << " nrCellLeafGID_:" << nrCellLeafGID_ );
02311 SUNDANCE_MSG3(verb() , " nrVertexLeafLID_:" << nrVertexLeafLID_ << " nrEdgeLeafLID_:" << nrEdgeLeafLID_
02312 << " nrFaceLeafLID_:" <<nrFaceLeafLID_ << " nrCellLeafLID_:" << nrCellLeafLID_);
02313 SUNDANCE_MSG3(verb() , " vertexLIDToLeafMapping_: " << vertexLIDToLeafMapping_);
02314 SUNDANCE_MSG3(verb() , "HNMesh3D::createLeafNumbering , DONE");
02315 }