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 "SundanceHomogeneousDOFMap.hpp"
00032 #include "SundanceCellFilter.hpp"
00033 #include "SundanceMaximalCellFilter.hpp"
00034 #include "PlayaMPIContainerComm.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "PlayaTabs.hpp"
00037 #include "Teuchos_Time.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042 using Playa::MPIComm;
00043 using Playa::MPIContainerComm;
00044
00045
00046 static Time& dofLookupTimer()
00047 {
00048 static RCP<Time> rtn
00049 = TimeMonitor::getNewTimer("unbatched dof lookup");
00050 return *rtn;
00051 }
00052
00053 static Time& dofBatchLookupTimer()
00054 {
00055 static RCP<Time> rtn
00056 = TimeMonitor::getNewTimer("batched dof lookup");
00057 return *rtn;
00058 }
00059
00060 HomogeneousDOFMap::HomogeneousDOFMap(const Mesh& mesh,
00061 const BasisFamily& basis,
00062 int numFuncs,
00063 int setupVerb)
00064 : DOFMapBase(mesh, setupVerb),
00065 dim_(mesh.spatialDim()),
00066 dofs_(mesh.spatialDim()+1),
00067 maximalDofs_(),
00068 haveMaximalDofs_(false),
00069 localNodePtrs_(mesh.spatialDim()+1),
00070 nNodesPerCell_(mesh.spatialDim()+1),
00071 totalNNodesPerCell_(mesh.spatialDim()+1, 0),
00072 numFacets_(mesh.spatialDim()+1),
00073 originalFacetOrientation_(2),
00074 basisIsContinuous_(false)
00075 {
00076 verbosity() = DOFMapBase::classVerbosity();
00077
00078 CellFilter maximalCells = new MaximalCellFilter();
00079 cellSets().append(maximalCells.getCells(mesh));
00080 cellDimOnCellSets().append(mesh.spatialDim());
00081
00082 allocate(mesh, basis, numFuncs);
00083 initMap();
00084 }
00085
00086
00087 HomogeneousDOFMap::HomogeneousDOFMap(const Mesh& mesh,
00088 const BasisFamily& basis,
00089 const Array<CellFilter>& subregions,
00090 int numFuncs)
00091 : DOFMapBase(mesh),
00092 dim_(mesh.spatialDim()),
00093 dofs_(mesh.spatialDim()+1),
00094 maximalDofs_(),
00095 haveMaximalDofs_(false),
00096 localNodePtrs_(mesh.spatialDim()+1),
00097 nNodesPerCell_(mesh.spatialDim()+1),
00098 totalNNodesPerCell_(mesh.spatialDim()+1, 0),
00099 numFacets_(mesh.spatialDim()+1),
00100 originalFacetOrientation_(2),
00101 basisIsContinuous_(false)
00102 {
00103 verbosity() = DOFMapBase::classVerbosity();
00104
00105 for (int r=0; r<subregions.size(); r++)
00106 {
00107 cellSets().append(subregions[r].getCells(mesh));
00108 cellDimOnCellSets().append(subregions[r].dimension(mesh));
00109 }
00110
00111 allocate(mesh, basis, numFuncs);
00112 initMap();
00113 }
00114
00115
00116 void HomogeneousDOFMap::allocate(const Mesh& mesh,
00117 const BasisFamily& basis,
00118 int numFuncs)
00119 {
00120 Tabs tab;
00121 SUNDANCE_MSG1(setupVerb(), tab << "allocating DOF map for nFuncs=" << numFuncs);
00122 Array<int> fid(numFuncs);
00123 for (int f=0; f<numFuncs; f++) fid[f] = f;
00124 funcIDOnCellSets().append(fid);
00125
00126
00127
00128 for (int d=0; d<=dim_; d++)
00129 {
00130 Tabs tab1;
00131 SUNDANCE_MSG2(setupVerb(), tab1 << "allocating d=" << d);
00132
00133
00134 numFacets_[d].resize(d);
00135 for (int fd=0; fd<d; fd++) numFacets_[d][fd]=mesh.numFacets(d, 0, fd);
00136 SUNDANCE_MSG3(setupVerb(), tab1 << "num facets for dimension " << d << " is "
00137 << numFacets_[d]);
00138
00139
00140
00141 basis.ptr()->getLocalDOFs(mesh.cellType(d), localNodePtrs_[d]);
00142
00143
00144 SUNDANCE_MSG3(setupVerb(), tab1 << "node ptrs for dimension " << d << " are "
00145 << localNodePtrs_[d]);
00146
00147
00148
00149 if (localNodePtrs_[d][d].size() > 0)
00150 {
00151 nNodesPerCell_[d] = localNodePtrs_[d][d][0].size();
00152 }
00153 else
00154 {
00155 nNodesPerCell_[d] = 0;
00156 }
00157 SUNDANCE_MSG3(setupVerb(), tab1 <<
00158 "num nodes for dimension " << d << " is "
00159 << nNodesPerCell_[d]);
00160
00161 totalNNodesPerCell_[d] = nNodesPerCell_[d];
00162 for (int dd=0; dd<d; dd++)
00163 {
00164 totalNNodesPerCell_[d] += numFacets_[d][dd]*nNodesPerCell_[dd];
00165 }
00166
00167
00168 if (nNodesPerCell_[d] > 0)
00169 {
00170 dofs_[d].resize(mesh.numCells(d));
00171 }
00172 else
00173 {
00174 dofs_[d].resize(0);
00175 }
00176
00177 if (d > 0 && d < dim_) originalFacetOrientation_[d-1].resize(mesh.numCells(d));
00178
00179
00180
00181 if (d < dim_ && nNodesPerCell_[d] > 0) basisIsContinuous_ = true;
00182
00183
00184
00185
00186 int numCells = dofs_[d].size();
00187 for (int c=0; c<numCells; c++)
00188 {
00189 dofs_[d][c].resize(funcIDList().size() * nNodesPerCell_[d]);
00190
00191 for (int i=0; i<dofs_[d][c].size(); i++)
00192 {
00193 dofs_[d][c][i] = uninitializedVal();
00194 }
00195 }
00196 }
00197 SUNDANCE_MSG1(setupVerb(), tab << "done allocating DOF map");
00198 }
00199
00200 void HomogeneousDOFMap::initMap()
00201 {
00202 Tabs tab;
00203 SUNDANCE_MSG1(setupVerb(), tab << "initializing DOF map");
00204
00205 int nextDOF = 0;
00206
00207
00208
00209
00210 Array<Array<Array<int> > > remoteCells(mesh().spatialDim()+1,
00211 mesh().comm().getNProc());
00212
00213 for (int r=0; r<numCellSets(); r++)
00214 {
00215
00216
00217
00218
00219
00220
00221
00222
00223 CellSet cells = cellSet(r);
00224 CellIterator iter;
00225 for (iter=cells.begin(); iter != cells.end(); iter++)
00226 {
00227
00228 int cellLID = *iter;
00229 int owner;
00230
00231 if (nNodesPerCell_[dim_] > 0)
00232 {
00233
00234
00235
00236 if (isRemote(dim_, cellLID, owner))
00237 {
00238 int dummy=0;
00239 int cellGID = mesh().mapLIDToGID(dim_, cellLID);
00240 remoteCells[dim_][owner].append(cellGID);
00241 setDOFs(dim_, cellLID, dummy);
00242 }
00243 else
00244
00245 {
00246 setDOFs(dim_, cellLID, nextDOF);
00247 }
00248 }
00249
00250
00251
00252
00253 if (basisIsContinuous_)
00254 {
00255 for (int d=0; d<dim_; d++)
00256 {
00257 if (nNodesPerCell_[d] > 0)
00258 {
00259 int nf = numFacets_[dim_][d];
00260 Array<int> facetLID(nf);
00261 Array<int> facetOrientations(nf);
00262
00263 mesh().getFacetArray(dim_, cellLID, d,
00264 facetLID, facetOrientations);
00265
00266 for (int f=0; f<nf; f++)
00267 {
00268
00269
00270 if (!hasBeenAssigned(d, facetLID[f]))
00271 {
00272
00273 if (isRemote(d, facetLID[f], owner))
00274 {
00275 int dummy=0;
00276 int facetGID = mesh().mapLIDToGID(d, facetLID[f]);
00277 remoteCells[d][owner].append(facetGID);
00278 setDOFs(d, facetLID[f], dummy);
00279 }
00280 else
00281 {
00282
00283 setDOFs(d, facetLID[f], nextDOF);
00284
00285 if (d > 0)
00286 {
00287 originalFacetOrientation_[d-1][facetLID[f]]
00288 = facetOrientations[f];
00289 }
00290 }
00291 }
00292 }
00293 }
00294 }
00295 }
00296 }
00297 }
00298
00299
00300
00301
00302 if (mesh().comm().getNProc() > 1)
00303 {
00304 for (int d=0; d<=dim_; d++)
00305 {
00306 if (nNodesPerCell_[d] > 0)
00307 {
00308 computeOffsets(d, nextDOF);
00309 shareDOFs(d, remoteCells[d]);
00310 }
00311 }
00312 }
00313 else
00314 {
00315 setLowestLocalDOF(0);
00316 setNumLocalDOFs(nextDOF);
00317 setTotalNumDOFs(nextDOF);
00318 }
00319 SUNDANCE_MSG1(setupVerb(), tab << "done initializing DOF map");
00320 }
00321
00322 void HomogeneousDOFMap::shareDOFs(int cellDim,
00323 const Array<Array<int> >& outgoingCellRequests)
00324 {
00325 int np = mesh().comm().getNProc();
00326 int rank = mesh().comm().getRank();
00327
00328 if (np == 1) return;
00329
00330 TEUCHOS_TEST_FOR_EXCEPTION(np != outgoingCellRequests.size(), std::runtime_error,
00331 "incorrect size of outgoingCellRequests array. Size is "
00332 << outgoingCellRequests.size() << ", should be np=" << np);
00333
00334
00335
00336 Array<Array<int> > incomingCellRequests;
00337 Array<Array<int> > outgoingDOFs(np);
00338 Array<Array<int> > incomingDOFs;
00339
00340 SUNDANCE_MSG2(setupVerb(),
00341 "p=" << mesh().comm().getRank()
00342 << "synchronizing DOFs for cells of dimension " << cellDim);
00343 SUNDANCE_MSG2(setupVerb(),
00344 "p=" << mesh().comm().getRank()
00345 << " sending cell reqs d=" << cellDim << " GID=" << outgoingCellRequests);
00346
00347
00348 MPIContainerComm<int>::allToAll(outgoingCellRequests,
00349 incomingCellRequests,
00350 mesh().comm());
00351
00352
00353
00354 SUNDANCE_MSG2(setupVerb(),
00355 "p=" << rank
00356 << "recvd DOF requests for cells of dimension " << cellDim
00357 << " GID=" << incomingCellRequests);
00358
00359
00360
00361 for (int p=0; p<np; p++)
00362 {
00363 if (p==rank) continue;
00364 const Array<int>& requestsFromProc = incomingCellRequests[p];
00365 int nReq = requestsFromProc.size();
00366 int blockSize = 1;
00367 if (cellDim > 0 && cellDim < dim_)
00368 {
00369 outgoingDOFs[p].resize(2 * nReq);
00370 blockSize = 2;
00371 }
00372 else
00373 {
00374 outgoingDOFs[p].resize(nReq);
00375 }
00376 SUNDANCE_MSG3(setupVerb(),
00377 "p=" << mesh().comm().getRank() << " recv'd from proc=" << p
00378 << " reqs for DOFs for cells " << requestsFromProc);
00379 for (int c=0; c<nReq; c++)
00380 {
00381 int GID = requestsFromProc[c];
00382 SUNDANCE_MSG3(setupVerb(),
00383 "p=" << rank
00384 << " processing cell with d=" << cellDim
00385 << " GID=" << GID);
00386 int LID = mesh().mapGIDToLID(cellDim, GID);
00387 SUNDANCE_MSG3(setupVerb(),
00388 "p=" << rank
00389 << " LID=" << LID << " dofs=" << dofs_[cellDim]);
00390 outgoingDOFs[p][blockSize*c] = dofs_[cellDim][LID][0];
00391 if (cellDim > 0 && cellDim < dim_)
00392 {
00393 outgoingDOFs[p][blockSize*c+1] = originalFacetOrientation_[cellDim-1][LID];
00394 }
00395 SUNDANCE_MSG3(setupVerb(),
00396 "p=" << rank
00397 << " done processing cell with GID=" << GID);
00398 }
00399 }
00400
00401
00402
00403 SUNDANCE_MSG2(setupVerb(),
00404 "p=" << mesh().comm().getRank()
00405 << "answering DOF requests for cells of dimension " << cellDim);
00406
00407
00408 MPIContainerComm<int>::allToAll(outgoingDOFs,
00409 incomingDOFs,
00410 mesh().comm());
00411
00412 SUNDANCE_MSG2(setupVerb(),
00413 "p=" << mesh().comm().getRank()
00414 << "communicated DOF answers for cells of dimension " << cellDim);
00415
00416
00417
00418 for (int p=0; p<mesh().comm().getNProc(); p++)
00419 {
00420 if (p==mesh().comm().getRank()) continue;
00421 const Array<int>& dofsFromProc = incomingDOFs[p];
00422 int blockSize = 1;
00423 if (cellDim > 0 && cellDim < dim_)
00424 {
00425 blockSize = 2;
00426 }
00427
00428 for (int c=0; c<dofsFromProc.size()/blockSize; c++)
00429 {
00430 int cellGID = outgoingCellRequests[p][c];
00431 int cellLID = mesh().mapGIDToLID(cellDim, cellGID);
00432 int dof = dofsFromProc[blockSize*c];
00433 setDOFs(cellDim, cellLID, dof, true);
00434 if (cellDim > 0 && cellDim < dim_)
00435 {
00436 originalFacetOrientation_[cellDim-1][cellLID] = dofsFromProc[blockSize*c+1];
00437 }
00438 }
00439 }
00440
00441 }
00442
00443
00444
00445 void HomogeneousDOFMap::setDOFs(int cellDim, int cellLID,
00446 int& nextDOF,
00447 bool isRemote)
00448 {
00449 Tabs tab;
00450 SUNDANCE_MSG3(setupVerb(), tab << "setting DOFs for " << cellDim << "-cell " << cellLID);
00451 Array<int>& cellDOFs = dofs_[cellDim][cellLID];
00452
00453 int nn = nNodesPerCell_[cellDim];
00454 int nf = funcIDList().size();
00455
00456 for (int i=0; i<nn; i++)
00457 {
00458 for (int f=0; f<nf; f++)
00459 {
00460 int k = nextDOF++;
00461 cellDOFs[funcIDList()[f] + nf*i] = k;
00462 if (isRemote) addGhostIndex(k);
00463 }
00464 }
00465 }
00466
00467
00468
00469 void HomogeneousDOFMap::getDOFsForCellBatch(int cellDim,
00470 const Array<int>& cellLID,
00471 Array<int>& dofs,
00472 int& nNodes) const
00473 {
00474 TimeMonitor timer(dofBatchLookupTimer());
00475
00476 Tabs tab;
00477 SUNDANCE_MSG3(setupVerb(),
00478 tab << "getDOFsForCellBatch(): cellDim=" << cellDim
00479 << " cellLID=" << cellLID);
00480
00481 if (cellLID.size()==0) return;
00482
00483 if (cellDim == dim_)
00484 {
00485 if (!haveMaximalDofs_)
00486 {
00487 buildMaximalDofTable();
00488 }
00489
00490 SUNDANCE_MSG4(setupVerb(), tab << "getting dofs for maximal cells");
00491
00492 nNodes = totalNNodesPerCell_[cellDim];
00493 int nCells = cellLID.size();
00494 int nTotalCells = mesh().numCells(cellDim);
00495 int nf = funcIDList().size();
00496 dofs.resize(totalNNodesPerCell_[cellDim] * cellLID.size() * nf);
00497
00498 for (int c=0; c<cellLID.size(); c++)
00499 {
00500 Tabs tab1;
00501 SUNDANCE_MSG4(setupVerb(), tab1 << "cell=" << c);
00502 for (int fid=0; fid<nf; fid++)
00503 {
00504 Tabs tab2;
00505 SUNDANCE_MSG4(setupVerb(), tab2 << "f= " << fid);
00506 for (int n=0; n<nNodes; n++)
00507 {
00508 Tabs tab3;
00509 dofs[(fid*nCells + c)*nNodes + n]
00510 = maximalDofs_[(fid*nTotalCells + cellLID[c])*nNodes + n];
00511 SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << n << " dof="
00512 << dofs[(fid*nCells + c)*nNodes + n]);
00513 }
00514
00515 }
00516 }
00517 }
00518 else
00519 {
00520 int nf = funcIDList().size();
00521 int nCells = cellLID.size();
00522 dofs.resize(totalNNodesPerCell_[cellDim] * cellLID.size() * nf);
00523
00524 Tabs tab1;
00525 SUNDANCE_MSG4(setupVerb(), tab1 << "getting dofs for non-maximal cells");
00526
00527 static Array<Array<int> > facetLID(3);
00528 static Array<Array<int> > facetOrientations(3);
00529 static Array<int> numFacets(3);
00530
00531 for (int d=0; d<cellDim; d++)
00532 {
00533 numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00534 mesh().getFacetLIDs(cellDim, cellLID, d, facetLID[d],
00535 facetOrientations[d]);
00536 }
00537
00538
00539 int nInteriorNodes = localNodePtrs_[cellDim][cellDim][0].size();
00540
00541 nNodes = totalNNodesPerCell_[cellDim];
00542
00543 for (int c=0; c<cellLID.size(); c++)
00544 {
00545 Tabs tab2;
00546 SUNDANCE_MSG4(setupVerb(), tab2 << "cell=" << c);
00547 SUNDANCE_MSG4(setupVerb(), tab2 << "doing interior nodes");
00548
00549
00550
00551 if (cellDim==0 || nInteriorNodes <= 1)
00552 {
00553 for (int m=0; m<nInteriorNodes; m++)
00554 {
00555 Tabs tab3;
00556 int ptr = localNodePtrs_[cellDim][cellDim][0][m];
00557 for (int f=0; f<nf; f++)
00558 {
00559 dofs[(f*nCells + c)*nNodes+ptr]
00560 = dofs_[cellDim][cellLID[c]][f + nf*m];
00561 SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << m << " f=" << f
00562 << " dof="
00563 << dofs[(f*nCells + c)*nNodes+ptr]);
00564 }
00565 }
00566 }
00567 else
00568 {
00569 int sign = originalFacetOrientation_[cellDim-1][cellLID[c]];
00570 for (int m=0; m<nInteriorNodes; m++)
00571 {
00572 int n = m;
00573 if (sign < 0) n = nInteriorNodes - 1 - m;
00574 Tabs tab3;
00575 int ptr = localNodePtrs_[cellDim][cellDim][0][m];
00576 for (int f=0; f<nf; f++)
00577 {
00578 dofs[(f*nCells + c)*nNodes+ptr]
00579 = dofs_[cellDim][cellLID[c]][f + nf*n];
00580 SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << m << " f=" << f
00581 << " dof="
00582 << dofs[(f*nCells + c)*nNodes+ptr]);
00583 }
00584 }
00585 }
00586
00587 SUNDANCE_MSG4(setupVerb(), tab1 << "doing facet nodes");
00588
00589 for (int d=0; d<cellDim; d++)
00590 {
00591 Tabs tab3;
00592 SUNDANCE_MSG4(setupVerb(), tab3 << "d= " << d);
00593 for (int f=0; f<numFacets[d]; f++)
00594 {
00595 Tabs tab4;
00596 int facetID = facetLID[d][c*numFacets[d]+f];
00597
00598 SUNDANCE_MSG4(setupVerb(), tab4 << "f= " << f << " facetLID = " << facetID);
00599 int nFacetNodes = localNodePtrs_[cellDim][d][f].size();
00600 if (d==0 || nFacetNodes <= 1)
00601 {
00602 Tabs tab5;
00603 for (int n=0; n<nFacetNodes; n++)
00604 {
00605 SUNDANCE_MSG3(setupVerb(),
00606 tab5 << "n=" << n);
00607 int ptr = localNodePtrs_[cellDim][d][f][n];
00608 SUNDANCE_MSG3(setupVerb(),
00609 tab5 << "local ptr=" << ptr);
00610 for (int funcID=0; funcID<nf; funcID++)
00611 {
00612 SUNDANCE_MSG3(setupVerb(),
00613 tab5 << "found dof="
00614 << dofs_[d][facetID][funcID + nf*n]);
00615 dofs[(funcID*nCells + c)*nNodes+ptr]
00616 = dofs_[d][facetID][funcID + nf*n];
00617 }
00618 }
00619 }
00620 else
00621 {
00622 Tabs tab5;
00623 int facetOrientation
00624 = facetOrientations[d][c*numFacets[d]+f]
00625 * originalFacetOrientation_[d-1][facetID];
00626 for (int m=0; m<nFacetNodes; m++)
00627 {
00628 int n = m;
00629 if (facetOrientation<0) n = nFacetNodes-1-m;
00630 SUNDANCE_MSG3(setupVerb(),
00631 tab5 << "n=" << n);
00632 int ptr = localNodePtrs_[cellDim][d][f][m];
00633 SUNDANCE_MSG3(setupVerb(),
00634 tab5 << "local ptr=" << ptr);
00635 for (int funcID=0; funcID<nf; funcID++)
00636 {
00637 SUNDANCE_MSG3(setupVerb(),
00638 tab5 << "found dof="
00639 << dofs_[d][facetID][funcID + nf*n]);
00640 dofs[(funcID*nCells + c)*nNodes+ptr]
00641 = dofs_[d][facetID][funcID + nf*n];
00642 }
00643 }
00644 }
00645 }
00646 }
00647 }
00648 }
00649 }
00650
00651 void HomogeneousDOFMap::buildMaximalDofTable() const
00652 {
00653 Tabs tab;
00654 int cellDim = dim_;
00655 int nf = funcIDList().size();
00656 int nCells = mesh().numCells(dim_);
00657
00658 SUNDANCE_MSG2(setupVerb(), tab << "building dofs for maximal cells");
00659
00660 SUNDANCE_MSG3(setupVerb(), tab << "nf=" << nf
00661 << " total nNodes=" << totalNNodesPerCell_[cellDim]);
00662
00663 static Array<Array<int> > facetLID(3);
00664 static Array<Array<int> > facetOrientations(3);
00665 static Array<int> numFacets(3);
00666
00667 Array<int> cellLID(nCells);
00668
00669 for (int c=0; c<cellLID.size(); c++) cellLID[c]=c;
00670
00671 for (int d=0; d<cellDim; d++)
00672 {
00673 numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00674 mesh().getFacetLIDs(cellDim, cellLID, d, facetLID[d], facetOrientations[d]);
00675 }
00676
00677
00678 int nInteriorNodes = localNodePtrs_[cellDim][cellDim][0].size();
00679
00680
00681 int nNodes = totalNNodesPerCell_[cellDim];
00682
00683 maximalDofs_.resize(nCells*nf*nNodes);
00684
00685 for (int c=0; c<nCells; c++)
00686 {
00687 Tabs tab1;
00688 SUNDANCE_MSG4(setupVerb(), tab1 << "working on cell=" << c << " LID=" << cellLID[c]);
00689
00690
00691 SUNDANCE_MSG4(setupVerb(), tab1 << "doing interior nodes");
00692 for (int n=0; n<nInteriorNodes; n++)
00693 {
00694 int ptr = localNodePtrs_[cellDim][cellDim][0][n];
00695 for (int f=0; f<nf; f++)
00696 {
00697 maximalDofs_[(f*nCells + c)*nNodes+ptr]
00698 = dofs_[cellDim][c][f + nf*n];
00699 }
00700 }
00701
00702 SUNDANCE_MSG4(setupVerb(), tab1 << "doing facet nodes");
00703
00704 for (int d=0; d<cellDim; d++)
00705 {
00706 Tabs tab2;
00707 SUNDANCE_MSG4(setupVerb(), tab2 << "facet dim=" << d);
00708 for (int f=0; f<numFacets[d]; f++)
00709 {
00710 Tabs tab3;
00711 int facetID = facetLID[d][c*numFacets[d]+f];
00712 SUNDANCE_MSG4(setupVerb(), tab2 << "f=" << f << " facetLID=" << facetID);
00713 int nFacetNodes = localNodePtrs_[cellDim][d][f].size();
00714 if (d == 0 || nFacetNodes <= 1)
00715 {
00716 for (int n=0; n<nFacetNodes; n++)
00717 {
00718 Tabs tab4;
00719 SUNDANCE_MSG3(setupVerb(),
00720 tab4 << "n=" << n);
00721 int ptr = localNodePtrs_[cellDim][d][f][n];
00722 for (int funcID=0; funcID<nf; funcID++)
00723 {
00724 SUNDANCE_MSG3(setupVerb(),
00725 tab4 << "found dof="
00726 << dofs_[d][facetID][funcID + nf*n])
00727 maximalDofs_[(funcID*nCells + c)*nNodes+ptr]
00728 = dofs_[d][facetID][funcID + nf*n];
00729 }
00730 }
00731 }
00732 else
00733 {
00734 Tabs tab4;
00735 int facetOrientation = facetOrientations[d][c*numFacets[d]+f]
00736 * originalFacetOrientation_[d-1][facetID];
00737 for (int m=0; m<nFacetNodes; m++)
00738 {
00739 Tabs tab4;
00740 int n = m;
00741 if (facetOrientation<0) n = nFacetNodes-1-m;
00742 SUNDANCE_MSG3(setupVerb(),
00743 tab4 << "n=" << n);
00744 int ptr = localNodePtrs_[cellDim][d][f][m];
00745 for (int funcID=0; funcID<nf; funcID++)
00746 {
00747 SUNDANCE_MSG3(setupVerb(),
00748 tab4 << "found dof="
00749 << dofs_[d][facetID][funcID + nf*n])
00750 maximalDofs_[(funcID*nCells + c)*nNodes+ptr]
00751 = dofs_[d][facetID][funcID + nf*n];
00752 }
00753 }
00754 }
00755 }
00756 }
00757 }
00758 haveMaximalDofs_ = true;
00759 }
00760
00761
00762
00763
00764
00765 void HomogeneousDOFMap::computeOffsets(int dim, int localCount)
00766 {
00767 if (verbosity() > 2)
00768 {
00769 comm().synchronize();
00770 comm().synchronize();
00771 comm().synchronize();
00772 comm().synchronize();
00773 }
00774 SUNDANCE_MSG2(setupVerb(),
00775 "p=" << mesh().comm().getRank()
00776 << " sharing offsets for DOF numbering for dim=" << dim);
00777
00778 SUNDANCE_MSG2(setupVerb(),
00779 "p=" << mesh().comm().getRank()
00780 << " I have " << localCount << " cells");
00781
00782 Array<int> dofOffsets;
00783 int totalDOFCount;
00784 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00785 mesh().comm());
00786 int myOffset = dofOffsets[mesh().comm().getRank()];
00787
00788 SUNDANCE_MSG2(setupVerb(),
00789 "p=" << mesh().comm().getRank()
00790 << " back from MPI accumulate");
00791
00792 if (verbosity() > 2)
00793 {
00794 comm().synchronize();
00795 comm().synchronize();
00796 comm().synchronize();
00797 comm().synchronize();
00798 }
00799
00800 for (int c=0; c<dofs_[dim].size(); c++)
00801 {
00802 if (hasBeenAssigned(dim, c))
00803 {
00804 for (int n=0; n<dofs_[dim][c].size(); n++)
00805 {
00806 dofs_[dim][c][n] += myOffset;
00807 }
00808 }
00809 }
00810
00811 setLowestLocalDOF(myOffset);
00812 setNumLocalDOFs(localCount);
00813 setTotalNumDOFs(totalDOFCount);
00814
00815 SUNDANCE_MSG2(setupVerb(),
00816 "p=" << mesh().comm().getRank()
00817 << " done sharing offsets for DOF numbering for dim=" << dim);
00818 if (verbosity() > 2)
00819 {
00820 comm().synchronize();
00821 comm().synchronize();
00822 comm().synchronize();
00823 comm().synchronize();
00824 }
00825
00826 }
00827
00828
00829
00830 void HomogeneousDOFMap::print(std::ostream& os) const
00831 {
00832 int myRank = mesh().comm().getRank();
00833
00834 Tabs tabs;
00835
00836 os << "DOFS = " << dofs_ << std::endl;
00837
00838 for (int p=0; p<mesh().comm().getNProc(); p++)
00839 {
00840 mesh().comm().synchronize();
00841 mesh().comm().synchronize();
00842 if (p == myRank)
00843 {
00844 os << tabs <<
00845 "========= DOFMap on proc p=" << p << " =============" << std::endl;
00846 for (int d=dim_; d>=0; d--)
00847 {
00848 Tabs tabs1;
00849 os << tabs1 << "dimension = " << d << std::endl;
00850 for (int c=0; c<mesh().numCells(d); c++)
00851 {
00852 Tabs tabs2;
00853 os << tabs2 << "Cell LID=" << c << " GID="
00854 << mesh().mapLIDToGID(d, c) << std::endl;
00855 for (int f=0; f<funcIDList().size(); f++)
00856 {
00857 Tabs tabs3;
00858 Array<int> dofs;
00859 getDOFsForCell(d, c, funcIDList()[f], dofs);
00860 os << tabs3 << "f=" << funcIDList()[f] << " "
00861 << dofs << std::endl;
00862 if (false)
00863 {
00864 os << tabs3 << "{";
00865 for (int i=0; i<dofs.size(); i++)
00866 {
00867 if (i != 0) os << ", ";
00868 if (isLocalDOF(dofs[i])) os << "L";
00869 else os << "R";
00870 }
00871 os << "}" << std::endl;
00872 }
00873 }
00874 }
00875 }
00876 }
00877 mesh().comm().synchronize();
00878 mesh().comm().synchronize();
00879 }
00880 }
00881
00882