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 "SundanceMesh.hpp"
00032 #include "PlayaExceptions.hpp"
00033 #include "SundanceVerboseFieldWriter.hpp"
00034 #include "SundanceFieldWriter.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "SundanceSet.hpp"
00037 #include "SundanceMap.hpp"
00038 #include "PlayaTabs.hpp"
00039 #include "PlayaMPIContainerComm.hpp"
00040
00041 using namespace Teuchos;
00042 using namespace Sundance;
00043 using Playa::Handle;
00044 using Playa::Handleable;
00045 using Playa::MPIComm;
00046 using Playa::MPIContainerComm;
00047
00048
00049 IncrementallyCreatableMesh* Mesh::creatableMesh()
00050 {
00051 IncrementallyCreatableMesh* mci
00052 = dynamic_cast<IncrementallyCreatableMesh*>(ptr().get());
00053 TEUCHOS_TEST_FOR_EXCEPTION(mci==0, std::runtime_error,
00054 "Mesh::creatableMesh() could not convert mesh to "
00055 "a type deriving from IncrementallyCreatableMesh");
00056
00057 return mci;
00058 }
00059
00060
00061 void Mesh::dump(const std::string& filename) const
00062 {
00063 FieldWriter w = new VerboseFieldWriter(filename);
00064 w.addMesh(*this);
00065 w.write();
00066 }
00067
00068 bool Mesh::checkConsistency(const std::string& filename) const
00069 {
00070 std::string f = filename;
00071
00072 if (comm().getNProc() > 1)
00073 {
00074 f = f + "." + Teuchos::toString(comm().getRank());
00075 }
00076 std::ofstream os(f.c_str());
00077 return checkConsistency(os);
00078 }
00079
00080 bool Mesh::checkConsistency(std::ostream& os) const
00081 {
00082
00083 if (comm().getNProc()==1)
00084 {
00085 os << "Mesh is serial, thus it is consistent across processors" << std::endl;
00086 return true;
00087 }
00088
00089
00090 if (spatialDim() >= 2)
00091 const_cast<Mesh*>(this)->assignIntermediateCellGIDs(1);
00092 if (spatialDim() >= 3)
00093 const_cast<Mesh*>(this)->assignIntermediateCellGIDs(2);
00094 bool isOK = checkVertexConsistency(os);
00095 for (int d=1; d<=spatialDim(); d++)
00096 {
00097 isOK = checkCellConsistency(os, d) && isOK;
00098 }
00099
00100 return isOK;
00101 }
00102
00103
00104 bool Mesh::checkCellConsistency(std::ostream& os, int dim) const
00105 {
00106 TEUCHOS_TEST_FOR_EXCEPTION(dim==0, std::runtime_error,
00107 "Mesh::checkCellConsistency called for d=0. "
00108 "checkVertexConsistency() should be called instead");
00109
00110 int myRank = comm().getRank();
00111 int nProcs = comm().getNProc();
00112 Array<int> nFacets(dim);
00113
00114 bool elemsAreOK = true;
00115
00116
00117
00118
00119 int dataSize = 2;
00120 for (int d=0; d<dim; d++)
00121 {
00122 nFacets[d] = numFacets(dim, 0, d);
00123 dataSize += 2*nFacets[d];
00124 if (d > 0) dataSize += nFacets[d]*numFacets(d, 0, 0);
00125 }
00126
00127 int nCells = numCells(dim);
00128 Array<int> elemData(dataSize*nCells);
00129
00130 for (int c=0; c<nCells; c++)
00131 {
00132 elemData[dataSize*c] = mapLIDToGID(dim, c);
00133 elemData[dataSize*c + 1] = ownerProcID(dim, c);
00134 int base = dataSize*c + 2;
00135 int k=0;
00136 for (int d=0; d<dim; d++)
00137 {
00138 for (int f=0; f<nFacets[d]; f++)
00139 {
00140 int orientation;
00141 int lid = facetLID(dim, c, d, f, orientation);
00142 int facetGID = mapLIDToGID(d, lid);
00143 int facetOwner = ownerProcID(d, lid);
00144 elemData[base + k++] = facetGID;
00145 elemData[base + k++] = facetOwner;
00146 if (d > 0)
00147 {
00148 int nNodes = numFacets(d, lid, 0);
00149 for (int v=0; v<nNodes; v++)
00150 {
00151 int nodeLID = facetLID(d, lid, 0, v, orientation);
00152 int nodeGID = mapLIDToGID(0, nodeLID);
00153 elemData[base + k++] = nodeGID;
00154 }
00155 }
00156 }
00157 }
00158 }
00159
00160
00161
00162 Array<Array<int> > outData(nProcs);
00163 for (int p=0; p<nProcs; p++)
00164 {
00165 outData[p] = elemData;
00166 }
00167
00168 Array<Array<int> > inData(nProcs);
00169
00170 MPIContainerComm<int>::allToAll(outData, inData, comm());
00171
00172 for (int p=0; p<nProcs; p++)
00173 {
00174 Tabs tab1;
00175 if (p==myRank) continue;
00176
00177 os << tab1 << "p=" << myRank << " testing " << dim
00178 << "-cells from remote p=" << p << std::endl;
00179
00180 const Array<int>& remoteData = inData[p];
00181 int nRemote = remoteData.size()/dataSize;
00182
00183 for (int c=0; c<nRemote; c++)
00184 {
00185 Tabs tab2;
00186 int cellGID = remoteData[dataSize*c];
00187 int cellOwner = remoteData[dataSize*c+1];
00188 int cellLID = -1;
00189 bool elemIsOK = checkRemoteEntity(os, p, dim, cellGID,
00190 cellOwner, false, cellLID);
00191 if (cellLID >= 0 && elemIsOK)
00192 {
00193 int base = dataSize*c + 2;
00194 int k = 0;
00195 for (int d=0; d<dim; d++)
00196 {
00197 Tabs tab3;
00198 int dir;
00199 os << tab3 << "checking " << d << "-facets" << std::endl;
00200
00201
00202
00203
00204
00205
00206
00207 Sundance::Map<int, int> remoteFacetToOwnerMap;
00208 Sundance::Map<int, int> localFacetToOwnerMap;
00209 Sundance::Map<int, Set<int> > remoteFacetToNodesMap;
00210 Sundance::Map<int, Set<int> > localFacetToNodesMap;
00211
00212 for (int f=0; f<nFacets[d]; f++)
00213 {
00214
00215 int lfLID = facetLID(dim, cellLID, d, f, dir);
00216 int lfGID = mapLIDToGID(d, lfLID);
00217 int lfOwner = ownerProcID(d, lfLID);
00218 localFacetToOwnerMap.put(lfGID, lfOwner);
00219 int rfGID = remoteData[base + k++];
00220 int rfOwner = remoteData[base + k++];
00221 remoteFacetToOwnerMap.put(rfGID, rfOwner);
00222 if (d > 0)
00223 {
00224 int numNodes = numFacets(d, 0, 0);
00225 Sundance::Set<int> localNodes;
00226 Sundance::Set<int> remoteNodes;
00227 for (int v=0; v<numNodes; v++)
00228 {
00229 int nodeLID = facetLID(d, lfLID, 0, v, dir);
00230 int nodeGID = mapLIDToGID(0, nodeLID);
00231 localNodes.put(nodeGID);
00232 int remoteNodeGID = remoteData[base + k++];
00233 remoteNodes.put(remoteNodeGID);
00234 }
00235 localFacetToNodesMap.put(lfGID, localNodes);
00236 remoteFacetToNodesMap.put(rfGID, remoteNodes);
00237 }
00238 }
00239
00240
00241 for (Sundance::Map<int,int>::const_iterator
00242 rf=remoteFacetToOwnerMap.begin(),
00243 lf=localFacetToOwnerMap.begin();
00244 rf != remoteFacetToOwnerMap.end();
00245 rf++, lf++)
00246 {
00247 int lfGID = lf->first;
00248 int lfOwner = lf->second;
00249 int rfGID = lf->first;
00250 int rfOwner = lf->second;
00251 Tabs tab4;
00252 os << tab4 << " local facet GID=" << lfGID
00253 << " remote GID=" << rfGID
00254 << " local Own=" << lfOwner
00255 << " remote Own=" << rfOwner << std::endl;
00256 elemIsOK = testIdentity(os, rfGID, lfGID,
00257 "facet GIDs") && elemIsOK;
00258 elemIsOK = testIdentity(os, rfOwner,
00259 lfOwner,
00260 "facet owners") && elemIsOK;
00261
00262 if (d > 0)
00263 {
00264 Tabs tab5;
00265
00266 const Set<int>& localNodes =
00267 localFacetToNodesMap.get(lfGID);
00268 const Set<int>& remoteNodes =
00269 remoteFacetToNodesMap.get(rfGID);
00270 elemIsOK = testIdentity(os, localNodes.elements(),
00271 remoteNodes.elements(),
00272 "facet nodes") && elemIsOK;
00273 os << tab5 << "facet node GIDs: local="
00274 << localNodes << " remote="
00275 << remoteNodes << std::endl;
00276 }
00277 }
00278 }
00279 }
00280
00281 elemsAreOK = elemsAreOK && elemIsOK;
00282 }
00283 }
00284
00285 return elemsAreOK;
00286 }
00287
00288 bool Mesh::testIdentity(std::ostream& os, int a, int b, const std::string& msg) const
00289 {
00290 Tabs tab;
00291 if (a != b)
00292 {
00293 os << tab << "CONFLICT in " << msg << std::endl;
00294 return false;
00295 }
00296 return true;
00297 }
00298
00299 bool Mesh::testIdentity(std::ostream& os,
00300 const Array<int>& a,
00301 const Array<int>& b, const std::string& msg) const
00302 {
00303 Tabs tab;
00304 bool ok = true;
00305 for (int i=0; i<a.size(); i++)
00306 {
00307 if (a[i] != b[i]) ok = false;
00308 }
00309 if (!ok)
00310 {
00311 os << tab << "CONFLICT in " << msg << std::endl;
00312 }
00313 return ok;
00314 }
00315
00316 bool Mesh::checkRemoteEntity(std::ostream& os, int p, int dim, int gid, int owner,
00317 bool mustExist, int& lid) const
00318 {
00319 Tabs tab;
00320 bool isOK = true;
00321
00322 int myRank = comm().getRank();
00323 lid = -1;
00324
00325 if (hasGID(dim, gid))
00326 {
00327 lid = mapGIDToLID(dim, gid);
00328 os << tab << "p=" << myRank << " got "
00329 << dim << "-cell GID=" << gid << " from proc=" << p
00330 << ", is LID=" << lid << " locally" << std::endl;
00331 }
00332 else
00333 {
00334 os << tab << "p=" << myRank << " got " << dim << "-cell GID=" << gid
00335 << " from proc=" << p << ", doesn't exist locally" << std::endl;
00336 if (mustExist) isOK = false;
00337 }
00338
00339 if (lid >= 0)
00340 {
00341 int localOwner = ownerProcID(dim, lid);
00342 os << tab << dim << "-cell GID=" << gid
00343 << " is locally LID=" << lid << " and owned by "
00344 << localOwner << std::endl;
00345 if (localOwner != owner)
00346 {
00347 os << tab << "OWNERSHIP CONFLICT: local p=" << myRank
00348 << " thinks " << dim << "-cell GID=" << gid << " is owned by "
00349 << localOwner << " but remote proc=" << p
00350 << " says it's owned by " << owner << std::endl;
00351 isOK = false;
00352 }
00353 }
00354
00355 return isOK;
00356
00357 }
00358
00359
00360
00361 bool Mesh::checkVertexConsistency(std::ostream& os) const
00362 {
00363 int dim = spatialDim();
00364 int myRank = comm().getRank();
00365 int nProcs = comm().getNProc();
00366 std::string pHdr = "p=" + Teuchos::toString(myRank);
00367
00368 Out::println(pHdr + " testing consistency of vertices");
00369
00370 int dataSize = 2;
00371 Array<int> vertData(dataSize*numCells(0));
00372 Array<double> vertPositions(numCells(0)*dim);
00373 for (int v=0; v<numCells(0); v++)
00374 {
00375 vertData[dataSize*v] = mapLIDToGID(0, v);
00376 vertData[dataSize*v + 1] = ownerProcID(0, v);
00377 Point x = nodePosition(v);
00378 for (int j=0; j<dim; j++)
00379 {
00380 vertPositions[dim*v + j] = x[j];
00381 }
00382 }
00383
00384
00385
00386 Array<Array<int> > outData(nProcs);
00387 Array<Array<double> > outPos(nProcs);
00388 for (int p=0; p<nProcs; p++)
00389 {
00390 outData[p] = vertData;
00391 outPos[p] = vertPositions;
00392 }
00393
00394 Array<Array<int> > inData(nProcs);
00395 Array<Array<double> > inPos(nProcs);
00396
00397 MPIContainerComm<int>::allToAll(outData,
00398 inData,
00399 comm());
00400
00401 MPIContainerComm<double>::allToAll(outPos, inPos, comm());
00402
00403 double tol = 1.0e-15;
00404
00405 bool ok = true;
00406 bool allVertsOK = true;
00407
00408
00409 for (int p=0; p<nProcs; p++)
00410 {
00411 Tabs tab1;
00412 if (p==myRank) continue;
00413
00414 os << tab1 << "p=" << myRank << " testing vertices from remote p="
00415 << p << std::endl;
00416
00417 int nVerts = inData[p].size()/dataSize;
00418
00419 for (int v=0; v<nVerts; v++)
00420 {
00421 Tabs tab2;
00422 int vertGID = inData[p][v*dataSize];
00423 int vertOwner = inData[p][v*dataSize+1];
00424 int vertLID = -1;
00425 bool vertIsOK = checkRemoteEntity(os, p, 0, vertGID, vertOwner,
00426 false, vertLID);
00427
00428 if (vertLID >= 0)
00429 {
00430 Point localX = nodePosition(vertLID);
00431
00432 Point remoteX = localX;
00433 for (int j=0; j<dim; j++) remoteX[j] = inPos[p][dim*v + j];
00434 double dx = remoteX.distance(localX);
00435 if (dx > tol)
00436 {
00437 os << tab2 << "POSITION CONFLICT: local p=" << myRank
00438 << " thinks GID=" << vertGID << " is at xLocal="
00439 << localX << " but remote proc=" << p
00440 << " says it's at xRemote" << remoteX << std::endl;
00441 os << tab2 << "distance = " << dx << std::endl;
00442 vertIsOK = false;
00443 }
00444 }
00445 allVertsOK = vertIsOK && allVertsOK;
00446 if (vertIsOK)
00447 {
00448 os << tab2 << "p=" << myRank
00449 << " says vertex GID=" << vertGID << " is OK" << std::endl;
00450 }
00451 }
00452 }
00453
00454 if (allVertsOK)
00455 {
00456 os << "p=" << myRank << " found all vertex data is OK" << std::endl;
00457 }
00458 else
00459 {
00460 os << "p=" << myRank << " detected conflicts in vertex data" << std::endl;
00461 }
00462
00463 ok = allVertsOK && ok;
00464
00465
00466
00467 return ok;
00468 }
00469
