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