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 "SundanceHomogeneousDOFMap.hpp"
00032 #include "SundanceCellFilter.hpp"
00033 #include "SundanceMaximalCellFilter.hpp"
00034 #include "PlayaMPIContainerComm.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "PlayaTabs.hpp"
00037 #include "Teuchos_Time.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 static Time& dofLookupTimer() 
00047 {
00048   static RCP<Time> rtn 
00049     = TimeMonitor::getNewTimer("unbatched dof lookup"); 
00050   return *rtn;
00051 }
00052 
00053 static Time& dofBatchLookupTimer() 
00054 {
00055   static RCP<Time> rtn 
00056     = TimeMonitor::getNewTimer("batched dof lookup"); 
00057   return *rtn;
00058 }
00059 
00060 HomogeneousDOFMap::HomogeneousDOFMap(const Mesh& mesh, 
00061   const BasisFamily& basis,
00062   int numFuncs, 
00063   int setupVerb)
00064   : DOFMapBase(mesh, setupVerb), 
00065     dim_(mesh.spatialDim()),
00066     dofs_(mesh.spatialDim()+1),
00067     maximalDofs_(),
00068     haveMaximalDofs_(false),
00069     localNodePtrs_(mesh.spatialDim()+1),
00070     nNodesPerCell_(mesh.spatialDim()+1),
00071     totalNNodesPerCell_(mesh.spatialDim()+1, 0),
00072     numFacets_(mesh.spatialDim()+1),
00073     originalFacetOrientation_(2),
00074     basisIsContinuous_(false)
00075 {
00076   verbosity() = DOFMapBase::classVerbosity();
00077   
00078   CellFilter maximalCells = new MaximalCellFilter();
00079   cellSets().append(maximalCells.getCells(mesh));
00080   cellDimOnCellSets().append(mesh.spatialDim());
00081 
00082   allocate(mesh, basis, numFuncs);
00083   initMap();
00084 }
00085 
00086 
00087 HomogeneousDOFMap::HomogeneousDOFMap(const Mesh& mesh, 
00088                                      const BasisFamily& basis,
00089                                      const Array<CellFilter>& subregions,
00090                                      int numFuncs)
00091   : DOFMapBase(mesh), 
00092     dim_(mesh.spatialDim()),
00093     dofs_(mesh.spatialDim()+1),
00094     maximalDofs_(),
00095     haveMaximalDofs_(false),
00096     localNodePtrs_(mesh.spatialDim()+1),
00097     nNodesPerCell_(mesh.spatialDim()+1),
00098     totalNNodesPerCell_(mesh.spatialDim()+1, 0),
00099     numFacets_(mesh.spatialDim()+1),
00100     originalFacetOrientation_(2),
00101     basisIsContinuous_(false)
00102 {
00103   verbosity() = DOFMapBase::classVerbosity();
00104   
00105   for (int r=0; r<subregions.size(); r++)
00106     {
00107       cellSets().append(subregions[r].getCells(mesh));
00108       cellDimOnCellSets().append(subregions[r].dimension(mesh));
00109     }
00110 
00111   allocate(mesh, basis, numFuncs);
00112   initMap();
00113 }
00114 
00115 
00116 void HomogeneousDOFMap::allocate(const Mesh& mesh, 
00117                                  const BasisFamily& basis,
00118                                  int numFuncs)
00119 {
00120   Tabs tab;
00121   SUNDANCE_MSG1(setupVerb(), tab << "allocating DOF map for nFuncs=" << numFuncs);
00122   Array<int> fid(numFuncs);
00123   for (int f=0; f<numFuncs; f++) fid[f] = f;
00124   funcIDOnCellSets().append(fid);
00125 
00126 
00127   
00128   for (int d=0; d<=dim_; d++)
00129     {
00130       Tabs tab1;
00131       SUNDANCE_MSG2(setupVerb(), tab1 << "allocating d=" << d);
00132       /* record the number of facets for each cell type so we're
00133        * not making a bunch of mesh calls */
00134       numFacets_[d].resize(d);
00135       for (int fd=0; fd<d; fd++) numFacets_[d][fd]=mesh.numFacets(d, 0, fd);
00136       SUNDANCE_MSG3(setupVerb(), tab1 << "num facets for dimension " << d << " is " 
00137                          << numFacets_[d]);
00138           
00139       /* look up the node pointer for this cell and for all of its
00140        * facets */
00141       basis.ptr()->getLocalDOFs(mesh.cellType(d), localNodePtrs_[d]);
00142 
00143 
00144       SUNDANCE_MSG3(setupVerb(), tab1 << "node ptrs for dimension " << d << " are " 
00145                          << localNodePtrs_[d]);
00146 
00147       /* with the node pointers in hand, we can work out the number
00148        * of nodes per cell in this dimension */
00149       if (localNodePtrs_[d][d].size() > 0) 
00150         {
00151           nNodesPerCell_[d] = localNodePtrs_[d][d][0].size();
00152         }
00153       else
00154         {
00155           nNodesPerCell_[d] = 0;
00156         }
00157       SUNDANCE_MSG3(setupVerb(), tab1 << 
00158                          "num nodes for dimension " << d << " is " 
00159                          << nNodesPerCell_[d]);
00160 
00161       totalNNodesPerCell_[d] = nNodesPerCell_[d];
00162       for (int dd=0; dd<d; dd++) 
00163         {
00164           totalNNodesPerCell_[d] += numFacets_[d][dd]*nNodesPerCell_[dd];
00165         }
00166 
00167       /* we know from the mesh the number of cells in this dimension */
00168       if (nNodesPerCell_[d] > 0)
00169         {
00170           dofs_[d].resize(mesh.numCells(d));
00171         }
00172       else
00173         {
00174           dofs_[d].resize(0);
00175         }
00176 
00177       if (d > 0 && d < dim_) originalFacetOrientation_[d-1].resize(mesh.numCells(d));
00178 
00179       /* If any nodes are associated with the facets, then we know we have
00180        * a continuous basis function */
00181       if (d < dim_ && nNodesPerCell_[d] > 0) basisIsContinuous_ = true;
00182 
00183 
00184       /* now that we know the number of nodes per cell for this dimension,
00185        * we can allocate space for the DOFs in this dimension */
00186       int numCells = dofs_[d].size();
00187       for (int c=0; c<numCells; c++)
00188         {
00189           dofs_[d][c].resize(funcIDList().size() * nNodesPerCell_[d]);
00190           /* set everything to uninitializedVal() */
00191           for (int i=0; i<dofs_[d][c].size(); i++) 
00192             {
00193               dofs_[d][c][i] = uninitializedVal();
00194             }
00195         }
00196     }
00197   SUNDANCE_MSG1(setupVerb(), tab << "done allocating DOF map");
00198 }
00199 
00200 void HomogeneousDOFMap::initMap()
00201 {
00202   Tabs tab;
00203   SUNDANCE_MSG1(setupVerb(), tab << "initializing DOF map");
00204   /* start the DOF count at zero. */
00205   int nextDOF = 0;
00206 
00207   /* Space in which to keep a list of remote cells needed by each processor
00208    * for each dimension. The first index is dimension, the second proc, the
00209    * third cell number. */
00210   Array<Array<Array<int> > > remoteCells(mesh().spatialDim()+1,
00211                                          mesh().comm().getNProc());
00212   
00213   for (int r=0; r<numCellSets(); r++)
00214     {
00215       /* Loop over maximal cells in the order specified by the cell iterator.
00216        * This might be reordered relative to the mesh. 
00217        *
00218        * At each maximal cell, we'll run through the facets and 
00219        * assign DOFs. That will take somewhat more work, but gives much 
00220        * better cache locality for the matrix because all the DOFs for
00221        * each maximal element and its facets are grouped together. */
00222 
00223       CellSet cells = cellSet(r);
00224       CellIterator iter;
00225       for (iter=cells.begin(); iter != cells.end(); iter++)
00226         {
00227           /* first assign any DOFs associated with the maximal cell */
00228           int cellLID = *iter;
00229           int owner;
00230       
00231           if (nNodesPerCell_[dim_] > 0)
00232             {
00233               /* if the maximal cell is owned by another processor,
00234                * put it in the list of cells for which we need to request 
00235                * DOF information from another processor */
00236               if (isRemote(dim_, cellLID, owner))
00237                 {
00238                   int dummy=0;
00239                   int cellGID = mesh().mapLIDToGID(dim_, cellLID);
00240                   remoteCells[dim_][owner].append(cellGID); 
00241                   setDOFs(dim_, cellLID, dummy);
00242                 }
00243               else /* the cell is locally owned, so we can 
00244                     * set its DOF numbers now */
00245                 {
00246                   setDOFs(dim_, cellLID, nextDOF);
00247                 }
00248             }
00249 
00250           /* Now assign any DOFs associated with the facets. */
00251           /* We can skip this step if the basis is discontinuous at element
00252            * boundaries, because the facets will own no nodes */
00253           if (basisIsContinuous_)
00254             {
00255               for (int d=0; d<dim_; d++)
00256                 {
00257                   if (nNodesPerCell_[d] > 0)
00258                     {
00259                       int nf = numFacets_[dim_][d];
00260                       Array<int> facetLID(nf);
00261                       Array<int> facetOrientations(nf);
00262                       /* look up the LIDs of the facets */
00263                       mesh().getFacetArray(dim_, cellLID, d, 
00264                                            facetLID, facetOrientations);
00265                       /* for each facet, process its DOFs */
00266                       for (int f=0; f<nf; f++)
00267                         {
00268                           /* if the facet's DOFs have been assigned already,
00269                            * we're done */
00270                           if (!hasBeenAssigned(d, facetLID[f]))
00271                             {
00272                               /* the facet may be owned by another processor */
00273                               if (isRemote(d, facetLID[f], owner))
00274                                 {
00275                                   int dummy=0;
00276                                   int facetGID = mesh().mapLIDToGID(d, facetLID[f]);
00277                                   remoteCells[d][owner].append(facetGID);
00278                                   setDOFs(d, facetLID[f], dummy);
00279                                 }
00280                               else /* we can assign a DOF locally */
00281                                 {
00282                                   /* assign DOF */
00283                                   setDOFs(d, facetLID[f], nextDOF);
00284                                   /* record the orientation wrt the maximal cell */
00285                                   if (d > 0) 
00286                                     {
00287                                       originalFacetOrientation_[d-1][facetLID[f]] 
00288                                         = facetOrientations[f];
00289                                     }
00290                                 }
00291                             }
00292                         }
00293                     }
00294                 }
00295             }
00296         }
00297     }
00298   /* Done with first pass, in which we have assigned DOFs for all
00299    * local processors. We now have to share DOF information between
00300    * processors */
00301 
00302   if (mesh().comm().getNProc() > 1)
00303     {
00304       for (int d=0; d<=dim_; d++)
00305         {
00306           if (nNodesPerCell_[d] > 0)
00307             {
00308               computeOffsets(d, nextDOF);
00309               shareDOFs(d, remoteCells[d]);
00310             }
00311         }
00312     }
00313   else
00314     {
00315       setLowestLocalDOF(0);
00316       setNumLocalDOFs(nextDOF);
00317       setTotalNumDOFs(nextDOF);
00318     }
00319   SUNDANCE_MSG1(setupVerb(), tab << "done initializing DOF map");
00320 }
00321 
00322 void HomogeneousDOFMap::shareDOFs(int cellDim,
00323                                   const Array<Array<int> >& outgoingCellRequests)
00324 {
00325   int np = mesh().comm().getNProc();
00326   int rank = mesh().comm().getRank();
00327 
00328   if (np == 1) return;
00329 
00330   TEUCHOS_TEST_FOR_EXCEPTION(np != outgoingCellRequests.size(), std::runtime_error,
00331                      "incorrect size of outgoingCellRequests array. Size is "
00332                      << outgoingCellRequests.size() << ", should be np=" << np);
00333 
00334 
00335 
00336   Array<Array<int> > incomingCellRequests;
00337   Array<Array<int> > outgoingDOFs(np);
00338   Array<Array<int> > incomingDOFs;
00339 
00340   SUNDANCE_MSG2(setupVerb(),  
00341                "p=" << mesh().comm().getRank()
00342                << "synchronizing DOFs for cells of dimension " << cellDim);
00343   SUNDANCE_MSG2(setupVerb(),  
00344                "p=" << mesh().comm().getRank()
00345                << " sending cell reqs d=" << cellDim << " GID=" << outgoingCellRequests);
00346 
00347   /* share the cell requests */
00348   MPIContainerComm<int>::allToAll(outgoingCellRequests, 
00349                                   incomingCellRequests,
00350                                   mesh().comm());
00351 
00352   
00353 
00354   SUNDANCE_MSG2(setupVerb(),  
00355                "p=" << rank
00356                << "recvd DOF requests for cells of dimension " << cellDim
00357                << " GID=" << incomingCellRequests);
00358 
00359   /* get orientations and DOF numbers for the first node of every cell that's been 
00360    * requested by someone else */
00361   for (int p=0; p<np; p++)
00362     {
00363       if (p==rank) continue;
00364       const Array<int>& requestsFromProc = incomingCellRequests[p];
00365       int nReq = requestsFromProc.size();
00366       int blockSize = 1;
00367       if (cellDim > 0 && cellDim < dim_)
00368         {
00369           outgoingDOFs[p].resize(2 * nReq);
00370           blockSize = 2;
00371         }
00372       else
00373         {
00374           outgoingDOFs[p].resize(nReq);
00375         }
00376       SUNDANCE_MSG3(setupVerb(),  
00377                    "p=" << mesh().comm().getRank() << " recv'd from proc=" << p
00378                    << " reqs for DOFs for cells " << requestsFromProc);
00379       for (int c=0; c<nReq; c++)
00380         {
00381           int GID = requestsFromProc[c];
00382           SUNDANCE_MSG3(setupVerb(),  
00383                        "p=" << rank
00384                        << " processing cell with d=" << cellDim 
00385                        << " GID=" << GID);
00386           int LID = mesh().mapGIDToLID(cellDim, GID);
00387           SUNDANCE_MSG3(setupVerb(),  
00388                        "p=" << rank
00389                        << " LID=" << LID << " dofs=" << dofs_[cellDim]);
00390           outgoingDOFs[p][blockSize*c] = dofs_[cellDim][LID][0];
00391           if (cellDim > 0 && cellDim < dim_) 
00392             {
00393               outgoingDOFs[p][blockSize*c+1] = originalFacetOrientation_[cellDim-1][LID];
00394             }
00395           SUNDANCE_MSG3(setupVerb(),  
00396                        "p=" << rank
00397                        << " done processing cell with GID=" << GID);
00398         }
00399     }
00400  
00401   
00402 
00403   SUNDANCE_MSG2(setupVerb(),  
00404                "p=" << mesh().comm().getRank()
00405                << "answering DOF requests for cells of dimension " << cellDim);
00406 
00407   /* share the DOF numbers */
00408   MPIContainerComm<int>::allToAll(outgoingDOFs,
00409                                   incomingDOFs,
00410                                   mesh().comm());
00411 
00412   SUNDANCE_MSG2(setupVerb(),  
00413                "p=" << mesh().comm().getRank()
00414                << "communicated DOF answers for cells of dimension " << cellDim);
00415 
00416   
00417   /* now assign the DOFs from the other procs */
00418   for (int p=0; p<mesh().comm().getNProc(); p++)
00419     {
00420       if (p==mesh().comm().getRank()) continue;
00421       const Array<int>& dofsFromProc = incomingDOFs[p];
00422       int blockSize = 1;
00423       if (cellDim > 0 && cellDim < dim_)
00424         {
00425           blockSize = 2;
00426         }
00427 
00428       for (int c=0; c<dofsFromProc.size()/blockSize; c++)
00429         {
00430           int cellGID = outgoingCellRequests[p][c];
00431           int cellLID = mesh().mapGIDToLID(cellDim, cellGID);
00432           int dof = dofsFromProc[blockSize*c];
00433           setDOFs(cellDim, cellLID, dof, true);
00434           if (cellDim > 0 && cellDim < dim_) 
00435             {
00436               originalFacetOrientation_[cellDim-1][cellLID] = dofsFromProc[blockSize*c+1];
00437             }
00438         }
00439     }
00440   
00441 }
00442 
00443 
00444 
00445 void HomogeneousDOFMap::setDOFs(int cellDim, int cellLID, 
00446                                 int& nextDOF,
00447                                 bool isRemote)
00448 {
00449   Tabs tab;
00450   SUNDANCE_MSG3(setupVerb(), tab << "setting DOFs for " << cellDim << "-cell " << cellLID);
00451   Array<int>& cellDOFs = dofs_[cellDim][cellLID];
00452   
00453   int nn = nNodesPerCell_[cellDim];
00454   int nf = funcIDList().size();
00455 
00456   for (int i=0; i<nn; i++)
00457     {
00458       for (int f=0; f<nf; f++)
00459         {
00460           int k = nextDOF++;
00461           cellDOFs[funcIDList()[f] + nf*i] = k;
00462           if (isRemote) addGhostIndex(k);
00463         }
00464     }
00465 }
00466 
00467 
00468 
00469 void HomogeneousDOFMap::getDOFsForCellBatch(int cellDim, 
00470                                             const Array<int>& cellLID,
00471                                             Array<int>& dofs,
00472                                             int& nNodes) const 
00473 {
00474   TimeMonitor timer(dofBatchLookupTimer());
00475 
00476   Tabs tab;
00477   SUNDANCE_MSG3(setupVerb(), 
00478                tab << "getDOFsForCellBatch(): cellDim=" << cellDim
00479                << " cellLID=" << cellLID);
00480 
00481   if (cellLID.size()==0) return;
00482 
00483   if (cellDim == dim_)
00484     {
00485       if (!haveMaximalDofs_) 
00486         {
00487           buildMaximalDofTable();
00488         }
00489 
00490       SUNDANCE_MSG4(setupVerb(), tab << "getting dofs for maximal cells");
00491 
00492       nNodes = totalNNodesPerCell_[cellDim];
00493       int nCells = cellLID.size();
00494       int nTotalCells = mesh().numCells(cellDim);
00495       int nf = funcIDList().size();
00496       dofs.resize(totalNNodesPerCell_[cellDim] * cellLID.size() * nf);
00497       
00498       for (int c=0; c<cellLID.size(); c++)
00499         {
00500           Tabs tab1;
00501           SUNDANCE_MSG4(setupVerb(), tab1 << "cell=" << c);
00502           for (int fid=0; fid<nf; fid++)
00503             {
00504               Tabs tab2;
00505               SUNDANCE_MSG4(setupVerb(), tab2 << "f= " << fid);
00506               for (int n=0; n<nNodes; n++) 
00507                 {
00508                   Tabs tab3;
00509                   dofs[(fid*nCells + c)*nNodes + n] 
00510                     = maximalDofs_[(fid*nTotalCells + cellLID[c])*nNodes + n];
00511                   SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << n << " dof=" 
00512                                         << dofs[(fid*nCells + c)*nNodes + n]);
00513                 }
00514               
00515             }
00516         }
00517     }
00518   else
00519     {
00520       int nf = funcIDList().size();
00521       int nCells = cellLID.size();
00522       dofs.resize(totalNNodesPerCell_[cellDim] * cellLID.size() * nf);
00523 
00524       Tabs tab1;
00525       SUNDANCE_MSG4(setupVerb(), tab1 << "getting dofs for non-maximal cells");
00526   
00527       static Array<Array<int> > facetLID(3);
00528       static Array<Array<int> > facetOrientations(3);
00529       static Array<int> numFacets(3);
00530 
00531       for (int d=0; d<cellDim; d++) 
00532         {
00533           numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00534           mesh().getFacetLIDs(cellDim, cellLID, d, facetLID[d], 
00535                               facetOrientations[d]);
00536         }
00537 
00538   
00539       int nInteriorNodes = localNodePtrs_[cellDim][cellDim][0].size();
00540   
00541       nNodes = totalNNodesPerCell_[cellDim];
00542       
00543       for (int c=0; c<cellLID.size(); c++)
00544         {
00545           Tabs tab2;
00546           SUNDANCE_MSG4(setupVerb(), tab2 << "cell=" << c);
00547           SUNDANCE_MSG4(setupVerb(), tab2 << "doing interior nodes");
00548 
00549           /* first get the DOFs for the nodes associated with 
00550            * the cell's interior */
00551           if (cellDim==0 || nInteriorNodes <= 1) /* orientation-independent */
00552             {
00553               for (int m=0; m<nInteriorNodes; m++)
00554                 {
00555                   Tabs tab3;
00556                   int ptr = localNodePtrs_[cellDim][cellDim][0][m];
00557                   for (int f=0; f<nf; f++)
00558                     {
00559                       dofs[(f*nCells + c)*nNodes+ptr] 
00560                         = dofs_[cellDim][cellLID[c]][f + nf*m];
00561                       SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << m << " f=" << f 
00562                                             << " dof=" 
00563                                             <<  dofs[(f*nCells + c)*nNodes+ptr]);
00564                     }
00565                 }
00566             }
00567           else
00568             {
00569               int sign = originalFacetOrientation_[cellDim-1][cellLID[c]];
00570               for (int m=0; m<nInteriorNodes; m++)
00571                 {
00572                   int n = m;
00573                   if (sign < 0) n = nInteriorNodes - 1 - m;
00574                   Tabs tab3;
00575                   int ptr = localNodePtrs_[cellDim][cellDim][0][m];
00576                   for (int f=0; f<nf; f++)
00577                     {
00578                       dofs[(f*nCells + c)*nNodes+ptr] 
00579                         = dofs_[cellDim][cellLID[c]][f + nf*n];
00580                       SUNDANCE_MSG4(setupVerb(), tab3 << "n=" << m << " f=" << f 
00581                                             << " dof=" 
00582                                             <<  dofs[(f*nCells + c)*nNodes+ptr]);
00583                     }
00584                 }
00585             }
00586 
00587           SUNDANCE_MSG4(setupVerb(), tab1 << "doing facet nodes");
00588           /* now get the DOFs for the nodes on the facets */
00589           for (int d=0; d<cellDim; d++)
00590             {
00591               Tabs tab3;
00592               SUNDANCE_MSG4(setupVerb(), tab3 << "d= " << d);
00593               for (int f=0; f<numFacets[d]; f++)
00594                 {
00595                   Tabs tab4;
00596                   int facetID = facetLID[d][c*numFacets[d]+f];
00597                   
00598                   SUNDANCE_MSG4(setupVerb(), tab4 << "f= " << f << " facetLID = " << facetID);
00599                   int nFacetNodes = localNodePtrs_[cellDim][d][f].size();
00600                   if (d==0 || nFacetNodes <= 1) /* orientation-independent */
00601                     {
00602                       Tabs tab5;
00603                       for (int n=0; n<nFacetNodes; n++)
00604                         {
00605                           SUNDANCE_MSG3(setupVerb(), 
00606                                        tab5 << "n=" << n);
00607                           int ptr = localNodePtrs_[cellDim][d][f][n];
00608                           SUNDANCE_MSG3(setupVerb(), 
00609                                        tab5 << "local ptr=" << ptr);
00610                           for (int funcID=0; funcID<nf; funcID++)
00611                             {
00612                               SUNDANCE_MSG3(setupVerb(), 
00613                                            tab5 << "found dof=" 
00614                                            << dofs_[d][facetID][funcID + nf*n]);
00615                               dofs[(funcID*nCells + c)*nNodes+ptr] 
00616                                 = dofs_[d][facetID][funcID + nf*n];
00617                             }
00618                         }
00619                     }
00620                   else /* orientation-dependent */
00621                     {
00622                       Tabs tab5;
00623                       int facetOrientation 
00624                         = facetOrientations[d][c*numFacets[d]+f] 
00625                         * originalFacetOrientation_[d-1][facetID];
00626                       for (int m=0; m<nFacetNodes; m++)
00627                         {
00628                           int n = m;
00629                           if (facetOrientation<0) n = nFacetNodes-1-m;
00630                           SUNDANCE_MSG3(setupVerb(), 
00631                                        tab5 << "n=" << n);
00632                           int ptr = localNodePtrs_[cellDim][d][f][m];
00633                           SUNDANCE_MSG3(setupVerb(), 
00634                                        tab5 << "local ptr=" << ptr);
00635                           for (int funcID=0; funcID<nf; funcID++)
00636                             {
00637                               SUNDANCE_MSG3(setupVerb(), 
00638                                            tab5 << "found dof=" 
00639                                            << dofs_[d][facetID][funcID + nf*n]);
00640                               dofs[(funcID*nCells + c)*nNodes+ptr] 
00641                                 = dofs_[d][facetID][funcID + nf*n];
00642                             }
00643                         }
00644                     }
00645                 }
00646             }
00647         }
00648     }
00649 }    
00650 
00651 void HomogeneousDOFMap::buildMaximalDofTable() const
00652 {
00653   Tabs tab;
00654   int cellDim = dim_;
00655   int nf = funcIDList().size();
00656   int nCells = mesh().numCells(dim_);
00657 
00658   SUNDANCE_MSG2(setupVerb(), tab << "building dofs for maximal cells");
00659 
00660   SUNDANCE_MSG3(setupVerb(), tab << "nf=" << nf
00661                << " total nNodes=" << totalNNodesPerCell_[cellDim]);
00662   
00663   static Array<Array<int> > facetLID(3);
00664   static Array<Array<int> > facetOrientations(3);
00665   static Array<int> numFacets(3);
00666 
00667   Array<int> cellLID(nCells);
00668 
00669   for (int c=0; c<cellLID.size(); c++) cellLID[c]=c;
00670   
00671   for (int d=0; d<cellDim; d++) 
00672     {
00673       numFacets[d] = mesh().numFacets(cellDim, cellLID[0], d);
00674       mesh().getFacetLIDs(cellDim, cellLID, d, facetLID[d], facetOrientations[d]);
00675     }
00676 
00677   
00678   int nInteriorNodes = localNodePtrs_[cellDim][cellDim][0].size();
00679   
00680 
00681   int nNodes = totalNNodesPerCell_[cellDim];
00682 
00683   maximalDofs_.resize(nCells*nf*nNodes);
00684 
00685   for (int c=0; c<nCells; c++)
00686     {
00687       Tabs tab1;
00688       SUNDANCE_MSG4(setupVerb(), tab1 << "working on cell=" << c << " LID=" << cellLID[c]);
00689       /* first get the DOFs for the nodes associated with 
00690        * the cell's interior */
00691       SUNDANCE_MSG4(setupVerb(), tab1 << "doing interior nodes");
00692       for (int n=0; n<nInteriorNodes; n++)
00693         {
00694           int ptr = localNodePtrs_[cellDim][cellDim][0][n];
00695           for (int f=0; f<nf; f++)
00696             {
00697               maximalDofs_[(f*nCells + c)*nNodes+ptr] 
00698                 = dofs_[cellDim][c][f + nf*n];
00699             }
00700         }
00701       
00702       SUNDANCE_MSG4(setupVerb(), tab1 << "doing facet nodes");
00703       /* now get the DOFs for the nodes on the facets */
00704       for (int d=0; d<cellDim; d++)
00705         {
00706           Tabs tab2;
00707           SUNDANCE_MSG4(setupVerb(), tab2 << "facet dim=" << d);
00708           for (int f=0; f<numFacets[d]; f++)
00709             {
00710               Tabs tab3;
00711               int facetID = facetLID[d][c*numFacets[d]+f];
00712               SUNDANCE_MSG4(setupVerb(), tab2 << "f=" << f << " facetLID=" << facetID);
00713               int nFacetNodes = localNodePtrs_[cellDim][d][f].size();
00714               if (d == 0 || nFacetNodes <= 1) /* orientation-independent */
00715                 {
00716                   for (int n=0; n<nFacetNodes; n++)
00717                     {
00718                       Tabs tab4;
00719                       SUNDANCE_MSG3(setupVerb(), 
00720                                    tab4 << "n=" << n);
00721                       int ptr = localNodePtrs_[cellDim][d][f][n];
00722                       for (int funcID=0; funcID<nf; funcID++)
00723                         {
00724                           SUNDANCE_MSG3(setupVerb(), 
00725                                        tab4 << "found dof=" 
00726                                        << dofs_[d][facetID][funcID + nf*n])
00727                             maximalDofs_[(funcID*nCells + c)*nNodes+ptr] 
00728                             = dofs_[d][facetID][funcID + nf*n];
00729                         }
00730                     }
00731                 }
00732               else /* orientation-dependent */
00733                 {
00734                   Tabs tab4;
00735                   int facetOrientation = facetOrientations[d][c*numFacets[d]+f]
00736                     * originalFacetOrientation_[d-1][facetID];
00737                   for (int m=0; m<nFacetNodes; m++)
00738                     {
00739                       Tabs tab4;
00740                       int n = m;
00741                       if (facetOrientation<0) n = nFacetNodes-1-m;
00742                       SUNDANCE_MSG3(setupVerb(), 
00743                                    tab4 << "n=" << n);
00744                       int ptr = localNodePtrs_[cellDim][d][f][m];
00745                       for (int funcID=0; funcID<nf; funcID++)
00746                         {
00747                           SUNDANCE_MSG3(setupVerb(), 
00748                                        tab4 << "found dof=" 
00749                                        << dofs_[d][facetID][funcID + nf*n])
00750                             maximalDofs_[(funcID*nCells + c)*nNodes+ptr] 
00751                             = dofs_[d][facetID][funcID + nf*n];
00752                         }
00753                     }
00754                 }
00755             }
00756         }
00757     }
00758   haveMaximalDofs_ = true;
00759 }
00760 
00761 
00762 
00763 
00764 
00765 void HomogeneousDOFMap::computeOffsets(int dim, int localCount)
00766 {
00767   if (verbosity() > 2)
00768     {
00769       comm().synchronize();
00770       comm().synchronize();
00771       comm().synchronize();
00772       comm().synchronize();
00773     }
00774   SUNDANCE_MSG2(setupVerb(), 
00775                "p=" << mesh().comm().getRank()
00776                << " sharing offsets for DOF numbering for dim=" << dim);
00777 
00778   SUNDANCE_MSG2(setupVerb(), 
00779                "p=" << mesh().comm().getRank()
00780                << " I have " << localCount << " cells");
00781 
00782   Array<int> dofOffsets;
00783   int totalDOFCount;
00784   MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00785                                     mesh().comm());
00786   int myOffset = dofOffsets[mesh().comm().getRank()];
00787 
00788   SUNDANCE_MSG2(setupVerb(), 
00789                "p=" << mesh().comm().getRank()
00790                << " back from MPI accumulate");
00791 
00792   if (verbosity() > 2)
00793     {
00794       comm().synchronize();
00795       comm().synchronize();
00796       comm().synchronize();
00797       comm().synchronize();
00798     }
00799 
00800   for (int c=0; c<dofs_[dim].size(); c++)
00801     {
00802       if (hasBeenAssigned(dim, c))
00803         {
00804           for (int n=0; n<dofs_[dim][c].size(); n++) 
00805             {
00806               dofs_[dim][c][n] += myOffset;
00807             }
00808         }
00809     }
00810 
00811   setLowestLocalDOF(myOffset);
00812   setNumLocalDOFs(localCount);
00813   setTotalNumDOFs(totalDOFCount);
00814 
00815   SUNDANCE_MSG2(setupVerb(), 
00816                "p=" << mesh().comm().getRank() 
00817                << " done sharing offsets for DOF numbering for dim=" << dim);
00818   if (verbosity() > 2)
00819     {
00820       comm().synchronize();
00821       comm().synchronize();
00822       comm().synchronize();
00823       comm().synchronize();
00824     }
00825 
00826 }                           
00827 
00828 
00829 
00830 void HomogeneousDOFMap::print(std::ostream& os) const
00831 {
00832   int myRank = mesh().comm().getRank();
00833 
00834   Tabs tabs;
00835 
00836   os << "DOFS = " << dofs_ << std::endl;
00837 
00838   for (int p=0; p<mesh().comm().getNProc(); p++)
00839     {
00840       mesh().comm().synchronize();
00841       mesh().comm().synchronize();
00842       if (p == myRank)
00843         {
00844           os << tabs << 
00845             "========= DOFMap on proc p=" << p << " =============" << std::endl;
00846           for (int d=dim_; d>=0; d--)
00847             {
00848               Tabs tabs1;
00849               os << tabs1 << "dimension = " << d << std::endl;
00850               for (int c=0; c<mesh().numCells(d); c++)
00851                 {
00852                   Tabs tabs2;
00853                   os << tabs2 << "Cell LID=" << c << " GID=" 
00854                      << mesh().mapLIDToGID(d, c) << std::endl;
00855                   for (int f=0; f<funcIDList().size(); f++)
00856                     {
00857                       Tabs tabs3;
00858                       Array<int> dofs;
00859                       getDOFsForCell(d, c, funcIDList()[f], dofs);
00860                       os << tabs3 << "f=" << funcIDList()[f] << " " 
00861                          << dofs << std::endl;
00862                       if (false)
00863                         {
00864                           os << tabs3 << "{";
00865                           for (int i=0; i<dofs.size(); i++)
00866                             {
00867                               if (i != 0) os << ", ";
00868                               if (isLocalDOF(dofs[i])) os << "L";
00869                               else os << "R";
00870                             }
00871                           os << "}" << std::endl;
00872                         }
00873                     }
00874                 }
00875             }
00876         }
00877       mesh().comm().synchronize();
00878       mesh().comm().synchronize();
00879     }
00880 }
00881 
00882
SundanceHomogeneousDOFMap.cpp
