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 "SundanceMixedDOFMapHN.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 using namespace Sundance;
00042 using namespace Teuchos;
00043 using Playa::MPIComm;
00044 using Playa::MPIContainerComm;
00045 using Playa::MPIDataType;
00046 using Playa::MPIOp;
00047
00048 static Time& mixedHNDOFCtorTimer()
00049 {
00050 static RCP<Time> rtn
00051 = TimeMonitor::getNewTimer("mixed hanging-node DOF map init");
00052 return *rtn;
00053 }
00054 static Time& maxDOFBuildTimer()
00055 {
00056 static RCP<Time> rtn
00057 = TimeMonitor::getNewTimer("max-cell dof table init");
00058 return *rtn;
00059 }
00060
00061 MixedDOFMapHN::MixedDOFMapHN(const Mesh& mesh,
00062 const Array<RCP<BasisDOFTopologyBase> >& basis,
00063 const CellFilter& maxCells,
00064 int setupVerb)
00065 : HNDoFMapBaseHomogeneous(mesh, basis.size(), setupVerb),
00066 maxCells_(maxCells),
00067 dim_(mesh.spatialDim()),
00068 dofs_(mesh.spatialDim()+1),
00069 maximalDofs_(),
00070 haveMaximalDofs_(false),
00071 localNodePtrs_(),
00072 hasCellHanging_(),
00073 isElementHanging_(mesh.spatialDim()+1),
00074 HN_To_globalFacetsLID_(),
00075 HN_To_globalFacetsDim_(),
00076 HN_To_coeffs_(),
00077 maxCellLIDwithHN_to_TrafoMatrix_(),
00078 matrixStore_(),
00079 nNodesPerCell_(),
00080 totalNNodesPerCell_(),
00081 cellHasAnyDOFs_(dim_+1),
00082 numFacets_(mesh.spatialDim()+1),
00083 originalFacetOrientation_(2),
00084 hasBeenAssigned_(mesh.spatialDim()+1),
00085 structure_(),
00086 nFuncs_()
00087 {
00088 TimeMonitor timer(mixedHNDOFCtorTimer());
00089 Tabs tab;
00090
00091 SUNDANCE_MSG1(setupVerb, tab << "building mixed hanging node DOF map");
00092
00093 Sundance::Map<OrderedHandle<BasisDOFTopologyBase>, int> basisToChunkMap;
00094 Array<RCP<BasisDOFTopologyBase> > chunkBases;
00095 Array<Array<int> > chunkFuncs;
00096
00097 int chunk = 0;
00098 int nBasis = basis.size();
00099 basis_.resize(nBasis);
00100
00101 nPoints_ = mesh.numCells(0);
00102
00103 for (int i=0; i<nBasis; i++)
00104 {
00105 OrderedHandle<BasisDOFTopologyBase> bh = basis[i];
00106 if (!basisToChunkMap.containsKey(bh))
00107 {
00108 chunkBases.append(basis[i]);
00109 basisToChunkMap.put(bh, chunk);
00110 chunkFuncs.append(tuple(i));
00111 basis_[chunk] = rcp_dynamic_cast<BasisFamilyBase>(basis[i]);
00112 chunk++;
00113 }
00114 else
00115 {
00116 int b = basisToChunkMap.get(bh);
00117 chunkFuncs[b].append(i);
00118 }
00119 }
00120 basis_.resize(chunk);
00121
00122
00123 matrixStore_.init(chunk);
00124
00125 Tabs tab1;
00126 SUNDANCE_MSG1(setupVerb, tab1 << "basis to chunk map = " << basisToChunkMap);
00127
00128
00129 structure_ = rcp(new MapStructure(basis.size(), chunkBases, chunkFuncs));
00130
00131 nFuncs_.resize(chunkBases.size());
00132 for (int i=0; i<nFuncs_.size(); i++)
00133 nFuncs_[i] = chunkFuncs[i].size();
00134
00135 allocate(mesh);
00136
00137 initMap();
00138
00139 buildMaximalDofTable();
00140
00141
00142
00143
00144 SUNDANCE_MSG1(setupVerb, tab << "done building mixed HN DOF map");
00145
00146 }
00147
00148
00149 void MixedDOFMapHN::allocate(const Mesh& mesh)
00150 {
00151 Tabs tab;
00152
00153
00154 SUNDANCE_MSG1(setupVerb(),tab << "grouping like basis functions");
00155
00156
00157 localNodePtrs_.resize(nBasisChunks());
00158 nNodesPerCell_.resize(nBasisChunks());
00159 totalNNodesPerCell_.resize(nBasisChunks());
00160 nDofsPerCell_.resize(nBasisChunks());
00161 totalNDofsPerCell_.resize(nBasisChunks());
00162 maximalDofs_.resize(nBasisChunks());
00163
00164 for (int b=0; b<nBasisChunks(); b++)
00165 {
00166 localNodePtrs_[b].resize(mesh.spatialDim()+1);
00167 nNodesPerCell_[b].resize(mesh.spatialDim()+1);
00168 totalNNodesPerCell_[b].resize(mesh.spatialDim()+1);
00169 nDofsPerCell_[b].resize(mesh.spatialDim()+1);
00170 totalNDofsPerCell_[b].resize(mesh.spatialDim()+1);
00171 }
00172
00173
00174 SUNDANCE_MSG1(setupVerb(),tab << "working out DOF map node counts");
00175
00176 numFacets_.resize(dim_+1);
00177 isElementHanging_.resize(dim_+1);
00178 hasCellHanging_.resize(mesh.numCells(dim_),false);
00179
00180 for (int d=0; d<=dim_; d++)
00181 {
00182 Tabs tab1;
00183 SUNDANCE_MSG1(setupVerb(),tab << "allocating maps for cell dim=" << d);
00184
00185
00186 numFacets_[d].resize(d);
00187
00188 for (int fd=0; fd<d; fd++) numFacets_[d][fd]=mesh.numFacets(d, 0, fd);
00189 SUNDANCE_MSG1(setupVerb(),tab1 << "num facets for dimension " << d << " is "
00190 << numFacets_[d]);
00191
00192 cellHasAnyDOFs_[d] = false;
00193 dofs_[d].resize(nBasisChunks());
00194
00195 int numCells = mesh.numCells(d);
00196 hasBeenAssigned_[d].resize(numCells);
00197
00198
00199 isElementHanging_[d].resize(numCells,false);
00200 for (int c=0; c<numCells; c++) { isElementHanging_[d][c] = false; }
00201 if (d == dim_){
00202 for (int c=0; c<numCells; c++) { hasCellHanging_[c] = false; }
00203 }
00204
00205 for (int b=0; b<nBasisChunks(); b++)
00206 {
00207 Tabs tab2;
00208 SUNDANCE_MSG2(setupVerb(),tab1 << "basis chunk=" << b);
00209 SUNDANCE_MSG2(setupVerb(),tab2 << "basis=" << basis(b)->description());
00210 int nNodes = basis(b).ptr()->nReferenceDOFsWithFacets(mesh.cellType(dim_), mesh.cellType(d));
00211 SUNDANCE_MSG2(setupVerb(),tab2 << "nNodes per cell=" << nNodes);
00212 if (nNodes == 0)
00213 {
00214 nNodesPerCell_[b][d] = 0;
00215 nDofsPerCell_[b][d] = 0;
00216 }
00217 else
00218 {
00219
00220
00221 basis(b).ptr()->getReferenceDOFs(mesh.cellType(dim_),
00222 mesh.cellType(d),
00223 localNodePtrs_[b][d]);
00224
00225
00226 SUNDANCE_MSG3(setupVerb(),tab2 << "node ptrs are "
00227 << localNodePtrs_[b][d]);
00228
00229
00230
00231 if (localNodePtrs_[b][d][d].size() > 0)
00232 {
00233 nNodesPerCell_[b][d] = localNodePtrs_[b][d][d][0].size();
00234 if (nNodesPerCell_[b][d] > 0) cellHasAnyDOFs_[d] = true;
00235 }
00236 else
00237 {
00238 nNodesPerCell_[b][d] = 0;
00239 }
00240 nDofsPerCell_[b][d] = nNodesPerCell_[b][d] * nFuncs(b);
00241 }
00242
00243
00244
00245
00246
00247 if (nDofsPerCell_[b][d] > 0 && d > 0 && d < mesh.spatialDim())
00248 {
00249 Mesh& tmpMesh = const_cast<Mesh&>(mesh);
00250 tmpMesh.assignIntermediateCellGIDs(d);
00251 }
00252
00253 SUNDANCE_MSG3(setupVerb(),tab2 <<
00254 "num nodes is "
00255 << nNodesPerCell_[b][d]);
00256
00257 totalNNodesPerCell_[b][d] = nNodesPerCell_[b][d];
00258 for (int dd=0; dd<d; dd++)
00259 {
00260 totalNNodesPerCell_[b][d]
00261 += numFacets_[d][dd]*nNodesPerCell_[b][dd];
00262 }
00263 totalNDofsPerCell_[b][d] = totalNNodesPerCell_[b][d] * nFuncs(b);
00264 SUNDANCE_MSG3(setupVerb(),tab2 << "totalNDofsPerCell " << totalNDofsPerCell_[b][d]);
00265
00266 if (nNodes > 0)
00267 {
00268
00269 dofs_[d][b].resize(nDofsPerCell_[b][d] * numCells);
00270
00271
00272 int nDof = dofs_[d][b].size();
00273 Array<int>& dofs = dofs_[d][b];
00274 int marker = uninitializedVal();
00275 for (int i=0; i<nDof; i++)
00276 {
00277 dofs[i] = marker;
00278 }
00279 }
00280
00281 if (d==dim_)
00282 {
00283 maximalDofs_[b].resize(totalNDofsPerCell_[b][d]*numCells);
00284 }
00285 }
00286
00287
00288 if (d > 0 && d < dim_)
00289 {
00290 originalFacetOrientation_[d-1].resize(numCells);
00291 }
00292
00293 }
00294 SUNDANCE_MSG1(setupVerb(),tab << "done allocating DOF map");
00295 }
00296
00297 void MixedDOFMapHN::initMap()
00298 {
00299 Tabs tab;
00300 SUNDANCE_MSG1(setupVerb(),tab << "initializing DOF map");
00301
00302 int nextDOF = 0;
00303
00304
00305
00306
00307 Array<Array<Array<int> > > remoteCells(mesh().spatialDim()+1);
00308
00309 for (int d=0; d<remoteCells.size(); d++)
00310 remoteCells[d].resize(mesh().comm().getNProc());
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 CellSet cells = maxCells_.getCells(mesh());
00321 CellIterator iter;
00322 for (iter=cells.begin(); iter != cells.end(); iter++)
00323 {
00324
00325 int cellLID = *iter;
00326 int owner;
00327
00328 if (cellHasAnyDOFs_[dim_])
00329 {
00330
00331
00332
00333 if (isRemote(dim_, cellLID, owner))
00334 {
00335
00336 int cellGID = mesh().mapLIDToGID(dim_, cellLID);
00337 SUNDANCE_MSG4(setupVerb(),"proc=" << comm().getRank()
00338 << " thinks d-" << dim_
00339 << " cell GID=" << cellGID
00340 << " is owned remotely by p="
00341 << owner);
00342 remoteCells[dim_][owner].append(cellGID);
00343 }
00344 else
00345
00346 {
00347 for (int b=0; b<nBasisChunks(); b++)
00348 {
00349 setDOFs(b, dim_, cellLID, nextDOF);
00350 }
00351 }
00352 }
00353
00354
00355 for (int d=0; d<dim_; d++)
00356 {
00357 if (cellHasAnyDOFs_[d])
00358 {
00359 int nf = numFacets_[dim_][d];
00360 Array<int> facetLID(nf);
00361 Array<int> facetOrientations(nf);
00362
00363
00364 mesh().getFacetArray(dim_, cellLID, d,
00365 facetLID, facetOrientations);
00366
00367 for (int f=0; f<nf; f++)
00368 {
00369
00370
00371 if (!hasBeenAssigned(d, facetLID[f]))
00372 {
00373 markAsAssigned(d, facetLID[f]);
00374
00375 if ( isRemote(d, facetLID[f], owner) && (!mesh().isElementHangingNode(d,facetLID[f])))
00376 {
00377 int facetGID
00378 = mesh().mapLIDToGID(d, facetLID[f]);
00379 SUNDANCE_MSG4(setupVerb(),"proc=" << comm().getRank()
00380 << " thinks d-" << d
00381 << " cell GID=" << facetGID
00382 << " is owned remotely by p=" << owner);
00383 remoteCells[d][owner].append(facetGID);
00384 }
00385 else
00386 {
00387
00388 for (int b=0; b<nBasisChunks(); b++)
00389 {
00390 setDOFs(b, d, facetLID[f], nextDOF);
00391 }
00392
00393 if (d > 0)
00394 {
00395 originalFacetOrientation_[d-1][facetLID[f]]
00396 = facetOrientations[f];
00397 }
00398 }
00399 }
00400 }
00401 }
00402 }
00403 }
00404
00405
00406
00407
00408
00409
00410 int numLocalDOFs = nextDOF;
00411 if (mesh().comm().getNProc() > 1)
00412 {
00413 for (int d=0; d<=dim_; d++)
00414 {
00415 if (cellHasAnyDOFs_[d])
00416 {
00417 computeOffsets(d, numLocalDOFs);
00418 shareDOFs(d, remoteCells[d]);
00419 }
00420 }
00421 }
00422 else
00423 {
00424 setLowestLocalDOF(0);
00425 setNumLocalDOFs(numLocalDOFs);
00426 setTotalNumDOFs(numLocalDOFs);
00427 }
00428 SUNDANCE_MSG1(setupVerb(),tab << "done initializing DOF map , numLocalDOFs:" << numLocalDOFs);
00429 }
00430
00431
00432 void MixedDOFMapHN::setDOFs(int basisChunk, int cellDim, int cellLID,
00433 int& nextDOF, bool isRemote)
00434 {
00435 int nDofs = nDofsPerCell_[basisChunk][cellDim];
00436 if (nDofs==0) return;
00437 Tabs tab;
00438 SUNDANCE_MSG4(setupVerb(),tab << "setting DOFs for " << cellDim
00439 << "-cell " << cellLID << " , basisChunk:" << basisChunk <<" , nDofs:" <<nDofs << " , nextDOF:" << nextDOF );
00440
00441 int* ptr = getInitialDOFPtrForCell(cellDim, cellLID, basisChunk);
00442
00443
00444
00445
00446
00447 if (isRemote)
00448 {
00449 if (mesh().isElementHangingNode(cellDim, cellLID)){
00450 isElementHanging_[cellDim][cellLID] = true;
00451 for (int i=0; i<nDofs; i++) {
00452 ptr[i] = -1;
00453
00454 }
00455 }
00456 else
00457 {
00458 for (int i=0; i<nDofs; i++, nextDOF++){
00459 ptr[i] = nextDOF;
00460 addGhostIndex(nextDOF);
00461 }
00462 }
00463 }
00464 else
00465 {
00466 if (mesh().isElementHangingNode(cellDim, cellLID)){
00467 isElementHanging_[cellDim][cellLID] = true;
00468 for (int i=0; i<nDofs; i++) {
00469
00470
00471
00472
00473 ptr[i] = -1;
00474 }
00475 }
00476 else
00477 {
00478 for (int i=0; i<nDofs; i++,nextDOF++) {
00479
00480
00481
00482 ptr[i] = nextDOF;
00483 }
00484 }
00485 }
00486 }
00487
00488
00489
00490
00491
00492 void MixedDOFMapHN::shareDOFs(int cellDim,
00493 const Array<Array<int> >& outgoingCellRequests)
00494 {
00495 int np = mesh().comm().getNProc();
00496 int rank = mesh().comm().getRank();
00497
00498 Array<Array<int> > incomingCellRequests;
00499 Array<Array<int> > outgoingDOFs(np);
00500 Array<Array<int> > incomingDOFs;
00501
00502 SUNDANCE_MSG2(setupVerb(),
00503 "p=" << mesh().comm().getRank()
00504 << "synchronizing DOFs for cells of dimension " << cellDim);
00505 SUNDANCE_MSG2(setupVerb(),
00506 "p=" << mesh().comm().getRank()
00507 << " sending cell reqs d=" << cellDim << " GID=" << outgoingCellRequests);
00508
00509
00510 MPIContainerComm<int>::allToAll(outgoingCellRequests,
00511 incomingCellRequests,
00512 mesh().comm());
00513
00514
00515
00516
00517
00518 int blockSize = 0;
00519 bool sendOrientation = false;
00520 for (int b=0; b<nBasisChunks(); b++)
00521 {
00522 int nDofs = nDofsPerCell_[b][cellDim];
00523 if (nDofs > 0) blockSize++;
00524 if (nDofs > 1 && cellDim > 0 && cellDim < dim_) sendOrientation = true;
00525 }
00526 blockSize += sendOrientation;
00527
00528 SUNDANCE_MSG2(setupVerb(),
00529 "p=" << rank
00530 << "recvd DOF requests for cells of dimension " << cellDim
00531 << " GID=" << incomingCellRequests);
00532
00533
00534
00535 for (int p=0; p<np; p++)
00536 {
00537 if (p==rank) continue;
00538 const Array<int>& requestsFromProc = incomingCellRequests[p];
00539 int nReq = requestsFromProc.size();
00540
00541 SUNDANCE_MSG4(setupVerb(),"p=" << mesh().comm().getRank()
00542 << " recv'd from proc=" << p
00543 << " reqs for DOFs for cells "
00544 << requestsFromProc);
00545
00546 outgoingDOFs[p].resize(nReq * blockSize);
00547
00548 for (int c=0; c<nReq; c++)
00549 {
00550 int GID = requestsFromProc[c];
00551 SUNDANCE_MSG3(setupVerb(),
00552 "p=" << rank
00553 << " processing cell with d=" << cellDim
00554 << " GID=" << GID);
00555 int LID = mesh().mapGIDToLID(cellDim, GID);
00556 SUNDANCE_MSG3(setupVerb(),
00557 "p=" << rank
00558 << " LID=" << LID << " dofs=" << dofs_[cellDim]);
00559 int blockOffset = 0;
00560 for (int b=0; b<nBasisChunks(); b++)
00561 {
00562 if (nDofsPerCell_[b][cellDim] == 0) continue;
00563 outgoingDOFs[p][blockSize*c+blockOffset]
00564 = getInitialDOFForCell(cellDim, LID, b);
00565 blockOffset++;
00566 }
00567 if (sendOrientation)
00568 {
00569 outgoingDOFs[p][blockSize*(c+1) - 1]
00570 = originalFacetOrientation_[cellDim-1][LID];
00571 }
00572 SUNDANCE_MSG3(setupVerb(),
00573 "p=" << rank
00574 << " done processing cell with GID=" << GID);
00575 }
00576 }
00577
00578
00579 SUNDANCE_MSG2(setupVerb(),
00580 "p=" << mesh().comm().getRank()
00581 << "answering DOF requests for cells of dimension " << cellDim);
00582
00583
00584 MPIContainerComm<int>::allToAll(outgoingDOFs,
00585 incomingDOFs,
00586 mesh().comm());
00587
00588 SUNDANCE_MSG2(setupVerb(),
00589 "p=" << mesh().comm().getRank()
00590 << "communicated DOF answers for cells of dimension " << cellDim);
00591
00592
00593
00594
00595 for (int p=0; p<mesh().comm().getNProc(); p++)
00596 {
00597 if (p==mesh().comm().getRank()) continue;
00598 const Array<int>& dofsFromProc = incomingDOFs[p];
00599 int numCells = dofsFromProc.size()/blockSize;
00600 for (int c=0; c<numCells; c++)
00601 {
00602 int cellGID = outgoingCellRequests[p][c];
00603 int cellLID = mesh().mapGIDToLID(cellDim, cellGID);
00604 int blockOffset = 0;
00605 for (int b=0; b<nBasisChunks(); b++)
00606 {
00607 if (nDofsPerCell_[b][cellDim] == 0) continue;
00608 int dof = dofsFromProc[blockSize*c+blockOffset];
00609 setDOFs(b, cellDim, cellLID, dof, true);
00610 blockOffset++;
00611 }
00612 if (sendOrientation)
00613 {
00614 originalFacetOrientation_[cellDim-1][cellLID]
00615 = dofsFromProc[blockSize*(c+1)-1];
00616 }
00617 }
00618 }
00619
00620 }
00621
00622
00623 RCP<const MapStructure> MixedDOFMapHN
00624 ::getDOFsForCellBatch(int cellDim,
00625 const Array<int>& cellLID,
00626 const Set<int>& requestedFuncSet,
00627 Array<Array<int> >& dofs,
00628 Array<int>& nNodes,
00629 int verbosity) const
00630 {
00631 TimeMonitor timer(batchedDofLookupTimer());
00632
00633 Tabs tab;
00634
00635 SUNDANCE_MSG2(verbosity,
00636 tab << "getDOFsForCellBatch(): cellDim=" << cellDim
00637 << " cellLID=" << cellLID);
00638
00639 dofs.resize(nBasisChunks());
00640 nNodes.resize(nBasisChunks());
00641
00642 int nCells = cellLID.size();
00643
00644 if (cellDim == dim_)
00645 {
00646 Tabs tab1;
00647
00648 SUNDANCE_MSG4(verbosity,tab1 << "getting dofs for maximal cells : " << cellLID );
00649
00650 for (int b=0; b<nBasisChunks(); b++)
00651 {
00652 nNodes[b] = totalNNodesPerCell_[b][cellDim];
00653 dofs[b].resize(nNodes[b]*nFuncs(b)*nCells);
00654 int dofsPerCell = nFuncs(b)*nNodes[b];
00655 Array<int>& chunkDofs = dofs[b];
00656
00657
00658
00659 for (int c=0; c<nCells; c++)
00660 {
00661 for (int i=0; i<dofsPerCell; i++)
00662 {
00663 chunkDofs[c*dofsPerCell + i]
00664 = maximalDofs_[b][cellLID[c]*dofsPerCell+i];
00665 }
00666 }
00667 SUNDANCE_MSG3(verbosity,tab1 << "dofs for chunk-" << b << " :" << chunkDofs);
00668 }
00669 }
00670 else
00671 {
00672 Tabs tab1;
00673 SUNDANCE_MSG3(verbosity,tab1 << "getting dofs for non-maximal cells");
00674
00675 static Array<Array<int> > facetLID(3);
00676 static Array<Array<int> > facetOrientations(3);
00677 static Array<int> numFacets(3);
00678
00679 for (int d=0; d<cellDim; d++)
00680 {
00681 numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00682 mesh().getFacetLIDs(cellDim, cellLID, d, facetLID[d],
00683 facetOrientations[d]);
00684 }
00685
00686 for (int b=0; b<nBasisChunks(); b++)
00687 {
00688 nNodes[b] = totalNNodesPerCell_[b][cellDim];
00689 dofs[b].resize(nNodes[b]*nFuncs(b)*nCells);
00690 int dofsPerCell = nFuncs(b)*nNodes[b];
00691
00692 Array<int>& toPtr = dofs[b];
00693 int nf = nFuncs(b);
00694
00695 for (int c=0; c<nCells; c++)
00696 {
00697 Tabs tab2;
00698 SUNDANCE_MSG3(verbosity,tab2 << "cell=" << c << " ,cellLID[c]:" << cellLID[c]);
00699 int offset = dofsPerCell*c;
00700
00701
00702
00703
00704 SUNDANCE_MSG3(verbosity,tab2 << "doing interior nodes");
00705 int nInteriorNodes = nNodesPerCell_[b][cellDim];
00706
00707 if (nInteriorNodes > 0)
00708 {
00709 if (cellDim==0 || nInteriorNodes <= 1)
00710 {
00711
00712 for (int func=0; func<nf; func++)
00713 {
00714 for (int n=0; n<nInteriorNodes; n++)
00715 {
00716 int ptr = localNodePtrs_[b][cellDim][cellDim][0][n];
00717 toPtr[offset + func*nNodes[b] + ptr]
00718 = dofs_[cellDim][b][cellLID[c]*nDofsPerCell_[b][cellDim]+func*nNodesPerCell_[b][cellDim]+n];
00719
00720 }
00721 }
00722 }
00723 else
00724 {
00725 int sign = originalFacetOrientation_[cellDim-1][cellLID[c]];
00726 int nInteriorNodes = localNodePtrs_[b][cellDim][cellDim][0].size();
00727
00728 for (int func=0; func<nf; func++)
00729 {
00730 for (int m=0; m<nInteriorNodes; m++)
00731 {
00732 int n = m;
00733 if (sign<0) n = nInteriorNodes-1-m;
00734 int ptr = localNodePtrs_[b][cellDim][cellDim][0][m];
00735 toPtr[offset + func*nNodes[b] + ptr] = dofs_[cellDim][b][cellLID[c]*nDofsPerCell_[b][cellDim]+func*nNodesPerCell_[b][cellDim]+n];
00736
00737 }
00738 }
00739 }
00740 }
00741
00742
00743 for (int d=0; d<cellDim; d++)
00744 {
00745 Tabs tab2;
00746 SUNDANCE_MSG3(verbosity,tab2 << "facet dim=" << d);
00747 if (nNodesPerCell_[b][d] == 0) continue;
00748 for (int f=0; f<numFacets[d]; f++)
00749 {
00750 Tabs tab3;
00751 int facetID = facetLID[d][c*numFacets[d]+f];
00752 SUNDANCE_MSG3(verbosity,tab2 << "f=" << f << " facetLID=" << facetID);
00753 if (localNodePtrs_[b][cellDim].size()==0) continue;
00754 int nFacetNodes = localNodePtrs_[b][cellDim][d][f].size();
00755
00756 int* toPtr1 = &(dofs[b][dofsPerCell*c]);
00757 const int* nodePtr = &(localNodePtrs_[b][cellDim][d][f][0]);
00758
00759
00760
00761 for (int func=0; func<nf; func++)
00762 {
00763 if (d == 0 || nFacetNodes <= 1)
00764 {
00765 for (int n=0; n<nFacetNodes; n++)
00766 {
00767 int ptr = nodePtr[n];
00768 toPtr1[func*nNodes[b] + ptr]
00769 = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
00770
00771
00772
00773 if (toPtr1[func*nNodes[b] + ptr] < 0 )
00774 {
00775 int fIndex = 1 , tmp1 = 1;
00776
00777
00778
00779 int maxCell = mesh().maxCofacetLID( d , facetID , 0 , fIndex);
00780
00781
00782
00783
00784
00785 Array<int> facetsLIDs;
00786 mesh().returnParentFacets( maxCell , d , facetsLIDs , tmp1 );
00787
00788
00789
00790
00791
00792 toPtr1[func*nNodes[b] + ptr]
00793 = dofs_[d][b][facetsLIDs[fIndex]*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
00794 }
00795 }
00796 }
00797 else
00798 {
00799 int facetOrientation = facetOrientations[d][c*numFacets[d]+f]
00800 * originalFacetOrientation_[d-1][facetID];
00801 for (int m=0; m<nFacetNodes; m++)
00802 {
00803 int n = m;
00804 if (facetOrientation<0) n = nFacetNodes-1-m;
00805 int ptr = nodePtr[n];
00806 toPtr1[func*nNodes[b]+ptr]
00807 = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
00808
00809
00810 if (toPtr1[func*nNodes[b]+ptr] < 0)
00811 {
00812 int fIndex = 1 , tmp1 = 1;
00813
00814
00815
00816 int maxCell = mesh().maxCofacetLID( d , facetID , 0 , fIndex);
00817 Array<int> facetsLIDs;
00818 mesh().returnParentFacets( maxCell , d , facetsLIDs , tmp1 );
00819
00820
00821
00822 toPtr1[func*nNodes[b] + ptr]
00823 = dofs_[d][b][facetsLIDs[fIndex]*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
00824 }
00825 }
00826 }
00827 }
00828 }
00829 }
00830 }
00831 }
00832
00833 }
00834 return structure_;
00835 }
00836
00837 void MixedDOFMapHN::buildMaximalDofTable()
00838 {
00839 TimeMonitor timer(maxDOFBuildTimer());
00840 Tabs tab;
00841 int cellDim = dim_;
00842 int nCells = mesh().numCells(dim_);
00843
00844 SUNDANCE_MSG1(setupVerb(),tab << "building dofs for maximal cells");
00845
00846 Array<Array<int> > facetLID(3);
00847 Array<Array<int> > facetOrientations(3);
00848 Array<int> numFacets(3);
00849
00850 Array<int> cellLID(nCells);
00851
00852 for (int c=0; c<nCells; c++) cellLID[c]=c;
00853
00854 for (int d=0; d<cellDim; d++)
00855 {
00856 numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00857 mesh().getFacetLIDs(cellDim, cellLID, d,
00858 facetLID[d], facetOrientations[d]);
00859 }
00860
00861 Array<int> nInteriorNodes(nBasisChunks());
00862 Array<int> nNodes(nBasisChunks());
00863 for (int b = 0; b<nBasisChunks(); b++)
00864 {
00865
00866 if (localNodePtrs_[b][cellDim].size() != 0)
00867 {
00868 nInteriorNodes[b] = localNodePtrs_[b][cellDim][cellDim][0].size();
00869 }
00870 else
00871 {
00872 nInteriorNodes[b] = 0;
00873 }
00874 nNodes[b] = totalNNodesPerCell_[b][cellDim];
00875 }
00876
00877
00878
00879
00880 Array< Array<Array<int> > >globalDoFs_Cell(nBasisChunks());
00881
00882 Array<Array<double> > trafoMatrix(nBasisChunks());
00883
00884 Array<int> localDoFs;
00885 Array<int> parentFacetDim;
00886 Array<int> parentFacetIndex;
00887 Array<int> parentFacetNode;
00888 Array<int> retFacets;
00889 Array<double> coefs;
00890
00891
00892 for (int c=0; c<nCells; c++)
00893 {
00894
00895
00896 for (int jj = 0 ; jj < numFacets[0] ; jj++){
00897 if (mesh().isElementHangingNode(0,facetLID[0][ c*numFacets[0] + jj])){
00898
00899
00900 hasCellHanging_[cellLID[c]] = true;
00901 break;
00902 }
00903 }
00904
00905 Tabs tab1;
00906 SUNDANCE_MSG3(setupVerb(),tab1 << "working on cell=" << c
00907 << " LID=" << cellLID[c]);
00908
00909
00910 SUNDANCE_MSG3(setupVerb(),tab1 << "doing interior nodes");
00911 for (int b=0; b<nBasisChunks(); b++)
00912 {
00913
00914 if (hasCellHanging_[cellLID[c]]){
00915
00916 globalDoFs_Cell[b].resize( nFuncs(b));
00917 trafoMatrix[b].resize(nNodes[b]*nNodes[b], 0.0);
00918 for (int jj = 0; jj < nNodes[b]*nNodes[b] ; jj++) trafoMatrix[b][jj] = 0.0;
00919 for (int func=0; func<nFuncs(b); func++)
00920 {
00921 globalDoFs_Cell[b][func].resize(nNodes[b]);
00922 for (int kk = 0 ; kk < nNodes[b] ; kk++ ){
00923 globalDoFs_Cell[b][func][kk] = -1;
00924 }
00925 }
00926 }
00927
00928 SUNDANCE_MSG3(setupVerb(),tab1 << "basis chunk=" << b);
00929 if (nInteriorNodes[b]>0)
00930 {
00931 SUNDANCE_MSG3(setupVerb(),tab1<< "nInteriorNodes = " <<nInteriorNodes[b]);
00932 int* toPtr = &(maximalDofs_[b][nNodes[b]*nFuncs(b)*cellLID[c]]);
00933 int nf = nFuncs(b);
00934 for (int func=0; func<nf; func++)
00935 {
00936 for (int n=0; n<nInteriorNodes[b]; n++)
00937 {
00938
00939
00940
00941
00942
00943
00944 int ptr = localNodePtrs_[b][cellDim][cellDim][0][n];
00945
00946 if (hasCellHanging_[cellLID[c]])
00947 {
00948 SUNDANCE_MSG1(setupVerb(), tab1 << "Cell has HDoF cellLID[c]=" << cellLID[c] );
00949 int dof_tmp = dofs_[cellDim][b][cellLID[c]*nDofsPerCell_[b][cellDim]+func*nNodesPerCell_[b][cellDim]+n];
00950
00951 SUNDANCE_MSG1(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , globalDoFs_Cell[b][func]:"<<globalDoFs_Cell[b][func]<<" dof=" << dof_tmp);
00952 globalDoFs_Cell[b][func][ptr] = dof_tmp;
00953 SUNDANCE_MSG1(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , ptr=" << ptr <<" dof=" << dof_tmp);
00954 if (func == 0){
00955 trafoMatrix[b][nNodes[b]*ptr + ptr] = 1.0;
00956 }
00957 }
00958
00959 toPtr[func*nNodes[b] + ptr] =
00960 dofs_[cellDim][b][cellLID[c]*nDofsPerCell_[b][cellDim]+func*nNodesPerCell_[b][cellDim]+n];
00961 SUNDANCE_MSG1(setupVerb(), tab1 << " dof=" << dofs_[cellDim][b][cellLID[c]*nDofsPerCell_[b][cellDim]+func*nNodesPerCell_[b][cellDim]+n]);
00962 }
00963 }
00964 }
00965 }
00966
00967 SUNDANCE_MSG3(setupVerb(),tab1 << "doing facet nodes");
00968
00969 for (int d=0; d<cellDim; d++)
00970 {
00971 Tabs tab2;
00972 SUNDANCE_MSG3(setupVerb(),tab2 << "facet dim=" << d);
00973
00974 for (int f=0; f<numFacets[d]; f++)
00975 {
00976 Tabs tab3;
00977 int facetID = facetLID[d][c*numFacets[d]+f];
00978 SUNDANCE_MSG3(setupVerb(),tab2 << "f=" << f << " facetLID=" << facetID);
00979
00980 for (int b=0; b<nBasisChunks(); b++)
00981 {
00982 Tabs tab4;
00983 SUNDANCE_MSG3(setupVerb(),tab4 << "basis chunk=" << b);
00984 SUNDANCE_MSG3(setupVerb(),tab4 << "num nodes per cell=" << nNodesPerCell_[b][d]);
00985 int nf = nFuncs(b);
00986 if (nDofsPerCell_[b][d]==0) continue;
00987 int nFacetNodes = 0;
00988 if (localNodePtrs_[b][cellDim].size()!=0)
00989 nFacetNodes = localNodePtrs_[b][cellDim][d][f].size();
00990 if (nFacetNodes == 0) continue;
00991
00992 SUNDANCE_MSG3(setupVerb(),tab4 << "dof table size=" << maximalDofs_[b].size());
00993
00994 int* toPtr = &(maximalDofs_[b][nNodes[b]*nFuncs(b)*cellLID[c]]);
00995
00996 const int* nodePtr = &(localNodePtrs_[b][cellDim][d][f][0]);
00997
00998 for (int func=0; func<nf; func++)
00999 {
01000 Tabs tab5;
01001 SUNDANCE_MSG3(setupVerb(),tab5 << "func=" << func);
01002 if (d == 0 || nFacetNodes <= 1)
01003 {
01004 Tabs tab6;
01005 for (int n=0; n<nFacetNodes; n++)
01006 {
01007 SUNDANCE_MSG3(setupVerb(),tab5 << "n=" << n);
01008 int ptr = nodePtr[n];
01009 SUNDANCE_MSG3(setupVerb(),tab6 << "ptr=" << ptr);
01010
01011
01012
01013
01014 if (hasCellHanging_[cellLID[c]]){
01015 int parentCellLID = 0;
01016 int dof_tmp = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
01017 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , d="<<d<<", f="<<f<<", n="<<n<<", dof=" << dof_tmp);
01018 if (dof_tmp < 0){
01019 Array<int> constFacetLID;
01020 Array<int> constFacetDim;
01021
01022 basis_[b]->getConstrainsForHNDoF(
01023 mesh().indexInParent(cellLID[c]), cellDim ,
01024 mesh().maxChildren(), d , f , n,
01025 localDoFs, parentFacetDim,
01026 parentFacetIndex, parentFacetNode, coefs );
01027 SUNDANCE_MSG3(setupVerb(), tab1 << " After basis_[b]->getConstrainsForHNDoF:");
01028 constFacetLID.resize(localDoFs.size());
01029 constFacetDim.resize(localDoFs.size());
01030 for (int indexi = 0 ; indexi < localDoFs.size() ; indexi++)
01031 {
01032
01033 int ptr1 = localNodePtrs_[b][cellDim][parentFacetDim[indexi]]
01034 [parentFacetIndex[indexi]][parentFacetNode[indexi]];
01035
01036 mesh().returnParentFacets(cellLID[c] , parentFacetDim[indexi], retFacets , parentCellLID );
01037 int parFacetLID = retFacets[parentFacetIndex[indexi]];
01038 int parFacetNode = parentFacetNode[indexi];
01039 SUNDANCE_MSG3(setupVerb(), tab1 << "returnParentFacets , parFacetLID=" << parFacetLID <<" parFacetNode=" << parFacetNode );
01040 SUNDANCE_MSG3(setupVerb(), tab1 << "returnParentFacets , retFacets=" << retFacets );
01041
01042 constFacetDim[indexi] = parentFacetDim[indexi];
01043 constFacetLID[indexi] = parFacetLID;
01044
01045 dof_tmp = dofs_[parentFacetDim[indexi]][b]
01046 [parFacetLID * nDofsPerCell_[b][parentFacetDim[indexi]]+func*nNodesPerCell_[b][parentFacetDim[indexi]]+parFacetNode];
01047
01048 globalDoFs_Cell[b][func][ptr1] = dof_tmp;
01049
01050 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , ptr=" << ptr << "ptr1=" << ptr1 <<" dof=" << dof_tmp);
01051 if (func == 0){
01052 trafoMatrix[b][nNodes[b]*ptr + ptr1] = coefs[indexi];
01053 }
01054 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , DONE Fill");
01055 }
01056
01057 if ( (d == 0) && (func == 0)) {
01058 HN_To_globalFacetsLID_.put( nPoints_*b + facetID , constFacetLID);
01059 HN_To_globalFacetsDim_.put( nPoints_*b + facetID , constFacetDim);
01060 HN_To_coeffs_.put( nPoints_*b + facetID , coefs);
01061 }
01062 }else {
01063
01064 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , ptr=" << ptr <<" dof=" << dof_tmp);
01065
01066
01067
01068
01069 globalDoFs_Cell[b][func][ptr] = dof_tmp;
01070 if (func == 0){
01071 trafoMatrix[b][nNodes[b]*ptr + ptr] = 1.0;
01072 }
01073 }
01074 } else {
01075 toPtr[func*nNodes[b] + ptr]
01076 = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
01077 SUNDANCE_MSG3(setupVerb(), tab1 << " dof=" << dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n]);
01078 }
01079 }
01080 }
01081 else
01082 {
01083 int facetOrientation = facetOrientations[d][c*numFacets[d]+f]
01084 * originalFacetOrientation_[d-1][facetID];
01085 for (int m=0; m<nFacetNodes; m++)
01086 {
01087 int n = m;
01088 if (facetOrientation<0) n = nFacetNodes-1-m;
01089 int ptr = nodePtr[m];
01090
01091
01092
01093
01094
01095 if (hasCellHanging_[cellLID[c]])
01096 {
01097 int parentCellLID = 0;
01098 int dof_tmp = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
01099 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , d="<<d<<", f="<<f<<", n="<<n<<", dof=" << dof_tmp);
01100 if (dof_tmp < 0){
01101 Array<int> constFacetLID;
01102 Array<int> constFacetDim;
01103
01104 basis_[b]->getConstrainsForHNDoF(
01105 mesh().indexInParent(cellLID[c]), cellDim ,
01106 mesh().maxChildren(), d , f , n,
01107 localDoFs, parentFacetDim,
01108 parentFacetIndex, parentFacetNode, coefs );
01109 SUNDANCE_MSG3(setupVerb(), tab1 << " After basis_[b]->getConstrainsForHNDoF:");
01110 constFacetLID.resize(localDoFs.size());
01111 constFacetDim.resize(localDoFs.size());
01112 for (int indexi = 0 ; indexi < localDoFs.size() ; indexi++)
01113 {
01114
01115 int ptr1 = localNodePtrs_[b][cellDim][parentFacetDim[indexi]]
01116 [parentFacetIndex[indexi]][parentFacetNode[indexi]];
01117
01118 mesh().returnParentFacets(cellLID[c] , parentFacetDim[indexi], retFacets , parentCellLID );
01119 int parFacetLID = retFacets[parentFacetIndex[indexi]];
01120 int parFacetNode = parentFacetNode[indexi];
01121 SUNDANCE_MSG3(setupVerb(), tab1 << "returnParentFacets , parFacetLID=" << parFacetLID <<" parFacetNode=" << parFacetNode );
01122 SUNDANCE_MSG3(setupVerb(), tab1 << "returnParentFacets , retFacets=" << retFacets );
01123
01124 constFacetDim[indexi] = parentFacetDim[indexi];
01125 constFacetLID[indexi] = parFacetLID;
01126
01127 dof_tmp = dofs_[parentFacetDim[indexi]][b]
01128 [parFacetLID * nDofsPerCell_[b][parentFacetDim[indexi]]+func*nNodesPerCell_[b][parentFacetDim[indexi]]+parFacetNode];
01129
01130 globalDoFs_Cell[b][func][ptr1] = dof_tmp;
01131
01132 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , ptr=" << ptr1 <<" dof=" << dof_tmp);
01133 if (func == 0){
01134 trafoMatrix[b][nNodes[b]*ptr + ptr1] = coefs[indexi];
01135 }
01136 }
01137
01138 if ( (d == 0) && (func == 0)) {
01139 HN_To_globalFacetsLID_.put( nPoints_*b + facetID , constFacetLID);
01140 HN_To_globalFacetsDim_.put( nPoints_*b + facetID , constFacetDim);
01141 HN_To_coeffs_.put( nPoints_*b + facetID , coefs);
01142 }
01143 }else{
01144
01145
01146 SUNDANCE_MSG3(setupVerb(), tab1 << "Cell has HDoF cellLID[c] , ptr=" << ptr <<" dof=" << dof_tmp);
01147 globalDoFs_Cell[b][func][ptr] = dof_tmp;
01148 if (func == 0){
01149 trafoMatrix[b][nNodes[b]*ptr + ptr] = 1.0;
01150 }
01151 }
01152 } else {
01153 toPtr[func*nNodes[b]+ptr]
01154 = dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n];
01155 SUNDANCE_MSG3(setupVerb(), tab1 << " dof=" << dofs_[d][b][facetID*nDofsPerCell_[b][d]+func*nNodesPerCell_[b][d]+n]);
01156 }
01157 }
01158 }
01159 }
01160 }
01161 }
01162 }
01163
01164
01165 if (hasCellHanging_[cellLID[c]]){
01166
01167 Array<int> matrixIndexes(nBasisChunks());
01168 for (int b=0; b<nBasisChunks(); b++){
01169 matrixIndexes[b] = matrixStore_.addMatrix( b , trafoMatrix[b] );
01170 }
01171 maxCellLIDwithHN_to_TrafoMatrix_.put( cellLID[c] , matrixIndexes );
01172
01173 for (int b=0; b<nBasisChunks(); b++)
01174 {
01175
01176 SUNDANCE_MSG3(setupVerb(), "trafoMatrix[b]=" << trafoMatrix[b]);
01177 for (int func=0; func<nFuncs(b); func++)
01178 {
01179
01180 SUNDANCE_MSG1(setupVerb(), "globalDoFs_Cell[b][func]=" << globalDoFs_Cell[b][func]);
01181 for (int jj=0 ; jj < nNodes[b] ; jj++){
01182
01183 maximalDofs_[b][(cellLID[c]*nFuncs(b) + func)*nNodes[b] + jj] = globalDoFs_Cell[b][func][jj];
01184 }
01185 }
01186 }
01187 }
01188 }
01189
01190 haveMaximalDofs_ = true;
01191
01192 SUNDANCE_MSG1(setupVerb(),tab << "done building dofs for maximal cells");
01193 }
01194
01195
01196 void MixedDOFMapHN::getTrafoMatrixForCell(
01197 int cellLID,
01198 int funcID,
01199 int& trafoMatrixSize,
01200 bool& doTransform,
01201 Array<double>& transfMatrix ) const {
01202
01203 trafoMatrixSize = 0;
01204
01205 if (maxCellLIDwithHN_to_TrafoMatrix_.containsKey(cellLID))
01206 {
01207
01208 const Array<int> &matrixIndexes = maxCellLIDwithHN_to_TrafoMatrix_.get( cellLID );
01209 matrixStore_.getMatrix( chunkForFuncID(funcID) , matrixIndexes[chunkForFuncID(funcID)] , transfMatrix );
01210
01211
01212
01213 trafoMatrixSize = sqrt((double) transfMatrix.size());
01214 SUNDANCE_MSG1(setupVerb(), "getTrafoMatrixForCell() cellLID:" << cellLID << ",funcID:" <<
01215 funcID << ",chunkForFuncID(funcID):" << chunkForFuncID(funcID) << ", trafoMatrixSize:" << trafoMatrixSize);
01216
01217 doTransform = true;
01218 }
01219 else
01220 {
01221 doTransform = false;
01222 }
01223 }
01224
01225 void MixedDOFMapHN::getTrafoMatrixForFacet(
01226 int cellDim,
01227 int cellLID,
01228 int facetIndex,
01229 int funcID,
01230 int& trafoMatrixSize,
01231 bool& doTransform,
01232 Array<double>& transfMatrix ) const {
01233
01234 int fIndex;
01235 int maxCellLID;
01236
01237 SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() cellDim :" << cellDim << ", cellLID:" << cellLID);
01238 maxCellLID = mesh().maxCofacetLID( cellDim, cellLID, 0 , fIndex);
01239 SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() testing :" << maxCellLID);
01240
01241
01242
01243 if (maxCellLIDwithHN_to_TrafoMatrix_.containsKey(maxCellLID))
01244 {
01245 const Array<int> &matrixIndexes = maxCellLIDwithHN_to_TrafoMatrix_.get( maxCellLID );
01246 matrixStore_.getMatrix( chunkForFuncID(funcID) , matrixIndexes[chunkForFuncID(funcID)] , transfMatrix );
01247 doTransform = true;
01248 }
01249 else
01250 {
01251 doTransform = false;
01252 }
01253 }
01254
01255
01256 void MixedDOFMapHN::getDOFsForHNCell(
01257 int cellDim,
01258 int cellLID,
01259 int funcID,
01260 Array<int>& dofs ,
01261 Array<double>& coefs ) const {
01262
01263
01264
01265
01266
01267 if ( (cellDim == 0) && ( HN_To_globalFacetsLID_.containsKey(nPoints_*chunkForFuncID(funcID) + cellLID)) )
01268 {
01269 Array<int> facetLIDs;
01270 Array<int> facetDim;
01271 SUNDANCE_MSG1(setupVerb(), "getDOFsForHNCell() cellDim:"<<cellDim<<" cellLID:"<<
01272 cellLID<<"funcID:" << funcID <<",chunkForFuncID(funcID)" << chunkForFuncID(funcID));
01273 facetLIDs = HN_To_globalFacetsLID_.get( nPoints_*chunkForFuncID(funcID) + cellLID );
01274 facetDim = HN_To_globalFacetsDim_.get( nPoints_*chunkForFuncID(funcID) + cellLID );
01275 dofs.resize(facetLIDs.size());
01276 int b = chunkForFuncID(funcID);
01277 int func = indexForFuncID(funcID);
01278
01279 for (int ind = 0 ; ind < facetLIDs.size() ; ind++){
01280
01281 dofs[ind] = dofs_[facetDim[ind]][b]
01282 [facetLIDs[ind]*nDofsPerCell_[b][facetDim[ind]]+func*nNodesPerCell_[b][facetDim[ind]] + 0 ];
01283 }
01284
01285 coefs = HN_To_coeffs_.get( nPoints_*chunkForFuncID(funcID) + cellLID );
01286 SUNDANCE_MSG1(setupVerb(), "b=" << b);
01287 SUNDANCE_MSG1(setupVerb(), "func=" << func);
01288 SUNDANCE_MSG1(setupVerb(), "facetLIDs=" << facetLIDs);
01289 SUNDANCE_MSG1(setupVerb(), "facetDim = " << facetDim);
01290 SUNDANCE_MSG1(setupVerb(), "dofs=" << dofs);
01291 SUNDANCE_MSG1(setupVerb(), "coefs = " << coefs);
01292 }
01293 }
01294
01295
01296 void MixedDOFMapHN::computeOffsets(int dim, int localCount)
01297 {
01298 if (setupVerb() > 2)
01299 {
01300 comm().synchronize();
01301 comm().synchronize();
01302 comm().synchronize();
01303 comm().synchronize();
01304 }
01305 SUNDANCE_MSG2(setupVerb(),
01306 "p=" << mesh().comm().getRank()
01307 << " sharing offsets for DOF numbering for dim=" << dim);
01308
01309 SUNDANCE_MSG2(setupVerb(),
01310 "p=" << mesh().comm().getRank()
01311 << " I have " << localCount << " cells");
01312
01313 Array<int> dofOffsets;
01314 int totalDOFCount;
01315 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
01316 mesh().comm());
01317 int myOffset = dofOffsets[mesh().comm().getRank()];
01318
01319 SUNDANCE_MSG2(setupVerb(),
01320 "p=" << mesh().comm().getRank()
01321 << " back from MPI accumulate");
01322
01323 if (setupVerb() > 2)
01324 {
01325 comm().synchronize();
01326 comm().synchronize();
01327 comm().synchronize();
01328 comm().synchronize();
01329 }
01330
01331 for (int chunk=0; chunk<dofs_[dim].size(); chunk++)
01332 {
01333 for (int n=0; n<dofs_[dim][chunk].size(); n++)
01334 {
01335
01336 if (dofs_[dim][chunk][n] >= 0) {
01337 dofs_[dim][chunk][n] += myOffset;
01338 }
01339 }
01340 }
01341
01342 setLowestLocalDOF(myOffset);
01343 setNumLocalDOFs(localCount);
01344 setTotalNumDOFs(totalDOFCount);
01345
01346 SUNDANCE_MSG2(setupVerb(),
01347 "p=" << mesh().comm().getRank()
01348 << " done sharing offsets for DOF numbering for dim=" << dim);
01349 if (setupVerb() > 2)
01350 {
01351 comm().synchronize();
01352 comm().synchronize();
01353 comm().synchronize();
01354 comm().synchronize();
01355 }
01356
01357 }
01358
01359
01360
01361 void MixedDOFMapHN::checkTable() const
01362 {
01363 int bad = 0;
01364 for (int d=0; d<dofs_.size(); d++)
01365 {
01366 for (int chunk=0; chunk<dofs_[d].size(); chunk++)
01367 {
01368 const Array<int>& dofs = dofs_[d][chunk];
01369 for (int n=0; n<dofs.size(); n++)
01370 {
01371 if (dofs[n] < 0) bad = 1;
01372 }
01373 }
01374 }
01375
01376 int anyBad = bad;
01377 comm().allReduce((void*) &bad, (void*) &anyBad, 1,
01378 MPIDataType::intType(), MPIOp::sumOp());
01379 TEUCHOS_TEST_FOR_EXCEPTION(anyBad > 0, std::runtime_error, "invalid DOF map");
01380 }