FEI Version of the Day
SNL_FEI_Structure.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include "fei_sstream.hpp"
00045 #include "fei_fstream.hpp"
00046 
00047 #include <limits>
00048 #include <cmath>
00049 #include <assert.h>
00050 
00051 #include "fei_defs.h"
00052 
00053 #include "fei_TemplateUtils.hpp"
00054 #include <fei_CommUtils.hpp>
00055 #include "snl_fei_Constraint.hpp"
00056 typedef snl_fei::Constraint<GlobalID> ConstraintType;
00057 
00058 #include "fei_NodeDescriptor.hpp"
00059 #include "fei_NodeCommMgr.hpp"
00060 #include "fei_NodeDatabase.hpp"
00061 
00062 #include "fei_SlaveVariable.hpp"
00063 
00064 #include "fei_BlockDescriptor.hpp"
00065 
00066 #include "snl_fei_PointBlockMap.hpp"
00067 #include "fei_ProcEqns.hpp"
00068 #include "fei_EqnBuffer.hpp"
00069 #include <fei_FillableMat.hpp>
00070 #include <fei_FillableVec.hpp>
00071 #include <fei_CSRMat.hpp>
00072 #include <fei_CSVec.hpp>
00073 #include "fei_EqnCommMgr.hpp"
00074 
00075 #include "fei_Lookup.hpp"
00076 #include "fei_ConnectivityTable.hpp"
00077 #include "snl_fei_Utils.hpp"
00078 #include "SNL_FEI_Structure.hpp"
00079 
00080 #undef fei_file
00081 #define fei_file "SNL_FEI_Structure.cpp"
00082 #include "fei_ErrMacros.hpp"
00083 
00084 //-----Constructor-------------------------------------------------------------
00085 SNL_FEI_Structure::SNL_FEI_Structure(MPI_Comm comm)
00086  : comm_(comm),
00087    localProc_(0),
00088    masterProc_(0),
00089    numProcs_(1),
00090    fieldIDs_(),
00091    fieldSizes_(),
00092    fieldDatabase_(new std::map<int,int>),
00093    fieldDofMap_(),
00094    workarray_(),
00095    blockIDs_(),
00096    blocks_(),
00097    connTables_(),
00098    nodeDatabase_(NULL),
00099    activeNodesInitialized_(false),
00100    globalNodeOffsets_(),
00101    globalEqnOffsets_(),
00102    globalBlkEqnOffsets_(),
00103    slaveVars_(NULL),
00104    slaveEqns_(NULL),
00105    slvEqnNumbers_(NULL),
00106    numSlvs_(0),
00107    lowestSlv_(0),
00108    highestSlv_(0),
00109    slaveMatrix_(NULL),
00110    globalNumNodesVanished_(),
00111    localVanishedNodeNumbers_(),
00112    nodeCommMgr_(NULL),
00113    eqnCommMgr_(NULL),
00114    slvCommMgr_(NULL),
00115    numGlobalEqns_(0),
00116    numLocalEqns_(0),
00117    localStartRow_(0),
00118    localEndRow_(0),
00119    numLocalNodalEqns_(0),
00120    numLocalElemDOF_(0),
00121    numLocalMultCRs_(0),
00122    reducedStartRow_(0),
00123    reducedEndRow_(0),
00124    numLocalReducedRows_(0),
00125    Kid_(NULL),
00126    Kdi_(NULL),
00127    Kdd_(NULL),
00128    tmpMat1_(),
00129    tmpMat2_(),
00130    reducedEqnCounter_(0),
00131    reducedRHSCounter_(0),
00132    rSlave_(),
00133    cSlave_(),
00134    work_nodePtrs_(),
00135    structureFinalized_(false),
00136    generateGraph_(true),
00137    sysMatIndices_(NULL),
00138    blockMatrix_(false),
00139    numGlobalEqnBlks_(0),
00140    numLocalEqnBlks_(0),
00141    numLocalReducedEqnBlks_(0),
00142    localBlkOffset_(0),
00143    localReducedBlkOffset_(0),
00144    globalMaxBlkSize_(0),
00145    firstLocalNodeNumber_(-1),
00146    numGlobalNodes_(0),
00147    sysBlkMatIndices_(NULL),
00148    matIndicesDestroyed_(false),
00149    workSpace_(),
00150    blkEqnMapper_(new snl_fei::PointBlockMap()),
00151    multCRs_(),
00152    penCRs_(),
00153    checkSharedNodes_(false),
00154    name_(),
00155    outputLevel_(0),
00156    debugOutput_(false),
00157    dbgPath_(),
00158    dbgOStreamPtr_(NULL),
00159    setDbgOutCalled_(true)
00160 {
00161   numProcs_ = 1, localProc_ = 0, masterProc_ = 0;
00162 
00163   localProc_ = fei::localProc(comm_);
00164   numProcs_ = fei::numProcs(comm_);
00165   masterProc_ = 0;
00166 
00167   slaveVars_ = new std::vector<SlaveVariable*>;
00168   slaveEqns_ = new EqnBuffer;
00169 
00170   nodeCommMgr_ = new NodeCommMgr(comm_, *this);
00171 
00172   eqnCommMgr_ = new EqnCommMgr(comm_);
00173   eqnCommMgr_->setNumRHSs(1);
00174 
00175   nodeDatabase_ = new NodeDatabase(fieldDatabase_, nodeCommMgr_);
00176 
00177   Kid_ = new fei::FillableMat;
00178   Kdi_ = new fei::FillableMat;
00179   Kdd_ = new fei::FillableMat;
00180 }
00181 
00182 //-----------------------------------------------------------------------------
00183 int SNL_FEI_Structure::parameters(int numParams, const char*const* paramStrings)
00184 {
00185   const char* param = NULL;
00186 
00187   param = snl_fei::getParamValue("outputLevel",numParams,paramStrings);
00188   if (param != NULL){
00189     std::string str(param);
00190     FEI_ISTRINGSTREAM isstr(str);
00191     isstr >> outputLevel_;
00192   }
00193 
00194   param = snl_fei::getParam("debugOutput",numParams,paramStrings);
00195   if (param != NULL){
00196     debugOutput_ = true;
00197   }
00198 
00199   param = snl_fei::getParam("debugOutputOff",numParams,paramStrings);
00200   if (param != NULL){
00201     debugOutput_ = false;
00202   }
00203 
00204   param = snl_fei::getParam("FEI_CHECK_SHARED_IDS", numParams,paramStrings);
00205   if (param != NULL){
00206     checkSharedNodes_ = true;
00207   }
00208 
00209   param = snl_fei::getParamValue("sharedNodeOwnership",
00210           numParams,paramStrings);
00211   if (param != NULL){
00212     if (!strcmp(param, "LowNumberedProc")) {
00213       nodeCommMgr_->setSharedOwnershipRule(NodeCommMgr::STRICTLY_LOW_PROC);
00214     }
00215     if (!strcmp(param, "ProcWithLocalElem")) {
00216       nodeCommMgr_->setSharedOwnershipRule(NodeCommMgr::PROC_WITH_LOCAL_ELEM);
00217     }
00218   }
00219 
00220   return(0);
00221 }
00222 
00223 //-----Destructor--------------------------------------------------------------
00224 SNL_FEI_Structure::~SNL_FEI_Structure()
00225 {
00226   for(size_t j=0; j<slaveVars_->size(); j++) {
00227     delete (*slaveVars_)[j];
00228   }
00229   delete slaveVars_;
00230 
00231   delete slaveEqns_;
00232   delete slaveMatrix_;
00233 
00234   destroyBlockRoster();
00235   destroyConnectivityTables();
00236 
00237   delete nodeCommMgr_;
00238   delete eqnCommMgr_;
00239   delete blkEqnMapper_;
00240 
00241   destroyMatIndices();
00242 
00243   deleteMultCRs();
00244 
00245   int numPCRs = penCRs_.size();
00246   if (numPCRs > 0) {
00247     fei::destroyValues(penCRs_);
00248     penCRs_.clear();
00249   }
00250 
00251   delete nodeDatabase_;
00252   delete fieldDatabase_;
00253 
00254   delete Kid_;
00255   delete Kdi_;
00256   delete Kdd_;
00257 }
00258 
00259 //------------------------------------------------------------------------------
00260 int SNL_FEI_Structure::setDbgOut(FEI_OSTREAM& ostr,
00261          const char* path, const char* feiName)
00262 {
00263   dbgOStreamPtr_ = &ostr;
00264   setDbgOutCalled_ = true;
00265   if (path != NULL) {
00266     dbgPath_ = path;
00267   }
00268   else {
00269     dbgPath_ = ".";
00270   }
00271   if (feiName != NULL) {
00272     name_ = feiName;
00273   }
00274 
00275   debugOutput_ = true;
00276   return(0);
00277 }
00278 
00279 //------------------------------------------------------------------------------
00280 void SNL_FEI_Structure::destroyBlockRoster()
00281 {
00282   for(size_t i=0; i<blockIDs_.size(); i++) delete blocks_[i];
00283   blocks_.resize(0);
00284 }
00285 
00286 //------------------------------------------------------------------------------
00287 void SNL_FEI_Structure::destroyConnectivityTables()
00288 {
00289   for(size_t i=0; i<blockIDs_.size(); i++) {
00290     delete connTables_[i];
00291   }
00292   connTables_.resize(0);
00293 }
00294 
00295 //------------------------------------------------------------------------------
00296 void SNL_FEI_Structure::destroyMatIndices()
00297 {
00298   if (matIndicesDestroyed_ == true) return;
00299 
00300   delete [] sysMatIndices_;
00301   sysMatIndices_ = NULL;
00302 
00303   delete [] sysBlkMatIndices_;
00304   sysBlkMatIndices_ = NULL;
00305 
00306   matIndicesDestroyed_ = true;
00307 }
00308 
00309 //------------------------------------------------------------------------------
00310 const int* SNL_FEI_Structure::getNumFieldsPerNode(GlobalID blockID)
00311 {
00312   BlockDescriptor* block = NULL;
00313   int err = getBlockDescriptor(blockID, block);
00314   if (err) return(NULL);
00315 
00316   return(block->fieldsPerNodePtr());
00317 }
00318 
00319 //------------------------------------------------------------------------------
00320 const int* const* SNL_FEI_Structure::getFieldIDsTable(GlobalID blockID)
00321 {
00322   BlockDescriptor* block = NULL;
00323   int err = getBlockDescriptor(blockID, block);
00324   if (err) return(NULL);
00325 
00326   return(block->fieldIDsTablePtr());
00327 }
00328 
00329 //------------------------------------------------------------------------------
00330 void SNL_FEI_Structure::getElemBlockInfo(GlobalID blockID,
00331                          int& interleaveStrategy, int& lumpingStrategy,
00332                          int& numElemDOF, int& numElements,
00333                          int& numNodesPerElem, int& numEqnsPerElem)
00334 {
00335   BlockDescriptor* block = NULL;
00336   int err = getBlockDescriptor(blockID, block);
00337   if (err) {
00338       interleaveStrategy = -1;
00339       lumpingStrategy = -1;
00340       numElemDOF = -1;
00341       numElements = -1;
00342       numNodesPerElem = -1;
00343       numEqnsPerElem = -1;
00344    }
00345 
00346   interleaveStrategy = block->getInterleaveStrategy();
00347   lumpingStrategy = block->getLumpingStrategy();
00348   numElemDOF = block->getNumElemDOFPerElement();
00349   numElements = block->getNumElements();
00350   numNodesPerElem = block->getNumNodesPerElement();
00351   numEqnsPerElem = block->getNumEqnsPerElement();
00352 }
00353 
00354 //------------------------------------------------------------------------------
00355 int SNL_FEI_Structure::getEqnNumber(int nodeNumber, int fieldID)
00356 {
00357   //this function is only used by clients of the Lookup interface, who are
00358   //expecting equations in the 'reduced' equation space.
00359 
00360   int eqnNumber;
00361 
00362   const NodeDescriptor* node = NULL;
00363   CHK_ERR( nodeDatabase_->getNodeWithNumber(nodeNumber, node) );
00364 
00365   bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
00366   if (!hasField) {
00367     fei::console_out() << "SNL_FEI_Structure::getEqnNumber: ERROR, node with nodeNumber "
00368    << nodeNumber << " doesn't have fieldID " << fieldID << FEI_ENDL;
00369     ERReturn(-1);
00370   }
00371 
00372   int reducedEqn = -1;
00373   bool isSlave = translateToReducedEqn(eqnNumber, reducedEqn);
00374   if (isSlave) return(-1);
00375 
00376   return(reducedEqn);
00377 }
00378 
00379 //------------------------------------------------------------------------------
00380 int SNL_FEI_Structure::getOwnerProcForEqn(int eqn)
00381 {
00382   if (eqn < 0) return(-1);
00383 
00384   int len = globalEqnOffsets_.size(); // len is numProcs+1...
00385 
00386   for(int i=0; i<len-1; i++) {
00387     if (eqn >= globalEqnOffsets_[i] && eqn < globalEqnOffsets_[i+1]) return(i);
00388   }
00389 
00390   return(-1);
00391 }
00392 
00393 //------------------------------------------------------------------------------
00394 int SNL_FEI_Structure::initFields(int numFields, 
00395          const int *fieldSizes, 
00396          const int *fieldIDs,
00397          const int *fieldTypes)
00398 {
00399   // store the incoming solution fields
00400   //
00401   if (debugOutput_) {
00402     FEI_OSTREAM& ostr = dbgOut();
00403     ostr << "FEI: initFields" << FEI_ENDL
00404    << "#num-fields" << FEI_ENDL << numFields << FEI_ENDL;
00405     int nf;
00406     ostr << "#field-sizes" << FEI_ENDL;
00407     for(nf=0; nf<numFields; nf++) {
00408       ostr <<fieldSizes[nf] << " ";
00409     }
00410     ostr << FEI_ENDL << "#field-ids" << FEI_ENDL;
00411     for(nf=0; nf<numFields; nf++) {
00412       ostr << fieldIDs[nf] << " ";
00413     }
00414     if (fieldTypes != NULL) {
00415       ostr << FEI_ENDL << "#field-types" << FEI_ENDL;
00416       for(nf=0; nf<numFields; nf++) {
00417         ostr << fieldTypes[nf] << " ";
00418       }
00419     }
00420     ostr<<FEI_ENDL;
00421   }
00422 
00423   for (int i=0; i<numFields; i++) {
00424     fieldDatabase_->insert(std::pair<int,int>(fieldIDs[i], fieldSizes[i]));
00425     int offs = fei::sortedListInsert(fieldIDs[i], fieldIDs_); 
00426     if (offs >= 0) fieldSizes_.insert(fieldSizes_.begin()+offs, fieldSizes[i]);
00427 
00428     if (fieldIDs[i] >= 0) {
00429       if (fieldTypes != NULL) {
00430         fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i], fieldTypes[i]);
00431       }
00432       else {
00433         fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i]);
00434       }
00435     }
00436   }
00437 
00438   return(FEI_SUCCESS);
00439 }
00440  
00441 //------------------------------------------------------------------------------
00442 int SNL_FEI_Structure::initElemBlock(GlobalID elemBlockID,
00443             int numElements,
00444             int numNodesPerElement,
00445             const int* numFieldsPerNode,
00446             const int* const* nodalFieldIDs,
00447             int numElemDofFieldsPerElement,
00448             const int* elemDofFieldIDs,
00449             int interleaveStrategy)
00450 {
00451   int nn, nf;
00452   if (debugOutput_) {
00453     FEI_OSTREAM& os = dbgOut();
00454     os << "FEI: initElemBlock" << FEI_ENDL << "#elemBlockID" << FEI_ENDL
00455        << (int)elemBlockID << FEI_ENDL;
00456     os << "#numElements"<< FEI_ENDL << numElements << FEI_ENDL;
00457     os << "#numNodesPerElement"<< FEI_ENDL <<numNodesPerElement << FEI_ENDL;
00458     os << "#numFieldsPerNode -- one entry per node " << FEI_ENDL;
00459     for(nn=0; nn<numNodesPerElement; nn++) os << numFieldsPerNode[nn]<<" ";
00460     os << FEI_ENDL << "#nodalFieldIDs -- one row per node" << FEI_ENDL;
00461     for(nn=0; nn<numNodesPerElement; ++nn) {
00462       for(nf=0; nf<numFieldsPerNode[nn]; ++nf) os << nodalFieldIDs[nn][nf] << " ";
00463       os << FEI_ENDL;
00464     }
00465     os << "#numElemDofFieldsPerElement" << FEI_ENDL
00466        << numElemDofFieldsPerElement<<FEI_ENDL;
00467     os << "#elemDofFieldIDs -- 'numElemDofFieldsPerElement' entries" << FEI_ENDL;
00468     for(nn=0; nn<numElemDofFieldsPerElement; ++nn) {
00469       os << elemDofFieldIDs[nn] << " ";
00470     }
00471     if (numElemDofFieldsPerElement > 0) os << FEI_ENDL;
00472     os << "#interleaveStrategy" << FEI_ENDL << interleaveStrategy << FEI_ENDL;
00473   }
00474    int j;
00475 
00476    CHK_ERR( addBlock(elemBlockID) );
00477    BlockDescriptor* block = NULL;
00478    CHK_ERR( getBlockDescriptor(elemBlockID, block) );
00479 
00480    CHK_ERR( block->setNumNodesPerElement(numNodesPerElement) );
00481    block->setNumElements(numElements);
00482    block->setElemDofFieldIDs(numElemDofFieldsPerElement, elemDofFieldIDs);
00483    block->setInterleaveStrategy(interleaveStrategy);
00484 
00485    int *fieldsPerNodePtr = block->fieldsPerNodePtr();
00486 
00487 // construct the list of nodal solution cardinalities for this block
00488 
00489    int numNodalEqns = 0;
00490    int countDOF;
00491    std::vector<int> distinctFields(0);
00492 
00493    for(j=0; j<numNodesPerElement; j++) {
00494 
00495       countDOF = 0;
00496       for(int k = 0; k < numFieldsPerNode[j]; k++) {
00497         fei::sortedListInsert(nodalFieldIDs[j][k], distinctFields);
00498 
00499          int fieldSize = getFieldSize(nodalFieldIDs[j][k]);
00500          if (fieldSize < 0) {
00501      fei::console_out() << "SNL_FEI_Structure::initElemBlock ERROR: fieldID " << 
00502        nodalFieldIDs[j][k] << " has negative size. " << FEI_ENDL;
00503      ERReturn(-1);
00504    }
00505          countDOF += fieldSize;
00506       }
00507 
00508       fieldsPerNodePtr[j] = numFieldsPerNode[j];
00509       numNodalEqns += countDOF;
00510    }
00511 
00512    block->setNumDistinctFields(distinctFields.size());
00513 
00514    int numElemDOFPerElement = 0;
00515    for(j=0; j<numElemDofFieldsPerElement; j++) {
00516       int fieldSize = getFieldSize(elemDofFieldIDs[j]);
00517       if (fieldSize < 0) {
00518   fei::console_out() << "SNL_FEI_Structure::initElemBlock ERROR: elemDoffieldID " << 
00519     elemDofFieldIDs[j] << " has negative size. " << FEI_ENDL;
00520   ERReturn(-1);
00521       }
00522       numElemDOFPerElement += fieldSize;
00523    }
00524 
00525    block->setNumElemDOFPerElement(numElemDOFPerElement);
00526    block->setNumEqnsPerElement(numNodalEqns + numElemDOFPerElement);
00527    block->setNumBlkEqnsPerElement(numNodesPerElement + numElemDOFPerElement);
00528 
00529 // cache a copy of the element fields array for later use...
00530 
00531    CHK_ERR( block->allocateFieldIDsTable() );
00532    int **fieldIDsTablePtr = block->fieldIDsTablePtr();
00533 
00534    for (j = 0; j < numNodesPerElement; j++) {
00535       for(int k = 0; k < numFieldsPerNode[j]; k++) {
00536          fieldIDsTablePtr[j][k] = nodalFieldIDs[j][k];
00537       }
00538    }
00539 
00540 //  create data structures for storage of element ID and topology info
00541 
00542    if (numElements > 0) {
00543       CHK_ERR( allocateBlockConnectivity(elemBlockID) );
00544    }
00545 
00546    return(FEI_SUCCESS);
00547 }
00548 
00549 //------------------------------------------------------------------------------
00550 int SNL_FEI_Structure::initElem(GlobalID elemBlockID,
00551              GlobalID elemID,
00552              const GlobalID* elemConn)
00553 {
00554   if (debugOutput_ && outputLevel_ > 2) {
00555     FEI_OSTREAM& os = dbgOut();
00556     os << "FEI: initElem"<<FEI_ENDL;
00557     os << "#blkID"<<FEI_ENDL<<(int)elemBlockID<<FEI_ENDL<<"#elmID"<<FEI_ENDL
00558        <<(int)elemID<< FEI_ENDL;
00559   }
00560 
00561   //first get the block-descriptor for this elemBlockID...
00562 
00563   BlockDescriptor* block = NULL;
00564   CHK_ERR( getBlockDescriptor(elemBlockID, block) );
00565 
00566   ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
00567 
00568   std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
00569   GlobalID* conn = &((*connTable.elem_conn_ids)[0]);
00570 
00571   int elemIndex = elemIDList.size();
00572   std::map<GlobalID,int>::iterator
00573     iter = elemIDList.find(elemID);
00574   
00575   bool redundantElem = false;
00576 
00577   if (iter != elemIDList.end()) {
00578     elemIndex = iter->second;
00579     redundantElem = true;
00580   }
00581   else {
00582     elemIDList.insert(std::make_pair(elemID,elemIndex));
00583   }
00584 
00585   int numNodes = block->getNumNodesPerElement();
00586 
00587   if (debugOutput_ && outputLevel_ > 2) {
00588     FEI_OSTREAM& os = dbgOut();
00589     os << "#n-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL<<"#nodeIDs"<<FEI_ENDL;
00590     for(int nn=0; nn<numNodes; nn++) os << (int)elemConn[nn] << " ";
00591     os << FEI_ENDL;
00592   }
00593 
00594   if (redundantElem) {
00595     //redundantElem == true means this element has been initialized before.
00596     //So we'll simply make sure the connectivity is the same as it was, and
00597     //if it isn't return -1.
00598 
00599     int offset = elemIndex*numNodes;
00600     for(int j=0; j<numNodes; j++) {
00601       if ( conn[offset+j] != elemConn[j]) {
00602   fei::console_out() << "SNL_FEI_Structure::initElem ERROR, elemID " << (int)elemID
00603        << " registered more than once, with differing connectivity."
00604        << FEI_ENDL;
00605   return(-1);
00606       }
00607     }
00608   }
00609   else {
00610     int offset = elemIndex*numNodes;
00611     for(int j = 0; j < numNodes; j++) {
00612       conn[offset+j] = elemConn[j];
00613     }
00614 
00615     CHK_ERR( nodeDatabase_->initNodeIDs(&(conn[offset]), numNodes) );
00616   }
00617 
00618   return(FEI_SUCCESS);
00619 }
00620 
00621 //------------------------------------------------------------------------------
00622 int SNL_FEI_Structure::initSlaveVariable(GlobalID slaveNodeID, 
00623            int slaveFieldID,
00624            int offsetIntoSlaveField,
00625            int numMasterNodes,
00626            const GlobalID* masterNodeIDs,
00627            const int* masterFieldIDs,
00628            const double* weights,
00629            double rhsValue)
00630 {
00631   if (debugOutput_) {
00632     FEI_OSTREAM& os = dbgOut();
00633     os << "FEI: initSlaveVariable" << FEI_ENDL;
00634     os << "#slvNodeID" << FEI_ENDL << (int)slaveNodeID << FEI_ENDL
00635        << "#slvFieldID"<< FEI_ENDL << slaveFieldID << FEI_ENDL
00636        << "#offsetIntoSlvField" << FEI_ENDL << offsetIntoSlaveField << FEI_ENDL
00637        << "#numMasterNodes" << FEI_ENDL << numMasterNodes << FEI_ENDL
00638        << "#masterNodeIDs" << FEI_ENDL;
00639     int nn;
00640     for(nn=0; nn<numMasterNodes; ++nn) {
00641       os <<(int)masterNodeIDs[nn]<<" ";
00642     }
00643     os << FEI_ENDL << "#masterFieldIDs" << FEI_ENDL;
00644     for(nn=0; nn<numMasterNodes; ++nn) {
00645       os <<masterFieldIDs[nn] << " ";
00646     }
00647     os << FEI_ENDL << "#field-sizes" << FEI_ENDL;
00648     for(nn=0; nn<numMasterNodes; ++nn) {
00649       int size = getFieldSize(masterFieldIDs[nn]);
00650       os << size << " ";
00651     }
00652     os << FEI_ENDL << "#weights" << FEI_ENDL;
00653     int offset = 0;
00654     for(nn=0; nn<numMasterNodes; ++nn) {
00655       int size = getFieldSize(masterFieldIDs[nn]);
00656       for(int j=0; j<size; ++j) {
00657   os << weights[offset++] << " ";
00658       }
00659     }
00660     os << FEI_ENDL << "#rhsValue" << FEI_ENDL << rhsValue << FEI_ENDL;
00661   }
00662 
00663   //Create and initialize a slave-variable object with the incoming data,
00664   //and add it to our list.
00665 
00666   SlaveVariable* svar = new SlaveVariable;
00667   svar->setNodeID(slaveNodeID);
00668   svar->setFieldID(slaveFieldID);
00669   svar->setFieldOffset(offsetIntoSlaveField);
00670 
00671   int woffset = 0;
00672 
00673   for(int i=0; i<numMasterNodes; i++) {
00674     svar->addMasterNodeID(masterNodeIDs[i]);
00675     svar->addMasterField(masterFieldIDs[i]);
00676     int fieldSize = getFieldSize(masterFieldIDs[i]);
00677     if (fieldSize < 0) ERReturn(-1);
00678 
00679     for(int j=0; j<fieldSize; j++) {
00680       svar->addWeight(weights[woffset++]);
00681     }
00682   }
00683 
00684   addSlaveVariable(svar);
00685 
00686   return(FEI_SUCCESS);
00687 }
00688 
00689 //------------------------------------------------------------------------------
00690 int SNL_FEI_Structure::deleteMultCRs()
00691 {
00692   int numMCRs = multCRs_.size();
00693   if (numMCRs > 0) {
00694     fei::destroyValues(multCRs_);
00695     multCRs_.clear();
00696   }
00697   return(0);
00698 }
00699 
00700 //------------------------------------------------------------------------------
00701 int SNL_FEI_Structure::initSharedNodes(int numSharedNodes,
00702               const GlobalID *sharedNodeIDs,  
00703               const int* numProcsPerNode, 
00704               const int *const *sharingProcIDs)
00705 {
00706   if (debugOutput_) {
00707     FEI_OSTREAM& os = dbgOut();
00708     os << "FEI: initSharedNodes"<<FEI_ENDL;
00709     os << "#n-nodes"<<FEI_ENDL<<numSharedNodes<<FEI_ENDL;
00710     os << "#num-procs-per-node"<<FEI_ENDL;
00711     int nn;
00712     for(nn=0; nn<numSharedNodes; ++nn) {
00713       os << numProcsPerNode[nn] << " ";
00714     }
00715     os << FEI_ENDL << "#following lines: nodeID   sharing-procs" << FEI_ENDL;
00716     for(nn=0; nn<numSharedNodes; nn++) {
00717       os << (int)sharedNodeIDs[nn]<<"   ";
00718       for(int np=0; np<numProcsPerNode[nn]; np++) os<<sharingProcIDs[nn][np]<<" ";
00719       os << FEI_ENDL;
00720     }
00721     os << "#end shared nodes"<<FEI_ENDL;
00722   }
00723 
00724   //  In this function we simply accumulate the incoming data into the
00725   //  NodeCommMgr object.
00726   //
00727   CHK_ERR( nodeCommMgr_->addSharedNodes(sharedNodeIDs, numSharedNodes,
00728           sharingProcIDs, numProcsPerNode) );
00729 
00730   if (activeNodesInitialized_) {
00731     //if active nodes have already been initialized, then we won't be
00732     //re-running the element-connectivities during the next call to
00733     //initComplete(), which means we won't have an opportunity to call the
00734     //NodeCommMgr::informLocal() method for nodes that reside locally. So we
00735     //need to check now for any locally residing shared nodes and make that
00736     //call. This is expensive, but it's a case that only arises when constraints
00737     //change between solves and nothing else changes. In general, constraints
00738     //are the only structural features allowed to change without requiring a
00739     //complete destruction and re-calculation of the FEI structure.
00740     //
00741     for(int i=0; i<numSharedNodes; ++i) {
00742       for(size_t nc=0; nc<connTables_.size(); ++nc) {
00743   if (connTables_[nc]->elem_conn_ids == NULL) continue;
00744   int len = connTables_[nc]->elem_conn_ids->size();
00745   if (len < 1) continue;
00746   GlobalID* conn_ids = &((*connTables_[nc]->elem_conn_ids)[0]);
00747   NodeDescriptor** nodes = &((*connTables_[nc]->elem_conn_ptrs)[0]);
00748 
00749   for(int j=0; j<len; ++j) {
00750     if (conn_ids[j] == sharedNodeIDs[i]) {
00751       CHK_ERR( nodeCommMgr_->informLocal(*(nodes[j])));
00752       break;
00753     }
00754   }
00755       }
00756     }
00757   }
00758 
00759   return(FEI_SUCCESS);
00760 }
00761 
00762 //------------------------------------------------------------------------------
00763 int SNL_FEI_Structure::initCRMult(int numCRNodes,
00764          const GlobalID* CRNodes,
00765          const int *CRFields,
00766          int& CRID)
00767 {
00768   if (debugOutput_) {
00769     FEI_OSTREAM& os = dbgOut();
00770     os << "FEI: initCRMult" << FEI_ENDL << "# numCRNodes: "<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
00771     os << "#CRNodes:"<<FEI_ENDL;
00772     int nn;
00773     for(nn=0; nn<numCRNodes; nn++) {
00774       os << (int)CRNodes[nn]<<" ";
00775     }
00776     os << FEI_ENDL<<"#fields:"<<FEI_ENDL;
00777     for(nn=0; nn<numCRNodes; nn++) {
00778       os << CRFields[nn]<<" ";
00779     }
00780     os<<FEI_ENDL;
00781   }
00782   //  tasks: add local constraint equations, determine sparsity
00783   //         patterns for the new constraint relation
00784   //
00785 
00786   CRID = localProc_*100000 + multCRs_.size();
00787   ConstraintType* multCRPtr = NULL;
00788   addCR(CRID, multCRPtr, multCRs_);
00789 
00790   ConstraintType& multCR = *multCRPtr;
00791 
00792   multCR.setConstraintID(CRID);
00793 
00794   std::vector<GlobalID>& CRNodeArray = multCR.getMasters();
00795   std::vector<int>& CRFieldArray = multCR.getMasterFieldIDs();
00796 
00797   for(int j = 0; j < numCRNodes; j++) {
00798     CRNodeArray.push_back(CRNodes[j]);
00799     CRFieldArray.push_back(CRFields[j]);
00800   }
00801 
00802   if (debugOutput_) dbgOut() << "#(output) CRID:"<<FEI_ENDL << CRID << FEI_ENDL;
00803 
00804   return(FEI_SUCCESS);
00805 }
00806 
00807 //------------------------------------------------------------------------------
00808 int SNL_FEI_Structure::initCRPen(int numCRNodes,
00809         const GlobalID* CRNodes,
00810         const int *CRFields,
00811         int& CRID)
00812 {
00813   if (debugOutput_) {
00814     FEI_OSTREAM& os = dbgOut();
00815     os << "initCRPen, numCRNodes: "<<numCRNodes<<FEI_ENDL;
00816     for(int nn=0; nn<numCRNodes; nn++) {
00817       os << "   crNodeID: "<<(int)CRNodes[nn]<<", field: "<<CRFields[nn]<<FEI_ENDL;
00818     }
00819   }
00820 
00821   CRID = localProc_*100000 + penCRs_.size();
00822 
00823   ConstraintType* penCRPtr = NULL;
00824   addCR(CRID, penCRPtr, penCRs_);
00825 
00826   ConstraintType& penCR = *penCRPtr;
00827 
00828   penCR.setConstraintID(CRID);
00829   penCR.setIsPenalty(true);
00830 
00831   std::vector<GlobalID>& CRNodesArray = penCR.getMasters();
00832 
00833   std::vector<int>& CRFieldArray = penCR.getMasterFieldIDs();
00834 
00835   for(int i = 0; i < numCRNodes; i++) {
00836     CRNodesArray.push_back(CRNodes[i]);
00837     CRFieldArray.push_back(CRFields[i]);
00838   }
00839 
00840   return(FEI_SUCCESS);
00841 }
00842 
00843 //------------------------------------------------------------------------------
00844 bool SNL_FEI_Structure::isInLocalElement(int nodeNumber)
00845 {
00846   //I'm going to make an assumption: if we're running in serial, then there
00847   //are no remote anythings, and the node in question must be in a local
00848   //element.
00849   if (numProcs_ < 2) return(true);
00850 
00851   const NodeDescriptor* node = NULL;
00852   int err = nodeDatabase_->getNodeWithNumber(nodeNumber, node);
00853   if (err != 0) return(false);
00854 
00855   GlobalID nodeID = node->getGlobalNodeID();
00856 
00857   //now let's find out if this node is a shared node.
00858   int shIndx = nodeCommMgr_->getSharedNodeIndex(nodeID);
00859   if (shIndx < 0) {
00860     //if shIndx < 0, then this isn't a shared node. For now, we will assume
00861     //that it must be a local node.
00862     return(true);
00863   }
00864 
00865   //If we reach this line, then the node is a shared node. Let's now ask if
00866   //it appears locally...
00867   std::vector<GlobalID>& localNodes = nodeCommMgr_->getLocalNodeIDs();
00868   int index = fei::binarySearch(nodeID, &localNodes[0], localNodes.size());
00869   if (index >= 0) return(true);
00870 
00871   //If we reach this line, then the node is shared, but doesn't appear in the
00872   //"localNodeIDs" held by NodeCommMgr. This means it is not in a local element,
00873   //so we'll return false.
00874   return(false);
00875 }
00876 
00877 //------------------------------------------------------------------------------
00878 int SNL_FEI_Structure::initComplete(bool generateGraph)
00879 {
00880   //This is the most significant function in SNL_FEI_Structure. This is where
00881   //the underlying matrix structure is calculated from all the data that has
00882   //been provided since construction. Calculating the equation space and
00883   //forming the matrix structure is a multi-step process, which proceeds like
00884   //this:
00885   //
00886   //1. finalizeActiveNodes()
00887   //  - makes sure that all shared nodes have been identified to the
00888   //    NodeDatabase, then has the NodeDatabase allocate its internal list of
00889   //    NodeDescriptor objects.
00890   //  - runs the element-connectivity tables and sets the nodal field and block
00891   //    info on the NodeDescriptors in the NodeDatabase. IMPORTANT: At this
00892   //    point, the nodeIDs in the connectivitiy tables are replaced with
00893   //    indices into the nodeDatabase to speed future lookups.
00894   //  - As the connectivity tables are being run, the NodeCommMgr is given the
00895   //    nodeID of each node that is connected to a local element, and the node's
00896   //    owner is set to the initial value of localProc_.
00897   //1a. finalizeNodeCommMgr()
00898   //  - NodeCommMgr::initComplete is called, at which time the ownership of all
00899   //    shared nodes is determined.
00900   //
00901   //2. Next, all local nodal degrees of freedom can be counted. This, together
00902   //  with the number of lagrange multipliers and element-centered DOFs are used
00903   //  by the function calcGlobalEqnInfo(), which determines global equation
00904   //  counts and offsets.
00905   //
00906   //3. setNodalEqnInfo() runs all locally active nodes and for each one that's
00907   //  locally owned, sets global equation-numbers on it, as well as its
00908   //  numNodalDOF and block-entry equation-number.
00909   //
00910   //4. setElemDOFEqnInfo() and setMultCREqnInfo() associate global equation-
00911   //  numbers with each of the lagrange multipliers and element-dofs.
00912   //
00913   //5. NodeDatabase::synchronize() is called, which associates nodeNumbers with
00914   //  each node, and also in turn calls NodeCommMgr::exchangeEqnInfo() to obtain
00915   //  equation-numbers and field/block info for remotely-owned shared nodes.
00916   //
00917   //6. setNumNodesAndEqnsPerBlock() then runs the active nodes again and counts
00918   //  active nodes and equations per element-block. This couldn't be done before
00919   //  because NodeCommMgr::exchangeEqnInfo() may have found out about new blocks
00920   //  that are associated with shared nodes on other processors.
00921   //
00922   //7. calculateSlaveEqns() associates equation-numbers with slave-nodes and
00923   //  master-nodes that have been identified locally, then gathers slave
00924   //  equation info from all processors onto all other processors, so that all
00925   //  processors know about all slave equations.
00926   //
00927   //8. initializeEqnCommMgr() runs the shared-nodes from NodeCommMgr and lets
00928   //  EqnCommMgr know which equations will be sent to other processors, and
00929   //  which will be received from other processors.
00930   //
00931   //9. translate localStartRow_ and localEndRow_ into the reduced equation
00932   //  space, where they are called reducedStartRow_ and reducedEndRow_, after
00933   //  which we know how many local equations there are not including any
00934   //  slave equations. At this point we can allocate the sysMatIndices_ array
00935   //  of arrays, which will hold the point-entry matrix structure.
00936   //
00937   //10. formMatrixStructure()
00938   //  - initElemBlockStructure() for each element in each element-block,
00939   //    obtains scatter-indices and makes the corresponding structurally
00940   //    symmetric contributions to the matrix structure. This process includes
00941   //    calculations associated with reducing out slave equations, and with
00942   //    putting off-processor contributions (element contributions with shared
00943   //    nodes) into the EqnCommMgr object.
00944   //  - initMultCRStructure() and initPenCRStructure() follow the same general
00945   //    routine as for the initElemBlockStructure() function.
00946   //  - EqnCommMgr::exchangeIndices() is called, which sends all shared
00947   //    contributions to the processors that own the corresponding equations,
00948   //    after which the receiving processors insert those contributions into
00949   //    the local matrix structure.
00950   //
00951   //11. initializeBlkEqnMapper() run all nodes, elem-dof, multiplier constraints
00952   //  and pass global point-equation-numbers with corresponding block-equation
00953   //  numbers to the snl_fei::PointBlockMap object which maintains a mapping
00954   //  between the point-entry and block-entry equation spaces.
00955   //
00956   //12. EqnCommMgr::exchangePtToBlkInfo() ensures that the point-to-block
00957   //  mapping data is globally consistent across all processors.
00958   //
00959   //13. The sysBlkMatIndices_ array is allocated, which is the array of
00960   //  arrays that will hold the block-entry matrix structure. This block-entry
00961   //  matrix structure is not filled, however, until later if/when the
00962   //  method getMatrixStructure is called specifically requesting block-entry
00963   //  structure data.
00964   //
00965   //14. Finally, if the debugOutputLevel_ is greater than 0, a file is written
00966   //  which contains a mapping from equations to nodes. This file is often very
00967   //  useful for post-mortem debugging operations, and is also necessary if the
00968   //  matrix produced from a serial run is to be compared with a matrix 
00969   //  produced from a similar run on multiple processors. The equation-to-node
00970   //  mapping is the only way to reconcile the equation-orderings between the
00971   //  two matrices.
00972   //
00973 
00974   if (debugOutput_) {
00975     FEI_OSTREAM& os = dbgOut();
00976     os << "FEI: initComplete" << FEI_ENDL;
00977   }
00978   // marshall all of the nodes we received in element connectivities into an
00979   // active node list, and set the fields, blocks and owning processors 
00980   // associated with those nodes.
00981 
00982   //IMPORTANT!!! Note that at this point, the nodeIDs in the
00983   //connectivitiy tables are being replaced with indices from the nodeDatabase.
00984 
00985   CHK_ERR( finalizeActiveNodes() );
00986 
00987   CHK_ERR( finalizeNodeCommMgr() );
00988 
00989   numLocalElemDOF_ = calcTotalNumElemDOF();
00990   numLocalMultCRs_ = calcNumMultCREqns();
00991 
00992   // ok, now the active equation calculations -- we need to count how many
00993   // local equations there are, then obtain the global equation info, then
00994   // set equation numbers on all of our active nodes.
00995 
00996   int numLocalNodes     = nodeDatabase_->countLocalNodeDescriptors(localProc_);
00997   numLocalNodalEqns_ = nodeDatabase_->countLocalNodalEqns(localProc_);
00998 
00999   numLocalEqns_ = numLocalNodalEqns_ + numLocalElemDOF_ + numLocalMultCRs_;
01000   numLocalEqnBlks_ = numLocalNodes   + numLocalElemDOF_ + numLocalMultCRs_;
01001 
01002   if (debugOutput_) {
01003     FEI_OSTREAM& os = dbgOut();
01004     os << "#   numMultCRs: " << numLocalMultCRs_ << ", numElemDOF: "
01005        << numLocalElemDOF_ << ", numLocalNodalEqns: " <<numLocalNodalEqns_
01006        << FEI_ENDL;
01007     os << "#   numLocalEqns_: " << numLocalEqns_ << FEI_ENDL;
01008   }
01009 
01010   calcGlobalEqnInfo(numLocalNodes, numLocalEqns_, numLocalEqnBlks_);
01011 
01012   CHK_ERR( setNodalEqnInfo() );
01013 
01014   //now we need to set equation numbers for the element-DOF's, and for the
01015   //lagrange-multipler constraints.
01016   //(exactly where to put these element DOF is an interesting issue... here, 
01017   //just toss them at the end of the nodal active eqns, which may or may not
01018   //be such a great choice.)
01019 
01020   setElemDOFEqnInfo();
01021 
01022   CHK_ERR( setMultCREqnInfo() );
01023 
01024   //now have the NodeDatabase run through the list of active nodes and
01025   //do whatever final initializations it needs to do. This function also
01026   //calls nodeCommMgr::exchangeEqnInfo.
01027   CHK_ERR( nodeDatabase_->synchronize(firstLocalNodeNumber_, localStartRow_,
01028               localProc_, comm_) );
01029 
01030   //now run the active nodes again and set the number of
01031   //active nodes and equations per block. We couldn't do this
01032   //before, because the nodeCommMgr::initComplete() and
01033   //nodeCommMgr::exchangeEqnInfo functions may both have found out about new
01034   //blocks that are associated with shared nodes.
01035 
01036   setNumNodesAndEqnsPerBlock();
01037 
01038   //at this point we'll run through any SlaveVariable records that were set,
01039   //and establish an EqnBuffer of slave equation numbers, and gather the slave
01040   //equations initialized on any other processors.
01041 
01042   slvCommMgr_ = new EqnCommMgr(comm_);
01043   slvCommMgr_->setNumRHSs(1);
01044 
01045   CHK_ERR( calculateSlaveEqns(comm_) );
01046 
01047   //From this point on, all computations involving equation numbers will assume
01048   //that those numbers are in the reduced equation space. So we'll start by
01049   //translating the contents of eqnCommMgr_ to the reduced space.
01050   //CHK_ERR( translateToReducedEqns(*eqnCommMgr_) );
01051 
01052   initializeEqnCommMgr();
01053 
01054   translateToReducedEqn(localEndRow_, reducedEndRow_);
01055   int startRow = localStartRow_;
01056   while(isSlaveEqn(startRow)) startRow++;
01057   translateToReducedEqn(startRow, reducedStartRow_);
01058   numLocalReducedRows_ = reducedEndRow_ - reducedStartRow_ + 1;
01059 
01060   generateGraph_ = generateGraph;
01061 
01062   if (generateGraph_) {
01063     sysMatIndices_ = new fei::ctg_set<int>[numLocalReducedRows_];
01064   }
01065 
01066   //Now we'll run through all the data structures that have been initialized
01067   //and form the matrix structure (lists of column-indices).
01068 
01069   CHK_ERR( formMatrixStructure() );
01070 
01071   CHK_ERR( initializeBlkEqnMapper() );
01072 
01073   if (globalMaxBlkSize_ > 1) {
01074     CHK_ERR( eqnCommMgr_->exchangePtToBlkInfo(*blkEqnMapper_) );
01075     if (numSlvs_ > 0) {
01076       try {
01077       CHK_ERR( slvCommMgr_->exchangeIndices() );
01078       }
01079       catch (std::runtime_error& exc) {
01080         fei::console_out() << exc.what() << FEI_ENDL;
01081         ERReturn(-1);
01082       }
01083 
01084       CHK_ERR( slvCommMgr_->exchangePtToBlkInfo(*blkEqnMapper_) );
01085     }
01086   }
01087 
01088   localReducedBlkOffset_ = ptEqnToBlkEqn(reducedStartRow_);
01089   int lastLocalReducedBlkEqn = ptEqnToBlkEqn(reducedEndRow_);
01090   numLocalReducedEqnBlks_ = lastLocalReducedBlkEqn - localReducedBlkOffset_ + 1;
01091 
01092   if (globalMaxBlkSize_ > 1 && generateGraph_) {
01093     sysBlkMatIndices_ = numLocalReducedEqnBlks_>0 ?
01094       new fei::ctg_set<int>[numLocalReducedEqnBlks_] : NULL;
01095   }
01096 
01097   delete slvCommMgr_;
01098 
01099   blockMatrix_ = true;
01100 
01101   structureFinalized_ = true;
01102   matIndicesDestroyed_ = false;
01103 
01104   //while we're here, let's write an equation -> nodeNumber map into a file
01105   //for possible debugging purposes... it will come in handy if we need to
01106   //compare matrices/vectors from different parallel runs.
01107   if (debugOutput_) writeEqn2NodeMap();
01108 
01109   if (debugOutput_) dbgOut() << "#leaving initComplete" << FEI_ENDL;
01110   return(FEI_SUCCESS);
01111 }
01112 
01113 //------------------------------------------------------------------------------
01114 int SNL_FEI_Structure::formMatrixStructure()
01115 {
01116   FEI_OSTREAM& os = dbgOut();
01117   if (debugOutput_) os << "#   formMatrixStructure" << FEI_ENDL;
01118 
01119   //do our local equation profile calculations (i.e., determine
01120   //how long each row is). use the element scatter arrays to determine 
01121   //the sparsity pattern
01122 
01123   CHK_ERR( initElemBlockStructure() );
01124 
01125   // next, handle the matrix structure imposed by the constraint relations...
01126   //
01127 
01128   CHK_ERR( initMultCRStructure() );
01129 
01130   // we also need to accomodate penalty constraints, so let's process them
01131   // now...
01132 
01133   CHK_ERR( initPenCRStructure() );
01134 
01135   //we've now processed all of the local data that produces matrix structure.
01136   //so now let's have the equation comm manager exchange indices with the
01137   //other processors so we can account for remote contributions to our
01138   //matrix structure.
01139   try {
01140   CHK_ERR( eqnCommMgr_->exchangeIndices(&os) );
01141   }
01142   catch(std::runtime_error& exc) {
01143     fei::console_out() << exc.what() << FEI_ENDL;
01144     ERReturn(-1);
01145   }
01146 
01147   //so now the remote contributions should be available, let's get them out
01148   //of the eqn comm mgr and put them into our local matrix structure.
01149 
01150   int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
01151   std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
01152   std::vector<fei::CSVec*>& recvEqns = eqnCommMgr_->localEqns();
01153   int i;
01154   if (debugOutput_) {
01155     os << "#     after eqnCommMgr_->exchangeIndices, numRecvEqns: "
01156        << numRecvEqns << FEI_ENDL;
01157   }
01158 
01159   for(i=0; i<numRecvEqns; i++) {
01160     int eqn = recvEqnNumbers[i];
01161     if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
01162       fei::console_out() << "SNL_FEI_Structure::initComplete: ERROR, recvEqn " << eqn
01163      << " out of range. (reducedStartRow_: " << reducedStartRow_
01164      << ", reducedEndRow_: " << reducedEndRow_ << ", localProc_: "
01165      << localProc_ << ")" << FEI_ENDL;
01166       return(-1);
01167     }
01168 
01169     for(size_t j=0; j<recvEqns[i]->size(); j++) {
01170       CHK_ERR( createMatrixPosition(eqn, recvEqns[i]->indices()[j],
01171             "frmMatStr") );
01172     }
01173   }
01174 
01175   if (debugOutput_) {
01176     os << "#  leaving formMatrixStructure" << FEI_ENDL;
01177   }
01178 
01179   return(0);
01180 }
01181 
01182 //------------------------------------------------------------------------------
01183 int SNL_FEI_Structure::initElemBlockStructure()
01184 {
01185   int numBlocks = getNumElemBlocks();
01186 
01187   FEI_OSTREAM& os = dbgOut();
01188   if (debugOutput_) {
01189     os << "#     initElemBlockStructure" << FEI_ENDL;
01190     os << "#        numElemBlocks: " << numBlocks << FEI_ENDL;
01191   }
01192 
01193   for (int bIndex = 0; bIndex < numBlocks; bIndex++) {
01194     BlockDescriptor* block = NULL;
01195     CHK_ERR( getBlockDescriptor_index(bIndex, block) );
01196 
01197     int numEqns = block->getNumEqnsPerElement();
01198     int interleave = block->getInterleaveStrategy();
01199     std::vector<int> scatterIndices(numEqns);
01200 
01201     GlobalID elemBlockID = block->getGlobalBlockID();
01202     ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
01203 
01204     std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
01205     block->setNumElements(elemIDList.size());
01206     int numBlockElems = block->getNumElements();
01207 
01208     //loop over all the elements, determining the elemental (both from nodes 
01209     //and from element DOF) contributions to the sparse matrix structure
01210     if (debugOutput_) {
01211       os << "#        block " << bIndex << ", numElems: " << numBlockElems<<FEI_ENDL;
01212     }
01213     for(int elemIndex = 0; elemIndex < numBlockElems; elemIndex++) {
01214 
01215       getScatterIndices_index(bIndex, elemIndex, interleave,
01216                               &scatterIndices[0]);
01217 
01218       //now, store the structure that will arise from this contribution,
01219       //after first performing any manipulations associated with slave
01220       //reduction.
01221 
01222       CHK_ERR( createSymmEqnStructure(scatterIndices) );
01223     }
01224   }
01225 
01226   if (reducedEqnCounter_ > 0) CHK_ERR( assembleReducedStructure() );
01227 
01228   if (debugOutput_) os << "#     leaving initElemBlockStructure" << FEI_ENDL;
01229   return(FEI_SUCCESS);
01230 }
01231 
01232 //------------------------------------------------------------------------------
01233 int SNL_FEI_Structure::getMatrixRowLengths(std::vector<int>& rowLengths)
01234 {
01235   if (!structureFinalized_) return(-1);
01236 
01237   rowLengths.resize(numLocalReducedRows_);
01238 
01239   for(int i=0; i<numLocalReducedRows_; i++) {
01240     rowLengths[i] = sysMatIndices_[i].size();
01241   }
01242   return(0);
01243 }
01244 
01245 //------------------------------------------------------------------------------
01246 int SNL_FEI_Structure::getMatrixStructure(int** indices,
01247            std::vector<int>& rowLengths)
01248 {
01249   if (!structureFinalized_) return(-1);
01250 
01251   rowLengths.resize(numLocalReducedRows_);
01252 
01253   for(int i=0; i<numLocalReducedRows_; i++) {
01254     rowLengths[i] = sysMatIndices_[i].size();
01255     fei::copySetToArray(sysMatIndices_[i], rowLengths[i], indices[i]);
01256   }
01257 
01258   return(0);
01259 }
01260 
01261 //------------------------------------------------------------------------------
01262 int SNL_FEI_Structure::getMatrixStructure(int** ptColIndices,
01263            std::vector<int>& ptRowLengths,
01264            int** blkColIndices,
01265             int* blkIndices_1D,
01266            std::vector<int>& blkRowLengths,
01267            std::vector<int>& numPtRowsPerBlkRow)
01268 {
01269   int err = getMatrixStructure(ptColIndices, ptRowLengths);
01270   if (err != 0) return(-1);
01271 
01272   if (globalMaxBlkSize_ == 1) {
01273     //No block-equations have more than one point-equation, so we'll just assign
01274     //the block-structure arrays to be the same as the point-structure arrays.
01275     int numRows = ptRowLengths.size();
01276     blkRowLengths.resize(numRows);
01277     numPtRowsPerBlkRow.assign(numRows, 1);
01278 
01279     blkRowLengths.resize(ptRowLengths.size());
01280     for(size_t ii=0; ii<ptRowLengths.size(); ++ii) {
01281       blkRowLengths[ii] = ptRowLengths[ii];
01282     }
01283 
01284     for(int i=0; i<numRows; i++) {
01285       blkColIndices[i] = ptColIndices[i];
01286     }
01287   }
01288   else {
01289     //There are block-equations with more than 1 point-equation, so let's
01290     //calculate the actual block-structure.
01291 
01292     std::map<int,int>* ptEqns = blkEqnMapper_->getPtEqns();
01293     int numPtEqns = ptEqns->size();
01294 
01295     std::map<int,int>::const_iterator
01296       pteq = ptEqns->begin(),
01297       pteq_end = ptEqns->end();
01298 
01299     int lastBlkRow = -1;
01300 
01301     for(int jj=0; jj<numPtEqns; ++jj, ++pteq) {
01302       int ptEqn = (*pteq).first;
01303       int localPtEqn = ptEqn - reducedStartRow_;
01304       if (localPtEqn < 0 || localPtEqn >= numLocalReducedRows_) continue;
01305 
01306       int rowLength = sysMatIndices_[localPtEqn].size();
01307 
01308       int blkRow = blkEqnMapper_->eqnToBlkEqn(ptEqn);
01309       if (blkRow < 0) {
01310   ERReturn(-1);
01311       }
01312 
01313       if (blkRow == lastBlkRow) continue;
01314 
01315       int localBlkRow = blkRow - localReducedBlkOffset_;
01316       if (localBlkRow < 0 || localBlkRow >= numLocalReducedEqnBlks_) {
01317   ERReturn(-1);
01318       }
01319 
01320       fei::ctg_set<int>& sysBlkIndices = sysBlkMatIndices_[localBlkRow];
01321 
01322       int* theseColIndices = ptColIndices[localPtEqn];
01323 
01324       for(int colj=0; colj<rowLength; colj++) {
01325   int blkCol = blkEqnMapper_->eqnToBlkEqn(theseColIndices[colj]);
01326   if (blkCol < 0) {
01327     fei::console_out() << localProc_
01328          <<"SNL_FEI_Structure::getMatrixStructure ERROR pt row "
01329          << ptEqn << ", pt col "
01330          << ptColIndices[localPtEqn][colj]
01331          << " doesn't have a corresponding block" << FEI_ENDL;
01332     blkCol = blkEqnMapper_->eqnToBlkEqn(theseColIndices[colj]);
01333     std::abort();
01334   }
01335 
01336   sysBlkIndices.insert2(blkCol);
01337       }
01338 
01339       lastBlkRow = blkRow;
01340     }
01341 
01342     //now we're ready to set the block-sizes...
01343 
01344     numPtRowsPerBlkRow.resize(numLocalReducedEqnBlks_);
01345     blkRowLengths.resize(numLocalReducedEqnBlks_);
01346     int offset = 0;
01347     for(int i=0; i<numLocalReducedEqnBlks_; i++) {
01348       blkRowLengths[i] = sysBlkMatIndices_[i].size();
01349       fei::copySetToArray(sysBlkMatIndices_[i],blkRowLengths[i],
01350                &(blkIndices_1D[offset]));
01351       offset += blkRowLengths[i];
01352 
01353       int blkEqn = localReducedBlkOffset_ + i;
01354       numPtRowsPerBlkRow[i] = blkEqnMapper_->getBlkEqnSize(blkEqn);
01355     }
01356   }
01357 
01358   return(0);
01359 }
01360 
01361 //------------------------------------------------------------------------------
01362 bool SNL_FEI_Structure::nodalEqnsAllSlaves(const NodeDescriptor* node,
01363              std::vector<int>& slaveEqns)
01364 {
01365   int numFields = node->getNumFields();
01366   const int* fieldIDs = node->getFieldIDList();
01367   const int* fieldEqns= node->getFieldEqnNumbers();
01368 
01369   for(int i=0; i<numFields; ++i) {
01370     int fieldSize = getFieldSize(fieldIDs[i]);
01371     for(int eqn=0; eqn<fieldSize; ++eqn) {
01372       int thisEqn = fieldEqns[i] + eqn;
01373       if (fei::binarySearch(thisEqn, slaveEqns) < 0) {
01374   //found a nodal eqn that's not a slave eqn, so return false.
01375   return(false);
01376       }
01377     }
01378   }
01379 
01380   return(true);
01381 }
01382 
01383 //------------------------------------------------------------------------------
01384 NodeDescriptor& SNL_FEI_Structure::findNodeDescriptor(GlobalID nodeID)
01385 {
01386 //
01387 //This function returns a NodeDescriptor reference if nodeID is an active node.
01388 //
01389   NodeDescriptor* node = NULL;
01390   int err = nodeDatabase_->getNodeWithID(nodeID, node);
01391 
01392   if (err != 0) {
01393     fei::console_out() << "ERROR, findNodeDescriptor unable to find node " << (int)nodeID
01394    << FEI_ENDL;
01395     std::abort();
01396   }
01397 
01398   return( *node );
01399 }
01400 
01401 //------------------------------------------------------------------------------
01402 NodeDescriptor* SNL_FEI_Structure::findNode(GlobalID nodeID)
01403 {
01404   NodeDescriptor* node = NULL;
01405   nodeDatabase_->getNodeWithID(nodeID, node);
01406   return node;
01407 }
01408 
01409 //------------------------------------------------------------------------------
01410 int SNL_FEI_Structure::initMultCRStructure() 
01411 {
01412   FEI_OSTREAM& os = dbgOut();
01413   if (debugOutput_) os << "#    initMultCRStructure" << FEI_ENDL;
01414   //
01415   //Private SNL_FEI_Structure function, to be called from initComplete.
01416   //
01417   //Records the system matrix structure arising from Lagrange Multiplier
01418   //Constraint Relations.
01419   //
01420   // since at initialization all we are being passed is the
01421   // general form of the constraint equations, we can't check to see if any
01422   // of the weight terms are zeros at this stage of the game.  Hence, we
01423   // have to reserve space for all the nodal weight vectors, even though
01424   // they might turn out to be zeros during the load step....
01425 
01426   std::map<GlobalID,ConstraintType*>::const_iterator
01427     cr_iter = multCRs_.begin(),
01428     cr_end  = multCRs_.end();
01429 
01430   while(cr_iter != cr_end) {
01431     ConstraintType& multCR = *((*cr_iter).second);
01432 
01433     int lenList = multCR.getMasters().size();
01434 
01435     std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
01436     GlobalID *CRNodePtr = &CRNode_vec[0];
01437     std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
01438     int* CRFieldPtr = &CRField_vec[0];
01439 
01440     int crEqn = multCR.getEqnNumber();
01441     int reducedCrEqn = 0;
01442     translateToReducedEqn(crEqn, reducedCrEqn);
01443 
01444     createMatrixPosition(reducedCrEqn, reducedCrEqn, "initMCRStr");
01445 
01446     for(int j=0; j<lenList; j++) {
01447       GlobalID nodeID = CRNodePtr[j];
01448       int fieldID = CRFieldPtr[j];
01449 
01450       NodeDescriptor *nodePtr = findNode(nodeID);
01451       if(nodePtr==NULL) continue;
01452       NodeDescriptor& node = *nodePtr;
01453 
01454       //first, store the column indices associated with this node, in
01455       //the constraint's equation. i.e., position (crEqn, node)
01456 
01457       storeNodalColumnIndices(crEqn, node, fieldID);
01458 
01459       //now we need to store the transpose of the above contribution,
01460       //i.e., position(s) (node, crEqn)
01461 
01462       if (node.getOwnerProc() == localProc_) {
01463   //if we own this node, we will simply store its equation
01464   //numbers in local rows (equations) of the matrix.
01465 
01466   storeNodalRowIndices(node, fieldID, crEqn);
01467       }
01468       else {
01469   //if we don't own it, then we need to register with the
01470   //eqn comm mgr that we'll be sending contributions to
01471   //column crEqn of the remote equations associated with this
01472   //node.
01473 
01474   storeNodalSendIndex(node, fieldID, crEqn);
01475       }
01476     }
01477     ++cr_iter;
01478   }
01479   return(FEI_SUCCESS);
01480 }
01481 
01482 //------------------------------------------------------------------------------
01483 int SNL_FEI_Structure::initPenCRStructure()
01484 {
01485   FEI_OSTREAM& os = dbgOut();
01486   if (debugOutput_) os << "#     initPenCRStructure" << FEI_ENDL;
01487 //
01488 //This function oversees the putting in of any matrix structure that results
01489 //from Penalty constraint relations.
01490 //
01491 // note that penalty constraints don't generate new equations
01492 // (unlike Lagrange constraints), but they do add terms to the system
01493 // stiffness matrix that couple all the nodes that contribute to the
01494 // penalty constraint.  In addition, each submatrix is defined by the pair
01495 // of nodes that creates its contribution, hence a submatrix can be defined
01496 // in terms of two weight vectors (of length p and q) instead of the
01497 // generally larger product matrix (of size pq)
01498 
01499 // the additional terms take the form of little submatrices that look a lot 
01500 // like an element stiffness and load matrix, where the nodes in the
01501 // constraint list take on the role of the nodes associated with an
01502 // element, and the individual matrix terms arise from the outer products
01503 // of the constraint nodal weight vectors 
01504 
01505 // rather than get elegant and treat this task as such an elemental energy
01506 // term, we'll use some brute force to construct these penalty contributions
01507 // on the fly, primarily to simplify -reading- this thing, so that the 
01508 // derivations in the annotated implementation document are more readily
01509 // followed...
01510   std::map<GlobalID,ConstraintType*>::const_iterator
01511     cr_iter = penCRs_.begin(),
01512     cr_end = penCRs_.end();
01513 
01514   while(cr_iter != cr_end) {
01515     ConstraintType& penCR = *((*cr_iter).second);
01516 
01517     int lenList = penCR.getMasters().size();
01518     std::vector<GlobalID>& CRNode_vec = penCR.getMasters();
01519     GlobalID* CRNodesPtr = &CRNode_vec[0];
01520 
01521     std::vector<int>& CRField_vec = penCR.getMasterFieldIDs();
01522     int* CRFieldPtr = &CRField_vec[0];
01523 
01524     // each constraint equation generates a set of nodal energy terms, so
01525     // we have to process a matrix of nodes for each constraint
01526 
01527     for(int i = 0; i < lenList; i++) {
01528       GlobalID iNodeID = CRNodesPtr[i];
01529       int iField = CRFieldPtr[i];
01530 
01531       NodeDescriptor* iNodePtr = findNode(iNodeID);
01532       if(!iNodePtr) continue;
01533       NodeDescriptor& iNode = *iNodePtr;
01534 
01535       for(int j = 0; j < lenList; j++) {
01536   GlobalID jNodeID = CRNodesPtr[j];
01537   int jField = CRFieldPtr[j];
01538 
01539   NodeDescriptor* jNodePtr = findNode(jNodeID);
01540         if(!jNodePtr) continue;
01541         NodeDescriptor &jNode = *jNodePtr;
01542 
01543   if (iNode.getOwnerProc() == localProc_) {
01544     //if iNode is local, we'll put equations into the local 
01545     //matrix structure.
01546 
01547     storeLocalNodeIndices(iNode, iField, jNode, jField);
01548   }
01549   else {
01550     //if iNode is remotely owned, we'll be registering equations
01551     //to send to its owning processor.
01552 
01553     storeNodalSendIndices(iNode, iField, jNode, jField);
01554   }
01555       }
01556     }   //   end i loop
01557     ++cr_iter;
01558    }   //   end while loop
01559    return(FEI_SUCCESS);
01560 }
01561  
01562 //------------------------------------------------------------------------------
01563 void SNL_FEI_Structure::storeNodalSendIndex(NodeDescriptor& node, int fieldID,
01564               int col)
01565 {
01566   //
01567   //This is a private SNL_FEI_Structure function. We can safely assume that it
01568   //will only be called with a node that is not locally owned.
01569   //
01570   //This function tells the eqn comm mgr that we'll be sending contributions
01571   //to column 'col' for the equations associated with 'fieldID', on 'node', on
01572   //node's owning processor.
01573   //
01574   int proc = node.getOwnerProc();
01575 
01576   int eqnNumber = -1;
01577   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
01578 
01579   int numEqns = getFieldSize(fieldID);
01580   if (numEqns < 1) {
01581     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01582    <<") with size "<<numEqns<<FEI_ENDL;
01583     voidERReturn;
01584   }
01585 
01586   for(int i=0; i<numEqns; i++) {
01587     eqnCommMgr_->addRemoteIndices(eqnNumber+i, proc, &col, 1);
01588   }
01589 }
01590 
01591 //------------------------------------------------------------------------------
01592 void SNL_FEI_Structure::storeNodalSendIndices(NodeDescriptor& iNode, int iField,
01593                                   NodeDescriptor& jNode, int jField){
01594 //
01595 //This function will register with the eqn comm mgr the equations associated
01596 //with iNode, field 'iField' having column indices that are the equations
01597 //associated with jNode, field 'jField', to be sent to the owner of iNode.
01598 //
01599    int proc = iNode.getOwnerProc();
01600 
01601    int iEqn = -1, jEqn = -1;
01602    if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
01603    if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
01604 
01605    int iNumParams = getFieldSize(iField);
01606    int jNumParams = getFieldSize(jField);
01607    if (iNumParams < 1 || jNumParams < 1) {
01608      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
01609     << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
01610     << jNumParams<<FEI_ENDL;
01611      voidERReturn;
01612    }
01613 
01614    for(int i=0; i<iNumParams; i++) {
01615       int eqn = iEqn + i;
01616 
01617       for(int j=0; j<jNumParams; j++) {
01618          int col = jEqn + j;
01619          eqnCommMgr_->addRemoteIndices(eqn, proc, &col, 1);
01620       }
01621    }
01622 }
01623 
01624 //------------------------------------------------------------------------------
01625 void SNL_FEI_Structure::storeLocalNodeIndices(NodeDescriptor& iNode, int iField,
01626                NodeDescriptor& jNode, int jField)
01627 {
01628 //
01629 //This function will add to the local matrix structure the equations associated
01630 //with iNode at iField, having column indices that are the equations associated
01631 //with jNode at jField.
01632 //
01633    int iEqn = -1, jEqn = -1;
01634    if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
01635    if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
01636 
01637    int iNumParams = getFieldSize(iField);
01638    int jNumParams = getFieldSize(jField);
01639    if (iNumParams < 1 || jNumParams < 1) {
01640      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
01641     << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
01642     << jNumParams<<FEI_ENDL;
01643      voidERReturn;
01644    }
01645 
01646    for(int i=0; i<iNumParams; i++) {
01647       int row = iEqn + i;
01648       int reducedRow = -1;
01649       bool isSlave = translateToReducedEqn(row, reducedRow);
01650       if (isSlave) continue;
01651       
01652       for(int j=0; j<jNumParams; j++) {
01653          int col = jEqn + j;
01654    int reducedCol = -1;
01655    isSlave = translateToReducedEqn(col, reducedCol);
01656    if (isSlave) continue;
01657 
01658          createMatrixPosition(reducedRow, reducedCol,
01659             "storeLocNdInd");
01660       }
01661    }
01662 }
01663 
01664 //------------------------------------------------------------------------------
01665 void SNL_FEI_Structure::storeNodalColumnIndices(int eqn, NodeDescriptor& node,
01666                  int fieldID)
01667 {
01668   //
01669   //This function stores the equation numbers associated with 'node' at
01670   //'fieldID' as column indices in row 'eqn' of the system matrix structure.
01671   //
01672   if ((localStartRow_ > eqn) || (eqn > localEndRow_)) {
01673     return;
01674   }
01675 
01676   int colNumber = -1;
01677   if (!node.getFieldEqnNumber(fieldID, colNumber)) voidERReturn;
01678 
01679   int numParams = getFieldSize(fieldID);
01680   if (numParams < 1) {
01681     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01682    <<") with size "<<numParams<<FEI_ENDL;
01683     voidERReturn;
01684   }
01685 
01686   int reducedEqn = -1;
01687   bool isSlave = translateToReducedEqn(eqn, reducedEqn);
01688   if (isSlave) return;
01689 
01690   for(int j=0; j<numParams; j++) {
01691     int reducedCol = -1;
01692     isSlave = translateToReducedEqn(colNumber+j, reducedCol);
01693     if (isSlave) continue;
01694 
01695     createMatrixPosition(reducedEqn, reducedCol,
01696        "storeNdColInd");
01697   }
01698 }
01699 
01700 //------------------------------------------------------------------------------
01701 void SNL_FEI_Structure::storeNodalRowIndices(NodeDescriptor& node,
01702               int fieldID, int eqn)
01703 {
01704   //
01705   //This function stores column 'eqn' in equation numbers associated with
01706   //'node' at 'fieldID' in the system matrix structure. 
01707   //
01708   int eqnNumber = -1;
01709   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
01710 
01711   int numParams = getFieldSize(fieldID);
01712   if (numParams < 1) {
01713     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01714    <<") with size "<<numParams<<FEI_ENDL;
01715     voidERReturn;
01716   }
01717 
01718   int reducedEqn = -1;
01719   bool isSlave = translateToReducedEqn(eqn, reducedEqn);
01720   if (isSlave) return;
01721 
01722   for(int j=0; j<numParams; j++) {
01723     int reducedRow = -1;
01724     isSlave = translateToReducedEqn(eqnNumber+j, reducedRow);
01725     if (isSlave) continue;
01726 
01727     createMatrixPosition(reducedRow, reducedEqn, "storeNdRowInd");
01728   }
01729 }
01730 
01731 //------------------------------------------------------------------------------
01732 int SNL_FEI_Structure::createSymmEqnStructure(std::vector<int>& scatterIndices)
01733 {
01734   //scatterIndices are both the rows and the columns for a structurally
01735   //symmetric 2-dimensional contribution to the matrix.
01736   //
01737 
01738   //if there aren't any slaves in the whole problem, store the scatter-indices
01739   //using a fast function (no translations-to-reduced-space) and return.
01740   if (numSlaveEquations() == 0) {
01741     CHK_ERR( storeElementScatterIndices_noSlaves(scatterIndices) );
01742     return(FEI_SUCCESS);
01743   }
01744 
01745   try {
01746 
01747   int len = scatterIndices.size();
01748   bool anySlaves = false;
01749   rSlave_.resize(len);
01750   for(int is=0; is<len; is++) { 
01751     rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
01752     if (rSlave_[is] == 1) anySlaves = true;
01753   }
01754 
01755   //if there aren't any slaves in this contribution, then just store it
01756   //and return
01757   if (!anySlaves) {
01758     CHK_ERR( storeElementScatterIndices(scatterIndices) );
01759     return(FEI_SUCCESS);
01760   }
01761 
01762   int* scatterPtr = &scatterIndices[0];
01763 
01764   workSpace_.resize(len);
01765   for(int j=0; j<len; j++) {
01766     translateToReducedEqn(scatterPtr[j], workSpace_[j]);
01767   }
01768   
01769   for(int i=0; i<len; i++) {
01770     int row = scatterPtr[i];
01771     if (rSlave_[i] == 1) {
01772       reducedEqnCounter_++;
01773       //
01774       //'row' is a slave equation, so add this row to Kdi_. But as we do that,
01775       //watch for columns that are slave equations and add them to Kid_.
01776       //
01777       for(int jj=0; jj<len; jj++) {
01778   int col = scatterPtr[jj];
01779 
01780   if (rSlave_[jj] == 1) {
01781     //'col' is a slave equation, so add this column to Kid_.
01782     for(int ii=0; ii<len; ii++) {
01783       int rowi = scatterPtr[ii];
01784 
01785       //only add the non-slave rows for this column.
01786       if (rSlave_[ii] == 1) continue;
01787 
01788       Kid_->createPosition(rowi, col);
01789     }
01790     //now skip to the next column 
01791     continue;
01792   }
01793 
01794   Kdi_->createPosition(row, col);
01795       }
01796 
01797       //finally, add all slave columns in this slave row to Kdd_.
01798       for(int kk=0; kk<len; kk++) {
01799   int colk = scatterPtr[kk];
01800 
01801   if (rSlave_[kk] != 1) continue;
01802 
01803   Kdd_->createPosition(row, colk);
01804       }
01805     }
01806     else {
01807       //this is not a slave row, so we will loop across it, creating matrix
01808       //positions for all non-slave columns in it.
01809       int reducedRow = -1;
01810       translateToReducedEqn(row, reducedRow);
01811 
01812       bool rowIsLocal = true;
01813       int owner = localProc_;
01814       if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
01815   rowIsLocal = false;
01816   owner = getOwnerProcForEqn(row);
01817       }
01818 
01819       for(int j=0; j<len; j++) {
01820   if (rSlave_[j] == 1) continue;
01821 
01822   int reducedCol = workSpace_[j];
01823 
01824   if (rowIsLocal) {
01825     CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
01826           "crtSymmEqnStr") );
01827   }
01828   else {
01829     eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
01830   }
01831       }
01832     }
01833   }
01834 
01835   if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
01836 
01837   }
01838   catch(std::runtime_error& exc) {
01839     fei::console_out() << exc.what() << FEI_ENDL;
01840     ERReturn(-1);
01841   }
01842 
01843   return(FEI_SUCCESS);
01844 }
01845 
01846 //------------------------------------------------------------------------------
01847 int SNL_FEI_Structure::createBlkSymmEqnStructure(std::vector<int>& scatterIndices)
01848 {
01849   //scatterIndices are both the rows and the columns for a structurally
01850   //symmetric 2-dimensional contribution to the matrix.
01851   //
01852 
01853   //if there aren't any slaves in the whole problem, store the scatter-indices
01854   //using a fast function (no translations-to-reduced-space) and return.
01855   if (numSlaveEquations() == 0) {
01856     return( storeElementScatterBlkIndices_noSlaves(scatterIndices) );
01857   }
01858 
01859   try {
01860 
01861   int len = scatterIndices.size();
01862   bool anySlaves = false;
01863   rSlave_.resize(len);
01864   for(int is=0; is<len; is++) { 
01865     rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
01866     if (rSlave_[is] == 1) anySlaves = true;
01867   }
01868 
01869   //if there aren't any slaves in this contribution, then just store it
01870   //and return
01871   if (!anySlaves) {
01872     CHK_ERR( storeElementScatterIndices(scatterIndices) );
01873     return(FEI_SUCCESS);
01874   }
01875 
01876   int* scatterPtr = &scatterIndices[0];
01877 
01878   workSpace_.resize(len);
01879   for(int j=0; j<len; j++) {
01880     translateToReducedEqn(scatterPtr[j], workSpace_[j]);
01881   }
01882   
01883   for(int i=0; i<len; i++) {
01884     int row = scatterPtr[i];
01885     if (rSlave_[i] == 1) {
01886       reducedEqnCounter_++;
01887       //
01888       //'row' is a slave equation, so add this row to Kdi_. But as we do that,
01889       //watch for columns that are slave equations and add them to Kid_.
01890       //
01891       for(int jj=0; jj<len; jj++) {
01892   int col = scatterPtr[jj];
01893 
01894   if (rSlave_[jj] == 1) {
01895     //'col' is a slave equation, so add this column to Kid_.
01896     for(int ii=0; ii<len; ii++) {
01897       int rowi = scatterPtr[ii];
01898 
01899       //only add the non-slave rows for this column.
01900       if (rSlave_[ii] == 1) continue;
01901 
01902       Kid_->createPosition(rowi, col);
01903     }
01904     //now skip to the next column 
01905     continue;
01906   }
01907 
01908   Kdi_->createPosition(row, col);
01909       }
01910 
01911       //finally, add all slave columns in this slave row to Kdd_.
01912       for(int kk=0; kk<len; kk++) {
01913   int colk = scatterPtr[kk];
01914 
01915   if (rSlave_[kk] != 1) continue;
01916 
01917   Kdd_->createPosition(row, colk);
01918       }
01919     }
01920     else {
01921       //this is not a slave row, so we will loop across it, creating matrix
01922       //positions for all non-slave columns in it.
01923       int reducedRow = -1;
01924       translateToReducedEqn(row, reducedRow);
01925 
01926       bool rowIsLocal = true;
01927       int owner = localProc_;
01928       if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
01929   rowIsLocal = false;
01930   owner = getOwnerProcForEqn(row);
01931       }
01932 
01933       for(int j=0; j<len; j++) {
01934   if (rSlave_[j] == 1) continue;
01935 
01936   int reducedCol = workSpace_[j];
01937 
01938   if (rowIsLocal) {
01939     CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
01940           "crtSymmEqnStr") );
01941   }
01942   else {
01943     eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
01944   }
01945       }
01946     }
01947   }
01948 
01949   if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
01950 
01951   }
01952   catch(std::runtime_error& exc) {
01953     fei::console_out() << exc.what() << FEI_ENDL;
01954     ERReturn(-1);
01955   }
01956 
01957   return(FEI_SUCCESS);
01958 }
01959 
01960 //------------------------------------------------------------------------------
01961 int SNL_FEI_Structure::storeElementScatterIndices(std::vector<int>& indices)
01962 {
01963 //This function takes a list of equation numbers, as input, and for each
01964 //one, if it is a local equation number, goes to that row of the sysMatIndices
01965 //structure and stores the whole list of equation numbers as column indices
01966 //in that row of the matrix structure. If the equation number is not local,
01967 //then it is given to the EqnCommMgr.
01968 //
01969    int numIndices = indices.size();
01970    int* indPtr = &indices[0];
01971 
01972    workSpace_.resize(numIndices);
01973    int* workPtr = &workSpace_[0];
01974    int reducedEqn = -1;
01975    for(size_t j=0; j<indices.size(); j++) {
01976      bool isSlave = translateToReducedEqn(indPtr[j], reducedEqn);
01977      if (isSlave) ERReturn(-1);
01978      workPtr[j] = reducedEqn;
01979    }
01980 
01981    for(int i=0; i<numIndices; i++) {
01982       int row = indPtr[i];
01983       int rrow = workPtr[i];
01984 
01985       if ((localStartRow_ > row) || (row > localEndRow_)) {
01986   
01987   int owner = getOwnerProcForEqn(row);
01988   eqnCommMgr_->addRemoteIndices(rrow, owner, workPtr, numIndices);
01989   continue;
01990       }
01991 
01992       CHK_ERR( createMatrixPositions(rrow, numIndices, workPtr,
01993              "storeElemScttrInd") );
01994    }
01995 
01996    return(FEI_SUCCESS);
01997 }
01998 
01999 //------------------------------------------------------------------------------
02000 int SNL_FEI_Structure::storeElementScatterIndices_noSlaves(std::vector<int>& indices)
02001 {
02002   //This function takes a list of equation numbers as input and for each one,
02003   //if it is a local equation number, goes to that row of the sysMatIndices
02004   //structure and stores the whole list of equation numbers as column indices
02005   //in that row of the matrix structure. If the equation number is not local,
02006   //then it is given to the EqnCommMgr.
02007   //
02008   int i, numIndices = indices.size();
02009   int* indPtr = &indices[0];
02010 
02011   for(i=0; i<numIndices; i++) {
02012     int row = indPtr[i];
02013     if (row < localStartRow_ || row > localEndRow_) {
02014       int owner = getOwnerProcForEqn(row);
02015       eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
02016       continue;
02017     }
02018 
02019     if (generateGraph_) {
02020       fei::ctg_set<int>& thisRow =
02021   sysMatIndices_[row - reducedStartRow_];
02022 
02023       for(int j=0; j<numIndices; ++j) {
02024   if (debugOutput_ && outputLevel_ > 2) {
02025     dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
02026        << row << "," << indPtr[j] << ")"<<FEI_ENDL;
02027   }
02028 
02029   thisRow.insert2(indPtr[j]);
02030       }
02031     }
02032   }
02033 
02034   return(FEI_SUCCESS);
02035 }
02036 
02037 //------------------------------------------------------------------------------
02038 int SNL_FEI_Structure::storeElementScatterBlkIndices_noSlaves(std::vector<int>& indices)
02039 {
02040   //This function takes a list of equation numbers as input and for each one,
02041   //if it is a local equation number, goes to that row of the sysMatIndices
02042   //structure and stores the whole list of equation numbers as column indices
02043   //in that row of the matrix structure. If the equation number is not local,
02044   //then it is given to the EqnCommMgr.
02045   //
02046   int i, numIndices = indices.size();
02047   int* indPtr = &indices[0];
02048 
02049   for(i=0; i<numIndices; i++) {
02050     int row = indPtr[i];
02051     if (row < localStartRow_ || row > localEndRow_) {
02052       int owner = getOwnerProcForEqn(row);
02053       eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
02054       continue;
02055     }
02056 
02057     if (generateGraph_) {
02058       fei::ctg_set<int>& thisRow =
02059   sysMatIndices_[row - reducedStartRow_];
02060 
02061       for(int j=0; j<numIndices; ++j) {
02062   if (debugOutput_ && outputLevel_ > 2) {
02063     dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
02064        << row << "," << indPtr[j] << ")"<<FEI_ENDL;
02065   }
02066 
02067   thisRow.insert2(indPtr[j]);
02068       }
02069     }
02070   }
02071 
02072   return(FEI_SUCCESS);
02073 }
02074 
02075 //------------------------------------------------------------------------------
02076 int SNL_FEI_Structure::assembleReducedStructure()
02077 {
02078   //This function performs the appropriate manipulations (matrix-matrix-products,
02079   //matrix-vector-products) on the Kid,Kdi,Kdd matrices and then assembles the
02080   //results into the global matrix structure. This function also adds any
02081   //resulting contributions to off-processor portions of the matrix structure to
02082   //the EqnCommMgr.
02083   //
02084   fei::FillableMat* D = getSlaveDependencies();
02085 
02086   csrD = *D;
02087   csrKid = *Kid_;
02088   csrKdi = *Kdi_;
02089   csrKdd = *Kdd_;
02090 
02091   //form tmpMat1_ = Kid_*D
02092   fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
02093 
02094   //form tmpMat2_ = D^T*Kdi_
02095   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
02096 
02097   CHK_ERR( translateMatToReducedEqns(tmpMat1_) );
02098   CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
02099   if (numProcs_ > 1) {
02100     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat1_, true) );
02101     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
02102   }
02103 
02104   CHK_ERR( createMatrixPositions(tmpMat1_) );
02105   CHK_ERR( createMatrixPositions(tmpMat2_) );
02106 
02107   //form tmpMat1_ = D^T*Kdd_
02108   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
02109 
02110   //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
02111   fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
02112 
02113 
02114   CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
02115   if (numProcs_ > 1) {
02116     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
02117   }
02118   CHK_ERR( createMatrixPositions(tmpMat2_) );
02119 
02120   Kdi_->clear();
02121   Kid_->clear();
02122   Kdd_->clear();
02123   reducedEqnCounter_ = 0;
02124 
02125   return(FEI_SUCCESS);
02126 }
02127 
02128 //------------------------------------------------------------------------------
02129 int SNL_FEI_Structure::translateToReducedEqns(EqnCommMgr& eqnCommMgr)
02130 {
02131   if (numSlvs_ < 1) return(0);
02132 
02133   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvProcEqns())) );
02134   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvEqns())) );
02135   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendProcEqns())) );
02136   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendEqns())) );
02137 
02138   return(0);
02139 }
02140 
02141 //------------------------------------------------------------------------------
02142 int SNL_FEI_Structure::translateToReducedEqns(EqnBuffer& eqnBuf)
02143 {
02144   int numEqns = eqnBuf.getNumEqns();
02145   int* eqnNumbers = &(eqnBuf.eqnNumbers()[0]);
02146   std::vector<fei::CSVec*>& eqnArray = eqnBuf.eqns();
02147   for(int i=0; i<numEqns; ++i) {
02148     int reducedEqn;
02149     translateToReducedEqn(eqnNumbers[i], reducedEqn);
02150     eqnNumbers[i] = reducedEqn;
02151 
02152     int* indicesPtr = &(eqnArray[i]->indices()[0]);
02153     int numIndices = eqnArray[i]->size();
02154     for(int j=0; j<numIndices; ++j) {
02155       translateToReducedEqn(indicesPtr[j], reducedEqn);
02156       indicesPtr[j] = reducedEqn;
02157     }
02158   }
02159 
02160   return(0);
02161 }
02162 
02163 //------------------------------------------------------------------------------
02164 int SNL_FEI_Structure::translateToReducedEqns(ProcEqns& procEqns)
02165 {
02166   std::vector<std::vector<int>*>& eqnTable = procEqns.procEqnNumbersPtr();
02167   int dim1 = eqnTable.size();
02168   for(int i=0; i<dim1; ++i) {
02169     int dim2 = eqnTable[i]->size();
02170     int* indicesPtr = &(*(eqnTable[i]))[0];
02171     for(int j=0; j<dim2; ++j) {
02172       int reducedEqn;
02173       translateToReducedEqn(indicesPtr[j], reducedEqn);
02174       indicesPtr[j] = reducedEqn;
02175     }
02176   }
02177 
02178   return(0);
02179 }
02180 
02181 //------------------------------------------------------------------------------
02182 int SNL_FEI_Structure::translateMatToReducedEqns(fei::CSRMat& mat)
02183 {
02184   //Given a matrix in unreduced numbering, convert all of its contents to the
02185   //"reduced" equation space. If any of the row-numbers or column-indices in
02186   //this matrix object are slave equations, they will not be referenced. In
02187   //this case, a positive warning-code will be returned.
02188 
02189   int reducedEqn = -1;
02190   int foundSlave = 0;
02191 
02192   fei::SparseRowGraph& srg = mat.getGraph();
02193 
02194   std::vector<int>& rowNumbers = srg.rowNumbers;
02195   for(size_t i=0; i<rowNumbers.size(); ++i) {
02196     bool isSlave = translateToReducedEqn(rowNumbers[i], reducedEqn);
02197     if (isSlave) foundSlave = 1;
02198     else rowNumbers[i] = reducedEqn;
02199   }
02200 
02201   std::vector<int>& colIndices = srg.packedColumnIndices;
02202   for(size_t i=0; i<colIndices.size(); ++i) {
02203     bool isSlave = translateToReducedEqn(colIndices[i], reducedEqn);
02204     if (isSlave) foundSlave = 1;
02205     else colIndices[i] = reducedEqn;
02206   }
02207 
02208   return(foundSlave);
02209 }
02210 
02211 //------------------------------------------------------------------------------
02212 int SNL_FEI_Structure::createMatrixPositions(fei::CSRMat& mat)
02213 {
02214   if (!generateGraph_) {
02215     return(0);
02216   }
02217 
02218   //This function must be called with a CSRMat object that already has its
02219   //contents numbered in "reduced" equation-numbers.
02220   //
02221   //This function has one simple task -- for each row,col pair stored in 'mat',
02222   //call 'createMatrixPosition' to make an entry in the global matrix structure
02223   //if it doesn't already exist.
02224   //
02225   const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
02226   const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
02227   const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
02228 
02229   for(size_t i=0; i<rowNumbers.size(); ++i) {
02230     int row = rowNumbers[i];
02231     int offset = rowOffsets[i];
02232     int rowlen = rowOffsets[i+1]-offset;
02233     const int* indices = &pckColInds[offset];
02234 
02235     for(int j=0; j<rowlen; j++) {
02236       CHK_ERR( createMatrixPosition(row, indices[j], "crtMatPos(CSRMat)") );
02237     }
02238   }
02239 
02240   return(0);
02241 }
02242 
02243 //------------------------------------------------------------------------------
02244 int SNL_FEI_Structure::createMatrixPosition(int row, int col,
02245               const char* callingFunction)
02246 {
02247   if (!generateGraph_) {
02248     return(0);
02249   }
02250 
02251   //
02252   //This function inserts the column index 'col' in the system matrix structure
02253   //in 'row', if it isn't already there.
02254   //
02255   //This is a private SNL_FEI_Structure function. These row/col numbers must be
02256   //global, 0-based equation numbers in the "reduced" equation space..
02257   //
02258 
02259   //if the equation isn't local, just return. (The assumption being that
02260   //this call is also being made on the processor that owns the equation.)
02261   if (row < reducedStartRow_ || row > reducedEndRow_) {
02262     if (debugOutput_) {
02263       dbgOut() << "         " << callingFunction << " passed " <<row<<","<<col
02264          << ", not local." << FEI_ENDL;
02265     }
02266 
02267     return(0);
02268   }
02269 
02270   if (debugOutput_ && outputLevel_ > 2) {
02271     dbgOut() << "#         " << callingFunction << " crtMatPos("
02272        << row << "," << col << ")"<<FEI_ENDL;
02273   }
02274 
02275   fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
02276 
02277   thisRow.insert2(col);
02278 
02279   return(0);
02280 }
02281 
02282 //------------------------------------------------------------------------------
02283 int SNL_FEI_Structure::createMatrixPositions(int row, int numCols, int* cols,
02284               const char* callingFunction)
02285 {
02286   if (!generateGraph_) {
02287     return(0);
02288   }
02289 
02290   //
02291   //This function inserts the column indices 'cols' in the system matrix structure
02292   //in 'row', if they aren't already there.
02293   //
02294   //This is a private SNL_FEI_Structure function. These row/col numbers must be
02295   //global, 0-based equation numbers in the "reduced" equation space..
02296   //
02297 
02298   //if the equation isn't local, just return. (The assumption being that
02299   //this call is also being made on the processor that owns the equation.)
02300   if (row < reducedStartRow_ || row > reducedEndRow_) {
02301     if (debugOutput_) {
02302       dbgOut() << "         " << callingFunction << " passed " <<row
02303          << ", not local." << FEI_ENDL;
02304     }
02305 
02306     return(0);
02307   }
02308 
02309   fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
02310 
02311   for(int i=0; i<numCols; i++) {
02312     if (debugOutput_ && outputLevel_ > 2) {
02313       dbgOut() << "#         " << callingFunction << " crtMatPoss("
02314          << row << "," << cols[i] << ")"<<FEI_ENDL;
02315     }
02316 
02317     thisRow.insert2(cols[i]);
02318   }
02319 
02320   return(0);
02321 }
02322 
02323 //------------------------------------------------------------------------------
02324 int SNL_FEI_Structure::initializeBlkEqnMapper()
02325 {
02326   //Run the nodes, elem-dof's and multiplier constraints, establishing the
02327   //point-equation to block-equation mappings.  Note that this data is 
02328   //initialized as non-reduced (i.e., not accounting for any slave equations)
02329   //equation numbers, then translated to the reduced space at the end of this
02330   //function.
02331 
02332   if (blkEqnMapper_ != NULL) delete blkEqnMapper_;
02333   blkEqnMapper_ = new snl_fei::PointBlockMap();
02334   snl_fei::PointBlockMap& blkEqnMapper = getBlkEqnMapper();
02335 
02336   //First, if there's only one (non-negative) fieldID and that fieldID's size
02337   //is 1, then we don't need to execute the rest of this function.
02338   int numFields = fieldDatabase_->size();
02339   const int* fieldIDs = getFieldIDsPtr();
02340   const int* fieldSizes = fieldIDs+numFields;
02341   int numNonNegativeFields = 0;
02342   int maxFieldSize = 0;
02343   for(int f=0; f<numFields; f++) {
02344     if (fieldIDs[f] >= 0) {
02345       numNonNegativeFields++;
02346       if (fieldSizes[f] > maxFieldSize) maxFieldSize = fieldSizes[f];
02347     }
02348   }
02349 
02350   if (numNonNegativeFields == 1 && maxFieldSize == 1) {
02351     globalMaxBlkSize_ = 1;
02352     blkEqnMapper.setPtEqualBlk();
02353     return(0);
02354   }
02355 
02356   int localVanishedNodeAdjustment = 0;
02357   for(int i=0; i<localProc_; ++i) {
02358     localVanishedNodeAdjustment += globalNumNodesVanished_[i];
02359   }
02360 
02361   int eqnNumber, blkEqnNumber;
02362 
02363   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02364   std::map<GlobalID,int>::iterator
02365     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02366 
02367   for(; iter!=iter_end; ++iter) {
02368     NodeDescriptor* node = NULL;
02369     nodeDatabase_->getNodeAtIndex(iter->second, node);
02370 
02371     // If the node doesn't exist, skip.
02372     if (node==NULL) continue;
02373 
02374     // If the node doesn't have any fields, skip
02375     int numFields  = node->getNumFields();
02376     if(numFields == 0) continue;
02377 
02378     int firstPtEqn = node->getFieldEqnNumbers()[0];
02379     int numNodalDOF = node->getNumNodalDOF();
02380     blkEqnNumber = node->getBlkEqnNumber() - localVanishedNodeAdjustment;
02381 
02382     int blkSize = numNodalDOF;
02383 
02384     for(int j=0; j<numNodalDOF; ++j) {
02385       bool isSlave = translateToReducedEqn(firstPtEqn+j, eqnNumber);
02386       if (isSlave) {
02387   --blkSize;
02388       }
02389       else {
02390   CHK_ERR( blkEqnMapper.setEqn(eqnNumber, blkEqnNumber) );
02391       }
02392     }
02393 
02394     if (blkSize > 0) {
02395       CHK_ERR( blkEqnMapper.setBlkEqnSize(blkEqnNumber, blkSize) );
02396     }
02397     else {
02398       ++localVanishedNodeAdjustment;
02399     }
02400   }
02401 
02402   //
02403   //Now the element-dofs...
02404   //
02405   int numBlocks = getNumElemBlocks();
02406   int numLocalNodes = globalNodeOffsets_[localProc_+1]-globalNodeOffsets_[localProc_];
02407   eqnNumber = localStartRow_ + numLocalNodalEqns_;
02408   blkEqnNumber = localBlkOffset_ + numLocalNodes;
02409 
02410   for(int i = 0; i < numBlocks; i++) {
02411     BlockDescriptor* block = NULL;
02412     CHK_ERR( getBlockDescriptor_index(i, block) );
02413 
02414     int numElemDOF = block->getNumElemDOFPerElement();
02415 
02416     if (numElemDOF > 0) {
02417       int numElems = block->getNumElements();
02418 
02419       for(int j=0; j<numElems; j++) {
02420 
02421   for(int ede=0; ede<numElemDOF; ede++) {
02422     blkEqnMapper.setEqn(eqnNumber+ede, blkEqnNumber+ede);
02423     blkEqnMapper.setBlkEqnSize(blkEqnNumber+ede, 1);
02424   }
02425 
02426   eqnNumber += numElemDOF;
02427   blkEqnNumber += numElemDOF;
02428       }
02429     }
02430   }
02431 
02432   //Now we'll set mappings for all local multiplier constraint-relations,
02433   //and also gather the equation numbers for other processors' multipliers
02434   //to record in our snl_fei::PointBlockMap object.
02435   //
02436   std::map<GlobalID,ConstraintType*>::const_iterator
02437     cr_iter = multCRs_.begin(),
02438     cr_end  = multCRs_.end();
02439 
02440   std::vector<int> localMultEqns;
02441   localMultEqns.reserve(multCRs_.size()*2);
02442   while(cr_iter != cr_end) {
02443     ConstraintType* multCR = (*cr_iter).second;
02444     eqnNumber = multCR->getEqnNumber();
02445     blkEqnNumber = multCR->getBlkEqnNumber();
02446 
02447     CHK_ERR( blkEqnMapper_->setEqn(eqnNumber, blkEqnNumber) );
02448     CHK_ERR( blkEqnMapper_->setBlkEqnSize(blkEqnNumber, 1) );
02449 
02450     localMultEqns.push_back(eqnNumber);
02451     localMultEqns.push_back(blkEqnNumber);
02452     ++cr_iter;
02453   }
02454 
02455   //Now gather and store the equation-numbers arising from other
02456   //processors' multiplier constraints. (We obviously only need to do this if
02457   //there is more than one processor, and we can skip it if the problem
02458   //only contains 1 scalar equation for each block-equation.)
02459 
02460   int localMaxBlkEqnSize = blkEqnMapper_->getMaxBlkEqnSize();
02461   globalMaxBlkSize_ = 0;
02462   CHK_ERR( fei::GlobalMax(comm_, localMaxBlkEqnSize, globalMaxBlkSize_) );
02463 
02464   blkEqnMapper_->setMaxBlkEqnSize(globalMaxBlkSize_);
02465 
02466   if (globalMaxBlkSize_ == 1) {
02467     //If the globalMax block-size is 1 then tell the blkEqnMapper, and it will
02468     //reclaim the memory it allocated in arrays, etc.
02469     blkEqnMapper_->setPtEqualBlk();
02470     return(0);
02471   }
02472 
02473   if (numProcs_ > 1) {
02474     std::vector<int> recvLengths, globalMultEqns;
02475     CHK_ERR(fei::Allgatherv(comm_, localMultEqns, recvLengths, globalMultEqns));
02476 
02477     int offset = 0;
02478     for(int p=0; p<numProcs_; p++) {
02479       if (p == localProc_) { offset += recvLengths[p]; continue; }
02480 
02481       for(int j=0; j<recvLengths[p]/2; j++) {
02482   blkEqnMapper_->setEqn(globalMultEqns[offset], globalMultEqns[offset+1]);
02483   blkEqnMapper_->setBlkEqnSize(globalMultEqns[offset+1], 1);
02484   offset += 2;
02485       }
02486     }
02487   }
02488 
02489   return(0);
02490 }
02491 
02492 //------------------------------------------------------------------------------
02493 int SNL_FEI_Structure::setMultCREqnInfo()
02494 {
02495   //We'll set equation-numbers on all local multiplier constraint-relations,
02496   //and also gather the equation numbers for other processors' multipliers
02497   //to record in our snl_fei::PointBlockMap object.
02498   //
02499   std::map<GlobalID,ConstraintType*>::const_iterator
02500     cr_iter = multCRs_.begin(),
02501     cr_end  = multCRs_.end();
02502 
02503   int eqnNumber = localStartRow_ + numLocalNodalEqns_ + numLocalElemDOF_;
02504   int blkEqnNum = localBlkOffset_ + numLocalEqnBlks_ - numLocalMultCRs_;
02505 
02506   while(cr_iter != cr_end) {
02507     ConstraintType* multCR = (*cr_iter).second;
02508 
02509     multCR->setEqnNumber(eqnNumber);
02510 
02511     multCR->setBlkEqnNumber(blkEqnNum);
02512 
02513     ++eqnNumber;
02514     ++blkEqnNum;
02515     ++cr_iter;
02516   }
02517 
02518   return(0);
02519 }
02520 
02521 //------------------------------------------------------------------------------
02522 int SNL_FEI_Structure::writeEqn2NodeMap()
02523 {
02524   FEI_OSTRINGSTREAM osstr;
02525   osstr << dbgPath_ << "/FEI" << name_ << "_nID_nNum_blkEq_ptEq."
02526   << numProcs_ << "." << localProc_;
02527 
02528   FEI_OFSTREAM e2nFile(osstr.str().c_str());
02529 
02530   if (!e2nFile) {
02531     fei::console_out() << "SNL_FEI_Structure::writeEqn2NodeMap: ERROR, couldn't open file "
02532          << osstr.str() << FEI_ENDL;
02533     ERReturn(-1);
02534   }
02535 
02536   std::map<GlobalID,ConstraintType*>::const_iterator
02537     cr_iter = multCRs_.begin(),
02538     cr_end  = multCRs_.end();
02539 
02540   while(cr_iter != cr_end) {
02541     ConstraintType* multCR = (*cr_iter).second;
02542     int eqnNumber = multCR->getEqnNumber();
02543     int blkEqnNumber = multCR->getBlkEqnNumber();
02544 
02545     e2nFile << "-999 -999 " << blkEqnNumber << " " << eqnNumber << FEI_ENDL;
02546     ++cr_iter;
02547   }
02548 
02549   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02550   std::map<GlobalID,int>::iterator
02551     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02552 
02553   for(; iter!=iter_end; ++iter) {
02554     NodeDescriptor* node = NULL;
02555     nodeDatabase_->getNodeAtIndex(iter->second, node);
02556 
02557     if (node==NULL || node->getOwnerProc() != localProc_) {
02558       continue;
02559     }
02560 
02561     const int* fieldIDList = node->getFieldIDList();
02562     const int* eqnNumbers = node->getFieldEqnNumbers();
02563     int nodeNum = node->getNodeNumber();
02564     int blkEqnNumber = node->getBlkEqnNumber();
02565     int numFields = node->getNumFields();
02566 
02567     for(int j=0; j<numFields; j++) {
02568       int numFieldParams = getFieldSize(fieldIDList[j]);
02569       assert( numFieldParams > 0 );
02570 
02571       for(int k=0; k<numFieldParams; k++) {
02572         e2nFile << (int)node->getGlobalNodeID() << " " << nodeNum
02573                 << " " << blkEqnNumber << " " << eqnNumbers[j]+k<<FEI_ENDL;
02574       }
02575     }
02576   }
02577 
02578   return(FEI_SUCCESS);
02579 }
02580 
02581 //------------------------------------------------------------------------------
02582 int SNL_FEI_Structure::calcTotalNumElemDOF()
02583 {
02584    int numBlocks = getNumElemBlocks();
02585    int totalNumElemDOF = 0;
02586 
02587    for(int i = 0; i < numBlocks; i++) {
02588      BlockDescriptor* block = NULL;
02589      CHK_ERR( getBlockDescriptor_index(i, block) );
02590       int numElemDOF = block->getNumElemDOFPerElement();
02591       if (numElemDOF > 0) {
02592          int numElems = block->getNumElements();
02593 
02594          totalNumElemDOF += numElems*numElemDOF;
02595       }
02596    }
02597 
02598    return(totalNumElemDOF);
02599 }
02600 
02601 //------------------------------------------------------------------------------
02602 int SNL_FEI_Structure::calcNumMultCREqns()
02603 {
02604   int totalNumMultCRs = 0;
02605   int numMCRecords = getNumMultConstRecords();
02606 
02607   totalNumMultCRs = numMCRecords;
02608 
02609   return( totalNumMultCRs );
02610 }
02611 
02612 //------------------------------------------------------------------------------
02613 void SNL_FEI_Structure::calcGlobalEqnInfo(int numLocallyOwnedNodes, 
02614            int numLocalEqns,
02615            int numLocalEqnBlks) 
02616 {
02617   FEI_OSTREAM& os = dbgOut();
02618   if (debugOutput_) os << "#  entering calcGlobalEqnInfo" << FEI_ENDL;
02619 
02620 //
02621 //This function does the inter-process communication necessary to obtain,
02622 //on each processor, the global number of equations, and the local starting
02623 //and ending equation numbers.
02624 //While we're here, we do the same thing for node-numbers.
02625 //
02626 
02627   //first, get each processor's local number of equations on the master proc.
02628 
02629   std::vector<int> numProcNodes(numProcs_);
02630   std::vector<int> numProcEqns(numProcs_);
02631   std::vector<int> numProcEqnBlks(numProcs_);
02632 
02633   if (numProcs_ > 1) {
02634 #ifndef FEI_SER
02635     int glist[3];
02636     std::vector<int> recvList(3*numProcs_);
02637 
02638     glist[0] = numLocallyOwnedNodes;
02639     glist[1] = numLocalEqns;
02640     glist[2] = numLocalEqnBlks;
02641     if (MPI_Gather(&glist[0], 3, MPI_INT, &recvList[0], 3,
02642        MPI_INT, masterProc_, comm_) != MPI_SUCCESS) {
02643       fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Gather" << FEI_ENDL;
02644       MPI_Abort(comm_, -1);
02645     }
02646 
02647     for(int p=0; p<numProcs_; p++) {
02648       numProcNodes[p]   = recvList[p*3];
02649       numProcEqns[p]    = recvList[p*3+1];
02650       numProcEqnBlks[p] = recvList[p*3+2];
02651     }
02652 #endif
02653   }
02654   else {
02655     numProcNodes[0] = numLocallyOwnedNodes;
02656     numProcEqns[0] = numLocalEqns;
02657     numProcEqnBlks[0] = numLocalEqnBlks;
02658   }
02659 
02660   //compute offsets for all processors (starting index for local equations)
02661 
02662   globalNodeOffsets_.resize(numProcs_+1);
02663   globalEqnOffsets_.resize(numProcs_+1);
02664   globalBlkEqnOffsets_.resize(numProcs_+1);
02665 
02666   globalNodeOffsets_[0] = 0; // 0-based node-numbers
02667   globalEqnOffsets_[0] = 0;    // we're going to start rows & cols at 0 (global)
02668   globalBlkEqnOffsets_[0] = 0;
02669 
02670   if (localProc_ == masterProc_) {
02671     for(int i=1;i<numProcs_;i++) {
02672       globalNodeOffsets_[i] = globalNodeOffsets_[i-1] +numProcNodes[i-1];
02673       globalEqnOffsets_[i] = globalEqnOffsets_[i-1] + numProcEqns[i-1];
02674       globalBlkEqnOffsets_[i] = globalBlkEqnOffsets_[i-1] +
02675                                   numProcEqnBlks[i-1];
02676     }
02677 
02678     globalNodeOffsets_[numProcs_] = globalNodeOffsets_[numProcs_-1] +
02679                                    numProcNodes[numProcs_-1];
02680     globalEqnOffsets_[numProcs_] = globalEqnOffsets_[numProcs_-1] +
02681                                   numProcEqns[numProcs_-1];
02682     globalBlkEqnOffsets_[numProcs_] = globalBlkEqnOffsets_[numProcs_-1] +
02683                                      numProcEqnBlks[numProcs_-1];
02684   }
02685 
02686   //now, broadcast vector of eqnOffsets to all processors
02687 
02688   if (numProcs_ > 1) {
02689 #ifndef FEI_SER
02690     std::vector<int> blist(3*(numProcs_+1));
02691 
02692     if (localProc_ == masterProc_) {
02693       int offset = 0;
02694       for(int i=0; i<numProcs_+1; i++) {
02695   blist[offset++] = globalNodeOffsets_[i];
02696   blist[offset++] = globalEqnOffsets_[i];
02697   blist[offset++] = globalBlkEqnOffsets_[i];
02698       }
02699     }
02700 
02701     if (MPI_Bcast(&blist[0], 3*(numProcs_+1), MPI_INT,
02702       masterProc_, comm_) != MPI_SUCCESS) {
02703       fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Bcast" << FEI_ENDL;
02704       MPI_Abort(comm_, -1);
02705     }
02706 
02707     if (localProc_ != masterProc_) {
02708       int offset = 0;
02709       for(int i=0; i<numProcs_+1; i++) {
02710   globalNodeOffsets_[i] = blist[offset++];
02711   globalEqnOffsets_[i] = blist[offset++];
02712   globalBlkEqnOffsets_[i] = blist[offset++];
02713       }
02714     }
02715 #endif
02716   }
02717 
02718   localStartRow_ = globalEqnOffsets_[localProc_];
02719   localEndRow_ = globalEqnOffsets_[localProc_+1]-1;
02720   numGlobalEqns_ = globalEqnOffsets_[numProcs_];
02721 
02722   firstLocalNodeNumber_ = globalNodeOffsets_[localProc_];
02723   numGlobalNodes_ = globalNodeOffsets_[numProcs_];
02724 
02725   localBlkOffset_ = globalBlkEqnOffsets_[localProc_];
02726   numGlobalEqnBlks_ = globalBlkEqnOffsets_[numProcs_];
02727 
02728   if (debugOutput_) {
02729     os << "#     localStartRow_: " << localStartRow_ << FEI_ENDL;
02730     os << "#     localEndRow_: "   << localEndRow_ << FEI_ENDL;
02731     os << "#     numGlobalEqns_: " << numGlobalEqns_ << FEI_ENDL;
02732     os << "#     numGlobalEqnBlks_: " << numGlobalEqnBlks_ << FEI_ENDL;
02733     os << "#     localBlkOffset_: " << localBlkOffset_ << FEI_ENDL;
02734     os << "#     firstLocalNodeNumber_: " << firstLocalNodeNumber_ << FEI_ENDL;
02735     size_t n;
02736     os << "#     numGlobalNodes_: " << numGlobalNodes_ << FEI_ENDL;
02737     os << "#     " << globalNodeOffsets_.size() << " globalNodeOffsets_: ";
02738     for(size_t nn=0; nn<globalNodeOffsets_.size(); nn++) 
02739       os << globalNodeOffsets_[nn] << " ";
02740     os << FEI_ENDL;
02741     os << "#     " << globalEqnOffsets_.size() << " globalEqnOffsets_: ";
02742     for(n=0; n<globalEqnOffsets_.size(); n++) 
02743       os << globalEqnOffsets_[n] << " ";
02744     os << FEI_ENDL;
02745     os << "#     " << globalBlkEqnOffsets_.size() << " globalBlkEqnOffsets_: ";
02746     for(n=0; n<globalBlkEqnOffsets_.size(); n++) 
02747       os << globalBlkEqnOffsets_[n] << " ";
02748     os << FEI_ENDL;
02749     os << "#  leaving calcGlobalEqnInfo" << FEI_ENDL;
02750   }
02751 }
02752 
02753 //------------------------------------------------------------------------------
02754 int SNL_FEI_Structure::setNodalEqnInfo()
02755 {
02756 //
02757 //Private SNL_FEI_Structure function, to be called in initComplete, after
02758 //'calcGlobalEqnInfo()'.
02759 //
02760 //This function sets global equation numbers on the nodes.
02761 //
02762 //The return value is an error-code.
02763 //
02764   int eqn = localStartRow_;
02765 
02766   //localBlkOffset_ is 0-based, and so is blkEqnNumber.
02767   int blkEqnNumber = localBlkOffset_;
02768 
02769   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02770   std::map<GlobalID,int>::iterator
02771     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02772 
02773   for(; iter!=iter_end; ++iter) {
02774     NodeDescriptor* node = NULL;
02775     nodeDatabase_->getNodeAtIndex(iter->second, node);
02776 
02777     if (node==NULL) continue;
02778 
02779     int numFields = node->getNumFields();
02780     const int* fieldIDs = node->getFieldIDList();
02781 
02782     if (node->getOwnerProc() == localProc_) {
02783       int numNodalDOF = 0;
02784 
02785       for(int j=0; j<numFields; j++) {
02786   node->setFieldEqnNumber(fieldIDs[j], eqn);
02787   int fieldSize = getFieldSize(fieldIDs[j]);
02788   eqn += fieldSize;
02789   numNodalDOF += fieldSize;
02790       }
02791 
02792       node->setNumNodalDOF(numNodalDOF);
02793       node->setBlkEqnNumber(blkEqnNumber++);
02794     }
02795   }
02796 
02797   return(0);
02798 }
02799 
02800 //------------------------------------------------------------------------------
02801 void SNL_FEI_Structure::setElemDOFEqnInfo()
02802 {
02803   //
02804   //Private SNL_FEI_Structure function, to be called from initComplete after
02805   //setBlkEqnInfo() has been called.
02806   //
02807   //Sets global equation numbers for all element-DOF.
02808   //
02809   int numBlocks = getNumElemBlocks();
02810 
02811   int eqnNumber = localStartRow_ + numLocalNodalEqns_;
02812 
02813   for(int i = 0; i < numBlocks; i++) {
02814     BlockDescriptor* block = NULL;
02815     int err = getBlockDescriptor_index(i, block);
02816     if (err) voidERReturn;
02817 
02818     int numElemDOF = block->getNumElemDOFPerElement();
02819 
02820     if (numElemDOF > 0) {
02821       std::vector<int>& elemDOFEqns = block->elemDOFEqnNumbers();
02822       std::map<GlobalID,int>& elemIDs = connTables_[i]->elemIDs;
02823 
02824       std::map<GlobalID,int>::const_iterator
02825         iter = elemIDs.begin(),
02826         iter_end = elemIDs.end();
02827 
02828       for(; iter != iter_end; ++iter) {
02829   elemDOFEqns[iter->second] = eqnNumber;
02830 
02831   eqnNumber += numElemDOF;
02832       }
02833     }
02834   }
02835 }
02836 
02837 //------------------------------------------------------------------------------
02838 int SNL_FEI_Structure::addBlock(GlobalID blockID) {
02839 //
02840 //Append a blockID to our (unsorted) list of blockIDs if it isn't
02841 //already present. Also, add a block-descriptor if blockID wasn't
02842 //already present, and a ConnectivityTable to go with it.
02843 //
02844   int insertPoint = -1;
02845   int found = fei::binarySearch(blockID, blockIDs_, insertPoint);
02846 
02847    if (found<0) {
02848       blockIDs_.insert(blockIDs_.begin()+insertPoint, blockID);
02849 
02850       BlockDescriptor* block = new BlockDescriptor;
02851       block->setGlobalBlockID(blockID);
02852 
02853       blocks_.insert(blocks_.begin()+insertPoint, block);
02854 
02855       ConnectivityTable* newConnTable = new ConnectivityTable;
02856       connTables_.insert(connTables_.begin()+insertPoint, newConnTable);
02857    }
02858 
02859    return(FEI_SUCCESS);
02860 }
02861 
02862 //------------------------------------------------------------------------------
02863 int SNL_FEI_Structure::getBlockDescriptor(GlobalID blockID,
02864                                          BlockDescriptor*& block)
02865 {
02866   int index = fei::binarySearch(blockID, blockIDs_);
02867 
02868    if (index < 0) {
02869       fei::console_out() << "SNL_FEI_Structure::getBlockDescriptor: ERROR, blockID "
02870            << (int)blockID << " not found." << FEI_ENDL;
02871       return(-1);
02872    }
02873 
02874    block = blocks_[index];
02875 
02876    return(FEI_SUCCESS);
02877 }
02878 
02879 //------------------------------------------------------------------------------
02880 int SNL_FEI_Structure::getIndexOfBlock(GlobalID blockID) const
02881 {
02882   int index = fei::binarySearch(blockID, blockIDs_);
02883   return(index);
02884 }
02885 
02886 //------------------------------------------------------------------------------
02887 int SNL_FEI_Structure::getBlockDescriptor_index(int index,
02888                  BlockDescriptor*& block)
02889 {
02890   if (index < 0 || index >= (int)blockIDs_.size()) ERReturn(FEI_FATAL_ERROR);
02891 
02892   block = blocks_[index];
02893 
02894   return(FEI_SUCCESS);
02895 }
02896 
02897 //------------------------------------------------------------------------------
02898 int SNL_FEI_Structure::allocateBlockConnectivity(GlobalID blockID) {
02899 
02900    int index = fei::binarySearch(blockID, blockIDs_);
02901 
02902    if (index < 0) {
02903       fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, blockID "
02904            << (int)blockID << " not found. Aborting." << FEI_ENDL;
02905       MPI_Abort(comm_, -1);
02906    }
02907 
02908    connTables_[index]->numNodesPerElem = 
02909                          blocks_[index]->getNumNodesPerElement();
02910 
02911    int numRows = blocks_[index]->getNumElements();
02912    int numCols = connTables_[index]->numNodesPerElem;
02913 
02914    if ((numRows <= 0) || (numCols <= 0)) {
02915       fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, either "
02916            << "numElems or numNodesPerElem not yet set for blockID "
02917            << (int)blockID << ". Aborting." << FEI_ENDL;
02918       MPI_Abort(comm_, -1);
02919    }
02920 
02921    connTables_[index]->elemNumbers.resize(numRows);
02922 
02923    connTables_[index]->elem_conn_ids = new std::vector<GlobalID>(numRows*numCols);
02924 
02925    return(FEI_SUCCESS);
02926 }
02927 
02928 //------------------------------------------------------------------------------
02929 ConnectivityTable& SNL_FEI_Structure::getBlockConnectivity(GlobalID blockID) {
02930 
02931    int index = fei::binarySearch(blockID, blockIDs_);
02932 
02933    if (index < 0) {
02934       fei::console_out() << "SNL_FEI_Structure::getBlockConnectivity: ERROR, blockID "
02935            << (int)blockID << " not found. Aborting." << FEI_ENDL;
02936       MPI_Abort(comm_, -1);
02937    }  
02938 
02939    return(*(connTables_[index]));
02940 }
02941 
02942 //------------------------------------------------------------------------------
02943 int SNL_FEI_Structure::finalizeActiveNodes()
02944 {
02945   FEI_OSTREAM& os = dbgOut();
02946   if (debugOutput_) {
02947     os << "#  finalizeActiveNodes" << FEI_ENDL;
02948   }
02949   //First, make sure that the shared-node IDs are in the nodeDatabase, since
02950   //they might not have been added yet...
02951 
02952   std::vector<GlobalID>& shNodeIDs = nodeCommMgr_->getSharedNodeIDs();
02953   for(unsigned i=0; i<shNodeIDs.size(); i++) {
02954     CHK_ERR( nodeDatabase_->initNodeID(shNodeIDs[i]) );
02955   }
02956 
02957   if (activeNodesInitialized_) return(0);
02958   //
02959   //Now run through the connectivity tables and set the nodal field and block
02960   //info. Also, replace the nodeID in the element-connectivity lists with an 
02961   //index from the NodeDatabase object. This will save us a lot of time when
02962   //looking up nodes later.
02963   //
02964   //One other thing we're going to do is give all elements an element-number.
02965   //This element-number will start at zero locally (meaning that it's not
02966   //globally unique).
02967   //
02968   int elemNumber = 0;
02969   int numBlocks = blockIDs_.size();
02970   for(int bIndex=0; bIndex<numBlocks; bIndex++) {
02971 
02972     GlobalID blockID = blocks_[bIndex]->getGlobalBlockID();
02973     ConnectivityTable& conn = *(connTables_[bIndex]);
02974 
02975     GlobalID* elemConn = NULL;
02976     NodeDescriptor** elemNodeDescPtrs = NULL;
02977 
02978     int numElems = conn.elemIDs.size();
02979     if (numElems > 0) {
02980       elemConn = &((*conn.elem_conn_ids)[0]);
02981       if (!activeNodesInitialized_) {
02982         int elemConnLen = conn.elem_conn_ids->size();
02983         conn.elem_conn_ptrs = new std::vector<NodeDescriptor*>(elemConnLen);
02984       }
02985       elemNodeDescPtrs = &((*conn.elem_conn_ptrs)[0]);
02986     }
02987     int nodesPerElem = conn.numNodesPerElem;
02988 
02989     int* fieldsPerNodePtr = blocks_[bIndex]->fieldsPerNodePtr();
02990     int** fieldIDsTablePtr = blocks_[bIndex]->fieldIDsTablePtr();
02991 
02992     //
02993     //Now we want to go through the connectivity table, and for each node,
02994     //add its fields and this block to the correct NodeDescriptor from the
02995     //nodeDatabase_.
02996     //(Note that the addField and addBlock functions only add the input if
02997     //it isn't already present in that NodeDescriptor.)
02998     //
02999     //We also need to set its owner proc to localProc_ as a default, and 
03000     //inform the nodeCommMgr that the node appears in local connectivities.
03001     //
03002 
03003     conn.elemNumbers.resize(numElems);
03004 
03005     int numDistinctFields = blocks_[bIndex]->getNumDistinctFields();
03006     int fieldID = fieldIDsTablePtr[0][0];
03007 
03008     int elem;
03009     for(elem=0; elem<numElems; elem++) {
03010       conn.elemNumbers[elem] = elemNumber++;
03011       GlobalID* elemNodes = &(elemConn[elem*nodesPerElem]);
03012       NodeDescriptor** elemNodePtrs = &(elemNodeDescPtrs[elem*nodesPerElem]);
03013 
03014       for(int n=0; n<nodesPerElem; n++) {
03015         NodeDescriptor* node = NULL;
03016         int index = nodeDatabase_->getNodeWithID( elemNodes[n], node );
03017         if (index < 0) {
03018           fei::console_out() << "ERROR in SNL_FEI_Structure::initializeActiveNodes, "
03019             << FEI_ENDL << "failed to find node "
03020             << (int)(elemNodes[n]) << FEI_ENDL;
03021         }
03022 
03023         if (numDistinctFields == 1) {
03024           node->addField(fieldID);
03025         }
03026         else {
03027           for(int i=0; i<fieldsPerNodePtr[n]; i++) {
03028             node->addField(fieldIDsTablePtr[n][i]);
03029           }
03030         }
03031 
03032         int blk_idx = getIndexOfBlock(blockID);
03033         if (blk_idx >= 0) {
03034           node->addBlockIndex(blk_idx);
03035         }
03036 
03037         //now store this NodeDescriptor pointer, for fast future lookups
03038         elemNodePtrs[n] = node;
03039 
03040         node->setOwnerProc(localProc_);
03041         if (numProcs_ > 1) nodeCommMgr_->informLocal(*node);
03042       }
03043     }
03044   }
03045 
03046   //
03047   //we also need to run through the penalty constraints and inform the
03048   //nodeCommMgr that the constrained nodes are local.
03049   //
03050 
03051   if (numProcs_ > 1) {
03052     std::map<GlobalID,ConstraintType*>::const_iterator
03053       cr_iter = getPenConstRecords().begin(),
03054       cr_end = getPenConstRecords().end();
03055 
03056     while(cr_iter != cr_end) {
03057       ConstraintType& cr = *((*cr_iter).second);
03058       std::vector<GlobalID>& nodeID_vec = cr.getMasters();
03059       GlobalID* nodeIDs = &nodeID_vec[0];
03060       int numNodes = cr.getMasters().size();
03061 
03062       NodeDescriptor* node = NULL;
03063       for(int k=0; k<numNodes; ++k) {
03064    int err = nodeDatabase_->getNodeWithID(nodeIDs[k], node) ;
03065          if ( err != -1 ) {
03066            nodeCommMgr_->informLocal(*node);
03067          }
03068   //CHK_ERR( nodeDatabase_->getNodeWithID(nodeIDs[k], node) );
03069   //CHK_ERR( nodeCommMgr_->informLocal(*node) );
03070       }
03071       ++cr_iter;
03072     }
03073   }
03074 
03075   activeNodesInitialized_ = true;
03076 
03077   if (debugOutput_) os << "#  leaving finalizeActiveNodes" << FEI_ENDL;
03078   return(FEI_SUCCESS);
03079 }
03080 
03081 //------------------------------------------------------------------------------
03082 int SNL_FEI_Structure::finalizeNodeCommMgr()
03083 {
03084   //call initComplete() on the nodeCommMgr, so that it can
03085   //finalize shared-node ownerships, etc.
03086 
03087   //request the safetyCheck if the user requested it via the
03088   //parameters function.
03089   bool safetyCheck = checkSharedNodes_;
03090 
03091   if (safetyCheck && localProc_==0 && numProcs_>1 && outputLevel_>0) {
03092     FEI_COUT << "FEI Info: A consistency-check of shared-node data will be "
03093    << "performed, which involves all-to-all communication. This check is "
03094    << "done only if explicitly requested by parameter "
03095          << "'FEI_CHECK_SHARED_IDS'."<<FEI_ENDL;
03096   }
03097 
03098   CHK_ERR( nodeCommMgr_->initComplete(*nodeDatabase_, safetyCheck) );
03099 
03100   if (debugOutput_) {
03101     FEI_OSTREAM& os = dbgOut();
03102     int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
03103     os << "#     numSharedNodes: " << numSharedNodes << FEI_ENDL;
03104     for(int ns=0; ns<numSharedNodes; ns++) {
03105       NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(ns);
03106       GlobalID nodeID = node.getGlobalNodeID();
03107       int proc = node.getOwnerProc();
03108       os << "#     shNodeID " << (int)nodeID << ", owned by proc "<<proc<<FEI_ENDL;
03109     }
03110   }
03111 
03112   return(0);
03113 }
03114 
03115 //------------------------------------------------------------------------------
03116 int SNL_FEI_Structure::setNumNodesAndEqnsPerBlock()
03117 {
03118   //This function will count how 
03119   //many active nodes there are for each block.
03120 
03121    int numBlocks = blockIDs_.size();
03122    std::vector<int> nodesPerBlock(numBlocks);
03123    std::vector<int> eqnsPerBlock(numBlocks);
03124    GlobalID* blockIDsPtr = &blockIDs_[0];
03125 
03126    int j;
03127    for(j=0; j<numBlocks; j++) {
03128       nodesPerBlock[j] = 0;
03129       eqnsPerBlock[j] = 0;
03130    }
03131 
03132    int numNodes = nodeDatabase_->getNumNodeDescriptors();
03133 
03134    for(j=0; j<numNodes; j++) {
03135      NodeDescriptor* node = NULL;
03136      nodeDatabase_->getNodeAtIndex(j, node);
03137      if (node==NULL) continue;
03138 
03139      const int* fieldIDList = node->getFieldIDList();
03140 
03141      int numFields = node->getNumFields();
03142 
03143      const std::vector<unsigned>& blkIndices = node->getBlockIndexList();
03144      for(std::vector<unsigned>::const_iterator b=blkIndices.begin(), bend=blkIndices.end();
03145          b!=bend; ++b) {
03146        nodesPerBlock[*b]++;
03147 
03148        for(int fld=0; fld<numFields; fld++) {
03149          if (blocks_[*b]->containsField(fieldIDList[fld])) {
03150            int fSize = getFieldSize(fieldIDList[fld]);
03151            assert(fSize >= 0);
03152            eqnsPerBlock[*b] += fSize;
03153          }
03154        }
03155      }
03156    }
03157 
03158    for(j=0; j<numBlocks; j++) {
03159      blocks_[j]->setNumActiveNodes(nodesPerBlock[j]);
03160    }
03161 
03162    //now we need to add the elem-dof to the eqn-count for each block,
03163    //and then set the total number of eqns on each block.
03164 
03165    for(j=0; j<numBlocks; j++) {
03166      eqnsPerBlock[j] += blocks_[j]->getNumElemDOFPerElement() *
03167                         blocks_[j]->getNumElements();
03168 
03169      blocks_[j]->setTotalNumEqns(eqnsPerBlock[j]);
03170    }
03171    return(0);
03172 }
03173 
03174 //------------------------------------------------------------------------------
03175 void SNL_FEI_Structure::initializeEqnCommMgr()
03176 {
03177   FEI_OSTREAM& os = dbgOut();
03178   if (debugOutput_) os << "#  initializeEqnCommMgr" << FEI_ENDL;
03179 
03180   //This function will take information from the (internal member) nodeCommMgr
03181   //and use it to tell the eqnCommMgr which equations we can expect other
03182   //processors to contribute to, and also which equations we'll be sending to
03183   //other processors.
03184   //
03185   //This function can not be called until after initComplete() has been called
03186   //on the nodeCommMgr.
03187   //
03188   int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
03189 
03190   for(int i=0; i<numSharedNodes; i++) {
03191     NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(i);
03192 
03193     int numFields = node.getNumFields();
03194     const int* nodeFieldIDList = node.getFieldIDList();
03195     const int* nodeEqnNums = node.getFieldEqnNumbers();
03196 
03197     int owner = node.getOwnerProc();
03198     if (owner == localProc_) {
03199       std::vector<int>& proc = nodeCommMgr_->getSharedNodeProcs(i);
03200       int numProcs = proc.size();
03201 
03202       for(int p=0; p<numProcs; p++) {
03203 
03204   if (proc[p] == localProc_) continue;
03205 
03206   for(int j=0; j<numFields; j++) {
03207     int numEqns = getFieldSize(nodeFieldIDList[j]);
03208     assert(numEqns >= 0);
03209 
03210     int eqn;
03211     for(eqn=0; eqn<numEqns; eqn++) {
03212       int reducedEqn = -1;
03213       bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn,
03214              reducedEqn);
03215       if (isSlave) continue;
03216 
03217       eqnCommMgr_->addLocalEqn(reducedEqn, proc[p]);
03218     }
03219   }
03220       }
03221     }
03222     else {
03223       for(int j=0; j<numFields; j++) {
03224   int numEqns = getFieldSize(nodeFieldIDList[j]);
03225   assert(numEqns >= 0);
03226   for(int eqn=0; eqn<numEqns; eqn++) {
03227     int reducedEqn = -1;
03228     bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn, reducedEqn);
03229     if (isSlave) continue;
03230 
03231     eqnCommMgr_->addRemoteIndices(reducedEqn, owner, NULL, 0);
03232   }
03233       }
03234     }
03235   }
03236 
03237   if (debugOutput_) os << "#  leaving initializeEqnCommMgr" << FEI_ENDL;
03238 }
03239 
03240 //------------------------------------------------------------------------------
03241 void SNL_FEI_Structure::getEqnInfo(int& numGlobalEqns, int& numLocalEqns,
03242                    int& localStartRow, int& localEndRow){
03243 
03244    numGlobalEqns = numGlobalEqns_;
03245    numLocalEqns = numLocalEqns_;
03246    localStartRow = localStartRow_;
03247    localEndRow = localEndRow_;
03248 }
03249 
03250 //------------------------------------------------------------------------------
03251 int SNL_FEI_Structure::getEqnNumbers(GlobalID ID,
03252                                      int idType, int fieldID,
03253                                      int& numEqns, int* eqnNumbers)
03254 {
03255   //for now, only allow node ID. allowance for element ID is coming soon!!!
03256   if (idType != FEI_NODE) ERReturn(-1);
03257 
03258   NodeDescriptor* node = NULL;
03259   CHK_ERR( nodeDatabase_->getNodeWithID(ID, node) );
03260 
03261   numEqns = getFieldSize(fieldID);
03262   if (numEqns < 0) ERReturn(-1);
03263 
03264   int nodeFieldEqnNumber = -1;
03265   if (!node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber)) {
03266     ERReturn(-1);
03267   }
03268 
03269   for(int i=0; i<numEqns; i++) eqnNumbers[i] = nodeFieldEqnNumber + i;
03270 
03271   return(FEI_SUCCESS);
03272 }
03273 
03274 //------------------------------------------------------------------------------
03275 int SNL_FEI_Structure::getEqnNumbers(int numIDs,
03276              const GlobalID* IDs,
03277              int idType, int fieldID,
03278              int& numEqns, int* eqnNumbers)
03279 {
03280     //for now, only allow node ID. allowance for element ID is coming soon!!!
03281     if (idType != FEI_NODE) ERReturn(-1);
03282 
03283     int fieldSize = getFieldSize(fieldID);
03284     if (fieldSize < 0) {
03285       ERReturn(-1);
03286     }
03287 
03288     int offset = 0;
03289     for(int i=0; i<numIDs; ++i) {
03290       NodeDescriptor* node = NULL;
03291 
03292       if ( nodeDatabase_->getNodeWithID(IDs[i], node) != 0 ) {
03293         // fei::console_out() << "SNL_FEI_Structure::getEqnNumbers: ERROR getting node " << IDs[i] << FEI_ENDL;
03294         for(int ii=0; ii<fieldSize; ii++) {
03295           eqnNumbers[offset++] = -1;
03296         }
03297         continue;
03298         // ERReturn(-1);
03299       }
03300 
03301       int nodeFieldEqnNumber = -1;
03302 
03303       if ( !node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber) ) {
03304         //if we didn't find a field eqn-number, set eqnNumbers to -1. These
03305         //positions will then be skipped by code that passes the data on down to
03306         //underlying solver objects.
03307         for(int ii=0; ii<fieldSize; ii++) {
03308           eqnNumbers[offset++] = -1;
03309         }
03310       }
03311       else {
03312         //else we did find valid eqn-numbers...
03313         for(int ii=0; ii<fieldSize; ii++) {
03314           eqnNumbers[offset++] = nodeFieldEqnNumber + ii;
03315         }
03316       }
03317     }
03318 
03319     numEqns = fieldSize*numIDs;
03320 
03321     return(FEI_SUCCESS);
03322 }
03323 
03324 //------------------------------------------------------------------------------
03325 int SNL_FEI_Structure::translateToReducedNodeNumber(int nodeNumber, int proc)
03326 {
03327   if (proc != localProc_) {
03328     return(nodeNumber - globalNumNodesVanished_[proc]);
03329   }
03330 
03331   int insertPoint = -1;
03332   int index =
03333     fei::binarySearch(nodeNumber, localVanishedNodeNumbers_, insertPoint);
03334 
03335   int localAdjustment = index < 0 ? insertPoint : index + 1;
03336 
03337   return(nodeNumber - localAdjustment - globalNumNodesVanished_[proc]);
03338 }
03339 
03340 //------------------------------------------------------------------------------
03341 void SNL_FEI_Structure::getEqnBlkInfo(int& numGlobalEqnBlks,
03342                                      int& numLocalEqnBlks,
03343                                      int& localBlkOffset) {
03344 
03345    numGlobalEqnBlks = numGlobalEqnBlks_;
03346    numLocalEqnBlks = numLocalEqnBlks_;
03347    localBlkOffset = localBlkOffset_;
03348 }
03349 
03350 //------------------------------------------------------------------------------
03351 int SNL_FEI_Structure::calculateSlaveEqns(MPI_Comm comm)
03352 {
03353   FEI_OSTREAM& os = dbgOut();
03354   if (debugOutput_) os << "#  calculateSlaveEqns" << FEI_ENDL;
03355 
03356   if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
03357   eqnCommMgr_ = new EqnCommMgr(comm_);
03358   eqnCommMgr_->setNumRHSs(1);
03359 
03360   std::vector<int> eqns;
03361   std::vector<int> mEqns;
03362   std::vector<double> mCoefs;
03363 
03364   for(size_t i=0; i<slaveVars_->size(); i++) {
03365     int numEqns;
03366     SlaveVariable* svar = (*slaveVars_)[i];
03367 
03368     eqns.resize( getFieldSize(svar->getFieldID()));
03369     CHK_ERR( getEqnNumbers(svar->getNodeID(), FEI_NODE, svar->getFieldID(),
03370          numEqns, &eqns[0]));
03371 
03372     int slaveEqn = eqns[svar->getFieldOffset()];
03373 
03374     const std::vector<GlobalID>* mNodes = svar->getMasterNodeIDs();
03375     const std::vector<int>* mFields = svar->getMasterFields();
03376     const std::vector<double>* mWeights = svar->getWeights();
03377     const std::vector<double>& mWeightsRef = *mWeights;
03378     int mwOffset = 0;
03379 
03380     for(size_t j=0; j<mNodes->size(); j++) {
03381       int mfSize = getFieldSize((*mFields)[j]);
03382 
03383       eqns.resize(mfSize);
03384       GlobalID mNodeID = (*mNodes)[j];
03385       int mFieldID = (*mFields)[j];
03386       CHK_ERR( getEqnNumbers( mNodeID, FEI_NODE, mFieldID,
03387             mfSize, &eqns[0]));
03388 
03389       double fei_eps = 1.e-49;
03390 
03391       for(int k=0; k<mfSize; k++) {
03392   if (std::abs(mWeightsRef[mwOffset++]) > fei_eps) {
03393     mEqns.push_back(eqns[k]);
03394     mCoefs.push_back(mWeightsRef[mwOffset-1]);
03395   }
03396       }
03397     }
03398 
03399     CHK_ERR( slaveEqns_->addEqn(slaveEqn, &mCoefs[0], &mEqns[0],
03400            mEqns.size(), false) );
03401     mEqns.resize(0);
03402     mCoefs.resize(0);
03403   }
03404 
03405 #ifndef FEI_SER
03406   int numLocalSlaves = slaveVars_->size();
03407   int globalMaxSlaves = 0;
03408   CHK_ERR( fei::GlobalMax(comm_, numLocalSlaves, globalMaxSlaves) );
03409 
03410   if (globalMaxSlaves > 0) {
03411     CHK_ERR( gatherSlaveEqns(comm, eqnCommMgr_, slaveEqns_) );
03412   }
03413 #endif
03414 
03415   globalNumNodesVanished_.resize(numProcs_+1, 0);
03416 
03417   slvEqnNumbers_ = &(slaveEqns_->eqnNumbers());
03418   numSlvs_ = slvEqnNumbers_->size();
03419   if (numSlvs_ > 0) {
03420     //first, let's remove any 'couplings' among the slave equations. A coupling
03421     //is where a slave depends on a master which is also a slave that depends on
03422     //other masters.
03423 
03424     if (debugOutput_) {
03425       os << "#  slave-equations:" << FEI_ENDL;
03426       os << *slaveEqns_;
03427       os << "#  leaving calculateSlaveEqns" << FEI_ENDL;
03428     }
03429 
03430     int levelsOfCoupling;
03431     CHK_ERR( removeCouplings(*slaveEqns_, levelsOfCoupling) );
03432 
03433     if (debugOutput_) {
03434       os << "#     SNL_FEI_Structure: " << levelsOfCoupling
03435    << " level(s) of coupling removed: " << FEI_ENDL;
03436     }
03437 
03438     lowestSlv_ = (*slvEqnNumbers_)[0];
03439     highestSlv_ = (*slvEqnNumbers_)[numSlvs_-1];
03440 
03441     //For each slave equation, we need to find out if we (the local proc) either
03442     //own or share the node from which that equation arises. If we own or share
03443     //a slave node, then we will need access to the solution for each of the
03444     //associated master equations, and all other processors that share the slave
03445     //will also need access to all of the master solutions.
03446     //So,
03447     // 1. for each slave node that we share but don't own, we need to add the
03448     //   master equations to the eqnCommMgr_ object using addRemoteIndices, if
03449     //   they're not local.
03450     // 2. for each slave node that we own, we need to add the master equations
03451     // to the eqnCommMgr_ object using addLocalEqn for each processor that
03452     // shares the slave node.
03453 
03454     std::vector<int>& slvEqns = *slvEqnNumbers_;
03455     std::vector<fei::CSVec*>& mstrEqns = slaveEqns_->eqns();
03456 
03457     //keep track of the number of locally owned nodes that vanish due to the
03458     //fact that all equations at that node are slave equations.
03459     int numLocalNodesVanished = 0;
03460 
03461     GlobalID lastNodeID = -1;
03462 
03463     for(int i=0; i<numSlvs_; i++) {
03464       const NodeDescriptor* node = NULL;
03465       int reducedSlaveEqn;
03466       translateToReducedEqn(slvEqns[i], reducedSlaveEqn);
03467       int numMasters = mstrEqns[i]->size();
03468 
03469       int err = nodeDatabase_->getNodeWithEqn(slvEqns[i], node);
03470       if (err != 0) {
03471         if (debugOutput_) {
03472           os << "#  no local node for slave eqn " << slvEqns[i] << FEI_ENDL;
03473         }
03474 
03475         continue;
03476       }
03477 
03478       if (node->getGlobalNodeID() != lastNodeID &&
03479           node->getOwnerProc() == localProc_) {
03480         if (nodalEqnsAllSlaves(node, slvEqns)) {
03481           numLocalNodesVanished++;
03482           if (fei::sortedListInsert(node->getNodeNumber(), localVanishedNodeNumbers_)
03483               == -2) {
03484             ERReturn(-1);
03485           }
03486         }
03487         lastNodeID = node->getGlobalNodeID();
03488       }
03489 
03490       GlobalID slvNodeID = node->getGlobalNodeID();
03491       int shIndex = nodeCommMgr_->getSharedNodeIndex(slvNodeID);
03492       if (shIndex < 0) {
03493         continue;
03494       }
03495 
03496       std::vector<int>& sharingProcs = nodeCommMgr_->getSharedNodeProcs(shIndex);
03497 
03498       for(int j=0; j<numMasters; j++) {
03499         int masterEquation = mstrEqns[i]->indices()[j];
03500         int owner = getOwnerProcForEqn( masterEquation );
03501 
03502         int reducedMasterEqn;
03503         translateToReducedEqn(masterEquation, reducedMasterEqn);
03504 
03505         if (owner == localProc_) {
03506           int numSharing = sharingProcs.size();
03507           for(int sp=0; sp<numSharing; sp++) {
03508             if (sharingProcs[sp] == localProc_) continue;
03509 
03510             if (debugOutput_) {
03511               os << "#     slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
03512                  << ", master eqn " << masterEquation << " eqnCommMgr_->addLocal "
03513                  << reducedMasterEqn << ", proc " << sharingProcs[sp] << FEI_ENDL;
03514             }
03515             eqnCommMgr_->addLocalEqn(reducedMasterEqn, sharingProcs[sp]);
03516       slvCommMgr_->addRemoteIndices(reducedSlaveEqn, sharingProcs[sp],
03517             &reducedMasterEqn, 1);
03518           }
03519         }
03520         else {
03521           if (debugOutput_) {
03522             os << "#     slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
03523                << ", master eqn " << masterEquation << " eqnCommMgr_->addRemote "
03524                << reducedMasterEqn << ", proc " << owner << FEI_ENDL;
03525           }
03526           eqnCommMgr_->addRemoteIndices(reducedMasterEqn, owner,
03527                                         &reducedMasterEqn, 1);
03528     slvCommMgr_->addLocalEqn(reducedSlaveEqn, owner);
03529         }
03530       }
03531     }
03532 
03533     std::vector<int> tmp(1), tmp2(numProcs_);
03534     tmp[0] = numLocalNodesVanished;
03535     CHK_ERR( fei::Allgatherv(comm_, tmp, tmp2, globalNumNodesVanished_) );
03536 
03537     if ((int)globalNumNodesVanished_.size() != numProcs_) {
03538       ERReturn(-1);
03539     }
03540 
03541     globalNumNodesVanished_.resize(numProcs_+1);
03542     int tmpNumVanished = 0;
03543     for(int proc=0; proc<numProcs_; ++proc) {
03544       int temporary = tmpNumVanished;
03545       tmpNumVanished += globalNumNodesVanished_[proc];
03546       globalNumNodesVanished_[proc] = temporary;
03547     }
03548     globalNumNodesVanished_[numProcs_] = tmpNumVanished;
03549   }
03550 
03551   if (slaveMatrix_ != NULL) delete slaveMatrix_;
03552   slaveMatrix_ = new fei::FillableMat(*slaveEqns_);
03553 
03554   if (debugOutput_) {
03555     os << "#  slave-equations:" << FEI_ENDL;
03556     os << *slaveEqns_;
03557     os << "#  leaving calculateSlaveEqns" << FEI_ENDL;
03558   }
03559 
03560   return(0);
03561 }
03562 
03563 //------------------------------------------------------------------------------
03564 int SNL_FEI_Structure::removeCouplings(EqnBuffer& eqnbuf, int& levelsOfCoupling)
03565 {
03566   std::vector<int>& eqnNumbers = eqnbuf.eqnNumbers();
03567   std::vector<fei::CSVec*>& eqns = eqnbuf.eqns();
03568 
03569   std::vector<double> tempCoefs;
03570   std::vector<int> tempIndices;
03571 
03572   levelsOfCoupling = 0;
03573   bool finished = false;
03574   while(!finished) {
03575     bool foundCoupling = false;
03576     for(size_t i=0; i<eqnNumbers.size(); i++) {
03577       int rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
03578 
03579       if (rowIndex == (int)i) {
03580   fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
03581        << " illegal master-slave constraint coupling. Eqn "
03582        << eqnNumbers[i] << " is both master and slave. " << FEI_ENDL;
03583   ERReturn(-1);
03584       }
03585 
03586       while(rowIndex >= 0) {
03587   foundCoupling = true;
03588   double coef = 0.0;
03589   CHK_ERR( eqnbuf.getCoefAndRemoveIndex( eqnNumbers[rowIndex],
03590                  eqnNumbers[i], coef) );
03591 
03592   std::vector<int>& indicesRef = eqns[i]->indices();
03593   std::vector<double>& coefsRef = eqns[i]->coefs();
03594 
03595   int len = indicesRef.size();
03596   tempCoefs.resize(len);
03597   tempIndices.resize(len);
03598 
03599   double* tempCoefsPtr = &tempCoefs[0];
03600   int* tempIndicesPtr = &tempIndices[0];
03601   double* coefsPtr = &coefsRef[0];
03602   int* indicesPtr = &indicesRef[0];
03603 
03604   for(int j=0; j<len; ++j) {
03605     tempIndicesPtr[j] = indicesPtr[j];
03606     tempCoefsPtr[j] = coef*coefsPtr[j];
03607   }
03608 
03609   CHK_ERR( eqnbuf.addEqn(eqnNumbers[rowIndex], tempCoefsPtr,
03610              tempIndicesPtr, len, true) );
03611 
03612   rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
03613       }
03614     }
03615     if (foundCoupling) ++levelsOfCoupling;
03616     else finished = true;
03617 
03618     if (levelsOfCoupling>1 && !finished) {
03619       fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
03620      << " too many (>1) levels of master-slave constraint coupling. "
03621      << "Hint: coupling is considered infinite if two slaves depend on "
03622      << "each other. This may or may not be the case here." << FEI_ENDL;
03623       ERReturn(-1);
03624     }
03625   }
03626 
03627   return(0);
03628 }
03629 
03630 #ifndef FEI_SER
03631 //------------------------------------------------------------------------------
03632 int SNL_FEI_Structure::gatherSlaveEqns(MPI_Comm comm,
03633                EqnCommMgr* eqnCommMgr,
03634                EqnBuffer* slaveEqns)
03635 {
03636   int numProcs = 0;
03637   if (MPI_Comm_size(comm, &numProcs) != MPI_SUCCESS) return(-1);
03638   if (numProcs == 1) return(0);
03639   int localProc;
03640   if (MPI_Comm_rank(comm, &localProc) != MPI_SUCCESS) return(-1);
03641 
03642   //We're going to send all of our slave equations to all other processors, and
03643   //receive the slave equations from all other processors.
03644   //So we'll first fill a ProcEqns object with all of our eqn/proc pairs,
03645   //then use the regular exchange functions from EqnCommMgr. (This may not be
03646   //the most efficient way to do it, but it involves the least amount of new
03647   //code.)
03648   ProcEqns localProcEqns, remoteProcEqns;
03649   std::vector<int>& slvEqnNums = slaveEqns->eqnNumbers();
03650   fei::CSVec** slvEqnsPtr = &(slaveEqns->eqns()[0]);
03651 
03652   for(size_t i=0; i<slvEqnNums.size(); i++) {
03653     for(int p=0; p<numProcs; p++) {
03654       if (p == localProc) continue;
03655 
03656       localProcEqns.addEqn(slvEqnNums[i], slvEqnsPtr[i]->size(), p);
03657     }
03658   }
03659 
03660   CHK_ERR( eqnCommMgr->mirrorProcEqns(localProcEqns, remoteProcEqns) );
03661 
03662   CHK_ERR( eqnCommMgr->mirrorProcEqnLengths(localProcEqns,
03663               remoteProcEqns));
03664 
03665   EqnBuffer remoteEqns;
03666   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm, &localProcEqns, slaveEqns,
03667             &remoteProcEqns, &remoteEqns,
03668             false) );
03669 
03670   //Now the remoteEqns buffer should hold the slave equations from all other
03671   //processors, so let's add them to the ones we already have.
03672   CHK_ERR( slaveEqns->addEqns(remoteEqns, false) );
03673 
03674   return(0);
03675 }
03676 #endif
03677 
03678 //------------------------------------------------------------------------------
03679 bool SNL_FEI_Structure::isSlaveEqn(int eqn)
03680 {
03681   if (numSlvs_ == 0) return(false);
03682 
03683   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03684   int insertPoint = -1;
03685   int foundOffset = fei::binarySearch(eqn, slvEqns, insertPoint);
03686 
03687   if (foundOffset >= 0) return(true);
03688   else return(false);
03689 }
03690 
03691 //------------------------------------------------------------------------------
03692 bool SNL_FEI_Structure::translateToReducedEqn(int eqn, int& reducedEqn)
03693 {
03694   if (numSlvs_ == 0) { reducedEqn = eqn; return(false); }
03695 
03696   if (eqn < lowestSlv_) {reducedEqn = eqn; return(false); }
03697   if (eqn > highestSlv_) {reducedEqn = eqn - numSlvs_; return(false); }
03698 
03699   int index = 0;
03700   int foundOffset = fei::binarySearch(eqn, *slvEqnNumbers_, index);
03701 
03702   bool isSlave = false;
03703 
03704   if (foundOffset < 0) {
03705     reducedEqn = eqn - index;
03706   }
03707   else {
03708     isSlave = true; reducedEqn = eqn - (foundOffset+1);
03709   }
03710 
03711   return(isSlave);
03712 }
03713 
03714 //------------------------------------------------------------------------------
03715 int SNL_FEI_Structure::translateFromReducedEqn(int reducedEqn)
03716 {
03717   int numSlvs = slaveEqns_->getNumEqns();
03718 
03719   if (numSlvs == 0) return(reducedEqn);
03720 
03721   const int* slvEqns = &(slaveEqns_->eqnNumbers()[0]);
03722 
03723   if (reducedEqn < slvEqns[0]) return(reducedEqn);
03724 
03725   int eqn = reducedEqn;
03726 
03727   for(int i=0; i<numSlvs; i++) {
03728     if (eqn < slvEqns[i]) return(eqn);
03729     eqn++;
03730   }
03731 
03732   return(eqn);
03733 }
03734 
03735 //------------------------------------------------------------------------------
03736 int SNL_FEI_Structure::getMasterEqnNumbers(int slaveEqn,
03737             std::vector<int>*& masterEqns)
03738 {
03739   if (slaveEqns_->getNumEqns() == 0) {
03740     masterEqns = NULL;
03741     return(0);
03742   }
03743 
03744   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03745   int index = 0;
03746   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03747 
03748   if (foundOffset >= 0) {
03749     masterEqns = &(slaveEqns_->eqns()[foundOffset]->indices());
03750   }
03751   else {
03752     masterEqns = NULL;
03753   }
03754 
03755   return(0);
03756 }
03757 
03758 //------------------------------------------------------------------------------
03759 int SNL_FEI_Structure::getMasterEqnCoefs(int slaveEqn,
03760           std::vector<double>*& masterCoefs)
03761 {
03762   if (slaveEqns_->getNumEqns() == 0) {
03763     masterCoefs = NULL;
03764     return(0);
03765   }
03766 
03767   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03768   int index = 0;
03769   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03770 
03771   if (foundOffset >= 0) {
03772     masterCoefs = &(slaveEqns_->eqns()[foundOffset]->coefs());
03773   }
03774   else {
03775     masterCoefs = NULL;
03776   }
03777 
03778   return(0);
03779 }
03780 
03781 //------------------------------------------------------------------------------
03782 int SNL_FEI_Structure::getMasterEqnRHS(int slaveEqn,
03783               double& rhsValue)
03784 {
03785   if (slaveEqns_->getNumEqns() == 0) {
03786     return(0);
03787   }
03788 
03789   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03790   int index = 0;
03791   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03792 
03793   if (foundOffset >= 0) {
03794     std::vector<double>* rhsCoefsPtr = (*(slaveEqns_->rhsCoefsPtr()))[foundOffset];
03795     rhsValue = (*rhsCoefsPtr)[0];
03796   }
03797 
03798   return(0);
03799 }
03800 
03801 //------------------------------------------------------------------------------
03802 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID, 
03803                                             int interleaveStrategy,
03804                                             int* scatterIndices)
03805 {
03806    int index = fei::binarySearch(blockID, blockIDs_);
03807 
03808    if (index < 0) {
03809       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
03810            << (int)blockID << " not found. Aborting." << FEI_ENDL;
03811       std::abort();
03812    }
03813 
03814    std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
03815 
03816    std::map<GlobalID,int>::const_iterator
03817      iter = elemIDs.find(elemID);
03818 
03819    if (iter == elemIDs.end()) {
03820       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: " 
03821            << (int)blockID << ", elemID "
03822            << (int)elemID << " not found. Aborting." << FEI_ENDL;
03823       std::abort();
03824    }
03825 
03826    int elemIndex = iter->second;
03827 
03828    getScatterIndices_index(index, elemIndex,
03829                            interleaveStrategy, scatterIndices);
03830 }
03831 
03832 //------------------------------------------------------------------------------
03833 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID, 
03834                int interleaveStrategy,
03835                int* scatterIndices,
03836                int* blkScatterIndices,
03837                int* blkSizes)
03838 {
03839    int index = fei::binarySearch(blockID, blockIDs_);
03840 
03841    if (index < 0) {
03842       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
03843            << (int)blockID << " not found. Aborting." << FEI_ENDL;
03844       std::abort();
03845    }
03846 
03847    std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
03848 
03849    std::map<GlobalID,int>::const_iterator
03850      iter = elemIDs.find(elemID);
03851 
03852    if (iter == elemIDs.end()) {
03853       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: " 
03854            << (int)blockID << ", elemID "
03855            << (int)elemID << " not found. Aborting." << FEI_ENDL;
03856       std::abort();
03857    }
03858 
03859    int elemIndex = iter->second;
03860 
03861    getScatterIndices_index(index, elemIndex,
03862                            interleaveStrategy, scatterIndices,
03863          blkScatterIndices, blkSizes);
03864 }
03865 
03866 //------------------------------------------------------------------------------
03867 int SNL_FEI_Structure::getBlkScatterIndices_index(int blockIndex,
03868               int elemIndex,
03869               int* scatterIndices)
03870 {
03871   BlockDescriptor& block = *(blocks_[blockIndex]);
03872   int numNodes = block.getNumNodesPerElement();
03873   work_nodePtrs_.resize(numNodes);
03874   NodeDescriptor** nodes = &work_nodePtrs_[0];
03875   int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03876   if (err) {
03877     fei::console_out() << "SNL_FEI_Structure::getBlkScatterIndices_index: ERROR getting"
03878    << " node descriptors." << FEI_ENDL;
03879     ERReturn(-1);
03880   }
03881 
03882   int offset = 0;
03883   return( getNodeBlkIndices(nodes, numNodes, scatterIndices, offset) );
03884 }
03885 
03886 //------------------------------------------------------------------------------
03887 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
03888             int interleaveStrategy,
03889             int* scatterIndices)
03890 {
03891 //On input, scatterIndices, is assumed to be allocated by the calling code,
03892 // and be of length the number of equations per element.
03893 //
03894    BlockDescriptor& block = *(blocks_[blockIndex]);
03895    int numNodes = block.getNumNodesPerElement();
03896    int* fieldsPerNode = block.fieldsPerNodePtr();
03897    int** fieldIDs = block.fieldIDsTablePtr();
03898 
03899    work_nodePtrs_.resize(numNodes);
03900    NodeDescriptor** nodes = &work_nodePtrs_[0];
03901 
03902    int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03903    if (err) {
03904       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
03905            << " node descriptors." << FEI_ENDL;
03906       std::abort();
03907    }
03908 
03909    int offset = 0;
03910    if (fieldDatabase_->size() == 1) {
03911      err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
03912             scatterIndices, offset);
03913      if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
03914    }
03915    else {
03916      switch (interleaveStrategy) {
03917      case 0:
03918        err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03919          scatterIndices, offset);
03920        if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
03921        break;
03922 
03923      case 1:
03924        err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03925           scatterIndices, offset);
03926        if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
03927        break;
03928 
03929      default:
03930        fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
03931        break;
03932      }
03933    }
03934 
03935    //now the element-DOF.
03936    int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
03937    std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
03938 
03939    for(int i=0; i<numElemDOF; i++) {
03940       scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
03941    }
03942 }
03943 
03944 //------------------------------------------------------------------------------
03945 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
03946             int interleaveStrategy,
03947             int* scatterIndices,
03948             int* blkScatterIndices,
03949             int* blkSizes)
03950 {
03951 //On input, scatterIndices, is assumed to be allocated by the calling code,
03952 // and be of length the number of equations per element.
03953 //
03954    BlockDescriptor& block = *(blocks_[blockIndex]);
03955    int numNodes = block.getNumNodesPerElement();
03956    int* fieldsPerNode = block.fieldsPerNodePtr();
03957    int** fieldIDs = block.fieldIDsTablePtr();
03958 
03959    work_nodePtrs_.resize(numNodes);
03960    NodeDescriptor** nodes = &work_nodePtrs_[0];
03961 
03962    int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03963    if (err) {
03964       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
03965            << " node descriptors." << FEI_ENDL;
03966       std::abort();
03967    }
03968 
03969    int offset = 0, blkOffset = 0;
03970    if (fieldDatabase_->size() == 1) {
03971      err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
03972          scatterIndices, offset,
03973          blkScatterIndices, blkSizes, blkOffset);
03974      if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
03975    }
03976    else {
03977      switch (interleaveStrategy) {
03978      case 0:
03979        err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03980          scatterIndices, offset,
03981          blkScatterIndices, blkSizes, blkOffset);
03982        if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
03983        break;
03984 
03985      case 1:
03986        err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03987           scatterIndices, offset);
03988        if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
03989        break;
03990 
03991      default:
03992        fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
03993        break;
03994      }
03995    }
03996 
03997    //now the element-DOF.
03998    int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
03999    std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
04000 
04001    for(int i=0; i<numElemDOF; i++) {
04002       scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
04003       if (interleaveStrategy == 0) {
04004   blkSizes[blkOffset] = 1;
04005   blkScatterIndices[blkOffset++] = elemDOFEqns[elemIndex] + i;
04006       }
04007    }
04008 }
04009 
04010 //------------------------------------------------------------------------------
04011 int SNL_FEI_Structure::getElemNodeDescriptors(int blockIndex, int elemIndex,
04012                                              NodeDescriptor** nodes)
04013 {
04014   //Private function, called by 'getScatterIndices_index'. We can safely
04015   //assume that blockIndex and elemIndex are valid in-range indices.
04016   //
04017   //This function's task is to obtain the NodeDescriptor objects, from the
04018   //nodeDatabase, that are connected to the specified element.
04019   //
04020   int numNodes = connTables_[blockIndex]->numNodesPerElem;
04021 
04022   if (activeNodesInitialized_) {
04023     NodeDescriptor** elemNodePtrs =
04024       &((*connTables_[blockIndex]->elem_conn_ptrs)[elemIndex*numNodes]);
04025     for(int i=0; i<numNodes; ++i) nodes[i] = elemNodePtrs[i];
04026   }
04027   else {
04028     const GlobalID* elemConn = 
04029       &((*connTables_[blockIndex]->elem_conn_ids)[elemIndex*numNodes]);
04030     for(int i=0; i<numNodes; ++i) {
04031       CHK_ERR( nodeDatabase_->getNodeWithID(elemConn[i], nodes[i]));
04032     }
04033   }
04034 
04035   return(FEI_SUCCESS);
04036 }
04037 
04038 //------------------------------------------------------------------------------
04039 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
04040                int numNodes,
04041                int fieldID,
04042                int* scatterIndices,
04043                int& offset)
04044 {
04045   int fieldSize = getFieldSize(fieldID);
04046 
04047   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04048     NodeDescriptor& node = *(nodes[nodeIndex]);
04049     const int* eqnNumbers = node.getFieldEqnNumbers();
04050     int eqn = eqnNumbers[0];
04051     scatterIndices[offset++] = eqn;
04052     if (fieldSize > 1) {
04053       for(int i=1; i<fieldSize; i++) {
04054   scatterIndices[offset++] = eqn+i;
04055       }
04056     }
04057   }
04058   return(0);
04059 }
04060 
04061 //------------------------------------------------------------------------------
04062 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
04063                int numNodes,
04064                int fieldID,
04065                int* scatterIndices,
04066                int& offset,
04067                int* blkScatterIndices,
04068                int* blkSizes,
04069                int& blkOffset)
04070 {
04071   int fieldSize = getFieldSize(fieldID);
04072 
04073   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04074     NodeDescriptor& node = *(nodes[nodeIndex]);
04075     const int* eqnNumbers = node.getFieldEqnNumbers();
04076     int eqn = eqnNumbers[0];
04077     scatterIndices[offset++] = eqn;
04078     if (fieldSize > 1) {
04079       for(int i=1; i<fieldSize; i++) {
04080   scatterIndices[offset++] = eqn+i;
04081       }
04082     }
04083     blkSizes[blkOffset] = node.getNumNodalDOF();
04084     blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
04085   }
04086   return(0);
04087 }
04088 
04089 //------------------------------------------------------------------------------
04090 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04091                                           int** fieldIDs, int* fieldsPerNode,
04092                                           int* scatterIndices, int& offset)
04093 {
04094   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04095 
04096       NodeDescriptor& node = *(nodes[nodeIndex]);
04097 
04098       const int* nodeFieldIDList = node.getFieldIDList();
04099       const int* nodeEqnNums = node.getFieldEqnNumbers();
04100       int numFields = node.getNumFields();
04101 
04102       int* fieldID_ind = fieldIDs[nodeIndex];
04103 
04104       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04105          int numEqns = getFieldSize(fieldID_ind[j]);
04106          assert(numEqns >= 0);
04107 
04108          int insert = -1;
04109          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04110                                              numFields, insert);
04111 
04112          if (nind >= 0) {
04113      int eqn = nodeEqnNums[nind];
04114 
04115      if (eqn < 0) {
04116        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04117       << (int)node.getGlobalNodeID()
04118       << ", field " << nodeFieldIDList[nind]
04119       << " has equation number " << eqn << FEI_ENDL;
04120        std::abort();
04121      }
04122 
04123      for(int jj=0; jj<numEqns; jj++) {
04124        scatterIndices[offset++] = eqn+jj;
04125      }
04126          }
04127          else {
04128      if (outputLevel_ > 2) {
04129        fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
04130       << " not found for node "
04131       << (int)(node.getGlobalNodeID()) << FEI_ENDL;
04132      }
04133          }
04134       }
04135    }
04136 
04137    return(FEI_SUCCESS);
04138 }
04139 
04140 //------------------------------------------------------------------------------
04141 int SNL_FEI_Structure::getNodeBlkIndices(NodeDescriptor** nodes,
04142            int numNodes,
04143            int* scatterIndices,
04144            int& offset)
04145 {
04146   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04147     NodeDescriptor* node = nodes[nodeIndex];
04148     scatterIndices[offset++] = node->getBlkEqnNumber();
04149   }
04150   return(0);
04151 }
04152 
04153 //------------------------------------------------------------------------------
04154 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04155              int** fieldIDs, int* fieldsPerNode,
04156              int* scatterIndices, int& offset,
04157              int* blkScatterIndices,
04158              int* blkSizes,
04159              int& blkOffset)
04160 {
04161   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04162 
04163       NodeDescriptor& node = *(nodes[nodeIndex]);
04164 
04165       const int* nodeFieldIDList = node.getFieldIDList();
04166       const int* nodeEqnNums = node.getFieldEqnNumbers();
04167       int numFields = node.getNumFields();
04168 
04169       blkSizes[blkOffset] = node.getNumNodalDOF();
04170       blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
04171 
04172       int* fieldID_ind = fieldIDs[nodeIndex];
04173 
04174       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04175          int numEqns = getFieldSize(fieldID_ind[j]);
04176          assert(numEqns >= 0);
04177 
04178          int insert = -1;
04179          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04180                                              numFields, insert);
04181 
04182          if (nind >= 0) {
04183      int eqn = nodeEqnNums[nind];
04184      if (eqn < 0) {
04185        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04186       << (int)node.getGlobalNodeID()
04187       << ", field " << nodeFieldIDList[nind]
04188       << " has equation number " << eqn << FEI_ENDL;
04189        std::abort();
04190      }
04191 
04192      for(int jj=0; jj<numEqns; jj++) {
04193        scatterIndices[offset++] = eqn+jj;
04194      }
04195          }
04196          else {
04197      if (outputLevel_ > 2) {
04198        fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
04199       << " not found for node "
04200       << (int)(node.getGlobalNodeID()) << FEI_ENDL;
04201      }
04202          }
04203       }
04204    }
04205 
04206    return(FEI_SUCCESS);
04207 }
04208 
04209 //------------------------------------------------------------------------------
04210 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04211                                           std::vector<int>* fieldIDs,
04212                                           std::vector<int>& fieldsPerNode,
04213                                           std::vector<int>& scatterIndices)
04214 {
04215    int offset = 0;
04216    scatterIndices.resize(0);
04217 
04218    for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04219 
04220       NodeDescriptor& node = *(nodes[nodeIndex]);
04221 
04222       const int* nodeFieldIDList = node.getFieldIDList();
04223       const int* nodeEqnNums = node.getFieldEqnNumbers();
04224       int numFields = node.getNumFields();
04225 
04226       int* fieldID_ind = &(fieldIDs[nodeIndex][0]);
04227 
04228       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04229          int numEqns = getFieldSize(fieldID_ind[j]);
04230          assert(numEqns >= 0);
04231 
04232          int insert = -1;
04233          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04234                                              numFields, insert);
04235 
04236          if (nind >= 0) {
04237      int eqn = nodeEqnNums[nind];
04238 
04239      if (eqn < 0) {
04240        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04241       << (int)node.getGlobalNodeID()
04242       << ", field " << nodeFieldIDList[nind]
04243       << " has equation number " << eqn << FEI_ENDL;
04244        MPI_Abort(comm_, -1);
04245      }
04246 
04247      scatterIndices.resize(offset+numEqns);
04248      int* inds = &scatterIndices[0];
04249 
04250      for(int jj=0; jj<numEqns; jj++) {
04251        inds[offset+jj] = eqn+jj;
04252      }
04253      offset += numEqns;
04254          }
04255          else {
04256      if (outputLevel_ > 2) {
04257        fei::console_out() << "WARNING, field ID " << fieldID_ind[j]
04258       << " not found for node "
04259       << (int)node.getGlobalNodeID() << FEI_ENDL;
04260            }
04261          }
04262       }
04263    }
04264 
04265    return(FEI_SUCCESS);
04266 }
04267 
04268 //------------------------------------------------------------------------------
04269 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
04270                                           int** fieldIDs, int* fieldsPerNode,
04271                                           int* scatterIndices, int& offset)
04272 {
04273   //In this function we want to group equations by field, but
04274   //in what order should the fields be?
04275   //Let's just run through the fieldIDs table, and add the fields to a
04276   //flat list, in the order we encounter them, but making sure no fieldID
04277   //gets added more than once.
04278 
04279   std::vector<int> fields;
04280   for(int ii=0; ii<numNodes; ii++) {
04281     for(int j=0; j<fieldsPerNode[ii]; j++) {
04282       if (std::find(fields.begin(), fields.end(), fieldIDs[ii][j]) == fields.end()) {
04283         fields.push_back(fieldIDs[ii][j]);
04284       }
04285     }
04286   }
04287 
04288   int* fieldsPtr = fields.size()>0 ? &fields[0] : NULL;
04289 
04290   //ok, we've got our flat list of fields, so let's proceed to get the
04291   //scatter indices.
04292 
04293   for(size_t i=0; i<fields.size(); i++) {
04294     int field = fieldsPtr[i];
04295 
04296     for(int nodeIndex = 0; nodeIndex < numNodes; ++nodeIndex) {
04297       int fidx = fei::searchList(field, fieldIDs[nodeIndex],
04298             fieldsPerNode[nodeIndex]);
04299       if (fidx < 0) {
04300   continue;
04301       }
04302 
04303       NodeDescriptor* node = nodes[nodeIndex];
04304 
04305       const int* nodeFieldIDList = node->getFieldIDList();
04306       const int* nodeEqnNums = node->getFieldEqnNumbers();
04307       int numFields = node->getNumFields();
04308 
04309       int numEqns = getFieldSize(field);
04310       assert(numEqns >= 0);
04311 
04312       int insert = -1;
04313       int nind = fei::binarySearch(field, nodeFieldIDList,
04314             numFields, insert);
04315 
04316       if (nind > -1) {
04317   for(int jj=0; jj<numEqns; ++jj) {
04318     scatterIndices[offset++] = nodeEqnNums[nind]+jj;
04319   }
04320       }
04321       else {
04322   ERReturn(-1);
04323       }
04324     }
04325   }
04326 
04327   return(FEI_SUCCESS);
04328 }
04329 
04330 //------------------------------------------------------------------------------
04331 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
04332                                           std::vector<int>* fieldIDs,
04333                                           std::vector<int>& fieldsPerNode,
04334                                           std::vector<int>& scatterIndices)
04335 {
04336    //In this function we want to group equations by field, but
04337    //in what order should the fields be?
04338    //Let's just run through the fieldIDs table, and add the fields to a
04339    //flat list, in the order we encounter them, but making sure no fieldID
04340    //gets added more than once.
04341 
04342    std::vector<int> fields;
04343    for(int ii=0; ii<numNodes; ii++) {
04344       for(int j=0; j<fieldsPerNode[ii]; j++) {
04345          std::vector<int>& fldIDs = fieldIDs[ii];
04346          if (std::find(fields.begin(), fields.end(), fldIDs[j]) == fields.end()) {
04347            fields.push_back(fldIDs[j]);
04348          }
04349       }
04350    }
04351 
04352    //ok, we've got our flat list of fields, so let's proceed to get the
04353    //scatter indices.
04354 
04355    int offset = 0;
04356    scatterIndices.resize(0);
04357 
04358    for(size_t i=0; i<fields.size(); i++) {
04359       for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04360 
04361          const int* nodeFieldIDList = nodes[nodeIndex]->getFieldIDList();
04362          const int* nodeEqnNums = nodes[nodeIndex]->getFieldEqnNumbers();
04363          int numFields = nodes[nodeIndex]->getNumFields();
04364 
04365          int numEqns = getFieldSize(fields[i]);
04366          assert(numEqns >= 0);
04367 
04368          int insert = -1;
04369          int nind = fei::binarySearch(fields[i], nodeFieldIDList,
04370                                              numFields, insert);
04371 
04372          if (nind >= 0) {
04373             for(int jj=0; jj<numEqns; jj++) {
04374                scatterIndices.push_back(nodeEqnNums[nind]+jj);
04375          offset++;
04376             }
04377          }
04378          else {
04379      if (outputLevel_ > 2) {
04380        fei::console_out() << "WARNING, field ID " << fields[i]
04381       << " not found for node "
04382       << (int)nodes[nodeIndex]->getGlobalNodeID() << FEI_ENDL;
04383      }
04384          }
04385       }
04386    }
04387 
04388    return(FEI_SUCCESS);
04389 }
04390 
04391 //------------------------------------------------------------------------------
04392 void SNL_FEI_Structure::addCR(int CRID,
04393            snl_fei::Constraint<GlobalID>*& cr,
04394          std::map<GlobalID,snl_fei::Constraint<GlobalID>* >& crDB)
04395 {
04396   std::map<GlobalID,snl_fei::Constraint<GlobalID>* >::iterator
04397     cr_iter = crDB.find(CRID);
04398 
04399   if (cr_iter == crDB.end()) {
04400     cr = new ConstraintType;
04401     crDB.insert(std::pair<GlobalID,snl_fei::Constraint<GlobalID>*>(CRID, cr));
04402   }
04403 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends