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