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