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 "SundanceMap.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceOrderedTuple.hpp"
00035 #include "SundanceNodalDOFMapHN.hpp"
00036 #include "SundanceLagrange.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039 
00040 
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043 using Playa::MPIComm;
00044 using Playa::MPIContainerComm;
00045 
00046 NodalDOFMapHN::NodalDOFMapHN(const Mesh& mesh,
00047   int nFuncs, 
00048   const CellFilter& maxCellFilter,
00049   int setupVerb)
00050   : HNDoFMapBaseHomogeneous(mesh, nFuncs, setupVerb),
00051     maxCellFilter_(maxCellFilter),
00052     dim_(mesh.spatialDim()),
00053     nFuncs_(nFuncs),
00054     nElems_(mesh.numCells(mesh.spatialDim())),
00055     nNodes_(mesh.numCells(0)),
00056     nNodesPerElem_(mesh.numFacets(mesh.spatialDim(),0,0)),
00057     elemDofs_(nElems_ * nFuncs * nNodesPerElem_),
00058     nodeDofs_(mesh.numCells(0)*nFuncs_, -2),
00059     structure_(rcp(new MapStructure(nFuncs_, rcp(new Lagrange(1))))),
00060     hasCellHanging_(nElems_),
00061     nodeIsHanging_(nElems_ * nFuncs * nNodesPerElem_),
00062     cellsWithHangingDoF_globalDoFs_(),
00063     cells_To_NodeLIDs_(),
00064     hangingNodeLID_to_NodesLIDs_(),
00065     hangindNodeLID_to_Coefs_(),
00066     matrixStore_(),
00067     facetLID_()
00068 {
00069   // init the matrix store with one chunck
00070   matrixStore_.init(1);
00071   init();
00072 }
00073 
00074 void NodalDOFMapHN::init()
00075 { 
00076   Tabs tab;
00077 
00078   SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing nodal DOF map nrFunc:"
00079       << nFuncs_ << "  nNodes_:" << nNodes_ << " nElems_:" <<nElems_);
00080 
00081   Array<Array<int> > remoteNodes(mesh().comm().getNProc());
00082   Array<int> hasProcessedCell(nNodes_, 0);
00083 
00084   /* start the DOF count at zero. */
00085   int nextDOF = 0;
00086   int owner;
00087 
00088   nFacets_ = mesh().numFacets(dim_, 0, 0);
00089   Array<int> elemLID(nElems_);
00090   Array<int> facetOrientations(nFacets_*nElems_);
00091 
00092   /* Assign node DOFs for locally owned 0-cells */
00093   CellSet maxCells = maxCellFilter_.getCells(mesh());
00094 
00095   int cc = 0;
00096   for (CellIterator iter=maxCells.begin(); iter != maxCells.end(); iter++, cc++)
00097     {
00098       int c = *iter;
00099       elemLID[cc] = c;
00100     }
00101   /* look up the LIDs of the facets (points)*/
00102   /* This is important so that we have locality in the DOF numbering */
00103   mesh().getFacetLIDs(dim_, elemLID, 0, facetLID_, facetOrientations);
00104   
00105   for (int c=0; c<nElems_; c++) {
00106     hasCellHanging_[c] = false;
00107       /* for each facet, process its DOFs */
00108       for (int f=0; f<nFacets_; f++) {
00109           /* if the facet's DOFs have been assigned already,
00110            * we're done */
00111           int fLID = facetLID_[c*nFacets_+f];
00112           SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() CellLID:"<< c <<"Try point LID:" << fLID << " facet:" << f);
00113           if (hasProcessedCell[fLID] == 0) {
00114               /* the facet may be owned by another processor */
00115             /* Do this only when that node is not HANGING! */
00116               if ( isRemote(0, fLID, owner) && (!mesh().isElementHangingNode(0,fLID)) ) {
00117                   int facetGID 
00118                     = mesh().mapLIDToGID(0, fLID);
00119                   remoteNodes[owner].append(facetGID);
00120                 }
00121               else {/* we can assign a DOF locally */
00122                   SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() Doing point LID:" << fLID << " facet:" << f);
00123                     // test if the node is not a hanging node
00124                     if ( mesh().isElementHangingNode(0,fLID) == false ){
00125                        /* assign DOFs , (for each function space) */
00126                        for (int i=0; i<nFuncs_; i++){
00127                           nodeDofs_[fLID*nFuncs_ + i] = nextDOF;
00128                           nextDOF++;
00129                        }
00130                     } else {
00131                        SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() Hanging node found LID:" << fLID);
00132                        hasCellHanging_[c] = true;
00133                        for (int i=0; i<nFuncs_; i++){
00134                          nodeDofs_[fLID*nFuncs_ + i] = -1; // this means that this is not golbal DoF
00135                        }
00136                        Array<int> pointsLIDs;
00137                        Array<int> facetIndex;
00138                        Array<double> coefs;
00139                        // get the global DoFs which contribute (what if that is not yet numbered?, then only the "fLIDs")
00140                        getPointLIDsForHN( fLID, f, c , pointsLIDs, coefs , facetIndex);
00141 
00142                      // also store the corresponding coefficients
00143                        hangingNodeLID_to_NodesLIDs_.put( fLID , pointsLIDs );
00144                        hangindNodeLID_to_Coefs_.put(fLID,coefs);
00145                     }
00146                 }
00147               hasProcessedCell[fLID] = 1;
00148             }
00149             // if this node is hanging then mark the cell as hanging
00150             if (mesh().isElementHangingNode(0,fLID) == true) {
00151                     hasCellHanging_[c] = true;
00152             }
00153        }
00154   }
00155   
00156   //SUNDANCE_MSG2(setupVerb(), "Before Communication: " << nodeDofs_);
00157 
00158   /* Compute offsets for each processor */
00159   int localCount = nextDOF;
00160   computeOffsets(localCount);
00161   
00162   /* Resolve remote DOF numbers */
00163   shareRemoteDOFs(remoteNodes);
00164 
00165   //SUNDANCE_MSG2(setupVerb(), "After Communication: " << nodeDofs_);
00166 
00167   /* Assign DOFs for elements */
00168   for (int c=0; c<nElems_; c++)
00169     {
00170       if (hasCellHanging_[c]){
00171         Array<int> HNDoFs(nFacets_);
00172         Array<double> transMatrix;
00173         Array<double> coefs;
00174 
00175         transMatrix.resize(nFacets_*nFacets_);
00176         for (int ii = 0 ; ii < nFacets_*nFacets_ ; ii++) transMatrix[ii] = 0.0;
00177 
00178         for (int f=0; f<nFacets_; f++)
00179         {
00180           int fLID = facetLID_[c*nFacets_+f];
00181         SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN cell:" << c << " facetLID:" << fLID
00182               << " hanging:"<< nodeDofs_[fLID*nFuncs_] << "  array:" << HNDoFs);
00183               if (nodeDofs_[fLID*nFuncs_] < 0)
00184               {
00185                   Array<int> pointsLIDs;
00186                   Array<int> facetIndex;
00187                   Array<double> coefs;
00188                   // get the composing points
00189                   getPointLIDsForHN( fLID, f, c , pointsLIDs, coefs , facetIndex);
00190 
00191                   for (int j=0 ; j < pointsLIDs.size() ; j++){
00192                     // look for tmpArray[j] in the existing set, put there coefs[j]
00193                     HNDoFs[facetIndex[j]] = pointsLIDs[j];
00194                     transMatrix[f*nFacets_  + facetIndex[j]] = coefs[j];
00195                   }
00196               }
00197               else
00198               {
00199                 // look for fLID in the actual set and put there 1.0, if not found then size+1
00200                 HNDoFs[f] = fLID;
00201                 transMatrix[f*nFacets_  + f] = 1.0;
00202               }
00203         }
00204         // store the matrix corresponding to this cell
00205         int matrixIndex = matrixStore_.addMatrix(0,transMatrix);
00206         maxCellLIDwithHN_to_TrafoMatrix_.put( c , matrixIndex );
00207 
00208         // store the point LID's which contribute to this cell
00209         cellsWithHangingDoF_globalDoFs_.put( c , HNDoFs );
00210 
00211         SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing cellLID:" << c << " array:" << HNDoFs);
00212         SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing cellLID:" << c << " Trafo array:" << transMatrix);
00213 
00214         // add the global DOFs to the array
00215         for (int f=0; f<nFacets_; f++)
00216         {
00217            int fLID = HNDoFs[f];
00218            for (int i=0; i<nFuncs_; i++) {
00219             elemDofs_[(c*nFuncs_+i)*nFacets_ + f] = nodeDofs_[fLID*nFuncs_ + i];
00220            }
00221 
00222         }
00223         SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN HN initializing cellLID:" << c << " elemDofs_:" << elemDofs_);
00224       }
00225       else {
00226     /* set the element DOFs given the dofs of the facets */
00227         for (int f=0; f<nFacets_; f++)
00228         {
00229           int fLID = facetLID_[c*nFacets_+f];
00230           for (int i=0; i<nFuncs_; i++)
00231           {
00232             elemDofs_[(c*nFuncs_+i)*nFacets_ + f] = nodeDofs_[fLID*nFuncs_ + i];
00233           }
00234         }
00235           SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN initializing cellLID:" << c << " elemDofs_:" << elemDofs_);
00236       }
00237     }
00238   SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN initializing DONE");
00239 
00240 }
00241 
00242 RCP<const MapStructure>
00243 NodalDOFMapHN::getDOFsForCellBatch(int cellDim,
00244   const Array<int>& cellLID,
00245   const Set<int>& requestedFuncSet,
00246   Array<Array<int> >& dofs,
00247   Array<int>& nNodes,
00248   int verbosity) const
00249 {
00250   TimeMonitor timer(batchedDofLookupTimer());
00251 
00252   Tabs tab;
00253   SUNDANCE_MSG2(verbosity,
00254                tab << "NodalDOFMapHN::getDOFsForCellBatch(): cellDim=" << cellDim
00255                << " requestedFuncSet:" << requestedFuncSet
00256                << " cellLID=" << cellLID);
00257 
00258 
00259   dofs.resize(1);
00260   nNodes.resize(1);
00261 
00262   int nCells = cellLID.size();
00263 
00264 
00265   if (cellDim == dim_)
00266     {
00267       int dofsPerElem = nFuncs_ * nNodesPerElem_;
00268       nNodes[0] = nNodesPerElem_;
00269       dofs[0].resize(nCells * dofsPerElem);
00270       Array<int>& dof0 = dofs[0];
00271       Array<int> tmpArray;
00272 
00273     tmpArray.resize(dofsPerElem);
00274 
00275       for (int c=0; c<nCells; c++)
00276         {
00277           // here , nothing to do ... in elemDofs_ should be the proper DoFs
00278             for (int i=0; i<dofsPerElem; i++) {
00279                dof0[c*dofsPerElem + i] = elemDofs_[cellLID[c]*dofsPerElem+i];
00280                tmpArray[i] = elemDofs_[cellLID[c]*dofsPerElem+i];
00281             }
00282             SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00283                 cellDim << " cellLID:" << c << " array:" << tmpArray);
00284         }
00285     }
00286   else if (cellDim == 0)
00287     { // point integrals are also mostly used for BCs
00288       nNodes[0] = 1;
00289       dofs[0].resize(nCells * nFuncs_);
00290       Array<int>& dof0 = dofs[0];
00291       Array<int> tmpArray;
00292     tmpArray.resize(nFuncs_);
00293 
00294       for (int c=0; c<nCells; c++)
00295         {
00296           for (int i=0; i<nFuncs_; i++)
00297             {
00298               dof0[c*nFuncs_ + i] = nodeDofs_[cellLID[c]*nFuncs_+i];
00299               tmpArray[i] = nodeDofs_[cellLID[c]*nFuncs_+i];
00300             }
00301           SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00302               cellDim << " pointLID:" << cellLID[c] << " array:" << tmpArray);
00303         }
00304     }
00305   else
00306     {
00307     // since edge or surface Integrals are mostly used for BCs, where are NO HN
00308       int nFacets = mesh().numFacets(cellDim, cellLID[0], 0);
00309       nNodes[0] = nFacets;
00310       int dofsPerCell = nFuncs_ * nNodes[0];
00311       dofs[0].resize(nCells * dofsPerCell);
00312       Array<int>& dof0 = dofs[0];
00313       Array<int> facetLIDs(nCells * nNodes[0]);
00314       Array<int> dummy(nCells * nNodes[0]);
00315       Array<int> tmpArray;
00316     tmpArray.resize(nFuncs_*nFacets);
00317 
00318       mesh().getFacetLIDs(cellDim, cellLID, 0, facetLIDs, dummy);
00319 
00320       for (int c=0; c<nCells; c++)
00321         {
00322           for (int f=0; f<nFacets; f++)
00323             {
00324               int facetCellLID = facetLIDs[c*nFacets+f];
00325               for (int i=0; i<nFuncs_; i++)
00326                 {
00327                   dof0[(c*nFuncs_+i)*nFacets + f]
00328                     = nodeDofs_[facetCellLID*nFuncs_+i];
00329                   // when we want to pass a hanging DoF , we correct the result
00330                   if ( nodeDofs_[facetCellLID*nFuncs_+i] < 0){
00331                     int tmp = 1;
00332                     // get any cell which has this element
00333                     int maxCell = mesh().maxCofacetLID( cellDim , cellLID[c] , 0 , tmp);
00334                       Array<int> facetsLIDs;
00335                       mesh().returnParentFacets( maxCell , 0 , facetsLIDs , tmp );
00336                       // get in the same way the point children from the child cell
00337                       Array<int> childFacets(mesh().numFacets(mesh().spatialDim(),0,0));
00338                       for (int jk = 0 ; jk < childFacets.size() ; jk++)
00339                         childFacets[jk] = mesh().facetLID( mesh().spatialDim() , maxCell, 0 , jk , tmp);
00340                       int fIndex = 0;
00341                       // get the correct point facet index, by comparing the child cell point facets to the HN LID
00342                       for (int jk = 0 ; jk < childFacets.size() ; jk++)
00343                         if ( childFacets[jk] == facetCellLID) { fIndex = jk; break;}
00344                       dof0[(c*nFuncs_+i)*nFacets + f] = nodeDofs_[facetsLIDs[fIndex]*nFuncs_ + i];
00345                   }
00346                   tmpArray[f*nFuncs_+i] = dof0[(c*nFuncs_+i)*nFacets + f];
00347                 }
00348             }
00349           SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00350               cellDim << " edge or face LID:" << cellLID[c] << " array:" << tmpArray);
00351         }
00352     }
00353   SUNDANCE_MSG2(verbosity,
00354                 tab << "NodalDOFMapHN::getDOFsForCellBatch(): DONE");
00355   return structure_;
00356 }
00357 
00358 
00359 void NodalDOFMapHN::getTrafoMatrixForCell(
00360       int cellLID,
00361       int funcID,
00362       int& trafoMatrixSize,
00363       bool& doTransform,
00364       Array<double>& transfMatrix ) const {
00365 
00366   trafoMatrixSize = nFacets_;
00367 
00368   if (cellsWithHangingDoF_globalDoFs_.containsKey(cellLID))
00369   {
00370     // return the transformation matrix from the Map
00371     int matrixIndex = maxCellLIDwithHN_to_TrafoMatrix_.get( cellLID ); // this should return a valid array
00372     matrixStore_.getMatrix( 0 , matrixIndex , transfMatrix );
00373     doTransform = true;
00374   }
00375   else  // no transformation needed, return false
00376   {
00377     doTransform = false;
00378   }
00379 }
00380 
00381 void NodalDOFMapHN::getTrafoMatrixForFacet(
00382       int cellDim,
00383       int cellLID,
00384       int facetIndex,
00385       int funcID,
00386       int& trafoMatrixSize,
00387       bool& doTransform,
00388       Array<double>& transfMatrix ) const {
00389 
00390   int fIndex;
00391   int maxCellLID;
00392   // here we ask for cofacet 0 , assuming that these are anyhow boundary facets
00393   SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() cellDim :" << cellDim << ", cellLID:" << cellLID);
00394   maxCellLID = mesh().maxCofacetLID( cellDim, cellLID, 0 , fIndex);
00395   SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() testing :" << maxCellLID);
00396 
00397   // todo: we might pre filter cases when this is not necessary
00398 
00399   if (cellsWithHangingDoF_globalDoFs_.containsKey(maxCellLID))
00400   {
00401     int matrixIndex    =    maxCellLIDwithHN_to_TrafoMatrix_.get( maxCellLID );
00402     matrixStore_.getMatrix( 0 , matrixIndex , transfMatrix );
00403     doTransform = true;
00404   }
00405   else  // no transformation needed, return false
00406   {
00407     doTransform = false;
00408   }
00409 }
00410 
00411 /** Function for nodal plotting */
00412 void NodalDOFMapHN::getDOFsForHNCell(
00413     int cellDim,
00414     int cellLID,
00415       int funcID,
00416       Array<int>& dofs ,
00417       Array<double>& coefs ) const {
00418 
00419   // treat the cells only of dimension zero
00420   if (  (cellDim == 0) && (hangingNodeLID_to_NodesLIDs_.containsKey(cellLID))  )
00421   {
00422     Array<int> pointLIDs;
00423     Array<int> pointCoefs;
00424     pointLIDs = hangingNodeLID_to_NodesLIDs_.get( cellLID );
00425     dofs.resize(pointLIDs.size());
00426     // return the DoFs belonging to the points and function which collaborate to the hanging node
00427     for (int ind = 0 ; ind < pointLIDs.size() ; ind++){
00428       dofs[ind] = nodeDofs_[ pointLIDs[ind] *nFuncs_ + funcID ]; // return the DoF corresp. to function ID
00429     }
00430     // get the coefficients
00431     coefs = hangindNodeLID_to_Coefs_.get( cellLID );
00432   }
00433 }
00434 
00435 /** This function is only for TRISECTION implemented */
00436 // we might use the method of BasisFunction->getConstraintForHN(), but this should be faster
00437 void NodalDOFMapHN::getPointLIDsForHN( int pointLID , int facetIndex ,
00438     int maxCellIndex ,Array<int>& glbLIDs, Array<double>& coefsArray , Array<int>& nodeIndex){
00439 
00440   Array<int> facets;
00441     int indexInParent , parentLID , facetCase = 0;
00442     double divRatio = 0.5;
00443   switch (dim_){
00444   case 2:
00445     // todo: eventually implement bisection as well (if it will be needed)
00446     glbLIDs.resize(2);
00447     nodeIndex.resize(2);
00448     indexInParent = mesh().indexInParent(maxCellIndex);
00449     switch (indexInParent){
00450       case 0: {facetCase = (facetIndex == 1)? 0 : 1;
00451               divRatio = (1.0/3.0);  break;}
00452       case 1: {facetCase = 0;
00453               divRatio = (facetIndex == 0)? (1.0/3.0) : (2.0/3.0);  break;}
00454       case 2: {facetCase = (facetIndex == 0)? 0 : 2;
00455                   divRatio = (facetIndex == 0)? (2.0/3.0) : (1.0/3.0);  break;}
00456       case 5: {facetCase = 1;
00457                   divRatio = (facetIndex == 0)? (1.0/3.0) : (2.0/3.0);  break;}
00458       case 3: {facetCase = 2;
00459                   divRatio = (facetIndex == 1)? (1.0/3.0) : (2.0/3.0);  break;}
00460       case 6: {facetCase = (facetIndex == 0)? 1 : 3;
00461                   divRatio = (facetIndex == 0)? (2.0/3.0) : (1.0/3.0);  break;}
00462       case 7: {facetCase = 3;
00463                   divRatio = (facetIndex == 2)? (1.0/3.0) : (2.0/3.0);  break;}
00464       case 8: {facetCase = (facetIndex == 1)? 2 : 3;
00465                   divRatio = (2.0/3.0);  break;}
00466     }
00467     mesh().returnParentFacets(maxCellIndex , 0 ,facets , parentLID );
00468     switch (facetCase){
00469        case 0: {glbLIDs[0] = facets[0]; glbLIDs[1] = facets[1];
00470                 nodeIndex[0] = 0;  nodeIndex[1] = 1;  break;}
00471        case 1: {glbLIDs[0] = facets[0]; glbLIDs[1] = facets[2];
00472                 nodeIndex[0] = 0;  nodeIndex[1] = 2;  break;}
00473        case 2: {glbLIDs[0] = facets[1]; glbLIDs[1] = facets[3];
00474                 nodeIndex[0] = 1;  nodeIndex[1] = 3;   break;}
00475        case 3: {glbLIDs[0] = facets[2]; glbLIDs[1] = facets[3];
00476                 nodeIndex[0] = 2;  nodeIndex[1] = 3;  break;}
00477     }
00478      coefsArray.resize(2); coefsArray[0] = 1 - divRatio; coefsArray[1] = divRatio;
00479      SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::getPointLIDsForHN() fc=" << facetCase << " R=" << divRatio
00480                    << " facetIndex:" << facetIndex <<   " indexInParent:" << indexInParent <<" glbLIDs:"
00481                    << glbLIDs << " maxCellIndex:" << maxCellIndex << " nodeIndex:" << nodeIndex);
00482     break;
00483   case 3:{
00484         // 3D hanging node
00485     double ind_x = 0 , ind_y = 0 , ind_z = 0;
00486     double f_x = 0.0 , f_y = 0.0 , f_z = 0.0;
00487     double xx = 0.0  , yy = 0.0  , zz = 0.0;
00488     double values[8];
00489     indexInParent = mesh().indexInParent(maxCellIndex);
00490     mesh().returnParentFacets(maxCellIndex , 0 ,facets , parentLID );
00491     //
00492     switch (indexInParent){
00493       case 0:  { ind_x = 0.0; ind_y = 0.0; ind_z = 0.0; break;}
00494       case 1:  { ind_x = 1.0; ind_y = 0.0; ind_z = 0.0; break;}
00495       case 2:  { ind_x = 2.0; ind_y = 0.0; ind_z = 0.0; break;}
00496       case 3:  { ind_x = 0.0; ind_y = 1.0; ind_z = 0.0; break;}
00497       case 4:  { ind_x = 1.0; ind_y = 1.0; ind_z = 0.0; break;}
00498       case 5:  { ind_x = 2.0; ind_y = 1.0; ind_z = 0.0; break;}
00499       case 6:  { ind_x = 0.0; ind_y = 2.0; ind_z = 0.0; break;}
00500       case 7:  { ind_x = 1.0; ind_y = 2.0; ind_z = 0.0; break;}
00501       case 8:  { ind_x = 2.0; ind_y = 2.0; ind_z = 0.0; break;}
00502       case 9:  { ind_x = 0.0; ind_y = 0.0; ind_z = 1.0; break;}
00503       case 10: { ind_x = 1.0; ind_y = 0.0; ind_z = 1.0; break;}
00504       case 11: { ind_x = 2.0; ind_y = 0.0; ind_z = 1.0; break;}
00505       case 12: { ind_x = 0.0; ind_y = 1.0; ind_z = 1.0; break;}
00506       case 14: { ind_x = 2.0; ind_y = 1.0; ind_z = 1.0; break;}
00507       case 15: { ind_x = 0.0; ind_y = 2.0; ind_z = 1.0; break;}
00508       case 16: { ind_x = 1.0; ind_y = 2.0; ind_z = 1.0; break;}
00509       case 17: { ind_x = 2.0; ind_y = 2.0; ind_z = 1.0; break;}
00510       case 18: { ind_x = 0.0; ind_y = 0.0; ind_z = 2.0; break;}
00511       case 19: { ind_x = 1.0; ind_y = 0.0; ind_z = 2.0; break;}
00512       case 20: { ind_x = 2.0; ind_y = 0.0; ind_z = 2.0; break;}
00513       case 21: { ind_x = 0.0; ind_y = 1.0; ind_z = 2.0; break;}
00514       case 22: { ind_x = 1.0; ind_y = 1.0; ind_z = 2.0; break;}
00515       case 23: { ind_x = 2.0; ind_y = 1.0; ind_z = 2.0; break;}
00516       case 24: { ind_x = 0.0; ind_y = 2.0; ind_z = 2.0; break;}
00517       case 25: { ind_x = 1.0; ind_y = 2.0; ind_z = 2.0; break;}
00518       case 26: { ind_x = 2.0; ind_y = 2.0; ind_z = 2.0; break;}
00519       default: {}// error this should not occur
00520     }
00521     switch (facetIndex){
00522       case 0:  { f_x = 0.0 , f_y = 0.0 , f_z = 0.0; break;}
00523       case 1:  { f_x = 1.0 , f_y = 0.0 , f_z = 0.0; break;}
00524       case 2:  { f_x = 0.0 , f_y = 1.0 , f_z = 0.0; break;}
00525       case 3:  { f_x = 1.0 , f_y = 1.0 , f_z = 0.0; break;}
00526       case 4:  { f_x = 0.0 , f_y = 0.0 , f_z = 1.0; break;}
00527       case 5:  { f_x = 1.0 , f_y = 0.0 , f_z = 1.0; break;}
00528       case 6:  { f_x = 0.0 , f_y = 1.0 , f_z = 1.0; break;}
00529       case 7:  { f_x = 1.0 , f_y = 1.0 , f_z = 1.0; break;}
00530       default: {}// error this should not occur
00531     }
00532     // evaluate the bilinear basis function at the given point
00533       xx = (ind_x + f_x)/3.0;
00534       yy = (ind_y + f_y)/3.0;
00535       zz = (ind_z + f_z)/3.0;
00536       values[0] = (1.0 - xx)*(1.0 - yy)*(1.0 - zz);
00537       values[1] = (xx)*(1.0 - yy)*(1.0 - zz);
00538       values[2] = (1.0 - xx)*(yy)*(1.0 - zz);
00539       values[3] = (xx)*(yy)*(1.0 - zz);
00540       values[4] = (1.0 - xx)*(1.0 - yy)*(zz);
00541       values[5] = (xx)*(1.0 - yy)*(zz);
00542       values[6] = (1.0 - xx)*(yy)*(zz);
00543       values[7] = (xx)*(yy)*(zz);
00544     // resize the array to zero and add
00545       glbLIDs.resize(0);
00546       nodeIndex.resize(0);
00547       coefsArray.resize(0);
00548       for (int ii = 0 ; ii < 8 ; ii++ ){
00549         // add the facet point if the basis is not zero
00550         if (fabs(values[ii]) > 1e-5){
00551           glbLIDs.append(facets[ii]);
00552           nodeIndex.append(ii);
00553           coefsArray.append(values[ii]);
00554         }
00555       }
00556       SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::getPointLIDsForHN()" << " facetIndex:" << facetIndex << " indexInParent:"
00557          << indexInParent <<" glbLIDs:" << glbLIDs << " maxCellIndex:" << maxCellIndex << " nodeIndex:" << nodeIndex);
00558     break;
00559       }
00560   }
00561 
00562 }
00563 
00564 
00565 
00566 void NodalDOFMapHN::computeOffsets(int localCount)
00567 {
00568   Array<int> dofOffsets;
00569   int totalDOFCount;
00570   int myOffset = 0;
00571   int np = mesh().comm().getNProc();
00572   if (np > 1)
00573     {
00574     SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::computeOffsets, localCount:" << localCount);
00575       MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00576                                         mesh().comm());
00577       myOffset = dofOffsets[mesh().comm().getRank()];
00578 
00579       SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::computeOffsets, offset:" << myOffset);
00580       int nDofs = nNodes_ * nFuncs_;
00581       for (int i=0; i<nDofs; i++)
00582         {
00583           if (nodeDofs_[i] >= 0) nodeDofs_[i] += myOffset;
00584         }
00585     }
00586   else
00587     {
00588       totalDOFCount = localCount;
00589     }
00590   
00591   setLowestLocalDOF(myOffset);
00592   setNumLocalDOFs(localCount);
00593   setTotalNumDOFs(totalDOFCount);
00594 }
00595 
00596 
00597 void NodalDOFMapHN::shareRemoteDOFs(const Array<Array<int> >& outgoingCellRequests)
00598 {
00599   int np = mesh().comm().getNProc();
00600   if (np==1) return;
00601   int rank = mesh().comm().getRank();
00602 
00603   Array<Array<int> > incomingCellRequests;
00604   Array<Array<int> > outgoingDOFs(np);
00605   Array<Array<int> > incomingDOFs;
00606 
00607   SUNDANCE_MSG2(setupVerb(),
00608                "p=" << mesh().comm().getRank()
00609                << "synchronizing DOFs for cells of dimension 0");
00610   SUNDANCE_MSG2(setupVerb(),
00611                "p=" << mesh().comm().getRank()
00612                << " sending cell reqs d=0, GID=" 
00613                << outgoingCellRequests);
00614 
00615   /* share the cell requests */
00616   MPIContainerComm<int>::allToAll(outgoingCellRequests, 
00617                                   incomingCellRequests,
00618                                   mesh().comm());
00619   
00620   /* get DOF numbers for the zeroth function index on every node that's been 
00621    * requested by someone else */
00622   for (int p=0; p<np; p++)
00623     {
00624       if (p==rank) continue;
00625       const Array<int>& requestsFromProc = incomingCellRequests[p];
00626       int nReq = requestsFromProc.size();
00627 
00628       SUNDANCE_MSG3(setupVerb(),"p=" << mesh().comm().getRank()
00629                             << " recv'd from proc=" << p
00630                             << " reqs for DOFs for cells " 
00631                             << requestsFromProc);
00632 
00633       outgoingDOFs[p].resize(nReq);
00634 
00635       for (int c=0; c<nReq; c++)
00636         {
00637           int GID = requestsFromProc[c];
00638           SUNDANCE_MSG2(setupVerb(),
00639                        "p=" << rank
00640                        << " processing zero-cell with GID=" << GID); 
00641           int LID = mesh().mapGIDToLID(0, GID);
00642           SUNDANCE_MSG2(setupVerb(),
00643                        "p=" << rank
00644                        << " LID=" << LID << " dofs=" 
00645                        << nodeDofs_[LID*nFuncs_]);
00646           outgoingDOFs[p][c] = nodeDofs_[LID*nFuncs_];
00647           SUNDANCE_MSG2(setupVerb(),
00648                        "p=" << rank
00649                        << " done processing cell with GID=" << GID);
00650         }
00651     }
00652 
00653   SUNDANCE_MSG2(setupVerb(),
00654                "p=" << mesh().comm().getRank()
00655                << "answering DOF requests for cells of dimension 0");
00656 
00657   /* share the DOF numbers */
00658   MPIContainerComm<int>::allToAll(outgoingDOFs,
00659                                   incomingDOFs,
00660                                   mesh().comm());
00661 
00662   SUNDANCE_MSG2(setupVerb(),
00663                "p=" << mesh().comm().getRank()
00664                << "communicated DOF answers for cells of dimension 0" );
00665 
00666   
00667   /* now assign the DOFs from the other procs */
00668 
00669   for (int p=0; p<mesh().comm().getNProc(); p++)
00670     {
00671       if (p==mesh().comm().getRank()) continue;
00672       const Array<int>& dofsFromProc = incomingDOFs[p];
00673       int numCells = dofsFromProc.size();
00674       for (int c=0; c<numCells; c++)
00675         {
00676           int cellGID = outgoingCellRequests[p][c];
00677           int cellLID = mesh().mapGIDToLID(0, cellGID);
00678           int dof = dofsFromProc[c];
00679           for (int i=0; i<nFuncs_; i++)
00680             {
00681               nodeDofs_[cellLID*nFuncs_ + i] = dof+i;
00682               addGhostIndex(dof+i);
00683             }
00684         }
00685     }
00686 }
SundanceNodalDOFMapHN.cpp
