Sundance 2.1

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