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 "SundanceMap.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceOrderedTuple.hpp"
00035 #include "SundanceInhomogeneousNodalDOFMap.hpp"
00036 #include "SundanceLagrange.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042 using Playa::MPIComm;
00043 using Playa::MPIContainerComm;
00044
00045
00046 InhomogeneousNodalDOFMap
00047 ::InhomogeneousNodalDOFMap(const Mesh& mesh,
00048 const Array<Map<Set<int>, CellFilter> >& funcSetToDomainMap,
00049 int setupVerb)
00050 : DOFMapBase(mesh, setupVerb),
00051 dim_(mesh.spatialDim()),
00052 basis_(rcp(new Lagrange(1))),
00053 nTotalFuncs_(),
00054 funcDomains_(),
00055 nodeDofs_(),
00056 elemDofs_(),
00057 nodeToFuncSetIndexMap_(mesh.numCells(0)),
00058 elemToFuncSetIndexMap_(mesh.numCells(dim_)),
00059 elemFuncSets_(),
00060 nodalFuncSets_(),
00061 nodeToOffsetMap_(mesh.numCells(0)),
00062 elemToOffsetMap_(mesh.numCells(0)),
00063 funcIndexWithinNodeFuncSet_(),
00064 elemStructure_(),
00065 nodeStructure_()
00066 {
00067 SUNDANCE_MSG1(setupVerb, "in InhomogeneousNodalDOFMap ctor");
00068 SUNDANCE_MSG4(setupVerb, "func set to domain map " << funcSetToDomainMap);
00069
00070
00071 Set<int> allFuncs;
00072 for (int d=dim_; d>=0; d--)
00073 {
00074 for (Map<Set<int>, CellFilter>::const_iterator
00075 i=funcSetToDomainMap[d].begin(); i!=funcSetToDomainMap[d].end(); i++)
00076 {
00077 allFuncs.merge(i->first);
00078 }
00079 }
00080
00081
00082 nTotalFuncs_ = allFuncs.size();
00083 SUNDANCE_MSG2(setupVerb, "found " << nTotalFuncs_ << " functions");
00084
00085
00086
00087 Array<CellFilter> subdomains;
00088 Array<CellFilter> maxSubdomains;
00089 Array<Array<int> > funcArrays;
00090 Array<Set<int> > funcSets;
00091
00092 for (int d=dim_; d>=0; d--)
00093 {
00094 for (Map<Set<int>, CellFilter>::const_iterator
00095 i=funcSetToDomainMap[d].begin(); i!=funcSetToDomainMap[d].end(); i++)
00096 {
00097 subdomains.append(i->second);
00098 if (d==dim_)
00099 {
00100 elemFuncSets_.append(i->first);
00101 maxSubdomains.append(i->second);
00102 }
00103
00104 Array<int> myFuncs = i->first.elements();
00105 funcArrays.append(myFuncs);
00106 funcSets.append(i->first);
00107 }
00108 }
00109
00110
00111
00112
00113
00114
00115 Array<Array<int> > cellLIDs(subdomains.size());
00116 Array<Array<int> > facetLIDs(subdomains.size());
00117 Array<Array<int> > facetOrientations(subdomains.size());
00118 Array<Set<int> > nodeToFuncSetMap(mesh.numCells(0));
00119
00120 for (int r=0; r<subdomains.size(); r++)
00121 {
00122 int d = subdomains[r].dimension(mesh);
00123 CellSet cells = subdomains[r].getCells(mesh);
00124 SUNDANCE_MSG4(setupVerb, "domain " << subdomains[r] << " has functions "
00125 << funcSets[r]);
00126
00127 for (CellIterator c=cells.begin(); c!=cells.end(); c++)
00128 {
00129 int cellLID = *c;
00130 cellLIDs[r].append(cellLID);
00131 }
00132 if (cellLIDs[r].size()==0) continue;
00133 int nFacets = mesh.numFacets(d, cellLIDs[r][0], 0);
00134 if (d==0)
00135 {
00136 for (int c=0; c<cellLIDs[r].size(); c++)
00137 {
00138 int fLID = cellLIDs[r][c];
00139 nodeToFuncSetMap[fLID].merge(funcSets[r]);
00140 }
00141 }
00142 else
00143 {
00144 Array<int>& facetLID = facetLIDs[r];
00145 facetLID.resize(nFacets*cellLIDs[r].size());
00146 facetOrientations[r].resize(nFacets*cellLIDs[r].size());
00147 mesh.getFacetLIDs(d, cellLIDs[r], 0, facetLID, facetOrientations[r]);
00148 for (int c=0; c<cellLIDs[r].size(); c++)
00149 {
00150 for (int f=0; f<nFacets; f++)
00151 {
00152 int fLID = facetLID[c*nFacets+f];
00153 nodeToFuncSetMap[fLID].merge(funcSets[r]);
00154 }
00155 }
00156 }
00157 }
00158
00159
00160
00161 Map<Set<int>, int> fsToFSIndexMap;
00162 Array<int> funcSetNodeCounts;
00163 int nNodes = mesh.numCells(0);
00164
00165 for (int n=0; n<nNodes; n++)
00166 {
00167 const Set<int>& f = nodeToFuncSetMap[n];
00168 int funcComboIndex;
00169 if (!fsToFSIndexMap.containsKey(f))
00170 {
00171 funcComboIndex = nodalFuncSets_.size();
00172 fsToFSIndexMap.put(f, funcComboIndex);
00173 nodalFuncSets_.append(f);
00174 funcSetNodeCounts.append(1);
00175 nodeToOffsetMap_[n] = 0;
00176 }
00177 else
00178 {
00179 funcComboIndex = fsToFSIndexMap.get(f);
00180 nodeToOffsetMap_[n] = f.size() * funcSetNodeCounts[funcComboIndex];
00181 funcSetNodeCounts[funcComboIndex]++;
00182 }
00183 nodeToFuncSetIndexMap_[n] = funcComboIndex;
00184 }
00185
00186 SUNDANCE_MSG2(setupVerb, "nodal func sets = " << nodalFuncSets_);
00187
00188 nodeDofs_.resize(nodalFuncSets_.size());
00189 nodeStructure_.resize(nodalFuncSets_.size());
00190 funcIndexWithinNodeFuncSet_.resize(nodalFuncSets_.size());
00191
00192 for (int f=0; f<nodalFuncSets_.size(); f++)
00193 {
00194 Array<int> funcs = nodalFuncSets_[f].elements();
00195 int nFuncs = funcs.size();
00196 funcIndexWithinNodeFuncSet_[f].resize(nTotalFuncs_, -1);
00197 for (int i=0; i<nFuncs; i++)
00198 {
00199 funcIndexWithinNodeFuncSet_[f][funcs[i]] = i;
00200 }
00201 nodeDofs_[f].resize(nFuncs * funcSetNodeCounts[f], -1);
00202 Array<RCP<BasisDOFTopologyBase> > nodeBases = tuple(basis_);
00203 Array<Array<int> > nodeFuncs = tuple(nodalFuncSets_[f].elements());
00204 nodeStructure_[f] = rcp(new MapStructure(nTotalFuncs_,
00205 nodeBases, nodeFuncs));
00206 }
00207
00208
00209
00210
00211 int nextDOF = 0;
00212 Array<int> hasProcessedNode(nNodes);
00213 Array<Array<int> > remoteNodes(mesh.comm().getNProc());
00214
00215 for (int r=0; r<subdomains.size(); r++)
00216 {
00217 int d = subdomains[r].dimension(mesh);
00218
00219 if (cellLIDs[r].size()==0) continue;
00220
00221 if (d==0)
00222 {
00223 for (int c=0; c<cellLIDs[r].size(); c++)
00224 {
00225 int fLID = cellLIDs[r][c];
00226 int funcSetIndex = nodeToFuncSetIndexMap_[fLID];
00227 int nFuncs = nodalFuncSets_[funcSetIndex].size();
00228 int dofOffset = nodeToOffsetMap_[fLID];
00229 if (!hasProcessedNode[fLID])
00230 {
00231 assignNode(fLID, funcSetIndex, dofOffset, nFuncs,
00232 remoteNodes, hasProcessedNode, nextDOF);
00233 }
00234 }
00235 }
00236 else
00237 {
00238 Array<int>& facetLID = facetLIDs[r];
00239 int nFacets = mesh.numFacets(d, cellLIDs[r][0], 0);
00240 for (int c=0; c<cellLIDs[r].size(); c++)
00241 {
00242 for (int f=0; f<nFacets; f++)
00243 {
00244 int fLID = facetLID[c*nFacets+f];
00245 int funcSetIndex = nodeToFuncSetIndexMap_[fLID];
00246 int dofOffset = nodeToOffsetMap_[fLID];
00247 int nFuncs = nodalFuncSets_[funcSetIndex].size();
00248 if (!hasProcessedNode[fLID])
00249 {
00250 assignNode(fLID, funcSetIndex, dofOffset, nFuncs,
00251 remoteNodes, hasProcessedNode, nextDOF);
00252 }
00253 }
00254 }
00255 }
00256 }
00257
00258
00259
00260 int localCount = nextDOF;
00261 computeOffsets(localCount);
00262
00263
00264 shareRemoteDOFs(remoteNodes);
00265
00266
00267 elemDofs_.resize(maxSubdomains.size());
00268 int count = 0;
00269 Array<Array<int> > elemFuncArrays;
00270
00271 funcDomains_.resize(nTotalFuncs_);
00272
00273 for (int r=0; r<subdomains.size(); r++)
00274 {
00275 int d = subdomains[r].dimension(mesh);
00276 if (d != dim_) continue;
00277 if (cellLIDs[r].size()==0) continue;
00278 int nFacets = mesh.numFacets(d, cellLIDs[r][0], 0);
00279 Array<int>& facetLID = facetLIDs[r];
00280 Array<Array<int> > dofs(1);
00281 dofs[0].resize(funcArrays[r].size() * cellLIDs[r].size() * nFacets);
00282 getFunctionDofs(d, cellLIDs[r], facetLID, funcArrays[r], dofs);
00283 elemDofs_[count] = dofs[0];
00284 elemFuncArrays.append(funcArrays[r]);
00285 Array<RCP<BasisDOFTopologyBase> > elemBases = tuple(basis_);
00286 Array<Array<int> > elemFuncs = tuple(elemFuncSets_[count].elements());
00287 elemStructure_.append(rcp(new MapStructure(nTotalFuncs_,
00288 elemBases, elemFuncs)));
00289 for (int c=0; c<cellLIDs[r].size(); c++)
00290 {
00291 elemToFuncSetIndexMap_[cellLIDs[r][c]] = count;
00292 }
00293 for (int i=0; i<elemFuncs[0].size(); i++)
00294 {
00295 funcDomains_[elemFuncs[0][i]] = funcDomains_[elemFuncs[0][i]] + subdomains[r];
00296 }
00297 count++;
00298 }
00299
00300 }
00301
00302
00303 void InhomogeneousNodalDOFMap::getFunctionDofs(int cellDim,
00304 const Array<int>& cellLID,
00305 const Array<int>& facetLID,
00306 const Array<int>& funcs,
00307 Array<Array<int> >& dofs) const
00308 {
00309 Array<int>& dofChunk = dofs[0];
00310 dofChunk.clear();
00311
00312 if (cellDim != 0)
00313 {
00314 int nFacets = mesh().numFacets(cellDim, cellLID[0], 0);
00315 dofChunk.resize(funcs.size() * cellLID.size() * nFacets, -1);
00316
00317
00318 for (int c=0; c<cellLID.size(); c++)
00319 {
00320 for (int f=0; f<nFacets; f++)
00321 {
00322 int fLID = facetLID[c*nFacets + f];
00323 int fci = nodeToFuncSetIndexMap_[fLID];
00324 int nodeOffset = nodeToOffsetMap_[fLID];
00325 for (int i=0; i<funcs.size(); i++)
00326 {
00327 int funcIndex = funcIndexWithinNodeFuncSet_[fci][funcs[i]];
00328 if (funcIndex >= 0)
00329 {
00330 dofChunk[(c*funcs.size()+i)*nFacets + f]
00331 = nodeDofs_[fci][nodeOffset+funcIndex];
00332
00333 }
00334 }
00335 }
00336 }
00337 }
00338 else
00339 {
00340 dofChunk.resize(funcs.size() * cellLID.size(), -1);
00341
00342 for (int c=0; c<cellLID.size(); c++)
00343 {
00344 int fci = nodeToFuncSetIndexMap_[cellLID[c]];
00345 int nodeOffset = nodeToOffsetMap_[cellLID[c]];
00346 for (int i=0; i<funcs.size(); i++)
00347 {
00348 int funcIndex = funcIndexWithinNodeFuncSet_[fci][funcs[i]];
00349 if (funcIndex >= 0)
00350 {
00351 dofChunk[c*funcs.size()+ i]
00352 = nodeDofs_[fci][nodeOffset+funcIndex];
00353 }
00354 }
00355 }
00356 }
00357 }
00358
00359 void InhomogeneousNodalDOFMap::assignNode(int fLID,
00360 int funcSetIndex,
00361 int dofOffset,
00362 int nFuncs,
00363 Array<Array<int> >& remoteNodes,
00364 Array<int>& hasProcessedCell,
00365 int& nextDOF)
00366 {
00367
00368 int owner;
00369 if (isRemote(0, fLID, owner))
00370 {
00371 int facetGID
00372 = mesh().mapLIDToGID(0, fLID);
00373 remoteNodes[owner].append(facetGID);
00374
00375 }
00376 else
00377 {
00378
00379 Array<int>& myDofs = nodeDofs_[funcSetIndex];
00380 for (int i=0; i<nFuncs; i++)
00381 {
00382 myDofs[dofOffset + i] = nextDOF;
00383 nextDOF++;
00384 }
00385 }
00386 hasProcessedCell[fLID] = 1;
00387 }
00388
00389 void InhomogeneousNodalDOFMap::computeOffsets(int localCount)
00390 {
00391 TEUCHOS_TEST_FOR_EXCEPTION(mesh().comm().getNProc() != 1,
00392 std::runtime_error,
00393 "parallel inhomogeneous DOF maps not yet supported");
00394
00395 int totalDOFCount = localCount;
00396 int myOffset = 0;
00397 setLowestLocalDOF(myOffset);
00398 setNumLocalDOFs(localCount);
00399 setTotalNumDOFs(totalDOFCount);
00400 }
00401
00402 void InhomogeneousNodalDOFMap::shareRemoteDOFs(const Array<Array<int> >& remoteNodes)
00403 {
00404 TEUCHOS_TEST_FOR_EXCEPTION(mesh().comm().getNProc() != 1,
00405 std::runtime_error,
00406 "parallel inhomogeneous DOF maps not yet supported");
00407 }
00408
00409 RCP<const Set<int> >
00410 InhomogeneousNodalDOFMap::allowedFuncsOnCellBatch(int cellDim,
00411 const Array<int>& cellLID) const
00412 {
00413 Set<int> rtn;
00414
00415 if (cellDim==0)
00416 {
00417 rtn = nodalFuncSets_[nodeToFuncSetIndexMap_[cellLID[0]]];
00418 for (int c=0; c<cellLID.size(); c++)
00419 {
00420 const Set<int>& s = nodalFuncSets_[nodeToFuncSetIndexMap_[cellLID[c]]];
00421 rtn = rtn.intersection(s);
00422 }
00423 }
00424 else if (cellDim==dim_)
00425 {
00426 rtn = elemFuncSets_[elemToFuncSetIndexMap_[cellLID[0]]];
00427 for (int c=0; c<cellLID.size(); c++)
00428 {
00429 const Set<int>& s = elemFuncSets_[elemToFuncSetIndexMap_[cellLID[c]]];
00430 rtn = rtn.intersection(s);
00431 }
00432 }
00433 else
00434 {
00435 Array<int> facetLID;
00436 Array<int> facetOrientations;
00437 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLID, facetOrientations);
00438 rtn = nodalFuncSets_[nodeToFuncSetIndexMap_[facetLID[0]]];
00439 for (int f=0; f<facetLID.size(); f++)
00440 {
00441 const Set<int>& s = nodalFuncSets_[nodeToFuncSetIndexMap_[facetLID[f]]];
00442 rtn = rtn.intersection(s);
00443 }
00444 }
00445 return rcp(new Set<int>(rtn));
00446 }
00447
00448
00449 RCP<const MapStructure>
00450 InhomogeneousNodalDOFMap::getDOFsForCellBatch(int cellDim,
00451 const Array<int>& cellLID,
00452 const Set<int>& requestedFuncSet,
00453 Array<Array<int> >& dofs,
00454 Array<int>& nNodes,
00455 int verb) const
00456 {
00457 TimeMonitor timer(batchedDofLookupTimer());
00458 Tabs tab0;
00459
00460 SUNDANCE_MSG2(verb, tab0 << "in InhomNodalDOFMap::getDOFsForCellBatch()");
00461
00462 if (cellDim==0)
00463 {
00464 Tabs tab1;
00465 SUNDANCE_MSG2(verb, tab1 << "cell dim = " << cellDim);
00466 bool isHomogeneous = true;
00467 int firstFuncSet = nodeToFuncSetIndexMap_[cellLID[0]];
00468 for (int c=0; c<cellLID.size(); c++)
00469 {
00470 if (nodeToFuncSetIndexMap_[cellLID[c]] != firstFuncSet)
00471 {
00472 isHomogeneous = false;
00473 break;
00474 }
00475 }
00476
00477 dofs.resize(1);
00478 nNodes.resize(1);
00479 nNodes[0] = 1;
00480 Array<int> dummyFacets;
00481
00482 if (isHomogeneous)
00483 {
00484 const Set<int>& funcSet = nodalFuncSets_[firstFuncSet];
00485 TEUCHOS_TEST_FOR_EXCEPT(requestedFuncSet.setDifference(funcSet).size() != 0);
00486 Array<int> funcs = funcSet.elements();
00487 getFunctionDofs(cellDim, cellLID, dummyFacets, funcs, dofs);
00488 return nodeStructure_[firstFuncSet];
00489 }
00490 else
00491 {
00492 Array<int> funcs = requestedFuncSet.elements();
00493 getFunctionDofs(cellDim, cellLID, dummyFacets, funcs, dofs);
00494 RCP<const MapStructure> rtn
00495 = rcp(new MapStructure(nTotalFuncs_, basis_, tuple(funcs)));
00496 return rtn;
00497 }
00498 }
00499 else if (cellDim==dim_)
00500 {
00501 Tabs tab1;
00502 SUNDANCE_MSG2(verb, tab1 << "cell dim = " << cellDim);
00503 bool isHomogeneous = true;
00504 int firstFuncSet = elemToFuncSetIndexMap_[cellLID[0]];
00505 SUNDANCE_MSG2(verb, tab1 << "first func set = " << firstFuncSet);
00506 for (int c=0; c<cellLID.size(); c++)
00507 {
00508 if (elemToFuncSetIndexMap_[cellLID[c]] != firstFuncSet)
00509 {
00510 isHomogeneous = false;
00511 break;
00512 }
00513 }
00514
00515 dofs.resize(1);
00516 nNodes.resize(1);
00517 nNodes[0] = mesh().numFacets(cellDim, cellLID[0], 0);
00518
00519 if (isHomogeneous)
00520 {
00521 Tabs tab2;
00522 Array<int> facetLID;
00523 Array<int> facetOrientations;
00524 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLID, facetOrientations);
00525 if (verb >= 2)
00526 {
00527 Out::os() << tab2 << "cellLID = " << cellLID << std::endl;
00528 Out::os() << tab2 << "facetLID = " << facetLID << std::endl;
00529 Out::os() << tab2 << "elem func sets = " << elemFuncSets_ << std::endl;
00530 }
00531 const Set<int>& funcSet = elemFuncSets_[firstFuncSet];
00532 TEUCHOS_TEST_FOR_EXCEPT(requestedFuncSet.setDifference(funcSet).size() != 0);
00533 Array<int> funcs = funcSet.elements();
00534 getFunctionDofs(cellDim, cellLID, facetLID, funcs, dofs);
00535 return elemStructure_[firstFuncSet];
00536 }
00537 else
00538 {
00539 Array<int> facetLID;
00540 Array<int> facetOrientations;
00541 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLID, facetOrientations);
00542 Array<int> funcs = requestedFuncSet.elements();
00543 getFunctionDofs(cellDim, cellLID, facetLID, funcs, dofs);
00544 RCP<const MapStructure> rtn
00545 = rcp(new MapStructure(nTotalFuncs_, basis_, tuple(funcs)));
00546 return rtn;
00547 }
00548 }
00549 else
00550 {
00551 Tabs tab1;
00552 SUNDANCE_MSG2(verb, tab1 << "cell dim = " << cellDim);
00553 Array<int> facetLID;
00554 Array<int> facetOrientations;
00555 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLID, facetOrientations);
00556 dofs.resize(1);
00557 nNodes.resize(1);
00558 nNodes[0] = mesh().numFacets(cellDim, cellLID[0], 0);
00559
00560 Array<int> funcs = requestedFuncSet.elements();
00561
00562 getFunctionDofs(cellDim, cellLID, facetLID, funcs, dofs);
00563 RCP<const MapStructure> rtn
00564 = rcp(new MapStructure(nTotalFuncs_, basis_, tuple(funcs)));
00565 return rtn;
00566 }
00567 }
00568
00569
00570
00571 void InhomogeneousNodalDOFMap::print(std::ostream& os) const
00572 {
00573 Tabs tab0;
00574 std::cout << tab0 << "dof map: " << std::endl;
00575 for (int d=0; d<=dim_; d++)
00576 {
00577 Tabs tab1;
00578 std::cout << tab1 << d << "-cells: " << std::endl;
00579 for (int n=0; n<mesh().numCells(d); n++)
00580 {
00581 Tabs tab2;
00582 std::cout << tab2 << "d=" << d << " cell=" << n << std::endl;
00583 RCP<const Set<int> > funcs = allowedFuncsOnCellBatch(d, tuple(n));
00584 for (Set<int>::const_iterator f=funcs->begin(); f!=funcs->end(); f++)
00585 {
00586 Tabs tab3;
00587 Array<int> dofs;
00588 getDOFsForCell(d, n, *f, dofs);
00589 std::cout << tab3 << " f=" << *f << " dofs=" << dofs << std::endl;
00590 }
00591 }
00592 }
00593 }
