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