00001 #include "SundanceSerialPartitionerBase.hpp"
00002 #include "SundanceMap.hpp"
00003 #include "SundanceBasicSimplicialMeshType.hpp"
00004
00005 using namespace Sundance;
00006 using namespace Sundance;
00007
00008 using namespace Teuchos;
00009 using namespace Sundance;
00010
00011 using std::max;
00012 using std::min;
00013 using std::endl;
00014 using std::cout;
00015 using std::cerr;
00016
00017
00018
00019 Set<int> SerialPartitionerBase::arrayToSet(const Array<int>& a) const
00020 {
00021 Set<int> rtn;
00022 for (int i=0; i<a.size(); i++)
00023 {
00024 rtn.put(a[i]);
00025 }
00026 return rtn;
00027 }
00028
00029 int SerialPartitionerBase::max(const Set<int>& s) const
00030 {
00031 return *(s.rbegin());
00032 }
00033
00034 void SerialPartitionerBase::getNeighbors(const Mesh& mesh,
00035 Array<Array<int> >& neighbors, int& nEdges) const
00036 {
00037 int dim = mesh.spatialDim();
00038 int nElems = mesh.numCells(dim);
00039 int nVerts = mesh.numCells(0);
00040
00041 neighbors.resize(nElems);
00042 nEdges = 0;
00043
00044 elemVerts_.resize(nElems);
00045 elemEdgewiseNbors_.resize(nElems);
00046 vertElems_.resize(nVerts);
00047
00048
00049 for (int c=0; c<nElems; c++)
00050 {
00051 Array<int> facetLID;
00052 Array<int> facetDir;
00053
00054 mesh.getFacetArray(dim, c, 0, facetLID, facetDir);
00055 elemVerts_[c] = arrayToSet(facetLID);
00056
00057
00058 facetLID.resize(0);
00059 facetDir.resize(0);
00060 mesh.getFacetArray(dim, c, 1, facetLID, facetDir);
00061 Set<int> adjElems;
00062 for (int f=0; f<facetLID.size(); f++)
00063 {
00064 Array<int> maxCofs;
00065 mesh.getCofacets(1, facetLID[f], dim, maxCofs);
00066 for (int cf=0; cf<maxCofs.size(); cf++)
00067 {
00068 if (maxCofs[cf] != c) adjElems.put(maxCofs[cf]);
00069 }
00070 }
00071 elemEdgewiseNbors_[c] = adjElems;
00072 }
00073
00074 for (int v=0; v<nVerts; v++)
00075 {
00076 Array<int> cofacetLID;
00077 mesh.getCofacets(0, v, dim, cofacetLID);
00078 vertElems_[v] = arrayToSet(cofacetLID);
00079 }
00080
00081
00082 for (int i=0; i<nElems; i++)
00083 {
00084 Set<int> allNbors;
00085 const Set<int>& ev = elemVerts_[i];
00086 for (Set<int>::const_iterator v=ev.begin(); v!=ev.end(); v++)
00087 {
00088 allNbors = allNbors.setUnion(vertElems_[*v]);
00089 }
00090 allNbors.erase(i);
00091
00092 Array<int> fullNbors;
00093 for (Set<int>::const_iterator j=allNbors.begin(); j!=allNbors.end(); j++)
00094 {
00095 Set<int> numCommonNodes = elemVerts_[i].intersection(elemVerts_[*j]);
00096 if ((int) numCommonNodes.size() == dim) fullNbors.append(*j);
00097 }
00098 nEdges += fullNbors.size();
00099 neighbors[i] = fullNbors;
00100 }
00101
00102 nEdges = nEdges/2;
00103 }
00104
00105
00106 void SerialPartitionerBase::getNodeAssignments(int nProc,
00107 const Array<int>& elemAssignments,
00108 Array<int>& nodeAssignments,
00109 Array<int>& nodeOwnerElems,
00110 Array<int>& nodesPerProc) const
00111 {
00112 nodesPerProc.resize(nProc);
00113 nodeAssignments.resize(vertElems_.size());
00114 nodeOwnerElems.resize(vertElems_.size());
00115
00116 for (int i=0; i<nodesPerProc.size(); i++) nodesPerProc[i] = 0;
00117
00118 for (int v=0; v<vertElems_.size(); v++)
00119 {
00120
00121 int ownerElem = max(vertElems_[v]);
00122 nodeOwnerElems[v] = ownerElem;
00123 nodeAssignments[v] = elemAssignments[ownerElem];
00124 nodesPerProc[nodeAssignments[v]]++;
00125 }
00126 }
00127
00128 void SerialPartitionerBase::getElemsPerProc(int nProc,
00129 const Array<int>& elemAssignments,
00130 Array<int>& elemsPerProc) const
00131 {
00132 elemsPerProc.resize(nProc);
00133 for (int i=0; i<elemsPerProc.size(); i++) elemsPerProc[i] = 0;
00134
00135 for (int e=0; e<elemAssignments.size(); e++)
00136 {
00137 elemsPerProc[elemAssignments[e]]++;
00138 }
00139 }
00140
00141 void SerialPartitionerBase::remapEntities(const Array<int>& assignments,
00142 int nProc,
00143 Array<int>& entityMap) const
00144 {
00145 Array<Array<int> > procBuckets(nProc);
00146 entityMap.resize(assignments.size());
00147
00148 for (int i=0; i<assignments.size(); i++)
00149 {
00150 procBuckets[assignments[i]].append(i);
00151 }
00152
00153 int g = 0;
00154
00155 for (int p=0; p<nProc; p++)
00156 {
00157 for (int i=0; i<procBuckets[p].size(); i++)
00158 {
00159 int lid = procBuckets[p][i];
00160 entityMap[lid] = g++;
00161 }
00162 }
00163 }
00164
00165
00166 void SerialPartitionerBase::getOffProcData(int p,
00167 const Array<int>& elemAssignments,
00168 const Array<int>& nodeAssignments,
00169 Set<int>& offProcNodes,
00170 Set<int>& offProcElems) const
00171 {
00172 Out::root() << "ignore ghosts = " << ignoreGhosts_ << endl;
00173
00174 for (int e=0; e<elemAssignments.size(); e++)
00175 {
00176 if (elemAssignments[e] != p) continue;
00177 for (Set<int>::const_iterator v=elemVerts_[e].begin(); v!=elemVerts_[e].end(); v++)
00178 {
00179 if (nodeAssignments[*v]!=p) offProcNodes.put(*v);
00180 }
00181 if (!ignoreGhosts_)
00182 {
00183 for (Set<int>::const_iterator
00184 n=elemEdgewiseNbors_[e].begin(); n!=elemEdgewiseNbors_[e].end(); n++)
00185 {
00186 if (elemAssignments[*n]!=p) offProcElems.put(*n);
00187 }
00188 }
00189 }
00190
00191 if (!ignoreGhosts_)
00192 {
00193
00194 for (int v=0; v<nodeAssignments.size(); v++)
00195 {
00196 if (nodeAssignments[v] != p) continue;
00197 const Set<int>& v2e = vertElems_[v];
00198 for (Set<int>::const_iterator e=v2e.begin(); e!=v2e.end(); e++)
00199 {
00200 if (elemAssignments[*e] != p) offProcElems.put(*e);
00201 }
00202 }
00203
00204
00205
00206 for (Set<int>::const_iterator e=offProcElems.begin(); e!=offProcElems.end(); e++)
00207 {
00208 for (Set<int>::const_iterator v=elemVerts_[*e].begin(); v!=elemVerts_[*e].end(); v++)
00209 {
00210 if (nodeAssignments[*v]!=p) offProcNodes.put(*v);
00211 }
00212 }
00213 }
00214 }
00215
00216
00217 Array<Mesh> SerialPartitionerBase::makeMeshParts(const Mesh& mesh, int np,
00218 Array<Sundance::Map<int, int> >& oldElemLIDToNewLIDMaps,
00219 Array<Sundance::Map<int, int> >& oldVertLIDToNewLIDMaps) const
00220 {
00221 int dim = mesh.spatialDim();
00222
00223 Array<int> elemAssignments;
00224 Array<int> nodeAssignments;
00225 Array<int> nodeOwnerElems;
00226 Array<int> nodesPerProc;
00227
00228 getAssignments(mesh, np, elemAssignments);
00229
00230 getNodeAssignments(np, elemAssignments, nodeAssignments, nodeOwnerElems,
00231 nodesPerProc);
00232
00233 Array<Array<int> > offProcNodes(np);
00234 Array<Array<int> > offProcElems(np);
00235
00236 Array<Set<int> > offProcNodeSet(np);
00237 Array<Set<int> > offProcElemSet(np);
00238
00239 Array<Array<int> > procNodes(np);
00240 Array<Array<int> > procElems(np);
00241
00242 for (int p=0; p<np; p++)
00243 {
00244 getOffProcData(p, elemAssignments, nodeAssignments,
00245 offProcNodeSet[p], offProcElemSet[p]);
00246 offProcNodes[p] = offProcNodeSet[p].elements();
00247 offProcElems[p] = offProcElemSet[p].elements();
00248 procElems[p].reserve(elemAssignments.size()/np);
00249 procNodes[p].reserve(nodeAssignments.size()/np);
00250 }
00251
00252 Array<int> remappedElems;
00253 Array<int> remappedNodes;
00254
00255 remapEntities(elemAssignments, np, remappedElems);
00256 remapEntities(nodeAssignments, np, remappedNodes);
00257
00258 for (int e=0; e<elemAssignments.size(); e++)
00259 {
00260 procElems[elemAssignments[e]].append(e);
00261 }
00262
00263 for (int n=0; n<nodeAssignments.size(); n++)
00264 {
00265 procNodes[nodeAssignments[n]].append(n);
00266 }
00267
00268
00269 Array<Mesh> rtn(np);
00270 oldVertLIDToNewLIDMaps.resize(np);
00271 oldElemLIDToNewLIDMaps.resize(np);
00272 for (int p=0; p<np; p++)
00273 {
00274 Sundance::Map<int, int>& oldVertLIDToNewLIDMap
00275 = oldVertLIDToNewLIDMaps[p];
00276 Sundance::Map<int, int>& oldElemLIDToNewLIDMap
00277 = oldElemLIDToNewLIDMaps[p];
00278 MeshType type = new BasicSimplicialMeshType();
00279 rtn[p] = type.createEmptyMesh(mesh.spatialDim(), MPIComm::world());
00280
00281 Set<int> unusedVertGID;
00282
00283
00284 for (int n=0; n<procNodes[p].size(); n++)
00285 {
00286 int oldLID = procNodes[p][n];
00287 int newGID = remappedNodes[oldLID];
00288 unusedVertGID.put(newGID);
00289 int newLID = rtn[p].addVertex(newGID, mesh.nodePosition(oldLID), p, 0);
00290 oldVertLIDToNewLIDMap.put(oldLID, newLID);
00291 }
00292
00293
00294 for (int n=0; n<offProcNodes[p].size(); n++)
00295 {
00296 int oldLID = offProcNodes[p][n];
00297 int nodeOwnerProc = nodeAssignments[oldLID];
00298 int newGID = remappedNodes[oldLID];
00299 unusedVertGID.put(newGID);
00300 int newLID = rtn[p].addVertex(newGID, mesh.nodePosition(oldLID), nodeOwnerProc, 0);
00301 oldVertLIDToNewLIDMap.put(oldLID, newLID);
00302 }
00303
00304
00305 for (int e=0; e<procElems[p].size(); e++)
00306 {
00307 int oldLID = procElems[p][e];
00308 int newGID = remappedElems[oldLID];
00309 Array<int> vertGIDs;
00310 Array<int> orientations;
00311 mesh.getFacetArray(dim, oldLID, 0, vertGIDs, orientations);
00312 for (int v=0; v<vertGIDs.size(); v++)
00313 {
00314 vertGIDs[v] = remappedNodes[vertGIDs[v]];
00315 if (unusedVertGID.contains(vertGIDs[v])) unusedVertGID.erase(newGID);
00316 }
00317 int newLID = rtn[p].addElement(newGID, vertGIDs, p, 1);
00318 oldElemLIDToNewLIDMap.put(oldLID, newLID);
00319 }
00320
00321
00322 for (int e=0; e<offProcElems[p].size(); e++)
00323 {
00324 int oldLID = offProcElems[p][e];
00325 int newGID = remappedElems[oldLID];
00326 int elemOwnerProc = elemAssignments[oldLID];
00327 Array<int> vertGIDs;
00328 Array<int> orientations;
00329 mesh.getFacetArray(dim, oldLID, 0, vertGIDs, orientations);
00330 for (int v=0; v<vertGIDs.size(); v++)
00331 {
00332 vertGIDs[v] = remappedNodes[vertGIDs[v]];
00333 if (unusedVertGID.contains(vertGIDs[v])) unusedVertGID.erase(newGID);
00334 }
00335 int newLID = rtn[p].addElement(newGID, vertGIDs, elemOwnerProc, 1);
00336 oldElemLIDToNewLIDMap.put(oldLID, newLID);
00337
00338
00339
00340 }
00341
00342
00343
00344 for (int d=1; d<dim; d++)
00345 {
00346 Set<int> labels = mesh.getAllLabelsForDimension(d);
00347 for (Set<int>::const_iterator i=labels.begin(); i!=labels.end(); i++)
00348 {
00349 Array<int> labeledCells;
00350 int label = *i;
00351 if (label == 0) continue;
00352 mesh.getLIDsForLabel(d, label, labeledCells);
00353 for (int c=0; c<labeledCells.size(); c++)
00354 {
00355 int lid = labeledCells[c];
00356 Array<int> cofacets;
00357 mesh.getCofacets(d, lid, dim, cofacets);
00358 for (int n=0; n<cofacets.size(); n++)
00359 {
00360 int cofLID = cofacets[n];
00361 if (elemAssignments[cofLID]==p
00362 || offProcElemSet[p].contains(cofLID))
00363 {
00364
00365
00366
00367
00368 Array<int> oldVerts;
00369 Array<int> newVerts;
00370 Array<int> orientation;
00371 mesh.getFacetArray(d, lid, 0, oldVerts, orientation);
00372 for (int v=0; v<oldVerts.size(); v++)
00373 {
00374 newVerts.append(remappedNodes[oldVerts[v]]);
00375 }
00376
00377
00378 Set<int> newVertSet = arrayToSet(newVerts);
00379
00380
00381 int newElemLID = oldElemLIDToNewLIDMap.get(cofLID);
00382
00383
00384 Array<int> submeshFacets;
00385 rtn[p].getFacetArray(dim, newElemLID, d,
00386 submeshFacets, orientation);
00387 int facetIndex = -1;
00388 for (int df=0; df<submeshFacets.size(); df++)
00389 {
00390
00391 int facetLID = submeshFacets[df];
00392 Array<int> verts;
00393 rtn[p].getFacetArray(d, facetLID, 0, verts, orientation);
00394 Array<int> vertGID(verts.size());
00395 for (int v=0; v<verts.size(); v++)
00396 {
00397 vertGID[v] = rtn[p].mapLIDToGID(0, verts[v]);
00398 }
00399 Set<int> subVertSet = arrayToSet(vertGID);
00400 if (subVertSet==newVertSet)
00401 {
00402 facetIndex = df;
00403 break;
00404 }
00405 }
00406 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex==-1, std::logic_error,
00407 "couldn't match new " << d << "-cell in submesh to old " << d
00408 << "cell. This should never happen");
00409
00410
00411
00412
00413 int o;
00414 int newFacetLID = rtn[p].facetLID(dim, newElemLID, d,
00415 facetIndex, o);
00416
00417 rtn[p].setLabel(d, newFacetLID, label);
00418 break;
00419
00420 }
00421 else
00422 {
00423 }
00424 }
00425 }
00426 }
00427 }
00428 }
00429
00430 return rtn;
00431 }
00432
00433