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 #include "SundanceBasicSimplicialMesh.hpp"
00032 #include "SundanceMeshType.hpp"
00033 #include "SundanceCellJacobianBatch.hpp"
00034 #include "SundanceMaximalCofacetBatch.hpp"
00035 #include "SundanceMeshSource.hpp"
00036 #include "SundanceDebug.hpp"
00037 #include "SundanceOut.hpp"
00038 #include "PlayaMPIContainerComm.hpp"
00039 #include "Teuchos_Time.hpp"
00040 #include "Teuchos_TimeMonitor.hpp"
00041 #include "SundanceObjectWithVerbosity.hpp"
00042 #include "SundanceCollectiveExceptionCheck.hpp"
00043
00044 #include <algorithm>
00045 #include <unistd.h>
00046
00047 using namespace Sundance;
00048 using namespace Teuchos;
00049 using Playa::MPIComm;
00050 using Playa::MPIContainerComm;
00051
00052 using std::endl;
00053
00054
00055 #ifdef BSM_LOW_LEVEL_TIMER
00056
00057 static Time& getAddElementTimer()
00058 {
00059 static RCP<Time> rtn
00060 = TimeMonitor::getNewTimer("BSM addElement");
00061 return *rtn;
00062 }
00063
00064 static Time& getFacetRegistrationTimer()
00065 {
00066 static RCP<Time> rtn
00067 = TimeMonitor::getNewTimer("BSM element facet registration");
00068 return *rtn;
00069 }
00070
00071
00072 static Time& getAddVertexTimer()
00073 {
00074 static RCP<Time> rtn
00075 = TimeMonitor::getNewTimer("BSM addVertex");
00076 return *rtn;
00077 }
00078
00079
00080 static Time& getAddEdgeTimer()
00081 {
00082 static RCP<Time> rtn
00083 = TimeMonitor::getNewTimer("BSM addEdge");
00084 return *rtn;
00085 }
00086
00087 static Time& getAddFaceTimer()
00088 {
00089 static RCP<Time> rtn
00090 = TimeMonitor::getNewTimer("BSM addFace");
00091 return *rtn;
00092 }
00093
00094 static Time& getLookupLIDTimer()
00095 {
00096 static RCP<Time> rtn
00097 = TimeMonitor::getNewTimer("BSM lookup LID");
00098 return *rtn;
00099 }
00100
00101 #endif
00102
00103 static Time& batchedFacetGrabTimer()
00104 {
00105 static RCP<Time> rtn
00106 = TimeMonitor::getNewTimer("batched facet grabbing");
00107 return *rtn;
00108 }
00109
00110 void sort(const Array<int>& in, Array<int>& rtn)
00111 {
00112 rtn.resize(in.size());
00113 const int* pIn = &(in[0]);
00114 int* pRtn = &(rtn[0]);
00115 for (int i=0; i<in.size(); i++) pRtn[i] = pIn[i];
00116
00117 for (int j=1; j<in.size(); j++)
00118 {
00119 int key = pRtn[j];
00120 int i=j-1;
00121 while (i>=0 && pRtn[i]>key)
00122 {
00123 pRtn[i+1]=pRtn[i];
00124 i--;
00125 }
00126 pRtn[i+1]=key;
00127 }
00128
00129 }
00130
00131
00132 static Time& getJacobianTimer()
00133 {
00134 static RCP<Time> rtn
00135 = TimeMonitor::getNewTimer("cell Jacobian grabbing");
00136 return *rtn;
00137 }
00138
00139
00140
00141
00142 BasicSimplicialMesh::BasicSimplicialMesh(int dim, const MPIComm& comm,
00143 const MeshEntityOrder& order)
00144 : IncrementallyCreatableMesh(dim, comm, order),
00145 numCells_(dim+1),
00146 points_(),
00147 edgeVerts_(2),
00148 faceVertLIDs_(dim),
00149 faceVertGIDs_(dim),
00150 faceEdges_(dim),
00151 faceEdgeSigns_(dim),
00152 elemVerts_(dim+1),
00153 elemEdges_(),
00154 elemEdgeSigns_(),
00155 elemFaces_(dim+1),
00156 elemFaceRotations_(dim+1),
00157 vertexSetToFaceIndexMap_(),
00158 edgeFaces_(),
00159 edgeCofacets_(),
00160 faceCofacets_(),
00161 vertEdges_(),
00162 vertFaces_(),
00163 vertCofacets_(),
00164 vertEdgePartners_(),
00165 LIDToGIDMap_(dim+1),
00166 GIDToLIDMap_(dim+1),
00167 labels_(dim+1),
00168 ownerProcID_(dim+1),
00169 faceVertGIDBase_(1),
00170 hasEdgeGIDs_(false),
00171 hasFaceGIDs_(false),
00172 neighbors_(),
00173 neighborsAreSynchronized_(false)
00174 {
00175 estimateNumVertices(1000);
00176 estimateNumElements(1000);
00177
00178
00179
00180
00181 faceVertGIDs_.resize(1);
00182 faceVertGIDBase_[0] = &(faceVertGIDs_.value(0,0));
00183 faceVertGIDs_.resize(0);
00184
00185
00186 if (spatialDim()==2)
00187 {
00188 elemEdges_.setTupleSize(3);
00189 elemEdgeSigns_.setTupleSize(3);
00190 }
00191 if (spatialDim()==3)
00192 {
00193 elemEdges_.setTupleSize(6);
00194 elemEdgeSigns_.setTupleSize(6);
00195 }
00196
00197 neighbors_.put(comm.getRank());
00198 }
00199
00200
00201 void BasicSimplicialMesh::getJacobians(int cellDim, const Array<int>& cellLID,
00202 CellJacobianBatch& jBatch) const
00203 {
00204 TimeMonitor timer(getJacobianTimer());
00205
00206 int flops = 0 ;
00207 int nCells = cellLID.size();
00208
00209 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00210 "cellDim=" << cellDim
00211 << " is not in expected range [0, " << spatialDim()
00212 << "]");
00213
00214 jBatch.resize(cellLID.size(), spatialDim(), cellDim);
00215
00216 if (cellDim < spatialDim())
00217 {
00218 switch(cellDim)
00219 {
00220 case 1:
00221 flops += 3*nCells;
00222 break;
00223 case 2:
00224
00225 flops += (4 + 10)*nCells;
00226 break;
00227 default:
00228 break;
00229 }
00230
00231 for (int i=0; i<nCells; i++)
00232 {
00233 int lid = cellLID[i];
00234 double* detJ = jBatch.detJ(i);
00235 switch(cellDim)
00236 {
00237 case 0:
00238 *detJ = 1.0;
00239 break;
00240 case 1:
00241 {
00242 int a = edgeVerts_.value(lid, 0);
00243 int b = edgeVerts_.value(lid, 1);
00244 const Point& pa = points_[a];
00245 const Point& pb = points_[b];
00246 Point dx = pb-pa;
00247 *detJ = sqrt(dx*dx);
00248 }
00249 break;
00250 case 2:
00251 {
00252 int a = faceVertLIDs_.value(lid, 0);
00253 int b = faceVertLIDs_.value(lid, 1);
00254 int c = faceVertLIDs_.value(lid, 2);
00255 const Point& pa = points_[a];
00256 const Point& pb = points_[b];
00257 const Point& pc = points_[c];
00258 Point directedArea = cross(pc-pa, pb-pa);
00259 *detJ = sqrt(directedArea*directedArea);
00260 }
00261 break;
00262 default:
00263 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00264 "cellDim=" << cellDim
00265 << " in BasicSimplicialMesh::getJacobians()");
00266 }
00267 }
00268 }
00269 else
00270 {
00271 Array<double> J(cellDim*cellDim);
00272
00273 flops += cellDim*cellDim*nCells;
00274
00275 for (int i=0; i<cellLID.size(); i++)
00276 {
00277 int lid = cellLID[i];
00278 double* J = jBatch.jVals(i);
00279 switch(cellDim)
00280 {
00281 case 0:
00282 J[0] = 1.0;
00283 break;
00284 case 1:
00285 {
00286 int a = elemVerts_.value(lid, 0);
00287 int b = elemVerts_.value(lid, 1);
00288 const Point& pa = points_[a];
00289 const Point& pb = points_[b];
00290 J[0] = fabs(pa[0]-pb[0]);
00291 }
00292 break;
00293 case 2:
00294 {
00295 int a = elemVerts_.value(lid, 0);
00296 int b = elemVerts_.value(lid, 1);
00297 int c = elemVerts_.value(lid, 2);
00298 const Point& pa = points_[a];
00299 const Point& pb = points_[b];
00300 const Point& pc = points_[c];
00301 J[0] = pb[0] - pa[0];
00302 J[1] = pc[0] - pa[0];
00303 J[2] = pb[1] - pa[1];
00304 J[3] = pc[1] - pa[1];
00305
00306 }
00307 break;
00308 case 3:
00309 {
00310 int a = elemVerts_.value(lid, 0);
00311 int b = elemVerts_.value(lid, 1);
00312 int c = elemVerts_.value(lid, 2);
00313 int d = elemVerts_.value(lid, 3);
00314 const Point& pa = points_[a];
00315 const Point& pb = points_[b];
00316 const Point& pc = points_[c];
00317 const Point& pd = points_[d];
00318 J[0] = pb[0] - pa[0];
00319 J[1] = pc[0] - pa[0];
00320 J[2] = pd[0] - pa[0];
00321 J[3] = pb[1] - pa[1];
00322 J[4] = pc[1] - pa[1];
00323 J[5] = pd[1] - pa[1];
00324 J[6] = pb[2] - pa[2];
00325 J[7] = pc[2] - pa[2];
00326 J[8] = pd[2] - pa[2];
00327 }
00328 break;
00329 default:
00330 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00331 "cellDim=" << cellDim
00332 << " in BasicSimplicialMesh::getJacobians()");
00333 }
00334 }
00335 }
00336 CellJacobianBatch::addFlops(flops);
00337 }
00338
00339
00340
00341 void BasicSimplicialMesh::getCellDiameters(int cellDim, const Array<int>& cellLID,
00342 Array<double>& cellDiameters) const
00343 {
00344 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00345 "cellDim=" << cellDim
00346 << " is not in expected range [0, " << spatialDim()
00347 << "]");
00348
00349 cellDiameters.resize(cellLID.size());
00350
00351 if (cellDim < spatialDim())
00352 {
00353 for (int i=0; i<cellLID.size(); i++)
00354 {
00355 int lid = cellLID[i];
00356 switch(cellDim)
00357 {
00358 case 0:
00359 cellDiameters[i] = 1.0;
00360 break;
00361 case 1:
00362 {
00363 int a = edgeVerts_.value(lid, 0);
00364 int b = edgeVerts_.value(lid, 1);
00365 const Point& pa = points_[a];
00366 const Point& pb = points_[b];
00367 Point dx = pb-pa;
00368 cellDiameters[i] = sqrt(dx*dx);
00369 }
00370 break;
00371 case 2:
00372 {
00373 int a = faceVertLIDs_.value(lid, 0);
00374 int b = faceVertLIDs_.value(lid, 1);
00375 int c = faceVertLIDs_.value(lid, 2);
00376 const Point& pa = points_[a];
00377 const Point& pb = points_[b];
00378 const Point& pc = points_[c];
00379 Point directedArea = cross(pc-pa, pb-pa);
00380 cellDiameters[i] = sqrt(directedArea*directedArea);
00381 }
00382 break;
00383 default:
00384 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00385 "cellDim=" << cellDim
00386 << " in BasicSimplicialMesh::getCellDiameters()");
00387 }
00388 }
00389 }
00390 else
00391 {
00392 for (int i=0; i<cellLID.size(); i++)
00393 {
00394 int lid = cellLID[i];
00395 switch(cellDim)
00396 {
00397 case 0:
00398 cellDiameters[i] = 1.0;
00399 break;
00400 case 1:
00401 {
00402 int a = elemVerts_.value(lid, 0);
00403 int b = elemVerts_.value(lid, 1);
00404 const Point& pa = points_[a];
00405 const Point& pb = points_[b];
00406 cellDiameters[i] = fabs(pa[0]-pb[0]);
00407 }
00408 break;
00409 case 2:
00410 {
00411 int a = elemVerts_.value(lid, 0);
00412 int b = elemVerts_.value(lid, 1);
00413 int c = elemVerts_.value(lid, 2);
00414 const Point& pa = points_[a];
00415 const Point& pb = points_[b];
00416 const Point& pc = points_[c];
00417 cellDiameters[i]
00418 = (pa.distance(pb) + pb.distance(pc) + pa.distance(pc))/3.0;
00419 }
00420 break;
00421 case 3:
00422 {
00423 int a = elemVerts_.value(lid, 0);
00424 int b = elemVerts_.value(lid, 1);
00425 int c = elemVerts_.value(lid, 2);
00426 int d = elemVerts_.value(lid, 3);
00427 const Point& pa = points_[a];
00428 const Point& pb = points_[b];
00429 const Point& pc = points_[c];
00430 const Point& pd = points_[d];
00431
00432 cellDiameters[i]
00433 = (pa.distance(pb) + pa.distance(pc) + pa.distance(pd)
00434 + pb.distance(pc) + pb.distance(pd)
00435 + pc.distance(pd))/6.0;
00436 }
00437 break;
00438 default:
00439 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00440 "cellDim=" << cellDim
00441 << " in BasicSimplicialMesh::getCellDiameters()");
00442 }
00443 }
00444 }
00445 }
00446
00447
00448 void BasicSimplicialMesh::pushForward(int cellDim, const Array<int>& cellLID,
00449 const Array<Point>& refQuadPts,
00450 Array<Point>& physQuadPts) const
00451 {
00452 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00453 "cellDim=" << cellDim
00454 << " is not in expected range [0, " << spatialDim()
00455 << "]");
00456 int flops = 0;
00457 int nQuad = refQuadPts.size();
00458 int nCells = cellLID.size();
00459 Array<double> J(cellDim*cellDim);
00460
00461 if (physQuadPts.size() > 0) physQuadPts.resize(0);
00462 physQuadPts.reserve(cellLID.size() * refQuadPts.size());
00463
00464 switch(cellDim)
00465 {
00466 case 1:
00467 flops += nCells * (1 + 2*nQuad);
00468 break;
00469 case 2:
00470 if (spatialDim()==2)
00471 {
00472 flops += nCells*(4 + 8*nQuad);
00473 }
00474 else
00475 {
00476 flops += 18*nCells*nQuad;
00477 }
00478 break;
00479 case 3:
00480 flops += 27*nCells*nQuad;
00481 break;
00482 default:
00483 break;
00484 }
00485
00486
00487 for (int i=0; i<cellLID.size(); i++)
00488 {
00489 int lid = cellLID[i];
00490 switch(cellDim)
00491 {
00492 case 0:
00493 physQuadPts.append(points_[lid]);
00494 break;
00495 case 1:
00496 {
00497 int a, b;
00498 if (spatialDim()==1)
00499 {
00500 a = elemVerts_.value(lid, 0);
00501 b = elemVerts_.value(lid, 1);
00502 }
00503 else
00504 {
00505 a = edgeVerts_.value(lid, 0);
00506 b = edgeVerts_.value(lid, 1);
00507 }
00508 const Point& pa = points_[a];
00509 const Point& pb = points_[b];
00510 Point dx = pb-pa;
00511 for (int q=0; q<nQuad; q++)
00512 {
00513 physQuadPts.append(pa + refQuadPts[q][0]*dx);
00514 }
00515 }
00516 break;
00517 case 2:
00518 {
00519 int a,b,c;
00520 if (spatialDim()==2)
00521 {
00522 a = elemVerts_.value(lid, 0);
00523 b = elemVerts_.value(lid, 1);
00524 c = elemVerts_.value(lid, 2);
00525 }
00526 else
00527 {
00528 a = faceVertLIDs_.value(lid, 0);
00529 b = faceVertLIDs_.value(lid, 1);
00530 c = faceVertLIDs_.value(lid, 2);
00531 }
00532 const Point& pa = points_[a];
00533 const Point& pb = points_[b];
00534 const Point& pc = points_[c];
00535
00536 if (spatialDim()==2)
00537 {
00538 J[0] = pb[0] - pa[0];
00539 J[1] = pc[0] - pa[0];
00540 J[2] = pb[1] - pa[1];
00541 J[3] = pc[1] - pa[1];
00542 for (int q=0; q<nQuad; q++)
00543 {
00544 physQuadPts.append(pa
00545 + Point(J[0]*refQuadPts[q][0] +J[1]*refQuadPts[q][1],
00546 J[2]*refQuadPts[q][0] +J[3]*refQuadPts[q][1]) );
00547 }
00548 }
00549 else
00550 {
00551 for (int q=0; q<nQuad; q++)
00552 {
00553 physQuadPts.append(pa
00554 + (pb-pa)*refQuadPts[q][0]
00555 + (pc-pa)*refQuadPts[q][1]);
00556 }
00557 }
00558
00559 }
00560 break;
00561 case 3:
00562 {
00563 int a = elemVerts_.value(lid, 0);
00564 int b = elemVerts_.value(lid, 1);
00565 int c = elemVerts_.value(lid, 2);
00566 int d = elemVerts_.value(lid, 3);
00567 const Point& pa = points_[a];
00568 const Point& pb = points_[b];
00569 const Point& pc = points_[c];
00570 const Point& pd = points_[d];
00571 J[0] = pb[0] - pa[0];
00572 J[1] = pc[0] - pa[0];
00573 J[2] = pd[0] - pa[0];
00574 J[3] = pb[1] - pa[1];
00575 J[4] = pc[1] - pa[1];
00576 J[5] = pd[1] - pa[1];
00577 J[6] = pb[2] - pa[2];
00578 J[7] = pc[2] - pa[2];
00579 J[8] = pd[2] - pa[2];
00580
00581 for (int q=0; q<nQuad; q++)
00582 {
00583 physQuadPts.append(pa
00584 + Point(J[0]*refQuadPts[q][0]
00585 + J[1]*refQuadPts[q][1]
00586 + J[2]*refQuadPts[q][2],
00587 J[3]*refQuadPts[q][0]
00588 + J[4]*refQuadPts[q][1]
00589 + J[5]*refQuadPts[q][2],
00590 J[6]*refQuadPts[q][0]
00591 + J[7]*refQuadPts[q][1]
00592 + J[8]*refQuadPts[q][2]));
00593
00594 }
00595
00596 }
00597 break;
00598 default:
00599 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00600 "in BasicSimplicialMesh::getJacobians()");
00601 }
00602 }
00603
00604 CellJacobianBatch::addFlops(flops);
00605 }
00606
00607 void BasicSimplicialMesh::estimateNumVertices(int nPts)
00608 {
00609 points_.reserve(nPts);
00610 vertCofacets_.reserve(nPts);
00611 vertEdges_.reserve(nPts);
00612 vertEdgePartners_.reserve(nPts);
00613
00614 ownerProcID_[0].reserve(nPts);
00615 LIDToGIDMap_[0].reserve(nPts);
00616 GIDToLIDMap_[0] = Hashtable<int,int>(nPts, 0.6);
00617 labels_[0].reserve(nPts);
00618 }
00619
00620 void BasicSimplicialMesh::estimateNumElements(int nElems)
00621 {
00622 int nEdges = 0;
00623 int nFaces = 0;
00624
00625 if (spatialDim()==3)
00626 {
00627 nFaces = 5*nElems;
00628 nEdges = 5*nElems;
00629 labels_[2].reserve(nFaces);
00630 }
00631 else if (spatialDim()==2)
00632 {
00633 nEdges = 3*nElems;
00634 }
00635
00636 labels_[1].reserve(nEdges);
00637 vertexSetToFaceIndexMap_ = Hashtable<VertexView, int>(nFaces);
00638
00639 edgeVerts_.reserve(nEdges);
00640 faceVertLIDs_.reserve(nFaces);
00641 faceVertGIDs_.reserve(nFaces);
00642 faceEdges_.reserve(nFaces);
00643 elemVerts_.reserve(nElems);
00644 elemEdges_.reserve(nElems);
00645 elemEdgeSigns_.reserve(nElems);
00646 elemFaces_.reserve(nElems);
00647 edgeCofacets_.reserve(nEdges);
00648 faceCofacets_.reserve(nFaces);
00649
00650 ownerProcID_[spatialDim()].reserve(nElems);
00651 LIDToGIDMap_[spatialDim()].reserve(nElems);
00652 GIDToLIDMap_[spatialDim()] = Hashtable<int,int>(nElems, 0.6);
00653 labels_[spatialDim()].reserve(nElems);
00654
00655
00656
00657 faceVertGIDs_.resize(1);
00658 faceVertGIDBase_[0] = &(faceVertGIDs_.value(0,0));
00659 faceVertGIDs_.resize(0);
00660 }
00661
00662
00663
00664 int BasicSimplicialMesh::numCells(int cellDim) const
00665 {
00666 return numCells_[cellDim];
00667 }
00668
00669 int BasicSimplicialMesh::ownerProcID(int cellDim, int cellLID) const
00670 {
00671 return ownerProcID_[cellDim][cellLID];
00672 }
00673
00674 int BasicSimplicialMesh::numFacets(int cellDim, int cellLID,
00675 int facetDim) const
00676 {
00677 if (cellDim==1)
00678 {
00679 return 2;
00680 }
00681 else if (cellDim==2)
00682 {
00683 return 3;
00684 }
00685 else
00686 {
00687 if (facetDim==0) return 4;
00688 if (facetDim==1) return 6;
00689 return 4;
00690 }
00691 }
00692
00693 void BasicSimplicialMesh::getFacetLIDs(int cellDim,
00694 const Array<int>& cellLID,
00695 int facetDim,
00696 Array<int>& facetLID,
00697 Array<int>& facetSign) const
00698 {
00699 TimeMonitor timer(batchedFacetGrabTimer());
00700
00701 int nf = numFacets(cellDim, cellLID[0], facetDim);
00702 facetLID.resize(cellLID.size() * nf);
00703 if (facetDim > 0) facetSign.resize(cellLID.size() * nf);
00704
00705
00706 if (facetDim==0)
00707 {
00708 if (cellDim == spatialDim())
00709 {
00710 int ptr=0;
00711 for (int c=0; c<cellLID.size(); c++)
00712 {
00713 const int* fPtr = &(elemVerts_.value(cellLID[c], 0));
00714 for (int f=0; f<nf; f++, ptr++)
00715 {
00716 facetLID[ptr] = fPtr[f];
00717 }
00718 }
00719 }
00720 else if (cellDim==1)
00721 {
00722 int ptr=0;
00723 for (int c=0; c<cellLID.size(); c++)
00724 {
00725 const int* fPtr = &(edgeVerts_.value(cellLID[c], 0));
00726 for (int f=0; f<nf; f++, ptr++)
00727 {
00728 facetLID[ptr] = fPtr[f];
00729 }
00730 }
00731 }
00732 else if (cellDim==2)
00733 {
00734 int ptr=0;
00735 for (int c=0; c<cellLID.size(); c++)
00736 {
00737 const int* fPtr = &(faceVertLIDs_.value(cellLID[c], 0));
00738 for (int f=0; f<nf; f++, ptr++)
00739 {
00740 facetLID[ptr] = fPtr[f];
00741 }
00742 }
00743 }
00744 }
00745 else if (facetDim==1)
00746 {
00747 int ptr=0;
00748 if (cellDim == spatialDim())
00749 {
00750 for (int c=0; c<cellLID.size(); c++)
00751 {
00752 const int* fPtr = &(elemEdges_.value(cellLID[c], 0));
00753
00754 for (int f=0; f<nf; f++, ptr++)
00755 {
00756 facetLID[ptr] = fPtr[f];
00757 facetSign[ptr] = 1;
00758 }
00759 }
00760 }
00761 else
00762 {
00763 for (int c=0; c<cellLID.size(); c++)
00764 {
00765 const int* fPtr = &(faceEdges_.value(cellLID[c], 0));
00766
00767 for (int f=0; f<nf; f++, ptr++)
00768 {
00769 facetLID[ptr] = fPtr[f];
00770
00771 }
00772 }
00773 }
00774 }
00775 else
00776 {
00777 int ptr=0;
00778 for (int c=0; c<cellLID.size(); c++)
00779 {
00780 const int* fPtr = &(elemFaces_.value(cellLID[c], 0));
00781
00782 for (int f=0; f<nf; f++, ptr++)
00783 {
00784 facetLID[ptr] = fPtr[f];
00785 facetSign[ptr] = 1;
00786 }
00787 }
00788 }
00789 }
00790
00791 int BasicSimplicialMesh::facetLID(int cellDim, int cellLID,
00792 int facetDim, int facetIndex, int& facetSign) const
00793 {
00794
00795 if (facetDim==0)
00796 {
00797 if (cellDim == spatialDim())
00798 {
00799 return elemVerts_.value(cellLID, facetIndex);
00800 }
00801 else if (cellDim==1) return edgeVerts_.value(cellLID, facetIndex);
00802 else if (cellDim==2) return faceVertLIDs_.value(cellLID, facetIndex);
00803 }
00804 if (facetDim==1)
00805 {
00806 if (cellDim==spatialDim())
00807 {
00808 facetSign = 1;
00809 return elemEdges_.value(cellLID, facetIndex);
00810 }
00811 else
00812 {
00813
00814 return faceEdges_.value(cellLID, facetIndex);
00815 }
00816 }
00817 else
00818 {
00819 facetSign = 1;
00820 return elemFaces_.value(cellLID, facetIndex);
00821 }
00822 }
00823
00824
00825 int BasicSimplicialMesh::numMaxCofacets(int cellDim, int cellLID) const
00826 {
00827
00828 if (cellDim==0)
00829 {
00830 return vertCofacets_[cellLID].length();
00831 }
00832 if (cellDim==1)
00833 {
00834 return edgeCofacets_[cellLID].length();
00835 }
00836 if (cellDim==2)
00837 {
00838 return faceCofacets_[cellLID].length();
00839 }
00840 return -1;
00841 }
00842
00843
00844
00845 int BasicSimplicialMesh::maxCofacetLID(int cellDim, int cellLID,
00846 int cofacetIndex,
00847 int& facetIndex) const
00848 {
00849 int rtn = -1;
00850 if (cellDim==0)
00851 {
00852 rtn = vertCofacets_[cellLID][cofacetIndex];
00853 }
00854 else if (cellDim==1)
00855 {
00856 rtn = edgeCofacets_[cellLID][cofacetIndex];
00857 }
00858 else if (cellDim==2)
00859 {
00860 rtn = faceCofacets_[cellLID][cofacetIndex];
00861 }
00862 else
00863 {
00864 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00865 "invalid cell dimension " << cellDim
00866 << " in request for cofacet");
00867 }
00868
00869
00870
00871 int dummy;
00872 for (int f=0; f<numFacets(spatialDim(), rtn, cellDim); f++)
00873 {
00874
00875 int c = facetLID(spatialDim(), rtn, cellDim, f, dummy);
00876
00877 if (cellLID==c)
00878 {
00879 facetIndex = f;
00880 return rtn;
00881 }
00882 }
00883 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "reverse pointer to facet not found"
00884 " in request for cofacet");
00885 return -1;
00886 }
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916 void BasicSimplicialMesh::getCofacets(int cellDim, int cellLID,
00917 int cofacetDim, Array<int>& cofacetLIDs) const
00918 {
00919
00920 TEUCHOS_TEST_FOR_EXCEPTION(cofacetDim > spatialDim() || cofacetDim < 0, std::runtime_error,
00921 "invalid cofacet dimension=" << cofacetDim);
00922 TEUCHOS_TEST_FOR_EXCEPTION( cofacetDim <= cellDim, std::runtime_error,
00923 "invalid cofacet dimension=" << cofacetDim
00924 << " for cell dim=" << cellDim);
00925 if (cofacetDim==spatialDim())
00926 {
00927 cofacetLIDs.resize(numMaxCofacets(cellDim, cellLID));
00928 int dum=0;
00929 for (int f=0; f<cofacetLIDs.size(); f++)
00930 {
00931 cofacetLIDs[f] = maxCofacetLID(cellDim, cellLID, f, dum);
00932 }
00933 }
00934 else
00935 {
00936 if (cellDim==0)
00937 {
00938 if (cofacetDim==1) cofacetLIDs = vertEdges_[cellLID];
00939 else if (cofacetDim==2) cofacetLIDs = vertFaces_[cellLID];
00940 else TEUCHOS_TEST_FOR_EXCEPT(true);
00941 }
00942 else if (cellDim==1)
00943 {
00944 if (cofacetDim==2) cofacetLIDs = edgeFaces_[cellLID];
00945 else TEUCHOS_TEST_FOR_EXCEPT(true);
00946 }
00947 else if (cellDim==2)
00948 {
00949
00950
00951 TEUCHOS_TEST_FOR_EXCEPT(true);
00952 }
00953 else
00954 {
00955 TEUCHOS_TEST_FOR_EXCEPT(true);
00956 }
00957 }
00958 }
00959
00960 int BasicSimplicialMesh::mapGIDToLID(int cellDim, int globalIndex) const
00961 {
00962 return GIDToLIDMap_[cellDim].get(globalIndex);
00963 }
00964
00965 bool BasicSimplicialMesh::hasGID(int cellDim, int globalIndex) const
00966 {
00967 return GIDToLIDMap_[cellDim].containsKey(globalIndex);
00968 }
00969
00970 int BasicSimplicialMesh::mapLIDToGID(int cellDim, int localIndex) const
00971 {
00972 return LIDToGIDMap_[cellDim][localIndex];
00973 }
00974
00975 CellType BasicSimplicialMesh::cellType(int cellDim) const
00976 {
00977
00978 switch(cellDim)
00979 {
00980 case 0:
00981 return PointCell;
00982 case 1:
00983 return LineCell;
00984 case 2:
00985 return TriangleCell;
00986 case 3:
00987 return TetCell;
00988 default:
00989 return NullCell;
00990 }
00991 }
00992
00993 int BasicSimplicialMesh::label(int cellDim,
00994 int cellLID) const
00995 {
00996 return labels_[cellDim][cellLID];
00997 }
00998
00999 void BasicSimplicialMesh::getLabels(int cellDim, const Array<int>& cellLID,
01000 Array<int>& labels) const
01001 {
01002 labels.resize(cellLID.size());
01003 const Array<int>& ld = labels_[cellDim];
01004 for (int i=0; i<cellLID.size(); i++)
01005 {
01006 labels[i] = ld[cellLID[i]];
01007 }
01008 }
01009
01010
01011 Set<int> BasicSimplicialMesh::getAllLabelsForDimension(int cellDim) const
01012 {
01013 Set<int> rtn;
01014
01015 const Array<int>& ld = labels_[cellDim];
01016 for (int i=0; i<ld.size(); i++)
01017 {
01018 rtn.put(ld[i]);
01019 }
01020
01021 return rtn;
01022 }
01023
01024 void BasicSimplicialMesh::getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const
01025 {
01026 cellLIDs.resize(0);
01027 const Array<int>& ld = labels_[cellDim];
01028 for (int i=0; i<ld.size(); i++)
01029 {
01030 if (ld[i] == label) cellLIDs.append(i);
01031 }
01032 }
01033
01034
01035
01036
01037 int BasicSimplicialMesh::addVertex(int globalIndex, const Point& x,
01038 int ownerProcID, int label)
01039 {
01040 #ifdef BSM_LOW_LEVEL_TIMER
01041 TimeMonitor timer(getAddVertexTimer());
01042 #endif
01043 int lid = points_.length();
01044 points_.append(x);
01045
01046 SUNDANCE_OUT(this->verb() > 4,
01047 "BSM added point " << x << " lid = " << lid);
01048
01049 numCells_[0]++;
01050
01051 LIDToGIDMap_[0].append(globalIndex);
01052 GIDToLIDMap_[0].put(globalIndex, lid);
01053
01054 if (comm().getRank() != ownerProcID) neighbors_.put(ownerProcID);
01055 ownerProcID_[0].append(ownerProcID);
01056 labels_[0].append(label);
01057
01058 vertCofacets_.resize(points_.length());
01059 vertEdges_.resize(points_.length());
01060 if (spatialDim() > 2) vertFaces_.resize(points_.length());
01061 vertEdgePartners_.resize(points_.length());
01062
01063 return lid;
01064 }
01065
01066 int BasicSimplicialMesh::addElement(int globalIndex,
01067 const Array<int>& vertGID,
01068 int ownerProcID, int label)
01069 {
01070 #ifdef BSM_LOW_LEVEL_TIMER
01071 TimeMonitor timer(getAddElementTimer());
01072 #endif
01073 SUNDANCE_VERB_HIGH("p=" << comm().getRank() << "adding element " << globalIndex
01074 << " with vertices:" << vertGID);
01075
01076 static Array<int> sortedVertGID;
01077 sort(vertGID, sortedVertGID);
01078
01079
01080
01081
01082
01083 int lid = elemVerts_.length();
01084 elemEdgeSigns_.resize(lid+1);
01085
01086 numCells_[spatialDim()]++;
01087
01088 LIDToGIDMap_[spatialDim()].append(globalIndex);
01089 GIDToLIDMap_[spatialDim()].put(globalIndex, lid);
01090 labels_[spatialDim()].append(label);
01091 ownerProcID_[spatialDim()].append(ownerProcID);
01092 if (comm().getRank() != ownerProcID) neighbors_.put(ownerProcID);
01093
01094
01095
01096 static Array<int> edges;
01097 static Array<int> faces;
01098 static Array<int> faceRotations;
01099 static Array<int> vertLID;
01100
01101
01102 {
01103 #ifdef BSM_LOW_LEVEL_TIMER
01104 TimeMonitor timer(getLookupLIDTimer());
01105 #endif
01106 vertLID.resize(vertGID.size());
01107 for (int i=0; i<vertGID.size(); i++)
01108 {
01109 vertLID[i] = GIDToLIDMap_[0].get(sortedVertGID[i]);
01110 }
01111 }
01112
01113
01114
01115
01116
01117
01118
01119
01120 if (spatialDim()==1)
01121 {
01122
01123 edges.resize(0);
01124 faces.resize(0);
01125
01126 vertCofacets_[vertLID[0]].append(lid);
01127 vertCofacets_[vertLID[1]].append(lid);
01128 }
01129 if (spatialDim()==2)
01130 {
01131
01132 edges.resize(3);
01133 faces.resize(0);
01134
01135
01136 edges[0] = addEdge(vertLID[1], vertLID[2], lid, globalIndex, 0);
01137 edges[1] = addEdge(vertLID[0], vertLID[2], lid, globalIndex, 1);
01138 edges[2] = addEdge(vertLID[0], vertLID[1], lid, globalIndex, 2);
01139
01140
01141 vertCofacets_[vertLID[0]].append(lid);
01142 vertCofacets_[vertLID[1]].append(lid);
01143 vertCofacets_[vertLID[2]].append(lid);
01144
01145 }
01146 else if (spatialDim()==3)
01147 {
01148
01149 edges.resize(6);
01150 faces.resize(4);
01151
01152
01153 edges[0] = addEdge(vertLID[2], vertLID[3], lid, globalIndex, 0);
01154 edges[1] = addEdge(vertLID[1], vertLID[3], lid, globalIndex, 1);
01155 edges[2] = addEdge(vertLID[1], vertLID[2], lid, globalIndex, 2);
01156 edges[3] = addEdge(vertLID[0], vertLID[3], lid, globalIndex, 3);
01157 edges[4] = addEdge(vertLID[0], vertLID[2], lid, globalIndex, 4);
01158 edges[5] = addEdge(vertLID[0], vertLID[1], lid, globalIndex, 5);
01159
01160
01161 vertCofacets_[vertLID[0]].append(lid);
01162 vertCofacets_[vertLID[1]].append(lid);
01163 vertCofacets_[vertLID[2]].append(lid);
01164 vertCofacets_[vertLID[3]].append(lid);
01165
01166
01167 static Array<int> tmpVertLID(3);
01168 static Array<int> tmpSortedVertGID(3);
01169 static Array<int> tmpEdges(4);
01170 #define BSM_OPT_CODE
01171 #ifndef BSM_OPT_CODE
01172 faces[0] = addFace(tuple(vertLID[1], vertLID[2], vertLID[3]),
01173 tuple(sortedVertGID[1], sortedVertGID[2], sortedVertGID[3]),
01174 tuple(edges[2], edges[0], edges[1]), lid, globalIndex);
01175
01176 faces[1] = addFace(tuple(vertLID[0], vertLID[2], vertLID[3]),
01177 tuple(sortedVertGID[0], sortedVertGID[2], sortedVertGID[3]),
01178 tuple(edges[3], edges[4], edges[0]), lid, globalIndex);
01179
01180 faces[2] = addFace(tuple(vertLID[0], vertLID[1], vertLID[3]),
01181 tuple(sortedVertGID[0], sortedVertGID[1], sortedVertGID[3]),
01182 tuple(edges[5], edges[1], edges[3]), lid, globalIndex);
01183
01184 faces[3] = addFace(tuple(vertLID[0], vertLID[1], vertLID[2]),
01185 tuple(sortedVertGID[0], sortedVertGID[1], sortedVertGID[2]),
01186 tuple(edges[5], edges[2], edges[4]), lid, globalIndex);
01187 #else
01188 tmpVertLID[0] = vertLID[1];
01189 tmpVertLID[1] = vertLID[2];
01190 tmpVertLID[2] = vertLID[3];
01191 tmpSortedVertGID[0] = sortedVertGID[1];
01192 tmpSortedVertGID[1] = sortedVertGID[2];
01193 tmpSortedVertGID[2] = sortedVertGID[3];
01194 tmpEdges[0] = edges[2];
01195 tmpEdges[1] = edges[0];
01196 tmpEdges[2] = edges[1];
01197 faces[0] = addFace(tmpVertLID, tmpSortedVertGID,
01198 tmpEdges, lid, globalIndex);
01199
01200 tmpVertLID[0] = vertLID[0];
01201 tmpVertLID[1] = vertLID[2];
01202 tmpVertLID[2] = vertLID[3];
01203 tmpSortedVertGID[0] = sortedVertGID[0];
01204 tmpSortedVertGID[1] = sortedVertGID[2];
01205 tmpSortedVertGID[2] = sortedVertGID[3];
01206 tmpEdges[0] = edges[3];
01207 tmpEdges[1] = edges[4];
01208 tmpEdges[2] = edges[0];
01209 faces[1] = addFace(tmpVertLID, tmpSortedVertGID,
01210 tmpEdges, lid, globalIndex);
01211
01212 tmpVertLID[0] = vertLID[0];
01213 tmpVertLID[1] = vertLID[1];
01214 tmpVertLID[2] = vertLID[3];
01215 tmpSortedVertGID[0] = sortedVertGID[0];
01216 tmpSortedVertGID[1] = sortedVertGID[1];
01217 tmpSortedVertGID[2] = sortedVertGID[3];
01218 tmpEdges[0] = edges[5];
01219 tmpEdges[1] = edges[1];
01220 tmpEdges[2] = edges[3];
01221 faces[2] = addFace(tmpVertLID, tmpSortedVertGID,
01222 tmpEdges, lid, globalIndex);
01223
01224 tmpVertLID[0] = vertLID[0];
01225 tmpVertLID[1] = vertLID[1];
01226 tmpVertLID[2] = vertLID[2];
01227 tmpSortedVertGID[0] = sortedVertGID[0];
01228 tmpSortedVertGID[1] = sortedVertGID[1];
01229 tmpSortedVertGID[2] = sortedVertGID[2];
01230 tmpEdges[0] = edges[5];
01231 tmpEdges[1] = edges[2];
01232 tmpEdges[2] = edges[4];
01233 faces[3] = addFace(tmpVertLID, tmpSortedVertGID,
01234 tmpEdges, lid, globalIndex);
01235
01236
01237 #endif
01238 }
01239
01240 {
01241 #ifdef BSM_LOW_LEVEL_TIMER
01242 TimeMonitor timer(getFacetRegistrationTimer());
01243 #endif
01244 elemVerts_.append(vertLID);
01245 if (edges.length() > 0) elemEdges_.append(edges);
01246
01247 if (faces.length() > 0)
01248 {
01249 elemFaces_.append(faces);
01250 elemFaceRotations_.append(faceRotations);
01251 }
01252
01253 }
01254
01255 for (int i=0; i<edges.length(); i++) edgeCofacets_[edges[i]].append(lid);
01256
01257
01258 for (int i=0; i<faces.length(); i++) faceCofacets_[faces[i]].append(lid);
01259
01260 return lid;
01261 }
01262
01263 int BasicSimplicialMesh::addFace(
01264 const Array<int>& vertLID,
01265 const Array<int>& vertGID,
01266 const Array<int>& edgeLID,
01267 int elemLID,
01268 int elemGID)
01269 {
01270 #ifdef BSM_LOW_LEVEL_TIMER
01271 TimeMonitor timer(getAddFaceTimer());
01272 #endif
01273
01274
01275 int lid = lookupFace(vertGID);
01276
01277 if (lid >= 0)
01278 {
01279 return lid;
01280 }
01281 else
01282 {
01283
01284 lid = faceVertLIDs_.length();
01285
01286
01287
01288 faceVertGIDs_.append(&(vertGID[0]), 3);
01289 faceVertLIDs_.append(&(vertLID[0]), 3);
01290 faceEdges_.append(&(edgeLID[0]), 3);
01291
01292 SUNDANCE_VERB_EXTREME("p=" << comm().getRank() << " adding face "
01293 << vertGID);
01294
01295
01296
01297
01298 int vert1Owner = ownerProcID_[0][vertLID[0]];
01299 int vert2Owner = ownerProcID_[0][vertLID[1]];
01300 int vert3Owner = ownerProcID_[0][vertLID[2]];
01301 int myRank = comm().getRank();
01302 if (vert1Owner==myRank && vert2Owner==myRank && vert3Owner==myRank)
01303 {
01304 ownerProcID_[2].append(myRank);
01305 }
01306 else ownerProcID_[2].append(-1);
01307
01308
01309 labels_[2].append(0);
01310
01311
01312
01313 faceVertGIDBase_[0] = &(faceVertGIDs_.value(0,0));
01314
01315
01316
01317
01318
01319 VertexView face(&(faceVertGIDBase_[0]), lid, 3);
01320
01321
01322
01323 vertexSetToFaceIndexMap_.put(face, lid);
01324
01325
01326 faceCofacets_.resize(lid+1);
01327
01328
01329 numCells_[spatialDim()-1]++;
01330
01331
01332 edgeFaces_[edgeLID[0]].append(lid);
01333 edgeFaces_[edgeLID[1]].append(lid);
01334 edgeFaces_[edgeLID[2]].append(lid);
01335
01336
01337 vertFaces_[vertLID[0]].append(lid);
01338 vertFaces_[vertLID[1]].append(lid);
01339 vertFaces_[vertLID[2]].append(lid);
01340
01341
01342 return lid;
01343 }
01344
01345 }
01346
01347
01348 int BasicSimplicialMesh::checkForExistingEdge(int vertLID1, int vertLID2)
01349 {
01350 const Array<int>& edgePartners = vertEdgePartners_[vertLID1];
01351 for (int i=0; i<edgePartners.length(); i++)
01352 {
01353 if (edgePartners[i] == vertLID2)
01354 {
01355 return vertEdges_[vertLID1][i];
01356 }
01357 }
01358 return -1;
01359 }
01360
01361 int BasicSimplicialMesh::lookupFace(const Array<int>& vertGID)
01362 {
01363
01364
01365 int* base = const_cast<int*>(&(vertGID[0]));
01366 int** basePtr = &base;
01367 VertexView inputFaceView(basePtr, 0, 3);
01368
01369 if (vertexSetToFaceIndexMap_.containsKey(inputFaceView))
01370 {
01371
01372 return vertexSetToFaceIndexMap_.get(inputFaceView);
01373 }
01374 else
01375 {
01376
01377 return -1;
01378 }
01379 }
01380
01381 int BasicSimplicialMesh::addEdge(int v1, int v2,
01382 int elemLID, int elemGID,
01383 int myFacetNumber)
01384 {
01385 #ifdef BSM_LOW_LEVEL_TIMER
01386 TimeMonitor timer(getAddEdgeTimer());
01387 #endif
01388
01389
01390
01391
01392 int lid = checkForExistingEdge(v1, v2);
01393
01394 if (lid >= 0)
01395 {
01396
01397 int g1 = LIDToGIDMap_[0][v1];
01398 int g2 = LIDToGIDMap_[0][v2];
01399 TEUCHOS_TEST_FOR_EXCEPT(g2 <= g1);
01400
01401
01402 return lid;
01403 }
01404 else
01405 {
01406
01407 lid = edgeVerts_.length();
01408
01409 edgeVerts_.append(tuple(v1,v2));
01410
01411
01412
01413 int vert1Owner = ownerProcID_[0][v1];
01414 int vert2Owner = ownerProcID_[0][v2];
01415 if (vert2Owner==comm().getRank() && vert1Owner==comm().getRank())
01416 ownerProcID_[1].append(vert1Owner);
01417 else ownerProcID_[1].append(-1);
01418
01419
01420 vertEdges_[v1].append(lid);
01421 vertEdgePartners_[v1].append(v2);
01422 vertEdges_[v2].append(lid);
01423 vertEdgePartners_[v2].append(v1);
01424
01425 edgeCofacets_.resize(lid+1);
01426 if (spatialDim() > 2) edgeFaces_.resize(lid+1);
01427
01428
01429
01430 labels_[1].append(0);
01431
01432
01433 numCells_[1]++;
01434 }
01435
01436 return lid;
01437 }
01438
01439
01440 void BasicSimplicialMesh::synchronizeGIDNumbering(int dim, int localCount)
01441 {
01442
01443 SUNDANCE_OUT(this->verb() > 4,
01444 "sharing offsets for GID numbering for dim=" << dim);
01445 SUNDANCE_OUT(this->verb() > 4,
01446 "I have " << localCount << " cells");
01447 Array<int> gidOffsets;
01448 int total;
01449 MPIContainerComm<int>::accumulate(localCount, gidOffsets, total, comm());
01450 int myOffset = gidOffsets[comm().getRank()];
01451
01452 SUNDANCE_OUT(this->verb() > 4,
01453 "back from MPI accumulate: offset for d=" << dim
01454 << " on p=" << comm().getRank() << " is " << myOffset);
01455
01456 for (int i=0; i<numCells(dim); i++)
01457 {
01458 if (LIDToGIDMap_[dim][i] == -1) continue;
01459 LIDToGIDMap_[dim][i] += myOffset;
01460 GIDToLIDMap_[dim].put(LIDToGIDMap_[dim][i], i);
01461 }
01462 SUNDANCE_OUT(this->verb() > 4,
01463 "done sharing offsets for GID numbering for dim=" << dim);
01464 }
01465
01466
01467 void BasicSimplicialMesh::assignIntermediateCellGIDs(int cellDim)
01468 {
01469 int myRank = comm().getRank();
01470
01471 SUNDANCE_VERB_MEDIUM("p=" << myRank << " assigning " << cellDim
01472 << "-cell owners");
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483 int numWaits = 1;
01484 if (staggerOutput())
01485 {
01486 SUNDANCE_VERB_LOW("p=" << myRank << " doing staggered output in "
01487 "assignIntermediateCellOwners()");
01488 numWaits = comm().getNProc();
01489 }
01490
01491
01492
01493
01494 Array<Array<int> > outgoingRequests(comm().getNProc());
01495 Array<Array<int> > incomingRequests(comm().getNProc());
01496
01497
01498 Array<Array<int> > outgoingRequestLIDs(comm().getNProc());
01499
01500
01501 Array<Array<int> > outgoingGIDs(comm().getNProc());
01502 Array<Array<int> > incomingGIDs(comm().getNProc());
01503
01504
01505 Array<int>& cellOwners = ownerProcID_[cellDim];
01506
01507
01508
01509 int localCount = 0;
01510 ArrayOfTuples<int>* cellVerts;
01511 if (cellDim==2) cellVerts = &(faceVertLIDs_);
01512 else cellVerts = &(edgeVerts_);
01513
01514 LIDToGIDMap_[cellDim].resize(numCells(cellDim));
01515
01516 SUNDANCE_OUT(this->verb() > 3,
01517 "starting loop over cells");
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538 try
01539 {
01540 resolveEdgeOwnership(cellDim);
01541 }
01542 catch(std::exception& e0)
01543 {
01544 SUNDANCE_TRACE_MSG(e0, "while resolving edge ownership");
01545 }
01546
01547
01548 for (int q=0; q<numWaits; q++)
01549 {
01550 if (numWaits > 1 && q != myRank)
01551 {
01552 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01553 "off-proc error detected on proc=" << myRank
01554 << " while computing GID resolution requests");
01555 continue;
01556 }
01557 try
01558 {
01559 for (int i=0; i<numCells(cellDim); i++)
01560 {
01561 SUNDANCE_OUT(this->verb() > 4,
01562 "p=" << myRank
01563 <<" cell " << i);
01564
01565
01566
01567
01568
01569
01570 int owner = cellOwners[i];
01571 if (owner==myRank)
01572 {
01573 SUNDANCE_VERB_EXTREME("p=" << myRank
01574 << " assigning GID=" << localCount
01575 << " to cell "
01576 << cellToStr(cellDim, i));
01577 LIDToGIDMap_[cellDim][i] = localCount++;
01578 }
01579 else
01580
01581
01582
01583
01584
01585
01586
01587
01588 {
01589 LIDToGIDMap_[cellDim][i] = -1;
01590 SUNDANCE_OUT(this->verb() > 4,
01591 "p=" << myRank << " adding cell LID=" << i
01592 << " to req list");
01593 outgoingRequestLIDs[owner].append(i);
01594 for (int v=0; v<=cellDim; v++)
01595 {
01596 int ptLID = cellVerts->value(i, v);
01597 int ptGID = LIDToGIDMap_[0][ptLID];
01598 SUNDANCE_OUT(this->verb() > 4,
01599 "p=" << myRank << " adding pt LID=" << ptLID
01600 << " GID=" << ptGID
01601 << " to req list");
01602 outgoingRequests[owner].append(ptGID);
01603 }
01604 SUNDANCE_VERB_EXTREME("p=" << myRank
01605 << " deferring GID assignment for cell "
01606 << cellToStr(cellDim, i));
01607 }
01608 }
01609 }
01610 catch(std::exception& e0)
01611 {
01612 reportFailure(comm());
01613 SUNDANCE_TRACE_MSG(e0, "while computing GID resolution requests");
01614 }
01615 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01616 "off-proc error detected on proc=" << myRank
01617 << " while computing GID resolution requests");
01618 }
01619
01620
01621
01622 try
01623 {
01624 synchronizeGIDNumbering(cellDim, localCount);
01625 }
01626 catch(std::exception& e0)
01627 {
01628 reportFailure(comm());
01629 SUNDANCE_TRACE_MSG(e0, "while computing GID offsets");
01630 }
01631 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01632 "off-proc error detected on proc=" << myRank
01633 << " while computing GID offsets");
01634
01635
01636 for (int q=0; q<numWaits; q++)
01637 {
01638 if (numWaits > 1) comm().synchronize();
01639 if (numWaits > 1 && q!=myRank) continue;
01640 SUNDANCE_OUT(this->verb() > 4,
01641 "p=" << myRank << "sending requests: "
01642 << outgoingRequests);
01643 }
01644
01645 try
01646 {
01647 MPIContainerComm<int>::allToAll(outgoingRequests,
01648 incomingRequests,
01649 comm());
01650 }
01651 catch(std::exception& e0)
01652 {
01653 reportFailure(comm());
01654 SUNDANCE_TRACE_MSG(e0, "while distributing GID resolution requests");
01655 }
01656 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01657 "off-proc error detected on proc=" << myRank
01658 << " while distributing GID resolution requests");
01659
01660
01661 for (int q=0; q<numWaits; q++)
01662 {
01663 if (numWaits > 1) comm().synchronize();
01664 if (numWaits > 1 && q!=myRank) continue;
01665 SUNDANCE_OUT(this->verb() > 4,
01666 "p=" << myRank << "recv'd requests: " << incomingRequests);
01667 }
01668
01669
01670 for (int q=0; q<numWaits; q++)
01671 {
01672 if (numWaits>1 && q != myRank)
01673 {
01674 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01675 "off-proc error detected on proc=" << myRank
01676 << " while computing edge/face GID responses");
01677 continue;
01678 }
01679 try
01680 {
01681
01682 for (int p=0; p<comm().getNProc(); p++)
01683 {
01684 const Array<int>& requestsFromProc = incomingRequests[p];
01685
01686 for (int c=0; c<requestsFromProc.size(); c+=(cellDim+1))
01687 {
01688 SUNDANCE_OUT(this->verb() > 4,
01689 "p=" << myRank << "processing request c="
01690 << c/(cellDim+1));
01691
01692
01693 int cellLID;
01694 if (cellDim==1)
01695 {
01696 int vertLID1 = mapGIDToLID(0, requestsFromProc[c]);
01697 int vertLID2 = mapGIDToLID(0, requestsFromProc[c+1]);
01698 SUNDANCE_VERB_HIGH("p=" << myRank
01699 << " edge vertex GIDs are "
01700 << requestsFromProc[c+1]
01701 << ", " << requestsFromProc[c]);
01702 cellLID = checkForExistingEdge(vertLID1, vertLID2);
01703 }
01704 else
01705 {
01706
01707 SUNDANCE_VERB_HIGH("p=" << myRank
01708 << " face vertex GIDs are "
01709 << requestsFromProc[c]
01710 << ", " << requestsFromProc[c+1]
01711 << ", " << requestsFromProc[c+2]);
01712 cellLID = lookupFace(tuple(requestsFromProc[c],
01713 requestsFromProc[c+1],
01714 requestsFromProc[c+2]));
01715 }
01716 SUNDANCE_VERB_HIGH("p=" << myRank << "cell LID is "
01717 << cellLID);
01718
01719 TEUCHOS_TEST_FOR_EXCEPTION(cellLID < 0, std::runtime_error,
01720 "unable to find requested cell "
01721 << cellStr(cellDim+1, &(requestsFromProc[c])));
01722
01723
01724
01725 int cellGID = mapLIDToGID(cellDim, cellLID);
01726 TEUCHOS_TEST_FOR_EXCEPTION(cellGID < 0, std::runtime_error,
01727 "proc " << myRank
01728 << " was asked by proc " << p
01729 << " to provide a GID for cell "
01730 << cellStr(cellDim+1, &(requestsFromProc[c]))
01731 << " with LID=" << cellLID
01732 << " but its GID is undefined");
01733
01734 outgoingGIDs[p].append(mapLIDToGID(cellDim, cellLID));
01735 }
01736 }
01737 }
01738 catch(std::exception& e0)
01739 {
01740 reportFailure(comm());
01741 SUNDANCE_TRACE_MSG(e0, "while computing edge/face GID responses");
01742 }
01743 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01744 "off-proc error detected "
01745 "on proc=" << myRank
01746 << " while computing edge/face GID responses");
01747 }
01748
01749
01750
01751 for (int q=0; q<numWaits; q++)
01752 {
01753 if (numWaits>1) comm().synchronize();
01754 if (numWaits>1 && q!=myRank) continue;
01755 SUNDANCE_OUT(this->verb() > 4,
01756 "p=" << myRank << "sending GIDs: " << outgoingGIDs);
01757 }
01758
01759
01760 try
01761 {
01762 MPIContainerComm<int>::allToAll(outgoingGIDs,
01763 incomingGIDs,
01764 comm());
01765 }
01766 catch(std::exception& e0)
01767 {
01768 reportFailure(comm());
01769 SUNDANCE_TRACE_MSG(e0, "while distributing edge/face GIDs");
01770 }
01771 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01772 "off-proc error detected on proc=" << myRank
01773 << " while distributing edge/face GIDs");
01774
01775
01776
01777
01778 for (int q=0; q<numWaits; q++)
01779 {
01780 if (numWaits > 1 && q != myRank)
01781 {
01782 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01783 "off-proc error detected on proc=" << myRank
01784 << " while setting remote edge/face GIDs");
01785 continue;
01786 }
01787 try
01788 {
01789 SUNDANCE_OUT(this->verb() > 4,
01790 "p=" << myRank << "recv'ing GIDs: " << incomingGIDs);
01791
01792
01793 for (int p=0; p<comm().getNProc(); p++)
01794 {
01795 const Array<int>& gidsFromProc = incomingGIDs[p];
01796 for (int c=0; c<gidsFromProc.size(); c++)
01797 {
01798 int cellLID = outgoingRequestLIDs[p][c];
01799 int cellGID = incomingGIDs[p][c];
01800 SUNDANCE_OUT(this->verb() > 4,
01801 "p=" << myRank <<
01802 " assigning GID=" << cellGID << " to LID="
01803 << cellLID);
01804 LIDToGIDMap_[cellDim][cellLID] = cellGID;
01805 GIDToLIDMap_[cellDim].put(cellGID, cellLID);
01806 }
01807 }
01808
01809
01810 SUNDANCE_OUT(this->verb() > 3,
01811 "p=" << myRank
01812 << "done synchronizing cells of dimension "
01813 << cellDim);
01814 }
01815 catch(std::exception& e0)
01816 {
01817 reportFailure(comm());
01818 SUNDANCE_TRACE_MSG(e0, " while setting remote edge/face GIDs");
01819 }
01820 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01821 "off-proc error detected on proc=" << myRank
01822 << " while setting remote edge/face GIDs");
01823 }
01824 if (cellDim==1) hasEdgeGIDs_ = true;
01825 else hasFaceGIDs_ = true;
01826 }
01827
01828
01829
01830
01831
01832
01833
01834 void BasicSimplicialMesh::synchronizeNeighborLists()
01835 {
01836 if (neighborsAreSynchronized_) return ;
01837
01838 Array<Array<int> > outgoingNeighborFlags(comm().getNProc());
01839 Array<Array<int> > incomingNeighborFlags(comm().getNProc());
01840
01841 int myRank = comm().getRank();
01842
01843 for (int p=0; p<comm().getNProc(); p++)
01844 {
01845 if (neighbors_.contains(p)) outgoingNeighborFlags[p].append(myRank);
01846 }
01847
01848 try
01849 {
01850 MPIContainerComm<int>::allToAll(outgoingNeighborFlags,
01851 incomingNeighborFlags,
01852 comm());
01853 }
01854 catch(std::exception& e0)
01855 {
01856 reportFailure(comm());
01857 SUNDANCE_TRACE_MSG(e0, "while distributing neighbor information");
01858 }
01859 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01860 "off-proc error detected on proc=" << myRank
01861 << " while distributing neighbor information");
01862
01863
01864
01865 for (int p=0; p<comm().getNProc(); p++)
01866 {
01867 if (incomingNeighborFlags[p].size() > 0) neighbors_.put(p);
01868 }
01869
01870 neighborsAreSynchronized_ = true;
01871 }
01872
01873
01874
01875
01876 void BasicSimplicialMesh::resolveEdgeOwnership(int cellDim)
01877 {
01878 int myRank = comm().getRank();
01879 SUNDANCE_VERB_LOW("p=" << myRank << " finding "
01880 << cellDim << "-cell owners");
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892 int numWaits = 1;
01893 if (staggerOutput())
01894 {
01895 SUNDANCE_VERB_LOW("p=" << myRank << " doing staggered output in "
01896 "assignIntermediateCellOwners()");
01897 numWaits = comm().getNProc();
01898 }
01899
01900 synchronizeNeighborLists();
01901
01902
01903 Array<int>& cellOwners = ownerProcID_[cellDim];
01904 ArrayOfTuples<int>* cellVerts;
01905 if (cellDim==2) cellVerts = &(faceVertLIDs_);
01906 else cellVerts = &(edgeVerts_);
01907
01908
01909 Array<int> neighbors = neighbors_.elements();
01910 SUNDANCE_VERB_LOW("p=" << myRank << " neighbors are " << neighbors);
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933 Array<int> reqIndexToEdgeLIDMap;
01934
01935
01936 Array<Array<int> > outgoingEdgeSpecifiers(comm().getNProc());
01937 Array<int> vertLID(cellDim+1);
01938 Array<int> vertGID(cellDim+1);
01939
01940 for (int q=0; q<numWaits; q++)
01941 {
01942 if (numWaits > 1 && q != myRank)
01943 {
01944 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
01945 "off-proc error detected on proc=" << myRank
01946 << " while determining edges to be resolved");
01947 continue;
01948 }
01949 try
01950 {
01951 for (int i=0; i<numCells(cellDim); i++)
01952 {
01953 int owner = cellOwners[i];
01954
01955
01956
01957
01958 if (owner==-1)
01959 {
01960
01961
01962 for (int v=0; v<=cellDim; v++)
01963 {
01964 vertLID[v] = cellVerts->value(i, v);
01965 vertGID[v] = LIDToGIDMap_[0][vertLID[v]];
01966 }
01967
01968 for (int n=0; n<neighbors.size(); n++)
01969 {
01970 for (int v=0; v<=cellDim; v++)
01971 {
01972 outgoingEdgeSpecifiers[neighbors[n]].append(vertGID[v]);
01973 }
01974 }
01975 SUNDANCE_VERB_HIGH("p=" << myRank
01976 << " is asking all neighbors to resolve edge "
01977 << cellToStr(cellDim, i));
01978
01979
01980
01981 reqIndexToEdgeLIDMap.append(i);
01982 }
01983 }
01984
01985 SUNDANCE_VERB_MEDIUM("p=" << myRank
01986 << " initially unresolved edge LIDs are "
01987 << reqIndexToEdgeLIDMap);
01988
01989 SUNDANCE_VERB_HIGH("p=" << myRank
01990 << " sending outgoing edge specs:") ;
01991 for (int p=0; p<comm().getNProc(); p++)
01992 {
01993 SUNDANCE_VERB_HIGH(outgoingEdgeSpecifiers[p].size()/(cellDim+1)
01994 << " edge reqs to proc "<< p) ;
01995 }
01996 SUNDANCE_VERB_HIGH("p=" << myRank << " outgoing edge specs are "
01997 << outgoingEdgeSpecifiers);
01998 }
01999 catch(std::exception& e0)
02000 {
02001 reportFailure(comm());
02002 SUNDANCE_TRACE_MSG(e0, "while determining edges to be resolved");
02003 }
02004 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02005 "off-proc error detected on proc=" << myRank
02006 << " while determining edges to be resolved");
02007 }
02008
02009
02010
02011
02012 Array<Array<int> > incomingEdgeSpecifiers(comm().getNProc());
02013
02014
02015 try
02016 {
02017 MPIContainerComm<int>::allToAll(outgoingEdgeSpecifiers,
02018 incomingEdgeSpecifiers,
02019 comm());
02020 }
02021 catch(std::exception& e0)
02022 {
02023 reportFailure(comm());
02024 SUNDANCE_TRACE_MSG(e0, "while distributing edge resolution reqs");
02025 }
02026 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02027 "off-proc error detected on proc=" << myRank
02028 << " while distributing edge resolution reqs");
02029
02030
02031 Array<Array<int> > outgoingEdgeContainers(comm().getNProc());
02032 for (int q=0; q<numWaits; q++)
02033 {
02034
02035 if (numWaits > 1 && q != myRank)
02036 {
02037 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02038 "off-proc error detected on proc=" << myRank
02039 << " while resolving edge ownership results");
02040 continue;
02041 }
02042
02043 try
02044 {
02045 SUNDANCE_VERB_HIGH("p=" << myRank << " incoming edge specs are "
02046 << incomingEdgeSpecifiers);
02047
02048
02049
02050
02051
02052
02053 for (int p=0; p<comm().getNProc(); p++)
02054 {
02055 const Array<int>& edgesFromProc = incomingEdgeSpecifiers[p];
02056
02057 int numVerts = cellDim+1;
02058 for (int c=0; c<edgesFromProc.size(); c+=numVerts)
02059 {
02060 SUNDANCE_VERB_LOW("p=" << myRank << " doing c=" << c/numVerts
02061 << " of " << edgesFromProc.size()/numVerts
02062 << " reqs from proc=" << p);
02063
02064
02065
02066
02067
02068
02069
02070
02071 int cellLID;
02072 SUNDANCE_VERB_LOW("p=" << myRank << " looking at "
02073 << cellStr(numVerts, &(edgesFromProc[c])));
02074
02075
02076
02077
02078
02079 int response = 0;
02080
02081 for (int v=0; v<numVerts; v++)
02082 {
02083 vertGID[v] = edgesFromProc[c+v];
02084 if (!GIDToLIDMap_[0].containsKey(vertGID[v]))
02085 {
02086 response = -1;
02087 break;
02088 }
02089 else
02090 {
02091 vertLID[v] = GIDToLIDMap_[0].get(vertGID[v]);
02092 }
02093 }
02094 if (response != -1)
02095 {
02096
02097
02098 if (cellDim==1)
02099 {
02100 cellLID = checkForExistingEdge(vertLID[0], vertLID[1]);
02101 }
02102 else
02103 {
02104 cellLID = lookupFace(tuple(vertGID[0], vertGID[1], vertGID[2]));
02105 }
02106
02107 if (cellLID==-1)
02108 {
02109 response = -1;
02110 }
02111 }
02112
02113 if (response == -1)
02114 {
02115 SUNDANCE_VERB_HIGH("p=" << myRank << " is unaware of cell "
02116 << cellStr(numVerts, &(vertGID[0])));
02117 }
02118 else
02119 {
02120 SUNDANCE_VERB_HIGH("p=" << myRank << " knows about cell "
02121 << cellStr(numVerts, &(vertGID[0])));
02122 if (ownerProcID_[cellDim][cellLID] != -1) response = 1;
02123 }
02124
02125 outgoingEdgeContainers[p].append(response);
02126 SUNDANCE_VERB_LOW("p=" << myRank << " did c=" << c/numVerts
02127 << " of " << edgesFromProc.size()/numVerts
02128 << " reqs from proc=" << p);
02129 }
02130 SUNDANCE_VERB_LOW("p=" << myRank << " done all reqs from proc "
02131 << p);
02132
02133 }
02134 SUNDANCE_VERB_LOW("p=" << myRank << " done processing edge reqs");
02135 SUNDANCE_VERB_HIGH("p=" << myRank
02136 << " outgoing edge containers are "
02137 << outgoingEdgeContainers);
02138 }
02139 catch(std::exception& e0)
02140 {
02141 reportFailure(comm());
02142 SUNDANCE_TRACE_MSG(e0, "while resolving edge ownership results");
02143 }
02144 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02145 "off-proc error detected on proc=" << myRank
02146 << " while resolving edge ownership results");
02147 }
02148
02149
02150
02151
02152
02153 Array<Array<int> > incomingEdgeContainers(comm().getNProc());
02154 try
02155 {
02156 MPIContainerComm<int>::allToAll(outgoingEdgeContainers,
02157 incomingEdgeContainers,
02158 comm());
02159 }
02160 catch(std::exception& e0)
02161 {
02162 reportFailure(comm());
02163 SUNDANCE_TRACE_MSG(e0, "while distributing edge resolution results");
02164 }
02165 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02166 "off-proc error detected on proc=" << myRank
02167 << " while distributing edge ownership results");
02168
02169
02170
02171
02172
02173
02174 for (int q=0; q<numWaits; q++)
02175 {
02176 if (numWaits > 1 && q != myRank)
02177 {
02178 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02179 "off-proc error detected on proc=" << myRank
02180 << " while finalizing edge ownership");
02181 continue;
02182 }
02183
02184 try
02185 {
02186 SUNDANCE_VERB_HIGH("p=" << myRank
02187 << " incoming edge containers are "
02188 << incomingEdgeContainers);
02189 Set<int> resolved;
02190 for (int p=0; p<comm().getNProc(); p++)
02191 {
02192 const Array<int>& edgeContainers = incomingEdgeContainers[p];
02193
02194 for (int c=0; c<edgeContainers.size(); c++)
02195 {
02196 int edgeLID = reqIndexToEdgeLIDMap[c];
02197
02198
02199
02200 if (edgeContainers[c] == 0 && !resolved.contains(edgeLID))
02201 {
02202
02203
02204 cellOwners[edgeLID] = p;
02205 }
02206
02207
02208
02209 else if (edgeContainers[c]==1)
02210 {
02211 cellOwners[edgeLID] = p;
02212 resolved.put(edgeLID);
02213 }
02214 }
02215 }
02216
02217
02218 for (int c=0; c<reqIndexToEdgeLIDMap.size(); c++)
02219 {
02220 int edgeLID = reqIndexToEdgeLIDMap[c];
02221 if (cellOwners[edgeLID] < 0) cellOwners[edgeLID] = myRank;
02222 SUNDANCE_VERB_EXTREME("p=" << myRank
02223 << " has decided cell "
02224 << cellToStr(cellDim, edgeLID)
02225 << " is owned by proc="
02226 << cellOwners[edgeLID]);
02227 }
02228 SUNDANCE_VERB_LOW("p=" << myRank << " done finalizing ownership");
02229 }
02230 catch(std::exception& e0)
02231 {
02232 reportFailure(comm());
02233 SUNDANCE_TRACE_MSG(e0, "while finalizing edge ownership");
02234 }
02235 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
02236 "off-proc error detected on proc=" << myRank
02237 << " while finalizing edge ownership");
02238 }
02239 }
02240
02241 string BasicSimplicialMesh::cellStr(int dim, const int* verts) const
02242 {
02243 std::string rtn="(";
02244 for (int i=0; i<dim; i++)
02245 {
02246 if (i!=0) rtn += ", ";
02247 rtn += Teuchos::toString(verts[i]);
02248 }
02249 rtn += ")";
02250 return rtn;
02251 }
02252
02253 string BasicSimplicialMesh::cellToStr(int cellDim, int cellLID) const
02254 {
02255 TeuchosOStringStream os;
02256
02257 if (cellDim > 0)
02258 {
02259 const ArrayOfTuples<int>* cellVerts;
02260 if (cellDim==spatialDim())
02261 {
02262 cellVerts = &(elemVerts_);
02263 }
02264 else
02265 {
02266 if (cellDim==2) cellVerts = &(faceVertLIDs_);
02267 else cellVerts = &(edgeVerts_);
02268 }
02269 os << "(";
02270 for (int j=0; j<cellVerts->tupleSize(); j++)
02271 {
02272 if (j!=0) os << ", ";
02273 os << mapLIDToGID(0, cellVerts->value(cellLID,j));
02274 }
02275 os << ")" << std::endl;
02276 }
02277 else
02278 {
02279 os << LIDToGIDMap_[0][cellLID] << std::endl;
02280 }
02281 return os.str();
02282 }
02283
02284 string BasicSimplicialMesh::printCells(int cellDim) const
02285 {
02286 TeuchosOStringStream os;
02287
02288 if (cellDim > 0)
02289 {
02290 const ArrayOfTuples<int>* cellVerts;
02291 if (cellDim==spatialDim())
02292 {
02293 cellVerts = &(elemVerts_);
02294 }
02295 else
02296 {
02297 if (cellDim==2) cellVerts = &(faceVertLIDs_);
02298 else cellVerts = &(edgeVerts_);
02299 }
02300
02301 os << "printing cells for processor " << comm().getRank() << std::endl;
02302 for (int i=0; i<cellVerts->length(); i++)
02303 {
02304 os << i << " (";
02305 for (int j=0; j<cellVerts->tupleSize(); j++)
02306 {
02307 if (j!=0) os << ", ";
02308 os << mapLIDToGID(0, cellVerts->value(i,j));
02309 }
02310 os << ")" << std::endl;
02311 }
02312 }
02313 else
02314 {
02315 os << LIDToGIDMap_[0] << std::endl;
02316 }
02317
02318 return os.str();
02319 }
02320
02321
02322