FEI Version of the Day
SNL_FEI_Structure.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include "fei_sstream.hpp"
00010 #include "fei_fstream.hpp"
00011 
00012 #include <limits>
00013 #include <cmath>
00014 #include <assert.h>
00015 
00016 #include "fei_defs.h"
00017 
00018 #include "fei_TemplateUtils.hpp"
00019 #include <fei_CommUtils.hpp>
00020 #include "snl_fei_Constraint.hpp"
00021 typedef snl_fei::Constraint<GlobalID> ConstraintType;
00022 
00023 #include "fei_NodeDescriptor.hpp"
00024 #include "fei_NodeCommMgr.hpp"
00025 #include "fei_NodeDatabase.hpp"
00026 
00027 #include "fei_SlaveVariable.hpp"
00028 
00029 #include "fei_BlockDescriptor.hpp"
00030 
00031 #include "snl_fei_PointBlockMap.hpp"
00032 #include "fei_ProcEqns.hpp"
00033 #include "fei_EqnBuffer.hpp"
00034 #include <fei_FillableMat.hpp>
00035 #include <fei_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::console_out() << "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::console_out() << "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::console_out() << "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::console_out() << "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::console_out() << 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::console_out() << 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::console_out() << "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::console_out() << 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::console_out() << "ERROR, findNodeDescriptor unable to find node " << (int)nodeID
01354    << FEI_ENDL;
01355     std::abort();
01356   }
01357 
01358   return( *node );
01359 }
01360 
01361 //------------------------------------------------------------------------------
01362 NodeDescriptor* SNL_FEI_Structure::findNode(GlobalID nodeID)
01363 {
01364   NodeDescriptor* node = NULL;
01365   nodeDatabase_->getNodeWithID(nodeID, node);
01366   return node;
01367 }
01368 
01369 //------------------------------------------------------------------------------
01370 int SNL_FEI_Structure::initMultCRStructure() 
01371 {
01372   FEI_OSTREAM& os = dbgOut();
01373   if (debugOutput_) os << "#    initMultCRStructure" << FEI_ENDL;
01374   //
01375   //Private SNL_FEI_Structure function, to be called from initComplete.
01376   //
01377   //Records the system matrix structure arising from Lagrange Multiplier
01378   //Constraint Relations.
01379   //
01380   // since at initialization all we are being passed is the
01381   // general form of the constraint equations, we can't check to see if any
01382   // of the weight terms are zeros at this stage of the game.  Hence, we
01383   // have to reserve space for all the nodal weight vectors, even though
01384   // they might turn out to be zeros during the load step....
01385 
01386   std::map<GlobalID,ConstraintType*>::const_iterator
01387     cr_iter = multCRs_.begin(),
01388     cr_end  = multCRs_.end();
01389 
01390   while(cr_iter != cr_end) {
01391     ConstraintType& multCR = *((*cr_iter).second);
01392 
01393     int lenList = multCR.getMasters().size();
01394 
01395     std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
01396     GlobalID *CRNodePtr = &CRNode_vec[0];
01397     std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
01398     int* CRFieldPtr = &CRField_vec[0];
01399 
01400     int crEqn = multCR.getEqnNumber();
01401     int reducedCrEqn = 0;
01402     translateToReducedEqn(crEqn, reducedCrEqn);
01403 
01404     createMatrixPosition(reducedCrEqn, reducedCrEqn, "initMCRStr");
01405 
01406     for(int j=0; j<lenList; j++) {
01407       GlobalID nodeID = CRNodePtr[j];
01408       int fieldID = CRFieldPtr[j];
01409 
01410       NodeDescriptor *nodePtr = findNode(nodeID);
01411       if(nodePtr==NULL) continue;
01412       NodeDescriptor& node = *nodePtr;
01413 
01414       //first, store the column indices associated with this node, in
01415       //the constraint's equation. i.e., position (crEqn, node)
01416 
01417       storeNodalColumnIndices(crEqn, node, fieldID);
01418 
01419       //now we need to store the transpose of the above contribution,
01420       //i.e., position(s) (node, crEqn)
01421 
01422       if (node.getOwnerProc() == localProc_) {
01423   //if we own this node, we will simply store its equation
01424   //numbers in local rows (equations) of the matrix.
01425 
01426   storeNodalRowIndices(node, fieldID, crEqn);
01427       }
01428       else {
01429   //if we don't own it, then we need to register with the
01430   //eqn comm mgr that we'll be sending contributions to
01431   //column crEqn of the remote equations associated with this
01432   //node.
01433 
01434   storeNodalSendIndex(node, fieldID, crEqn);
01435       }
01436     }
01437     ++cr_iter;
01438   }
01439   return(FEI_SUCCESS);
01440 }
01441 
01442 //------------------------------------------------------------------------------
01443 int SNL_FEI_Structure::initPenCRStructure()
01444 {
01445   FEI_OSTREAM& os = dbgOut();
01446   if (debugOutput_) os << "#     initPenCRStructure" << FEI_ENDL;
01447 //
01448 //This function oversees the putting in of any matrix structure that results
01449 //from Penalty constraint relations.
01450 //
01451 // note that penalty constraints don't generate new equations
01452 // (unlike Lagrange constraints), but they do add terms to the system
01453 // stiffness matrix that couple all the nodes that contribute to the
01454 // penalty constraint.  In addition, each submatrix is defined by the pair
01455 // of nodes that creates its contribution, hence a submatrix can be defined
01456 // in terms of two weight vectors (of length p and q) instead of the
01457 // generally larger product matrix (of size pq)
01458 
01459 // the additional terms take the form of little submatrices that look a lot 
01460 // like an element stiffness and load matrix, where the nodes in the
01461 // constraint list take on the role of the nodes associated with an
01462 // element, and the individual matrix terms arise from the outer products
01463 // of the constraint nodal weight vectors 
01464 
01465 // rather than get elegant and treat this task as such an elemental energy
01466 // term, we'll use some brute force to construct these penalty contributions
01467 // on the fly, primarily to simplify -reading- this thing, so that the 
01468 // derivations in the annotated implementation document are more readily
01469 // followed...
01470   std::map<GlobalID,ConstraintType*>::const_iterator
01471     cr_iter = penCRs_.begin(),
01472     cr_end = penCRs_.end();
01473 
01474   while(cr_iter != cr_end) {
01475     ConstraintType& penCR = *((*cr_iter).second);
01476 
01477     int lenList = penCR.getMasters().size();
01478     std::vector<GlobalID>& CRNode_vec = penCR.getMasters();
01479     GlobalID* CRNodesPtr = &CRNode_vec[0];
01480 
01481     std::vector<int>& CRField_vec = penCR.getMasterFieldIDs();
01482     int* CRFieldPtr = &CRField_vec[0];
01483 
01484     // each constraint equation generates a set of nodal energy terms, so
01485     // we have to process a matrix of nodes for each constraint
01486 
01487     for(int i = 0; i < lenList; i++) {
01488       GlobalID iNodeID = CRNodesPtr[i];
01489       int iField = CRFieldPtr[i];
01490 
01491       NodeDescriptor* iNodePtr = findNode(iNodeID);
01492       if(!iNodePtr) continue;
01493       NodeDescriptor& iNode = *iNodePtr;
01494 
01495       for(int j = 0; j < lenList; j++) {
01496   GlobalID jNodeID = CRNodesPtr[j];
01497   int jField = CRFieldPtr[j];
01498 
01499   NodeDescriptor* jNodePtr = findNode(jNodeID);
01500         if(!jNodePtr) continue;
01501         NodeDescriptor &jNode = *jNodePtr;
01502 
01503   if (iNode.getOwnerProc() == localProc_) {
01504     //if iNode is local, we'll put equations into the local 
01505     //matrix structure.
01506 
01507     storeLocalNodeIndices(iNode, iField, jNode, jField);
01508   }
01509   else {
01510     //if iNode is remotely owned, we'll be registering equations
01511     //to send to its owning processor.
01512 
01513     storeNodalSendIndices(iNode, iField, jNode, jField);
01514   }
01515       }
01516     }   //   end i loop
01517     ++cr_iter;
01518    }   //   end while loop
01519    return(FEI_SUCCESS);
01520 }
01521  
01522 //------------------------------------------------------------------------------
01523 void SNL_FEI_Structure::storeNodalSendIndex(NodeDescriptor& node, int fieldID,
01524               int col)
01525 {
01526   //
01527   //This is a private SNL_FEI_Structure function. We can safely assume that it
01528   //will only be called with a node that is not locally owned.
01529   //
01530   //This function tells the eqn comm mgr that we'll be sending contributions
01531   //to column 'col' for the equations associated with 'fieldID', on 'node', on
01532   //node's owning processor.
01533   //
01534   int proc = node.getOwnerProc();
01535 
01536   int eqnNumber = -1;
01537   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
01538 
01539   int numEqns = getFieldSize(fieldID);
01540   if (numEqns < 1) {
01541     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01542    <<") with size "<<numEqns<<FEI_ENDL;
01543     voidERReturn;
01544   }
01545 
01546   for(int i=0; i<numEqns; i++) {
01547     eqnCommMgr_->addRemoteIndices(eqnNumber+i, proc, &col, 1);
01548   }
01549 }
01550 
01551 //------------------------------------------------------------------------------
01552 void SNL_FEI_Structure::storeNodalSendIndices(NodeDescriptor& iNode, int iField,
01553                                   NodeDescriptor& jNode, int jField){
01554 //
01555 //This function will register with the eqn comm mgr the equations associated
01556 //with iNode, field 'iField' having column indices that are the equations
01557 //associated with jNode, field 'jField', to be sent to the owner of iNode.
01558 //
01559    int proc = iNode.getOwnerProc();
01560 
01561    int iEqn = -1, jEqn = -1;
01562    if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
01563    if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
01564 
01565    int iNumParams = getFieldSize(iField);
01566    int jNumParams = getFieldSize(jField);
01567    if (iNumParams < 1 || jNumParams < 1) {
01568      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
01569     << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
01570     << jNumParams<<FEI_ENDL;
01571      voidERReturn;
01572    }
01573 
01574    for(int i=0; i<iNumParams; i++) {
01575       int eqn = iEqn + i;
01576 
01577       for(int j=0; j<jNumParams; j++) {
01578          int col = jEqn + j;
01579          eqnCommMgr_->addRemoteIndices(eqn, proc, &col, 1);
01580       }
01581    }
01582 }
01583 
01584 //------------------------------------------------------------------------------
01585 void SNL_FEI_Structure::storeLocalNodeIndices(NodeDescriptor& iNode, int iField,
01586                NodeDescriptor& jNode, int jField)
01587 {
01588 //
01589 //This function will add to the local matrix structure the equations associated
01590 //with iNode at iField, having column indices that are the equations associated
01591 //with jNode at jField.
01592 //
01593    int iEqn = -1, jEqn = -1;
01594    if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
01595    if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
01596 
01597    int iNumParams = getFieldSize(iField);
01598    int jNumParams = getFieldSize(jField);
01599    if (iNumParams < 1 || jNumParams < 1) {
01600      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
01601     << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
01602     << jNumParams<<FEI_ENDL;
01603      voidERReturn;
01604    }
01605 
01606    for(int i=0; i<iNumParams; i++) {
01607       int row = iEqn + i;
01608       int reducedRow = -1;
01609       bool isSlave = translateToReducedEqn(row, reducedRow);
01610       if (isSlave) continue;
01611       
01612       for(int j=0; j<jNumParams; j++) {
01613          int col = jEqn + j;
01614    int reducedCol = -1;
01615    isSlave = translateToReducedEqn(col, reducedCol);
01616    if (isSlave) continue;
01617 
01618          createMatrixPosition(reducedRow, reducedCol,
01619             "storeLocNdInd");
01620       }
01621    }
01622 }
01623 
01624 //------------------------------------------------------------------------------
01625 void SNL_FEI_Structure::storeNodalColumnIndices(int eqn, NodeDescriptor& node,
01626                  int fieldID)
01627 {
01628   //
01629   //This function stores the equation numbers associated with 'node' at
01630   //'fieldID' as column indices in row 'eqn' of the system matrix structure.
01631   //
01632   if ((localStartRow_ > eqn) || (eqn > localEndRow_)) {
01633     return;
01634   }
01635 
01636   int colNumber = -1;
01637   if (!node.getFieldEqnNumber(fieldID, colNumber)) voidERReturn;
01638 
01639   int numParams = getFieldSize(fieldID);
01640   if (numParams < 1) {
01641     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01642    <<") with size "<<numParams<<FEI_ENDL;
01643     voidERReturn;
01644   }
01645 
01646   int reducedEqn = -1;
01647   bool isSlave = translateToReducedEqn(eqn, reducedEqn);
01648   if (isSlave) return;
01649 
01650   for(int j=0; j<numParams; j++) {
01651     int reducedCol = -1;
01652     isSlave = translateToReducedEqn(colNumber+j, reducedCol);
01653     if (isSlave) continue;
01654 
01655     createMatrixPosition(reducedEqn, reducedCol,
01656        "storeNdColInd");
01657   }
01658 }
01659 
01660 //------------------------------------------------------------------------------
01661 void SNL_FEI_Structure::storeNodalRowIndices(NodeDescriptor& node,
01662               int fieldID, int eqn)
01663 {
01664   //
01665   //This function stores column 'eqn' in equation numbers associated with
01666   //'node' at 'fieldID' in the system matrix structure. 
01667   //
01668   int eqnNumber = -1;
01669   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
01670 
01671   int numParams = getFieldSize(fieldID);
01672   if (numParams < 1) {
01673     fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
01674    <<") with size "<<numParams<<FEI_ENDL;
01675     voidERReturn;
01676   }
01677 
01678   int reducedEqn = -1;
01679   bool isSlave = translateToReducedEqn(eqn, reducedEqn);
01680   if (isSlave) return;
01681 
01682   for(int j=0; j<numParams; j++) {
01683     int reducedRow = -1;
01684     isSlave = translateToReducedEqn(eqnNumber+j, reducedRow);
01685     if (isSlave) continue;
01686 
01687     createMatrixPosition(reducedRow, reducedEqn, "storeNdRowInd");
01688   }
01689 }
01690 
01691 //------------------------------------------------------------------------------
01692 int SNL_FEI_Structure::createSymmEqnStructure(std::vector<int>& scatterIndices)
01693 {
01694   //scatterIndices are both the rows and the columns for a structurally
01695   //symmetric 2-dimensional contribution to the matrix.
01696   //
01697 
01698   //if there aren't any slaves in the whole problem, store the scatter-indices
01699   //using a fast function (no translations-to-reduced-space) and return.
01700   if (numSlaveEquations() == 0) {
01701     CHK_ERR( storeElementScatterIndices_noSlaves(scatterIndices) );
01702     return(FEI_SUCCESS);
01703   }
01704 
01705   try {
01706 
01707   int len = scatterIndices.size();
01708   bool anySlaves = false;
01709   rSlave_.resize(len);
01710   for(int is=0; is<len; is++) { 
01711     rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
01712     if (rSlave_[is] == 1) anySlaves = true;
01713   }
01714 
01715   //if there aren't any slaves in this contribution, then just store it
01716   //and return
01717   if (!anySlaves) {
01718     CHK_ERR( storeElementScatterIndices(scatterIndices) );
01719     return(FEI_SUCCESS);
01720   }
01721 
01722   int* scatterPtr = &scatterIndices[0];
01723 
01724   workSpace_.resize(len);
01725   for(int j=0; j<len; j++) {
01726     translateToReducedEqn(scatterPtr[j], workSpace_[j]);
01727   }
01728   
01729   for(int i=0; i<len; i++) {
01730     int row = scatterPtr[i];
01731     if (rSlave_[i] == 1) {
01732       reducedEqnCounter_++;
01733       //
01734       //'row' is a slave equation, so add this row to Kdi_. But as we do that,
01735       //watch for columns that are slave equations and add them to Kid_.
01736       //
01737       for(int jj=0; jj<len; jj++) {
01738   int col = scatterPtr[jj];
01739 
01740   if (rSlave_[jj] == 1) {
01741     //'col' is a slave equation, so add this column to Kid_.
01742     for(int ii=0; ii<len; ii++) {
01743       int rowi = scatterPtr[ii];
01744 
01745       //only add the non-slave rows for this column.
01746       if (rSlave_[ii] == 1) continue;
01747 
01748       Kid_->createPosition(rowi, col);
01749     }
01750     //now skip to the next column 
01751     continue;
01752   }
01753 
01754   Kdi_->createPosition(row, col);
01755       }
01756 
01757       //finally, add all slave columns in this slave row to Kdd_.
01758       for(int kk=0; kk<len; kk++) {
01759   int colk = scatterPtr[kk];
01760 
01761   if (rSlave_[kk] != 1) continue;
01762 
01763   Kdd_->createPosition(row, colk);
01764       }
01765     }
01766     else {
01767       //this is not a slave row, so we will loop across it, creating matrix
01768       //positions for all non-slave columns in it.
01769       int reducedRow = -1;
01770       translateToReducedEqn(row, reducedRow);
01771 
01772       bool rowIsLocal = true;
01773       int owner = localProc_;
01774       if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
01775   rowIsLocal = false;
01776   owner = getOwnerProcForEqn(row);
01777       }
01778 
01779       for(int j=0; j<len; j++) {
01780   if (rSlave_[j] == 1) continue;
01781 
01782   int reducedCol = workSpace_[j];
01783 
01784   if (rowIsLocal) {
01785     CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
01786           "crtSymmEqnStr") );
01787   }
01788   else {
01789     eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
01790   }
01791       }
01792     }
01793   }
01794 
01795   if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
01796 
01797   }
01798   catch(std::runtime_error& exc) {
01799     fei::console_out() << exc.what() << FEI_ENDL;
01800     ERReturn(-1);
01801   }
01802 
01803   return(FEI_SUCCESS);
01804 }
01805 
01806 //------------------------------------------------------------------------------
01807 int SNL_FEI_Structure::createBlkSymmEqnStructure(std::vector<int>& scatterIndices)
01808 {
01809   //scatterIndices are both the rows and the columns for a structurally
01810   //symmetric 2-dimensional contribution to the matrix.
01811   //
01812 
01813   //if there aren't any slaves in the whole problem, store the scatter-indices
01814   //using a fast function (no translations-to-reduced-space) and return.
01815   if (numSlaveEquations() == 0) {
01816     return( storeElementScatterBlkIndices_noSlaves(scatterIndices) );
01817   }
01818 
01819   try {
01820 
01821   int len = scatterIndices.size();
01822   bool anySlaves = false;
01823   rSlave_.resize(len);
01824   for(int is=0; is<len; is++) { 
01825     rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
01826     if (rSlave_[is] == 1) anySlaves = true;
01827   }
01828 
01829   //if there aren't any slaves in this contribution, then just store it
01830   //and return
01831   if (!anySlaves) {
01832     CHK_ERR( storeElementScatterIndices(scatterIndices) );
01833     return(FEI_SUCCESS);
01834   }
01835 
01836   int* scatterPtr = &scatterIndices[0];
01837 
01838   workSpace_.resize(len);
01839   for(int j=0; j<len; j++) {
01840     translateToReducedEqn(scatterPtr[j], workSpace_[j]);
01841   }
01842   
01843   for(int i=0; i<len; i++) {
01844     int row = scatterPtr[i];
01845     if (rSlave_[i] == 1) {
01846       reducedEqnCounter_++;
01847       //
01848       //'row' is a slave equation, so add this row to Kdi_. But as we do that,
01849       //watch for columns that are slave equations and add them to Kid_.
01850       //
01851       for(int jj=0; jj<len; jj++) {
01852   int col = scatterPtr[jj];
01853 
01854   if (rSlave_[jj] == 1) {
01855     //'col' is a slave equation, so add this column to Kid_.
01856     for(int ii=0; ii<len; ii++) {
01857       int rowi = scatterPtr[ii];
01858 
01859       //only add the non-slave rows for this column.
01860       if (rSlave_[ii] == 1) continue;
01861 
01862       Kid_->createPosition(rowi, col);
01863     }
01864     //now skip to the next column 
01865     continue;
01866   }
01867 
01868   Kdi_->createPosition(row, col);
01869       }
01870 
01871       //finally, add all slave columns in this slave row to Kdd_.
01872       for(int kk=0; kk<len; kk++) {
01873   int colk = scatterPtr[kk];
01874 
01875   if (rSlave_[kk] != 1) continue;
01876 
01877   Kdd_->createPosition(row, colk);
01878       }
01879     }
01880     else {
01881       //this is not a slave row, so we will loop across it, creating matrix
01882       //positions for all non-slave columns in it.
01883       int reducedRow = -1;
01884       translateToReducedEqn(row, reducedRow);
01885 
01886       bool rowIsLocal = true;
01887       int owner = localProc_;
01888       if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
01889   rowIsLocal = false;
01890   owner = getOwnerProcForEqn(row);
01891       }
01892 
01893       for(int j=0; j<len; j++) {
01894   if (rSlave_[j] == 1) continue;
01895 
01896   int reducedCol = workSpace_[j];
01897 
01898   if (rowIsLocal) {
01899     CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
01900           "crtSymmEqnStr") );
01901   }
01902   else {
01903     eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
01904   }
01905       }
01906     }
01907   }
01908 
01909   if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
01910 
01911   }
01912   catch(std::runtime_error& exc) {
01913     fei::console_out() << exc.what() << FEI_ENDL;
01914     ERReturn(-1);
01915   }
01916 
01917   return(FEI_SUCCESS);
01918 }
01919 
01920 //------------------------------------------------------------------------------
01921 int SNL_FEI_Structure::storeElementScatterIndices(std::vector<int>& indices)
01922 {
01923 //This function takes a list of equation numbers, as input, and for each
01924 //one, if it is a local equation number, goes to that row of the sysMatIndices
01925 //structure and stores the whole list of equation numbers as column indices
01926 //in that row of the matrix structure. If the equation number is not local,
01927 //then it is given to the EqnCommMgr.
01928 //
01929    int numIndices = indices.size();
01930    int* indPtr = &indices[0];
01931 
01932    workSpace_.resize(numIndices);
01933    int* workPtr = &workSpace_[0];
01934    int reducedEqn = -1;
01935    for(size_t j=0; j<indices.size(); j++) {
01936      bool isSlave = translateToReducedEqn(indPtr[j], reducedEqn);
01937      if (isSlave) ERReturn(-1);
01938      workPtr[j] = reducedEqn;
01939    }
01940 
01941    for(int i=0; i<numIndices; i++) {
01942       int row = indPtr[i];
01943       int rrow = workPtr[i];
01944 
01945       if ((localStartRow_ > row) || (row > localEndRow_)) {
01946   
01947   int owner = getOwnerProcForEqn(row);
01948   eqnCommMgr_->addRemoteIndices(rrow, owner, workPtr, numIndices);
01949   continue;
01950       }
01951 
01952       CHK_ERR( createMatrixPositions(rrow, numIndices, workPtr,
01953              "storeElemScttrInd") );
01954    }
01955 
01956    return(FEI_SUCCESS);
01957 }
01958 
01959 //------------------------------------------------------------------------------
01960 int SNL_FEI_Structure::storeElementScatterIndices_noSlaves(std::vector<int>& indices)
01961 {
01962   //This function takes a list of equation numbers as input and for each one,
01963   //if it is a local equation number, goes to that row of the sysMatIndices
01964   //structure and stores the whole list of equation numbers as column indices
01965   //in that row of the matrix structure. If the equation number is not local,
01966   //then it is given to the EqnCommMgr.
01967   //
01968   int i, numIndices = indices.size();
01969   int* indPtr = &indices[0];
01970 
01971   for(i=0; i<numIndices; i++) {
01972     int row = indPtr[i];
01973     if (row < localStartRow_ || row > localEndRow_) {
01974       int owner = getOwnerProcForEqn(row);
01975       eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
01976       continue;
01977     }
01978 
01979     if (generateGraph_) {
01980       fei::ctg_set<int>& thisRow =
01981   sysMatIndices_[row - reducedStartRow_];
01982 
01983       for(int j=0; j<numIndices; ++j) {
01984   if (debugOutput_ && outputLevel_ > 2) {
01985     dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
01986        << row << "," << indPtr[j] << ")"<<FEI_ENDL;
01987   }
01988 
01989   thisRow.insert2(indPtr[j]);
01990       }
01991     }
01992   }
01993 
01994   return(FEI_SUCCESS);
01995 }
01996 
01997 //------------------------------------------------------------------------------
01998 int SNL_FEI_Structure::storeElementScatterBlkIndices_noSlaves(std::vector<int>& indices)
01999 {
02000   //This function takes a list of equation numbers as input and for each one,
02001   //if it is a local equation number, goes to that row of the sysMatIndices
02002   //structure and stores the whole list of equation numbers as column indices
02003   //in that row of the matrix structure. If the equation number is not local,
02004   //then it is given to the EqnCommMgr.
02005   //
02006   int i, numIndices = indices.size();
02007   int* indPtr = &indices[0];
02008 
02009   for(i=0; i<numIndices; i++) {
02010     int row = indPtr[i];
02011     if (row < localStartRow_ || row > localEndRow_) {
02012       int owner = getOwnerProcForEqn(row);
02013       eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
02014       continue;
02015     }
02016 
02017     if (generateGraph_) {
02018       fei::ctg_set<int>& thisRow =
02019   sysMatIndices_[row - reducedStartRow_];
02020 
02021       for(int j=0; j<numIndices; ++j) {
02022   if (debugOutput_ && outputLevel_ > 2) {
02023     dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
02024        << row << "," << indPtr[j] << ")"<<FEI_ENDL;
02025   }
02026 
02027   thisRow.insert2(indPtr[j]);
02028       }
02029     }
02030   }
02031 
02032   return(FEI_SUCCESS);
02033 }
02034 
02035 //------------------------------------------------------------------------------
02036 int SNL_FEI_Structure::assembleReducedStructure()
02037 {
02038   //This function performs the appropriate manipulations (matrix-matrix-products,
02039   //matrix-vector-products) on the Kid,Kdi,Kdd matrices and then assembles the
02040   //results into the global matrix structure. This function also adds any
02041   //resulting contributions to off-processor portions of the matrix structure to
02042   //the EqnCommMgr.
02043   //
02044   fei::FillableMat* D = getSlaveDependencies();
02045 
02046   csrD = *D;
02047   csrKid = *Kid_;
02048   csrKdi = *Kdi_;
02049   csrKdd = *Kdd_;
02050 
02051   //form tmpMat1_ = Kid_*D
02052   fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
02053 
02054   //form tmpMat2_ = D^T*Kdi_
02055   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
02056 
02057   CHK_ERR( translateMatToReducedEqns(tmpMat1_) );
02058   CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
02059   if (numProcs_ > 1) {
02060     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat1_, true) );
02061     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
02062   }
02063 
02064   CHK_ERR( createMatrixPositions(tmpMat1_) );
02065   CHK_ERR( createMatrixPositions(tmpMat2_) );
02066 
02067   //form tmpMat1_ = D^T*Kdd_
02068   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
02069 
02070   //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
02071   fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
02072 
02073 
02074   CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
02075   if (numProcs_ > 1) {
02076     CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
02077   }
02078   CHK_ERR( createMatrixPositions(tmpMat2_) );
02079 
02080   Kdi_->clear();
02081   Kid_->clear();
02082   Kdd_->clear();
02083   reducedEqnCounter_ = 0;
02084 
02085   return(FEI_SUCCESS);
02086 }
02087 
02088 //------------------------------------------------------------------------------
02089 int SNL_FEI_Structure::translateToReducedEqns(EqnCommMgr& eqnCommMgr)
02090 {
02091   if (numSlvs_ < 1) return(0);
02092 
02093   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvProcEqns())) );
02094   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvEqns())) );
02095   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendProcEqns())) );
02096   CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendEqns())) );
02097 
02098   return(0);
02099 }
02100 
02101 //------------------------------------------------------------------------------
02102 int SNL_FEI_Structure::translateToReducedEqns(EqnBuffer& eqnBuf)
02103 {
02104   int numEqns = eqnBuf.getNumEqns();
02105   int* eqnNumbers = &(eqnBuf.eqnNumbers()[0]);
02106   std::vector<fei::CSVec*>& eqnArray = eqnBuf.eqns();
02107   for(int i=0; i<numEqns; ++i) {
02108     int reducedEqn;
02109     translateToReducedEqn(eqnNumbers[i], reducedEqn);
02110     eqnNumbers[i] = reducedEqn;
02111 
02112     int* indicesPtr = &(eqnArray[i]->indices()[0]);
02113     int numIndices = eqnArray[i]->size();
02114     for(int j=0; j<numIndices; ++j) {
02115       translateToReducedEqn(indicesPtr[j], reducedEqn);
02116       indicesPtr[j] = reducedEqn;
02117     }
02118   }
02119 
02120   return(0);
02121 }
02122 
02123 //------------------------------------------------------------------------------
02124 int SNL_FEI_Structure::translateToReducedEqns(ProcEqns& procEqns)
02125 {
02126   std::vector<std::vector<int>*>& eqnTable = procEqns.procEqnNumbersPtr();
02127   int dim1 = eqnTable.size();
02128   for(int i=0; i<dim1; ++i) {
02129     int dim2 = eqnTable[i]->size();
02130     int* indicesPtr = &(*(eqnTable[i]))[0];
02131     for(int j=0; j<dim2; ++j) {
02132       int reducedEqn;
02133       translateToReducedEqn(indicesPtr[j], reducedEqn);
02134       indicesPtr[j] = reducedEqn;
02135     }
02136   }
02137 
02138   return(0);
02139 }
02140 
02141 //------------------------------------------------------------------------------
02142 int SNL_FEI_Structure::translateMatToReducedEqns(fei::CSRMat& mat)
02143 {
02144   //Given a matrix in unreduced numbering, convert all of its contents to the
02145   //"reduced" equation space. If any of the row-numbers or column-indices in
02146   //this matrix object are slave equations, they will not be referenced. In
02147   //this case, a positive warning-code will be returned.
02148 
02149   int reducedEqn = -1;
02150   int foundSlave = 0;
02151 
02152   fei::SparseRowGraph& srg = mat.getGraph();
02153 
02154   std::vector<int>& rowNumbers = srg.rowNumbers;
02155   for(size_t i=0; i<rowNumbers.size(); ++i) {
02156     bool isSlave = translateToReducedEqn(rowNumbers[i], reducedEqn);
02157     if (isSlave) foundSlave = 1;
02158     else rowNumbers[i] = reducedEqn;
02159   }
02160 
02161   std::vector<int>& colIndices = srg.packedColumnIndices;
02162   for(size_t i=0; i<colIndices.size(); ++i) {
02163     bool isSlave = translateToReducedEqn(colIndices[i], reducedEqn);
02164     if (isSlave) foundSlave = 1;
02165     else colIndices[i] = reducedEqn;
02166   }
02167 
02168   return(foundSlave);
02169 }
02170 
02171 //------------------------------------------------------------------------------
02172 int SNL_FEI_Structure::createMatrixPositions(fei::CSRMat& mat)
02173 {
02174   if (!generateGraph_) {
02175     return(0);
02176   }
02177 
02178   //This function must be called with a CSRMat object that already has its
02179   //contents numbered in "reduced" equation-numbers.
02180   //
02181   //This function has one simple task -- for each row,col pair stored in 'mat',
02182   //call 'createMatrixPosition' to make an entry in the global matrix structure
02183   //if it doesn't already exist.
02184   //
02185   const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
02186   const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
02187   const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
02188 
02189   for(size_t i=0; i<rowNumbers.size(); ++i) {
02190     int row = rowNumbers[i];
02191     int offset = rowOffsets[i];
02192     int rowlen = rowOffsets[i+1]-offset;
02193     const int* indices = &pckColInds[offset];
02194 
02195     for(int j=0; j<rowlen; j++) {
02196       CHK_ERR( createMatrixPosition(row, indices[j], "crtMatPos(CSRMat)") );
02197     }
02198   }
02199 
02200   return(0);
02201 }
02202 
02203 //------------------------------------------------------------------------------
02204 int SNL_FEI_Structure::createMatrixPosition(int row, int col,
02205               const char* callingFunction)
02206 {
02207   if (!generateGraph_) {
02208     return(0);
02209   }
02210 
02211   //
02212   //This function inserts the column index 'col' in the system matrix structure
02213   //in 'row', if it isn't already there.
02214   //
02215   //This is a private SNL_FEI_Structure function. These row/col numbers must be
02216   //global, 0-based equation numbers in the "reduced" equation space..
02217   //
02218 
02219   //if the equation isn't local, just return. (The assumption being that
02220   //this call is also being made on the processor that owns the equation.)
02221   if (row < reducedStartRow_ || row > reducedEndRow_) {
02222     if (debugOutput_) {
02223       dbgOut() << "         " << callingFunction << " passed " <<row<<","<<col
02224          << ", not local." << FEI_ENDL;
02225     }
02226 
02227     return(0);
02228   }
02229 
02230   if (debugOutput_ && outputLevel_ > 2) {
02231     dbgOut() << "#         " << callingFunction << " crtMatPos("
02232        << row << "," << col << ")"<<FEI_ENDL;
02233   }
02234 
02235   fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
02236 
02237   thisRow.insert2(col);
02238 
02239   return(0);
02240 }
02241 
02242 //------------------------------------------------------------------------------
02243 int SNL_FEI_Structure::createMatrixPositions(int row, int numCols, int* cols,
02244               const char* callingFunction)
02245 {
02246   if (!generateGraph_) {
02247     return(0);
02248   }
02249 
02250   //
02251   //This function inserts the column indices 'cols' in the system matrix structure
02252   //in 'row', if they aren't already there.
02253   //
02254   //This is a private SNL_FEI_Structure function. These row/col numbers must be
02255   //global, 0-based equation numbers in the "reduced" equation space..
02256   //
02257 
02258   //if the equation isn't local, just return. (The assumption being that
02259   //this call is also being made on the processor that owns the equation.)
02260   if (row < reducedStartRow_ || row > reducedEndRow_) {
02261     if (debugOutput_) {
02262       dbgOut() << "         " << callingFunction << " passed " <<row
02263          << ", not local." << FEI_ENDL;
02264     }
02265 
02266     return(0);
02267   }
02268 
02269   fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
02270 
02271   for(int i=0; i<numCols; i++) {
02272     if (debugOutput_ && outputLevel_ > 2) {
02273       dbgOut() << "#         " << callingFunction << " crtMatPoss("
02274          << row << "," << cols[i] << ")"<<FEI_ENDL;
02275     }
02276 
02277     thisRow.insert2(cols[i]);
02278   }
02279 
02280   return(0);
02281 }
02282 
02283 //------------------------------------------------------------------------------
02284 int SNL_FEI_Structure::initializeBlkEqnMapper()
02285 {
02286   //Run the nodes, elem-dof's and multiplier constraints, establishing the
02287   //point-equation to block-equation mappings.  Note that this data is 
02288   //initialized as non-reduced (i.e., not accounting for any slave equations)
02289   //equation numbers, then translated to the reduced space at the end of this
02290   //function.
02291 
02292   if (blkEqnMapper_ != NULL) delete blkEqnMapper_;
02293   blkEqnMapper_ = new snl_fei::PointBlockMap();
02294   snl_fei::PointBlockMap& blkEqnMapper = getBlkEqnMapper();
02295 
02296   //First, if there's only one (non-negative) fieldID and that fieldID's size
02297   //is 1, then we don't need to execute the rest of this function.
02298   int numFields = fieldDatabase_->size();
02299   const int* fieldIDs = getFieldIDsPtr();
02300   const int* fieldSizes = fieldIDs+numFields;
02301   int numNonNegativeFields = 0;
02302   int maxFieldSize = 0;
02303   for(int f=0; f<numFields; f++) {
02304     if (fieldIDs[f] >= 0) {
02305       numNonNegativeFields++;
02306       if (fieldSizes[f] > maxFieldSize) maxFieldSize = fieldSizes[f];
02307     }
02308   }
02309 
02310   if (numNonNegativeFields == 1 && maxFieldSize == 1) {
02311     globalMaxBlkSize_ = 1;
02312     blkEqnMapper.setPtEqualBlk();
02313     return(0);
02314   }
02315 
02316   int localVanishedNodeAdjustment = 0;
02317   for(int i=0; i<localProc_; ++i) {
02318     localVanishedNodeAdjustment += globalNumNodesVanished_[i];
02319   }
02320 
02321   int eqnNumber, blkEqnNumber;
02322 
02323   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02324   std::map<GlobalID,int>::iterator
02325     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02326 
02327   for(; iter!=iter_end; ++iter) {
02328     NodeDescriptor* node = NULL;
02329     nodeDatabase_->getNodeAtIndex(iter->second, node);
02330 
02331     // If the node doesn't exist, skip.
02332     if (node==NULL) continue;
02333 
02334     // If the node doesn't have any fields, skip
02335     int numFields  = node->getNumFields();
02336     if(numFields == 0) continue;
02337 
02338     int firstPtEqn = node->getFieldEqnNumbers()[0];
02339     int numNodalDOF = node->getNumNodalDOF();
02340     blkEqnNumber = node->getBlkEqnNumber() - localVanishedNodeAdjustment;
02341 
02342     int blkSize = numNodalDOF;
02343 
02344     for(int j=0; j<numNodalDOF; ++j) {
02345       bool isSlave = translateToReducedEqn(firstPtEqn+j, eqnNumber);
02346       if (isSlave) {
02347   --blkSize;
02348       }
02349       else {
02350   CHK_ERR( blkEqnMapper.setEqn(eqnNumber, blkEqnNumber) );
02351       }
02352     }
02353 
02354     if (blkSize > 0) {
02355       CHK_ERR( blkEqnMapper.setBlkEqnSize(blkEqnNumber, blkSize) );
02356     }
02357     else {
02358       ++localVanishedNodeAdjustment;
02359     }
02360   }
02361 
02362   //
02363   //Now the element-dofs...
02364   //
02365   int numBlocks = getNumElemBlocks();
02366   int numLocalNodes = globalNodeOffsets_[localProc_+1]-globalNodeOffsets_[localProc_];
02367   eqnNumber = localStartRow_ + numLocalNodalEqns_;
02368   blkEqnNumber = localBlkOffset_ + numLocalNodes;
02369 
02370   for(int i = 0; i < numBlocks; i++) {
02371     BlockDescriptor* block = NULL;
02372     CHK_ERR( getBlockDescriptor_index(i, block) );
02373 
02374     int numElemDOF = block->getNumElemDOFPerElement();
02375 
02376     if (numElemDOF > 0) {
02377       int numElems = block->getNumElements();
02378 
02379       for(int j=0; j<numElems; j++) {
02380 
02381   for(int ede=0; ede<numElemDOF; ede++) {
02382     blkEqnMapper.setEqn(eqnNumber+ede, blkEqnNumber+ede);
02383     blkEqnMapper.setBlkEqnSize(blkEqnNumber+ede, 1);
02384   }
02385 
02386   eqnNumber += numElemDOF;
02387   blkEqnNumber += numElemDOF;
02388       }
02389     }
02390   }
02391 
02392   //Now we'll set mappings for all local multiplier constraint-relations,
02393   //and also gather the equation numbers for other processors' multipliers
02394   //to record in our snl_fei::PointBlockMap object.
02395   //
02396   std::map<GlobalID,ConstraintType*>::const_iterator
02397     cr_iter = multCRs_.begin(),
02398     cr_end  = multCRs_.end();
02399 
02400   std::vector<int> localMultEqns;
02401   localMultEqns.reserve(multCRs_.size()*2);
02402   while(cr_iter != cr_end) {
02403     ConstraintType* multCR = (*cr_iter).second;
02404     eqnNumber = multCR->getEqnNumber();
02405     blkEqnNumber = multCR->getBlkEqnNumber();
02406 
02407     CHK_ERR( blkEqnMapper_->setEqn(eqnNumber, blkEqnNumber) );
02408     CHK_ERR( blkEqnMapper_->setBlkEqnSize(blkEqnNumber, 1) );
02409 
02410     localMultEqns.push_back(eqnNumber);
02411     localMultEqns.push_back(blkEqnNumber);
02412     ++cr_iter;
02413   }
02414 
02415   //Now gather and store the equation-numbers arising from other
02416   //processors' multiplier constraints. (We obviously only need to do this if
02417   //there is more than one processor, and we can skip it if the problem
02418   //only contains 1 scalar equation for each block-equation.)
02419 
02420   int localMaxBlkEqnSize = blkEqnMapper_->getMaxBlkEqnSize();
02421   globalMaxBlkSize_ = 0;
02422   CHK_ERR( fei::GlobalMax(comm_, localMaxBlkEqnSize, globalMaxBlkSize_) );
02423 
02424   blkEqnMapper_->setMaxBlkEqnSize(globalMaxBlkSize_);
02425 
02426   if (globalMaxBlkSize_ == 1) {
02427     //If the globalMax block-size is 1 then tell the blkEqnMapper, and it will
02428     //reclaim the memory it allocated in arrays, etc.
02429     blkEqnMapper_->setPtEqualBlk();
02430     return(0);
02431   }
02432 
02433   if (numProcs_ > 1) {
02434     std::vector<int> recvLengths, globalMultEqns;
02435     CHK_ERR(fei::Allgatherv(comm_, localMultEqns, recvLengths, globalMultEqns));
02436 
02437     int offset = 0;
02438     for(int p=0; p<numProcs_; p++) {
02439       if (p == localProc_) { offset += recvLengths[p]; continue; }
02440 
02441       for(int j=0; j<recvLengths[p]/2; j++) {
02442   blkEqnMapper_->setEqn(globalMultEqns[offset], globalMultEqns[offset+1]);
02443   blkEqnMapper_->setBlkEqnSize(globalMultEqns[offset+1], 1);
02444   offset += 2;
02445       }
02446     }
02447   }
02448 
02449   return(0);
02450 }
02451 
02452 //------------------------------------------------------------------------------
02453 int SNL_FEI_Structure::setMultCREqnInfo()
02454 {
02455   //We'll set equation-numbers on all local multiplier constraint-relations,
02456   //and also gather the equation numbers for other processors' multipliers
02457   //to record in our snl_fei::PointBlockMap object.
02458   //
02459   std::map<GlobalID,ConstraintType*>::const_iterator
02460     cr_iter = multCRs_.begin(),
02461     cr_end  = multCRs_.end();
02462 
02463   int eqnNumber = localStartRow_ + numLocalNodalEqns_ + numLocalElemDOF_;
02464   int blkEqnNum = localBlkOffset_ + numLocalEqnBlks_ - numLocalMultCRs_;
02465 
02466   while(cr_iter != cr_end) {
02467     ConstraintType* multCR = (*cr_iter).second;
02468 
02469     multCR->setEqnNumber(eqnNumber);
02470 
02471     multCR->setBlkEqnNumber(blkEqnNum);
02472 
02473     ++eqnNumber;
02474     ++blkEqnNum;
02475     ++cr_iter;
02476   }
02477 
02478   return(0);
02479 }
02480 
02481 //------------------------------------------------------------------------------
02482 int SNL_FEI_Structure::writeEqn2NodeMap()
02483 {
02484   FEI_OSTRINGSTREAM osstr;
02485   osstr << dbgPath_ << "/FEI" << name_ << "_nID_nNum_blkEq_ptEq."
02486   << numProcs_ << "." << localProc_;
02487 
02488   FEI_OFSTREAM e2nFile(osstr.str().c_str());
02489 
02490   if (!e2nFile) {
02491     fei::console_out() << "SNL_FEI_Structure::writeEqn2NodeMap: ERROR, couldn't open file "
02492          << osstr.str() << FEI_ENDL;
02493     ERReturn(-1);
02494   }
02495 
02496   std::map<GlobalID,ConstraintType*>::const_iterator
02497     cr_iter = multCRs_.begin(),
02498     cr_end  = multCRs_.end();
02499 
02500   while(cr_iter != cr_end) {
02501     ConstraintType* multCR = (*cr_iter).second;
02502     int eqnNumber = multCR->getEqnNumber();
02503     int blkEqnNumber = multCR->getBlkEqnNumber();
02504 
02505     e2nFile << "-999 -999 " << blkEqnNumber << " " << eqnNumber << FEI_ENDL;
02506     ++cr_iter;
02507   }
02508 
02509   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02510   std::map<GlobalID,int>::iterator
02511     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02512 
02513   for(; iter!=iter_end; ++iter) {
02514     NodeDescriptor* node = NULL;
02515     nodeDatabase_->getNodeAtIndex(iter->second, node);
02516 
02517     if (node==NULL || node->getOwnerProc() != localProc_) {
02518       continue;
02519     }
02520 
02521     const int* fieldIDList = node->getFieldIDList();
02522     const int* eqnNumbers = node->getFieldEqnNumbers();
02523     int nodeNum = node->getNodeNumber();
02524     int blkEqnNumber = node->getBlkEqnNumber();
02525     int numFields = node->getNumFields();
02526 
02527     for(int j=0; j<numFields; j++) {
02528       int numFieldParams = getFieldSize(fieldIDList[j]);
02529       assert( numFieldParams > 0 );
02530 
02531       for(int k=0; k<numFieldParams; k++) {
02532         e2nFile << (int)node->getGlobalNodeID() << " " << nodeNum
02533                 << " " << blkEqnNumber << " " << eqnNumbers[j]+k<<FEI_ENDL;
02534       }
02535     }
02536   }
02537 
02538   return(FEI_SUCCESS);
02539 }
02540 
02541 //------------------------------------------------------------------------------
02542 int SNL_FEI_Structure::calcTotalNumElemDOF()
02543 {
02544    int numBlocks = getNumElemBlocks();
02545    int totalNumElemDOF = 0;
02546 
02547    for(int i = 0; i < numBlocks; i++) {
02548      BlockDescriptor* block = NULL;
02549      CHK_ERR( getBlockDescriptor_index(i, block) );
02550       int numElemDOF = block->getNumElemDOFPerElement();
02551       if (numElemDOF > 0) {
02552          int numElems = block->getNumElements();
02553 
02554          totalNumElemDOF += numElems*numElemDOF;
02555       }
02556    }
02557 
02558    return(totalNumElemDOF);
02559 }
02560 
02561 //------------------------------------------------------------------------------
02562 int SNL_FEI_Structure::calcNumMultCREqns()
02563 {
02564   int totalNumMultCRs = 0;
02565   int numMCRecords = getNumMultConstRecords();
02566 
02567   totalNumMultCRs = numMCRecords;
02568 
02569   return( totalNumMultCRs );
02570 }
02571 
02572 //------------------------------------------------------------------------------
02573 void SNL_FEI_Structure::calcGlobalEqnInfo(int numLocallyOwnedNodes, 
02574            int numLocalEqns,
02575            int numLocalEqnBlks) 
02576 {
02577   FEI_OSTREAM& os = dbgOut();
02578   if (debugOutput_) os << "#  entering calcGlobalEqnInfo" << FEI_ENDL;
02579 
02580 //
02581 //This function does the inter-process communication necessary to obtain,
02582 //on each processor, the global number of equations, and the local starting
02583 //and ending equation numbers.
02584 //While we're here, we do the same thing for node-numbers.
02585 //
02586 
02587   //first, get each processor's local number of equations on the master proc.
02588 
02589   std::vector<int> numProcNodes(numProcs_);
02590   std::vector<int> numProcEqns(numProcs_);
02591   std::vector<int> numProcEqnBlks(numProcs_);
02592 
02593   if (numProcs_ > 1) {
02594 #ifndef FEI_SER
02595     int glist[3];
02596     std::vector<int> recvList(3*numProcs_);
02597 
02598     glist[0] = numLocallyOwnedNodes;
02599     glist[1] = numLocalEqns;
02600     glist[2] = numLocalEqnBlks;
02601     if (MPI_Gather(&glist[0], 3, MPI_INT, &recvList[0], 3,
02602        MPI_INT, masterProc_, comm_) != MPI_SUCCESS) {
02603       fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Gather" << FEI_ENDL;
02604       MPI_Abort(comm_, -1);
02605     }
02606 
02607     for(int p=0; p<numProcs_; p++) {
02608       numProcNodes[p]   = recvList[p*3];
02609       numProcEqns[p]    = recvList[p*3+1];
02610       numProcEqnBlks[p] = recvList[p*3+2];
02611     }
02612 #endif
02613   }
02614   else {
02615     numProcNodes[0] = numLocallyOwnedNodes;
02616     numProcEqns[0] = numLocalEqns;
02617     numProcEqnBlks[0] = numLocalEqnBlks;
02618   }
02619 
02620   //compute offsets for all processors (starting index for local equations)
02621 
02622   globalNodeOffsets_.resize(numProcs_+1);
02623   globalEqnOffsets_.resize(numProcs_+1);
02624   globalBlkEqnOffsets_.resize(numProcs_+1);
02625 
02626   globalNodeOffsets_[0] = 0; // 0-based node-numbers
02627   globalEqnOffsets_[0] = 0;    // we're going to start rows & cols at 0 (global)
02628   globalBlkEqnOffsets_[0] = 0;
02629 
02630   if (localProc_ == masterProc_) {
02631     for(int i=1;i<numProcs_;i++) {
02632       globalNodeOffsets_[i] = globalNodeOffsets_[i-1] +numProcNodes[i-1];
02633       globalEqnOffsets_[i] = globalEqnOffsets_[i-1] + numProcEqns[i-1];
02634       globalBlkEqnOffsets_[i] = globalBlkEqnOffsets_[i-1] +
02635                                   numProcEqnBlks[i-1];
02636     }
02637 
02638     globalNodeOffsets_[numProcs_] = globalNodeOffsets_[numProcs_-1] +
02639                                    numProcNodes[numProcs_-1];
02640     globalEqnOffsets_[numProcs_] = globalEqnOffsets_[numProcs_-1] +
02641                                   numProcEqns[numProcs_-1];
02642     globalBlkEqnOffsets_[numProcs_] = globalBlkEqnOffsets_[numProcs_-1] +
02643                                      numProcEqnBlks[numProcs_-1];
02644   }
02645 
02646   //now, broadcast vector of eqnOffsets to all processors
02647 
02648   if (numProcs_ > 1) {
02649 #ifndef FEI_SER
02650     std::vector<int> blist(3*(numProcs_+1));
02651 
02652     if (localProc_ == masterProc_) {
02653       int offset = 0;
02654       for(int i=0; i<numProcs_+1; i++) {
02655   blist[offset++] = globalNodeOffsets_[i];
02656   blist[offset++] = globalEqnOffsets_[i];
02657   blist[offset++] = globalBlkEqnOffsets_[i];
02658       }
02659     }
02660 
02661     if (MPI_Bcast(&blist[0], 3*(numProcs_+1), MPI_INT,
02662       masterProc_, comm_) != MPI_SUCCESS) {
02663       fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Bcast" << FEI_ENDL;
02664       MPI_Abort(comm_, -1);
02665     }
02666 
02667     if (localProc_ != masterProc_) {
02668       int offset = 0;
02669       for(int i=0; i<numProcs_+1; i++) {
02670   globalNodeOffsets_[i] = blist[offset++];
02671   globalEqnOffsets_[i] = blist[offset++];
02672   globalBlkEqnOffsets_[i] = blist[offset++];
02673       }
02674     }
02675 #endif
02676   }
02677 
02678   localStartRow_ = globalEqnOffsets_[localProc_];
02679   localEndRow_ = globalEqnOffsets_[localProc_+1]-1;
02680   numGlobalEqns_ = globalEqnOffsets_[numProcs_];
02681 
02682   firstLocalNodeNumber_ = globalNodeOffsets_[localProc_];
02683   numGlobalNodes_ = globalNodeOffsets_[numProcs_];
02684 
02685   localBlkOffset_ = globalBlkEqnOffsets_[localProc_];
02686   numGlobalEqnBlks_ = globalBlkEqnOffsets_[numProcs_];
02687 
02688   if (debugOutput_) {
02689     os << "#     localStartRow_: " << localStartRow_ << FEI_ENDL;
02690     os << "#     localEndRow_: "   << localEndRow_ << FEI_ENDL;
02691     os << "#     numGlobalEqns_: " << numGlobalEqns_ << FEI_ENDL;
02692     os << "#     numGlobalEqnBlks_: " << numGlobalEqnBlks_ << FEI_ENDL;
02693     os << "#     localBlkOffset_: " << localBlkOffset_ << FEI_ENDL;
02694     os << "#     firstLocalNodeNumber_: " << firstLocalNodeNumber_ << FEI_ENDL;
02695     size_t n;
02696     os << "#     numGlobalNodes_: " << numGlobalNodes_ << FEI_ENDL;
02697     os << "#     " << globalNodeOffsets_.size() << " globalNodeOffsets_: ";
02698     for(size_t nn=0; nn<globalNodeOffsets_.size(); nn++) 
02699       os << globalNodeOffsets_[nn] << " ";
02700     os << FEI_ENDL;
02701     os << "#     " << globalEqnOffsets_.size() << " globalEqnOffsets_: ";
02702     for(n=0; n<globalEqnOffsets_.size(); n++) 
02703       os << globalEqnOffsets_[n] << " ";
02704     os << FEI_ENDL;
02705     os << "#     " << globalBlkEqnOffsets_.size() << " globalBlkEqnOffsets_: ";
02706     for(n=0; n<globalBlkEqnOffsets_.size(); n++) 
02707       os << globalBlkEqnOffsets_[n] << " ";
02708     os << FEI_ENDL;
02709     os << "#  leaving calcGlobalEqnInfo" << FEI_ENDL;
02710   }
02711 }
02712 
02713 //------------------------------------------------------------------------------
02714 int SNL_FEI_Structure::setNodalEqnInfo()
02715 {
02716 //
02717 //Private SNL_FEI_Structure function, to be called in initComplete, after
02718 //'calcGlobalEqnInfo()'.
02719 //
02720 //This function sets global equation numbers on the nodes.
02721 //
02722 //The return value is an error-code.
02723 //
02724   int eqn = localStartRow_;
02725 
02726   //localBlkOffset_ is 0-based, and so is blkEqnNumber.
02727   int blkEqnNumber = localBlkOffset_;
02728 
02729   std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
02730   std::map<GlobalID,int>::iterator
02731     iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
02732 
02733   for(; iter!=iter_end; ++iter) {
02734     NodeDescriptor* node = NULL;
02735     nodeDatabase_->getNodeAtIndex(iter->second, node);
02736 
02737     if (node==NULL) continue;
02738 
02739     int numFields = node->getNumFields();
02740     const int* fieldIDs = node->getFieldIDList();
02741 
02742     if (node->getOwnerProc() == localProc_) {
02743       int numNodalDOF = 0;
02744 
02745       for(int j=0; j<numFields; j++) {
02746   node->setFieldEqnNumber(fieldIDs[j], eqn);
02747   int fieldSize = getFieldSize(fieldIDs[j]);
02748   eqn += fieldSize;
02749   numNodalDOF += fieldSize;
02750       }
02751 
02752       node->setNumNodalDOF(numNodalDOF);
02753       node->setBlkEqnNumber(blkEqnNumber++);
02754     }
02755   }
02756 
02757   return(0);
02758 }
02759 
02760 //------------------------------------------------------------------------------
02761 void SNL_FEI_Structure::setElemDOFEqnInfo()
02762 {
02763   //
02764   //Private SNL_FEI_Structure function, to be called from initComplete after
02765   //setBlkEqnInfo() has been called.
02766   //
02767   //Sets global equation numbers for all element-DOF.
02768   //
02769   int numBlocks = getNumElemBlocks();
02770 
02771   int eqnNumber = localStartRow_ + numLocalNodalEqns_;
02772 
02773   for(int i = 0; i < numBlocks; i++) {
02774     BlockDescriptor* block = NULL;
02775     int err = getBlockDescriptor_index(i, block);
02776     if (err) voidERReturn;
02777 
02778     int numElemDOF = block->getNumElemDOFPerElement();
02779 
02780     if (numElemDOF > 0) {
02781       std::vector<int>& elemDOFEqns = block->elemDOFEqnNumbers();
02782       std::map<GlobalID,int>& elemIDs = connTables_[i]->elemIDs;
02783 
02784       std::map<GlobalID,int>::const_iterator
02785         iter = elemIDs.begin(),
02786         iter_end = elemIDs.end();
02787 
02788       for(; iter != iter_end; ++iter) {
02789   elemDOFEqns[iter->second] = eqnNumber;
02790 
02791   eqnNumber += numElemDOF;
02792       }
02793     }
02794   }
02795 }
02796 
02797 //------------------------------------------------------------------------------
02798 int SNL_FEI_Structure::addBlock(GlobalID blockID) {
02799 //
02800 //Append a blockID to our (unsorted) list of blockIDs if it isn't
02801 //already present. Also, add a block-descriptor if blockID wasn't
02802 //already present, and a ConnectivityTable to go with it.
02803 //
02804   int insertPoint = -1;
02805   int found = fei::binarySearch(blockID, blockIDs_, insertPoint);
02806 
02807    if (found<0) {
02808       blockIDs_.insert(blockIDs_.begin()+insertPoint, blockID);
02809 
02810       BlockDescriptor* block = new BlockDescriptor;
02811       block->setGlobalBlockID(blockID);
02812 
02813       blocks_.insert(blocks_.begin()+insertPoint, block);
02814 
02815       ConnectivityTable* newConnTable = new ConnectivityTable;
02816       connTables_.insert(connTables_.begin()+insertPoint, newConnTable);
02817    }
02818 
02819    return(FEI_SUCCESS);
02820 }
02821 
02822 //------------------------------------------------------------------------------
02823 int SNL_FEI_Structure::getBlockDescriptor(GlobalID blockID,
02824                                          BlockDescriptor*& block)
02825 {
02826   int index = fei::binarySearch(blockID, blockIDs_);
02827 
02828    if (index < 0) {
02829       fei::console_out() << "SNL_FEI_Structure::getBlockDescriptor: ERROR, blockID "
02830            << (int)blockID << " not found." << FEI_ENDL;
02831       return(-1);
02832    }
02833 
02834    block = blocks_[index];
02835 
02836    return(FEI_SUCCESS);
02837 }
02838 
02839 //------------------------------------------------------------------------------
02840 int SNL_FEI_Structure::getIndexOfBlock(GlobalID blockID)
02841 {
02842   int index = fei::binarySearch(blockID, blockIDs_);
02843   return(index);
02844 }
02845 
02846 //------------------------------------------------------------------------------
02847 int SNL_FEI_Structure::getBlockDescriptor_index(int index,
02848                  BlockDescriptor*& block)
02849 {
02850   if (index < 0 || index >= (int)blockIDs_.size()) ERReturn(FEI_FATAL_ERROR);
02851 
02852   block = blocks_[index];
02853 
02854   return(FEI_SUCCESS);
02855 }
02856 
02857 //------------------------------------------------------------------------------
02858 int SNL_FEI_Structure::allocateBlockConnectivity(GlobalID blockID) {
02859 
02860    int index = fei::binarySearch(blockID, blockIDs_);
02861 
02862    if (index < 0) {
02863       fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, blockID "
02864            << (int)blockID << " not found. Aborting." << FEI_ENDL;
02865       MPI_Abort(comm_, -1);
02866    }
02867 
02868    connTables_[index]->numNodesPerElem = 
02869                          blocks_[index]->getNumNodesPerElement();
02870 
02871    int numRows = blocks_[index]->getNumElements();
02872    int numCols = connTables_[index]->numNodesPerElem;
02873 
02874    if ((numRows <= 0) || (numCols <= 0)) {
02875       fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, either "
02876            << "numElems or numNodesPerElem not yet set for blockID "
02877            << (int)blockID << ". Aborting." << FEI_ENDL;
02878       MPI_Abort(comm_, -1);
02879    }
02880 
02881    connTables_[index]->elemNumbers.resize(numRows);
02882 
02883    connTables_[index]->elem_conn_ids = new std::vector<GlobalID>(numRows*numCols);
02884 
02885    return(FEI_SUCCESS);
02886 }
02887 
02888 //------------------------------------------------------------------------------
02889 ConnectivityTable& SNL_FEI_Structure::getBlockConnectivity(GlobalID blockID) {
02890 
02891    int index = fei::binarySearch(blockID, blockIDs_);
02892 
02893    if (index < 0) {
02894       fei::console_out() << "SNL_FEI_Structure::getBlockConnectivity: ERROR, blockID "
02895            << (int)blockID << " not found. Aborting." << FEI_ENDL;
02896       MPI_Abort(comm_, -1);
02897    }  
02898 
02899    return(*(connTables_[index]));
02900 }
02901 
02902 //------------------------------------------------------------------------------
02903 int SNL_FEI_Structure::finalizeActiveNodes()
02904 {
02905   FEI_OSTREAM& os = dbgOut();
02906   if (debugOutput_) {
02907     os << "#  finalizeActiveNodes" << FEI_ENDL;
02908   }
02909   //First, make sure that the shared-node IDs are in the nodeDatabase, since
02910   //they might not have been added yet...
02911 
02912   std::vector<GlobalID>& shNodeIDs = nodeCommMgr_->getSharedNodeIDs();
02913   for(unsigned i=0; i<shNodeIDs.size(); i++) {
02914     CHK_ERR( nodeDatabase_->initNodeID(shNodeIDs[i]) );
02915   }
02916 
02917   if (activeNodesInitialized_) return(0);
02918   //
02919   //Now run through the connectivity tables and set the nodal field and block
02920   //info. Also, replace the nodeID in the element-connectivity lists with an 
02921   //index from the NodeDatabase object. This will save us a lot of time when
02922   //looking up nodes later.
02923   //
02924   //One other thing we're going to do is give all elements an element-number.
02925   //This element-number will start at zero locally (meaning that it's not
02926   //globally unique).
02927   //
02928   int elemNumber = 0;
02929   int numBlocks = blockIDs_.size();
02930   for(int bIndex=0; bIndex<numBlocks; bIndex++) {
02931 
02932     GlobalID blockID = blocks_[bIndex]->getGlobalBlockID();
02933     ConnectivityTable& conn = *(connTables_[bIndex]);
02934 
02935     GlobalID* elemConn = NULL;
02936     NodeDescriptor** elemNodeDescPtrs = NULL;
02937 
02938     int numElems = conn.elemIDs.size();
02939     if (numElems > 0) {
02940       elemConn = &((*conn.elem_conn_ids)[0]);
02941       if (!activeNodesInitialized_) {
02942   int elemConnLen = conn.elem_conn_ids->size();
02943   conn.elem_conn_ptrs = new std::vector<NodeDescriptor*>(elemConnLen);
02944       }
02945       elemNodeDescPtrs = &((*conn.elem_conn_ptrs)[0]);
02946     }
02947     int nodesPerElem = conn.numNodesPerElem;
02948 
02949     int* fieldsPerNodePtr = blocks_[bIndex]->fieldsPerNodePtr();
02950     int** fieldIDsTablePtr = blocks_[bIndex]->fieldIDsTablePtr();
02951 
02952     //
02953     //Now we want to go through the connectivity table, and for each node,
02954     //add its fields and this block to the correct NodeDescriptor from the
02955     //nodeDatabase_.
02956     //(Note that the addField and addBlock functions only add the input if
02957     //it isn't already present in that NodeDescriptor.)
02958     //
02959     //We also need to set its owner proc to localProc_ as a default, and 
02960     //inform the nodeCommMgr that the node appears in local connectivities.
02961     //
02962 
02963     conn.elemNumbers.resize(numElems);
02964 
02965     int numDistinctFields = blocks_[bIndex]->getNumDistinctFields();
02966     int fieldID = fieldIDsTablePtr[0][0];
02967 
02968     int elem;
02969     for(elem=0; elem<numElems; elem++) {
02970       conn.elemNumbers[elem] = elemNumber++;
02971       GlobalID* elemNodes = &(elemConn[elem*nodesPerElem]);
02972       NodeDescriptor** elemNodePtrs = &(elemNodeDescPtrs[elem*nodesPerElem]);
02973 
02974       for(int n=0; n<nodesPerElem; n++) {
02975   NodeDescriptor* node = NULL;
02976   int index = nodeDatabase_->getNodeWithID( elemNodes[n], node );
02977   if (index < 0) {
02978     fei::console_out() << "ERROR in SNL_FEI_Structure::initializeActiveNodes, "
02979          << FEI_ENDL << "failed to find node "
02980          << (int)(elemNodes[n]) << FEI_ENDL;
02981   }
02982 
02983   if (numDistinctFields == 1) {
02984     node->addField(fieldID);
02985   }
02986   else {
02987     for(int i=0; i<fieldsPerNodePtr[n]; i++) {
02988       node->addField(fieldIDsTablePtr[n][i]);
02989     }
02990   }
02991 
02992   node->addBlock(blockID);
02993 
02994   //now store this NodeDescriptor pointer, for fast future lookups
02995   elemNodePtrs[n] = node;
02996 
02997   node->setOwnerProc(localProc_);
02998   if (numProcs_ > 1) nodeCommMgr_->informLocal(*node);
02999       }
03000     }
03001   }
03002 
03003   //
03004   //we also need to run through the penalty constraints and inform the
03005   //nodeCommMgr that the constrained nodes are local.
03006   //
03007 
03008   if (numProcs_ > 1) {
03009     std::map<GlobalID,ConstraintType*>::const_iterator
03010       cr_iter = getPenConstRecords().begin(),
03011       cr_end = getPenConstRecords().end();
03012 
03013     while(cr_iter != cr_end) {
03014       ConstraintType& cr = *((*cr_iter).second);
03015       std::vector<GlobalID>& nodeID_vec = cr.getMasters();
03016       GlobalID* nodeIDs = &nodeID_vec[0];
03017       int numNodes = cr.getMasters().size();
03018 
03019       NodeDescriptor* node = NULL;
03020       for(int k=0; k<numNodes; ++k) {
03021    int err = nodeDatabase_->getNodeWithID(nodeIDs[k], node) ;
03022          if ( err != -1 ) {
03023            nodeCommMgr_->informLocal(*node);
03024          }
03025   //CHK_ERR( nodeDatabase_->getNodeWithID(nodeIDs[k], node) );
03026   //CHK_ERR( nodeCommMgr_->informLocal(*node) );
03027       }
03028       ++cr_iter;
03029     }
03030   }
03031 
03032   activeNodesInitialized_ = true;
03033 
03034   if (debugOutput_) os << "#  leaving finalizeActiveNodes" << FEI_ENDL;
03035   return(FEI_SUCCESS);
03036 }
03037 
03038 //------------------------------------------------------------------------------
03039 int SNL_FEI_Structure::finalizeNodeCommMgr()
03040 {
03041   //call initComplete() on the nodeCommMgr, so that it can
03042   //finalize shared-node ownerships, etc.
03043 
03044   //request the safetyCheck if the user requested it via the
03045   //parameters function.
03046   bool safetyCheck = checkSharedNodes_;
03047 
03048   if (safetyCheck && localProc_==0 && numProcs_>1 && outputLevel_>0) {
03049     FEI_COUT << "FEI Info: A consistency-check of shared-node data will be "
03050    << "performed, which involves all-to-all communication. This check is "
03051    << "done only if explicitly requested by parameter "
03052          << "'FEI_CHECK_SHARED_IDS'."<<FEI_ENDL;
03053   }
03054 
03055   CHK_ERR( nodeCommMgr_->initComplete(*nodeDatabase_, safetyCheck) );
03056 
03057   if (debugOutput_) {
03058     FEI_OSTREAM& os = dbgOut();
03059     int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
03060     os << "#     numSharedNodes: " << numSharedNodes << FEI_ENDL;
03061     for(int ns=0; ns<numSharedNodes; ns++) {
03062       NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(ns);
03063       GlobalID nodeID = node.getGlobalNodeID();
03064       int proc = node.getOwnerProc();
03065       os << "#     shNodeID " << (int)nodeID << ", owned by proc "<<proc<<FEI_ENDL;
03066     }
03067   }
03068 
03069   return(0);
03070 }
03071 
03072 //------------------------------------------------------------------------------
03073 int SNL_FEI_Structure::setNumNodesAndEqnsPerBlock()
03074 {
03075   //This function will count how 
03076   //many active nodes there are for each block.
03077 
03078    int numBlocks = blockIDs_.size();
03079    std::vector<int> nodesPerBlock(numBlocks);
03080    std::vector<int> eqnsPerBlock(numBlocks);
03081    GlobalID* blockIDsPtr = &blockIDs_[0];
03082 
03083    int j;
03084    for(j=0; j<numBlocks; j++) {
03085       nodesPerBlock[j] = 0;
03086       eqnsPerBlock[j] = 0;
03087    }
03088 
03089    int numNodes = nodeDatabase_->getNumNodeDescriptors();
03090 
03091    for(j=0; j<numNodes; j++) {
03092      NodeDescriptor* node = NULL;
03093      nodeDatabase_->getNodeAtIndex(j, node);
03094      if (node==NULL) continue;
03095 
03096      const int* fieldIDList = node->getFieldIDList();
03097 
03098      int numFields = node->getNumFields();
03099 
03100      for(int blk=0; blk<numBlocks; blk++) {
03101        GlobalID blockID = blockIDsPtr[blk];
03102 
03103        if (node->containedInBlock(blockID)) {
03104    nodesPerBlock[blk]++;
03105        }
03106        else {
03107    continue;
03108        }
03109 
03110        for(int fld=0; fld<numFields; fld++) {
03111    if (blocks_[blk]->containsField(fieldIDList[fld])) {
03112      int fSize = getFieldSize(fieldIDList[fld]);
03113      assert(fSize >= 0);
03114      eqnsPerBlock[blk] += fSize;
03115    }
03116        }
03117      }
03118    }
03119 
03120    for(j=0; j<numBlocks; j++) {
03121      blocks_[j]->setNumActiveNodes(nodesPerBlock[j]);
03122    }
03123 
03124    //now we need to add the elem-dof to the eqn-count for each block,
03125    //and then set the total number of eqns on each block.
03126 
03127    for(j=0; j<numBlocks; j++) {
03128      eqnsPerBlock[j] += blocks_[j]->getNumElemDOFPerElement() *
03129                         blocks_[j]->getNumElements();
03130 
03131      blocks_[j]->setTotalNumEqns(eqnsPerBlock[j]);
03132    }
03133    return(0);
03134 }
03135 
03136 //------------------------------------------------------------------------------
03137 void SNL_FEI_Structure::initializeEqnCommMgr()
03138 {
03139   FEI_OSTREAM& os = dbgOut();
03140   if (debugOutput_) os << "#  initializeEqnCommMgr" << FEI_ENDL;
03141 
03142   //This function will take information from the (internal member) nodeCommMgr
03143   //and use it to tell the eqnCommMgr which equations we can expect other
03144   //processors to contribute to, and also which equations we'll be sending to
03145   //other processors.
03146   //
03147   //This function can not be called until after initComplete() has been called
03148   //on the nodeCommMgr.
03149   //
03150   int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
03151 
03152   for(int i=0; i<numSharedNodes; i++) {
03153     NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(i);
03154 
03155     int numFields = node.getNumFields();
03156     const int* nodeFieldIDList = node.getFieldIDList();
03157     const int* nodeEqnNums = node.getFieldEqnNumbers();
03158 
03159     int owner = node.getOwnerProc();
03160     if (owner == localProc_) {
03161       std::vector<int>& proc = nodeCommMgr_->getSharedNodeProcs(i);
03162       int numProcs = proc.size();
03163 
03164       for(int p=0; p<numProcs; p++) {
03165 
03166   if (proc[p] == localProc_) continue;
03167 
03168   for(int j=0; j<numFields; j++) {
03169     int numEqns = getFieldSize(nodeFieldIDList[j]);
03170     assert(numEqns >= 0);
03171 
03172     int eqn;
03173     for(eqn=0; eqn<numEqns; eqn++) {
03174       int reducedEqn = -1;
03175       bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn,
03176              reducedEqn);
03177       if (isSlave) continue;
03178 
03179       eqnCommMgr_->addLocalEqn(reducedEqn, proc[p]);
03180     }
03181   }
03182       }
03183     }
03184     else {
03185       for(int j=0; j<numFields; j++) {
03186   int numEqns = getFieldSize(nodeFieldIDList[j]);
03187   assert(numEqns >= 0);
03188   for(int eqn=0; eqn<numEqns; eqn++) {
03189     int reducedEqn = -1;
03190     bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn, reducedEqn);
03191     if (isSlave) continue;
03192 
03193     eqnCommMgr_->addRemoteIndices(reducedEqn, owner, NULL, 0);
03194   }
03195       }
03196     }
03197   }
03198 
03199   if (debugOutput_) os << "#  leaving initializeEqnCommMgr" << FEI_ENDL;
03200 }
03201 
03202 //------------------------------------------------------------------------------
03203 void SNL_FEI_Structure::getEqnInfo(int& numGlobalEqns, int& numLocalEqns,
03204                    int& localStartRow, int& localEndRow){
03205 
03206    numGlobalEqns = numGlobalEqns_;
03207    numLocalEqns = numLocalEqns_;
03208    localStartRow = localStartRow_;
03209    localEndRow = localEndRow_;
03210 }
03211 
03212 //------------------------------------------------------------------------------
03213 int SNL_FEI_Structure::getEqnNumbers(GlobalID ID,
03214                                      int idType, int fieldID,
03215                                      int& numEqns, int* eqnNumbers)
03216 {
03217   //for now, only allow node ID. allowance for element ID is coming soon!!!
03218   if (idType != FEI_NODE) ERReturn(-1);
03219 
03220   NodeDescriptor* node = NULL;
03221   CHK_ERR( nodeDatabase_->getNodeWithID(ID, node) );
03222 
03223   numEqns = getFieldSize(fieldID);
03224   if (numEqns < 0) ERReturn(-1);
03225 
03226   int nodeFieldEqnNumber = -1;
03227   if (!node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber)) {
03228     ERReturn(-1);
03229   }
03230 
03231   for(int i=0; i<numEqns; i++) eqnNumbers[i] = nodeFieldEqnNumber + i;
03232 
03233   return(FEI_SUCCESS);
03234 }
03235 
03236 //------------------------------------------------------------------------------
03237 int SNL_FEI_Structure::getEqnNumbers(int numIDs,
03238              const GlobalID* IDs,
03239              int idType, int fieldID,
03240              int& numEqns, int* eqnNumbers)
03241 {
03242     //for now, only allow node ID. allowance for element ID is coming soon!!!
03243     if (idType != FEI_NODE) ERReturn(-1);
03244 
03245     int fieldSize = getFieldSize(fieldID);
03246     if (fieldSize < 0) {
03247       ERReturn(-1);
03248     }
03249 
03250     int offset = 0;
03251     for(int i=0; i<numIDs; ++i) {
03252       NodeDescriptor* node = NULL;
03253 
03254       if ( nodeDatabase_->getNodeWithID(IDs[i], node) != 0 ) {
03255         // fei::console_out() << "SNL_FEI_Structure::getEqnNumbers: ERROR getting node " << IDs[i] << FEI_ENDL;
03256         for(int ii=0; ii<fieldSize; ii++) {
03257           eqnNumbers[offset++] = -1;
03258         }
03259         continue;
03260         // ERReturn(-1);
03261       }
03262 
03263       int nodeFieldEqnNumber = -1;
03264 
03265       if ( !node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber) ) {
03266         //if we didn't find a field eqn-number, set eqnNumbers to -1. These
03267         //positions will then be skipped by code that passes the data on down to
03268         //underlying solver objects.
03269         for(int ii=0; ii<fieldSize; ii++) {
03270           eqnNumbers[offset++] = -1;
03271         }
03272       }
03273       else {
03274         //else we did find valid eqn-numbers...
03275         for(int ii=0; ii<fieldSize; ii++) {
03276           eqnNumbers[offset++] = nodeFieldEqnNumber + ii;
03277         }
03278       }
03279     }
03280 
03281     numEqns = fieldSize*numIDs;
03282 
03283     return(FEI_SUCCESS);
03284 }
03285 
03286 //------------------------------------------------------------------------------
03287 int SNL_FEI_Structure::translateToReducedNodeNumber(int nodeNumber, int proc)
03288 {
03289   if (proc != localProc_) {
03290     return(nodeNumber - globalNumNodesVanished_[proc]);
03291   }
03292 
03293   int insertPoint = -1;
03294   int index =
03295     fei::binarySearch(nodeNumber, localVanishedNodeNumbers_, insertPoint);
03296 
03297   int localAdjustment = index < 0 ? insertPoint : index + 1;
03298 
03299   return(nodeNumber - localAdjustment - globalNumNodesVanished_[proc]);
03300 }
03301 
03302 //------------------------------------------------------------------------------
03303 void SNL_FEI_Structure::getEqnBlkInfo(int& numGlobalEqnBlks,
03304                                      int& numLocalEqnBlks,
03305                                      int& localBlkOffset) {
03306 
03307    numGlobalEqnBlks = numGlobalEqnBlks_;
03308    numLocalEqnBlks = numLocalEqnBlks_;
03309    localBlkOffset = localBlkOffset_;
03310 }
03311 
03312 //------------------------------------------------------------------------------
03313 int SNL_FEI_Structure::calculateSlaveEqns(MPI_Comm comm)
03314 {
03315   FEI_OSTREAM& os = dbgOut();
03316   if (debugOutput_) os << "#  calculateSlaveEqns" << FEI_ENDL;
03317 
03318   if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
03319   eqnCommMgr_ = new EqnCommMgr(comm_);
03320   eqnCommMgr_->setNumRHSs(1);
03321 
03322   std::vector<int> eqns;
03323   std::vector<int> mEqns;
03324   std::vector<double> mCoefs;
03325 
03326   for(size_t i=0; i<slaveVars_->size(); i++) {
03327     int numEqns;
03328     SlaveVariable* svar = (*slaveVars_)[i];
03329 
03330     eqns.resize( getFieldSize(svar->getFieldID()));
03331     CHK_ERR( getEqnNumbers(svar->getNodeID(), FEI_NODE, svar->getFieldID(),
03332          numEqns, &eqns[0]));
03333 
03334     int slaveEqn = eqns[svar->getFieldOffset()];
03335 
03336     const std::vector<GlobalID>* mNodes = svar->getMasterNodeIDs();
03337     const std::vector<int>* mFields = svar->getMasterFields();
03338     const std::vector<double>* mWeights = svar->getWeights();
03339     const std::vector<double>& mWeightsRef = *mWeights;
03340     int mwOffset = 0;
03341 
03342     for(size_t j=0; j<mNodes->size(); j++) {
03343       int mfSize = getFieldSize((*mFields)[j]);
03344 
03345       eqns.resize(mfSize);
03346       GlobalID mNodeID = (*mNodes)[j];
03347       int mFieldID = (*mFields)[j];
03348       CHK_ERR( getEqnNumbers( mNodeID, FEI_NODE, mFieldID,
03349             mfSize, &eqns[0]));
03350 
03351       double fei_eps = 1.e-49;
03352 
03353       for(int k=0; k<mfSize; k++) {
03354   if (std::abs(mWeightsRef[mwOffset++]) > fei_eps) {
03355     mEqns.push_back(eqns[k]);
03356     mCoefs.push_back(mWeightsRef[mwOffset-1]);
03357   }
03358       }
03359     }
03360 
03361     CHK_ERR( slaveEqns_->addEqn(slaveEqn, &mCoefs[0], &mEqns[0],
03362            mEqns.size(), false) );
03363     mEqns.resize(0);
03364     mCoefs.resize(0);
03365   }
03366 
03367 #ifndef FEI_SER
03368   int numLocalSlaves = slaveVars_->size();
03369   int globalMaxSlaves = 0;
03370   CHK_ERR( fei::GlobalMax(comm_, numLocalSlaves, globalMaxSlaves) );
03371 
03372   if (globalMaxSlaves > 0) {
03373     CHK_ERR( gatherSlaveEqns(comm, eqnCommMgr_, slaveEqns_) );
03374   }
03375 #endif
03376 
03377   globalNumNodesVanished_.resize(numProcs_+1, 0);
03378 
03379   slvEqnNumbers_ = &(slaveEqns_->eqnNumbers());
03380   numSlvs_ = slvEqnNumbers_->size();
03381   if (numSlvs_ > 0) {
03382     //first, let's remove any 'couplings' among the slave equations. A coupling
03383     //is where a slave depends on a master which is also a slave that depends on
03384     //other masters.
03385 
03386     if (debugOutput_) {
03387       os << "#  slave-equations:" << FEI_ENDL;
03388       os << *slaveEqns_;
03389       os << "#  leaving calculateSlaveEqns" << FEI_ENDL;
03390     }
03391 
03392     int levelsOfCoupling;
03393     CHK_ERR( removeCouplings(*slaveEqns_, levelsOfCoupling) );
03394 
03395     if (debugOutput_) {
03396       os << "#     SNL_FEI_Structure: " << levelsOfCoupling
03397    << " level(s) of coupling removed: " << FEI_ENDL;
03398     }
03399 
03400     lowestSlv_ = (*slvEqnNumbers_)[0];
03401     highestSlv_ = (*slvEqnNumbers_)[numSlvs_-1];
03402 
03403     //For each slave equation, we need to find out if we (the local proc) either
03404     //own or share the node from which that equation arises. If we own or share
03405     //a slave node, then we will need access to the solution for each of the
03406     //associated master equations, and all other processors that share the slave
03407     //will also need access to all of the master solutions.
03408     //So,
03409     // 1. for each slave node that we share but don't own, we need to add the
03410     //   master equations to the eqnCommMgr_ object using addRemoteIndices, if
03411     //   they're not local.
03412     // 2. for each slave node that we own, we need to add the master equations
03413     // to the eqnCommMgr_ object using addLocalEqn for each processor that
03414     // shares the slave node.
03415 
03416     std::vector<int>& slvEqns = *slvEqnNumbers_;
03417     std::vector<fei::CSVec*>& mstrEqns = slaveEqns_->eqns();
03418 
03419     //keep track of the number of locally owned nodes that vanish due to the
03420     //fact that all equations at that node are slave equations.
03421     int numLocalNodesVanished = 0;
03422 
03423     GlobalID lastNodeID = -1;
03424 
03425     for(int i=0; i<numSlvs_; i++) {
03426       const NodeDescriptor* node = NULL;
03427       int reducedSlaveEqn;
03428       translateToReducedEqn(slvEqns[i], reducedSlaveEqn);
03429       int numMasters = mstrEqns[i]->size();
03430 
03431       int err = nodeDatabase_->getNodeWithEqn(slvEqns[i], node);
03432       if (err != 0) {
03433         if (debugOutput_) {
03434           os << "#  no local node for slave eqn " << slvEqns[i] << FEI_ENDL;
03435         }
03436 
03437         continue;
03438       }
03439 
03440       if (node->getGlobalNodeID() != lastNodeID &&
03441           node->getOwnerProc() == localProc_) {
03442         if (nodalEqnsAllSlaves(node, slvEqns)) {
03443           numLocalNodesVanished++;
03444           if (fei::sortedListInsert(node->getNodeNumber(), localVanishedNodeNumbers_)
03445               == -2) {
03446             ERReturn(-1);
03447           }
03448         }
03449         lastNodeID = node->getGlobalNodeID();
03450       }
03451 
03452       GlobalID slvNodeID = node->getGlobalNodeID();
03453       int shIndex = nodeCommMgr_->getSharedNodeIndex(slvNodeID);
03454       if (shIndex < 0) {
03455         continue;
03456       }
03457 
03458       std::vector<int>& sharingProcs = nodeCommMgr_->getSharedNodeProcs(shIndex);
03459 
03460       for(int j=0; j<numMasters; j++) {
03461         int masterEquation = mstrEqns[i]->indices()[j];
03462         int owner = getOwnerProcForEqn( masterEquation );
03463 
03464         int reducedMasterEqn;
03465         translateToReducedEqn(masterEquation, reducedMasterEqn);
03466 
03467         if (owner == localProc_) {
03468           int numSharing = sharingProcs.size();
03469           for(int sp=0; sp<numSharing; sp++) {
03470             if (sharingProcs[sp] == localProc_) continue;
03471 
03472             if (debugOutput_) {
03473               os << "#     slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
03474                  << ", master eqn " << masterEquation << " eqnCommMgr_->addLocal "
03475                  << reducedMasterEqn << ", proc " << sharingProcs[sp] << FEI_ENDL;
03476             }
03477             eqnCommMgr_->addLocalEqn(reducedMasterEqn, sharingProcs[sp]);
03478       slvCommMgr_->addRemoteIndices(reducedSlaveEqn, sharingProcs[sp],
03479             &reducedMasterEqn, 1);
03480           }
03481         }
03482         else {
03483           if (debugOutput_) {
03484             os << "#     slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
03485                << ", master eqn " << masterEquation << " eqnCommMgr_->addRemote "
03486                << reducedMasterEqn << ", proc " << owner << FEI_ENDL;
03487           }
03488           eqnCommMgr_->addRemoteIndices(reducedMasterEqn, owner,
03489                                         &reducedMasterEqn, 1);
03490     slvCommMgr_->addLocalEqn(reducedSlaveEqn, owner);
03491         }
03492       }
03493     }
03494 
03495     std::vector<int> tmp(1), tmp2(numProcs_);
03496     tmp[0] = numLocalNodesVanished;
03497     CHK_ERR( fei::Allgatherv(comm_, tmp, tmp2, globalNumNodesVanished_) );
03498 
03499     if ((int)globalNumNodesVanished_.size() != numProcs_) {
03500       ERReturn(-1);
03501     }
03502 
03503     globalNumNodesVanished_.resize(numProcs_+1);
03504     int tmpNumVanished = 0;
03505     for(int proc=0; proc<numProcs_; ++proc) {
03506       int temporary = tmpNumVanished;
03507       tmpNumVanished += globalNumNodesVanished_[proc];
03508       globalNumNodesVanished_[proc] = temporary;
03509     }
03510     globalNumNodesVanished_[numProcs_] = tmpNumVanished;
03511   }
03512 
03513   if (slaveMatrix_ != NULL) delete slaveMatrix_;
03514   slaveMatrix_ = new fei::FillableMat(*slaveEqns_);
03515 
03516   if (debugOutput_) {
03517     os << "#  slave-equations:" << FEI_ENDL;
03518     os << *slaveEqns_;
03519     os << "#  leaving calculateSlaveEqns" << FEI_ENDL;
03520   }
03521 
03522   return(0);
03523 }
03524 
03525 //------------------------------------------------------------------------------
03526 int SNL_FEI_Structure::removeCouplings(EqnBuffer& eqnbuf, int& levelsOfCoupling)
03527 {
03528   std::vector<int>& eqnNumbers = eqnbuf.eqnNumbers();
03529   std::vector<fei::CSVec*>& eqns = eqnbuf.eqns();
03530 
03531   std::vector<double> tempCoefs;
03532   std::vector<int> tempIndices;
03533 
03534   levelsOfCoupling = 0;
03535   bool finished = false;
03536   while(!finished) {
03537     bool foundCoupling = false;
03538     for(size_t i=0; i<eqnNumbers.size(); i++) {
03539       int rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
03540 
03541       if (rowIndex == (int)i) {
03542   fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
03543        << " illegal master-slave constraint coupling. Eqn "
03544        << eqnNumbers[i] << " is both master and slave. " << FEI_ENDL;
03545   ERReturn(-1);
03546       }
03547 
03548       while(rowIndex >= 0) {
03549   foundCoupling = true;
03550   double coef = 0.0;
03551   CHK_ERR( eqnbuf.getCoefAndRemoveIndex( eqnNumbers[rowIndex],
03552                  eqnNumbers[i], coef) );
03553 
03554   std::vector<int>& indicesRef = eqns[i]->indices();
03555   std::vector<double>& coefsRef = eqns[i]->coefs();
03556 
03557   int len = indicesRef.size();
03558   tempCoefs.resize(len);
03559   tempIndices.resize(len);
03560 
03561   double* tempCoefsPtr = &tempCoefs[0];
03562   int* tempIndicesPtr = &tempIndices[0];
03563   double* coefsPtr = &coefsRef[0];
03564   int* indicesPtr = &indicesRef[0];
03565 
03566   for(int j=0; j<len; ++j) {
03567     tempIndicesPtr[j] = indicesPtr[j];
03568     tempCoefsPtr[j] = coef*coefsPtr[j];
03569   }
03570 
03571   CHK_ERR( eqnbuf.addEqn(eqnNumbers[rowIndex], tempCoefsPtr,
03572              tempIndicesPtr, len, true) );
03573 
03574   rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
03575       }
03576     }
03577     if (foundCoupling) ++levelsOfCoupling;
03578     else finished = true;
03579 
03580     if (levelsOfCoupling>1 && !finished) {
03581       fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
03582      << " too many (>1) levels of master-slave constraint coupling. "
03583      << "Hint: coupling is considered infinite if two slaves depend on "
03584      << "each other. This may or may not be the case here." << FEI_ENDL;
03585       ERReturn(-1);
03586     }
03587   }
03588 
03589   return(0);
03590 }
03591 
03592 #ifndef FEI_SER
03593 //------------------------------------------------------------------------------
03594 int SNL_FEI_Structure::gatherSlaveEqns(MPI_Comm comm,
03595                EqnCommMgr* eqnCommMgr,
03596                EqnBuffer* slaveEqns)
03597 {
03598   int numProcs = 0;
03599   if (MPI_Comm_size(comm, &numProcs) != MPI_SUCCESS) return(-1);
03600   if (numProcs == 1) return(0);
03601   int localProc;
03602   if (MPI_Comm_rank(comm, &localProc) != MPI_SUCCESS) return(-1);
03603 
03604   //We're going to send all of our slave equations to all other processors, and
03605   //receive the slave equations from all other processors.
03606   //So we'll first fill a ProcEqns object with all of our eqn/proc pairs,
03607   //then use the regular exchange functions from EqnCommMgr. (This may not be
03608   //the most efficient way to do it, but it involves the least amount of new
03609   //code.)
03610   ProcEqns localProcEqns, remoteProcEqns;
03611   std::vector<int>& slvEqnNums = slaveEqns->eqnNumbers();
03612   fei::CSVec** slvEqnsPtr = &(slaveEqns->eqns()[0]);
03613 
03614   for(size_t i=0; i<slvEqnNums.size(); i++) {
03615     for(int p=0; p<numProcs; p++) {
03616       if (p == localProc) continue;
03617 
03618       localProcEqns.addEqn(slvEqnNums[i], slvEqnsPtr[i]->size(), p);
03619     }
03620   }
03621 
03622   CHK_ERR( eqnCommMgr->mirrorProcEqns(localProcEqns, remoteProcEqns) );
03623 
03624   CHK_ERR( eqnCommMgr->mirrorProcEqnLengths(localProcEqns,
03625               remoteProcEqns));
03626 
03627   EqnBuffer remoteEqns;
03628   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm, &localProcEqns, slaveEqns,
03629             &remoteProcEqns, &remoteEqns,
03630             false) );
03631 
03632   //Now the remoteEqns buffer should hold the slave equations from all other
03633   //processors, so let's add them to the ones we already have.
03634   CHK_ERR( slaveEqns->addEqns(remoteEqns, false) );
03635 
03636   return(0);
03637 }
03638 #endif
03639 
03640 //------------------------------------------------------------------------------
03641 bool SNL_FEI_Structure::isSlaveEqn(int eqn)
03642 {
03643   if (numSlvs_ == 0) return(false);
03644 
03645   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03646   int insertPoint = -1;
03647   int foundOffset = fei::binarySearch(eqn, slvEqns, insertPoint);
03648 
03649   if (foundOffset >= 0) return(true);
03650   else return(false);
03651 }
03652 
03653 //------------------------------------------------------------------------------
03654 bool SNL_FEI_Structure::translateToReducedEqn(int eqn, int& reducedEqn)
03655 {
03656   if (numSlvs_ == 0) { reducedEqn = eqn; return(false); }
03657 
03658   if (eqn < lowestSlv_) {reducedEqn = eqn; return(false); }
03659   if (eqn > highestSlv_) {reducedEqn = eqn - numSlvs_; return(false); }
03660 
03661   int index = 0;
03662   int foundOffset = fei::binarySearch(eqn, *slvEqnNumbers_, index);
03663 
03664   bool isSlave = false;
03665 
03666   if (foundOffset < 0) {
03667     reducedEqn = eqn - index;
03668   }
03669   else {
03670     isSlave = true; reducedEqn = eqn - (foundOffset+1);
03671   }
03672 
03673   return(isSlave);
03674 }
03675 
03676 //------------------------------------------------------------------------------
03677 int SNL_FEI_Structure::translateFromReducedEqn(int reducedEqn)
03678 {
03679   int numSlvs = slaveEqns_->getNumEqns();
03680 
03681   if (numSlvs == 0) return(reducedEqn);
03682 
03683   const int* slvEqns = &(slaveEqns_->eqnNumbers()[0]);
03684 
03685   if (reducedEqn < slvEqns[0]) return(reducedEqn);
03686 
03687   int eqn = reducedEqn;
03688 
03689   for(int i=0; i<numSlvs; i++) {
03690     if (eqn < slvEqns[i]) return(eqn);
03691     eqn++;
03692   }
03693 
03694   return(eqn);
03695 }
03696 
03697 //------------------------------------------------------------------------------
03698 int SNL_FEI_Structure::getMasterEqnNumbers(int slaveEqn,
03699             std::vector<int>*& masterEqns)
03700 {
03701   if (slaveEqns_->getNumEqns() == 0) {
03702     masterEqns = NULL;
03703     return(0);
03704   }
03705 
03706   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03707   int index = 0;
03708   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03709 
03710   if (foundOffset >= 0) {
03711     masterEqns = &(slaveEqns_->eqns()[foundOffset]->indices());
03712   }
03713   else {
03714     masterEqns = NULL;
03715   }
03716 
03717   return(0);
03718 }
03719 
03720 //------------------------------------------------------------------------------
03721 int SNL_FEI_Structure::getMasterEqnCoefs(int slaveEqn,
03722           std::vector<double>*& masterCoefs)
03723 {
03724   if (slaveEqns_->getNumEqns() == 0) {
03725     masterCoefs = NULL;
03726     return(0);
03727   }
03728 
03729   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03730   int index = 0;
03731   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03732 
03733   if (foundOffset >= 0) {
03734     masterCoefs = &(slaveEqns_->eqns()[foundOffset]->coefs());
03735   }
03736   else {
03737     masterCoefs = NULL;
03738   }
03739 
03740   return(0);
03741 }
03742 
03743 //------------------------------------------------------------------------------
03744 int SNL_FEI_Structure::getMasterEqnRHS(int slaveEqn,
03745               double& rhsValue)
03746 {
03747   if (slaveEqns_->getNumEqns() == 0) {
03748     return(0);
03749   }
03750 
03751   std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
03752   int index = 0;
03753   int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
03754 
03755   if (foundOffset >= 0) {
03756     std::vector<double>* rhsCoefsPtr = (*(slaveEqns_->rhsCoefsPtr()))[foundOffset];
03757     rhsValue = (*rhsCoefsPtr)[0];
03758   }
03759 
03760   return(0);
03761 }
03762 
03763 //------------------------------------------------------------------------------
03764 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID, 
03765                                             int interleaveStrategy,
03766                                             int* scatterIndices)
03767 {
03768    int index = fei::binarySearch(blockID, blockIDs_);
03769 
03770    if (index < 0) {
03771       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
03772            << (int)blockID << " not found. Aborting." << FEI_ENDL;
03773       std::abort();
03774    }
03775 
03776    std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
03777 
03778    std::map<GlobalID,int>::const_iterator
03779      iter = elemIDs.find(elemID);
03780 
03781    if (iter == elemIDs.end()) {
03782       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: " 
03783            << (int)blockID << ", elemID "
03784            << (int)elemID << " not found. Aborting." << FEI_ENDL;
03785       std::abort();
03786    }
03787 
03788    int elemIndex = iter->second;
03789 
03790    getScatterIndices_index(index, elemIndex,
03791                            interleaveStrategy, scatterIndices);
03792 }
03793 
03794 //------------------------------------------------------------------------------
03795 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID, 
03796                int interleaveStrategy,
03797                int* scatterIndices,
03798                int* blkScatterIndices,
03799                int* blkSizes)
03800 {
03801    int index = fei::binarySearch(blockID, blockIDs_);
03802 
03803    if (index < 0) {
03804       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
03805            << (int)blockID << " not found. Aborting." << FEI_ENDL;
03806       std::abort();
03807    }
03808 
03809    std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
03810 
03811    std::map<GlobalID,int>::const_iterator
03812      iter = elemIDs.find(elemID);
03813 
03814    if (iter == elemIDs.end()) {
03815       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: " 
03816            << (int)blockID << ", elemID "
03817            << (int)elemID << " not found. Aborting." << FEI_ENDL;
03818       std::abort();
03819    }
03820 
03821    int elemIndex = iter->second;
03822 
03823    getScatterIndices_index(index, elemIndex,
03824                            interleaveStrategy, scatterIndices,
03825          blkScatterIndices, blkSizes);
03826 }
03827 
03828 //------------------------------------------------------------------------------
03829 int SNL_FEI_Structure::getBlkScatterIndices_index(int blockIndex,
03830               int elemIndex,
03831               int* scatterIndices)
03832 {
03833   BlockDescriptor& block = *(blocks_[blockIndex]);
03834   int numNodes = block.getNumNodesPerElement();
03835   work_nodePtrs_.resize(numNodes);
03836   NodeDescriptor** nodes = &work_nodePtrs_[0];
03837   int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03838   if (err) {
03839     fei::console_out() << "SNL_FEI_Structure::getBlkScatterIndices_index: ERROR getting"
03840    << " node descriptors." << FEI_ENDL;
03841     ERReturn(-1);
03842   }
03843 
03844   int offset = 0;
03845   return( getNodeBlkIndices(nodes, numNodes, scatterIndices, offset) );
03846 }
03847 
03848 //------------------------------------------------------------------------------
03849 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
03850             int interleaveStrategy,
03851             int* scatterIndices)
03852 {
03853 //On input, scatterIndices, is assumed to be allocated by the calling code,
03854 // and be of length the number of equations per element.
03855 //
03856    BlockDescriptor& block = *(blocks_[blockIndex]);
03857    int numNodes = block.getNumNodesPerElement();
03858    int* fieldsPerNode = block.fieldsPerNodePtr();
03859    int** fieldIDs = block.fieldIDsTablePtr();
03860 
03861    work_nodePtrs_.resize(numNodes);
03862    NodeDescriptor** nodes = &work_nodePtrs_[0];
03863 
03864    int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03865    if (err) {
03866       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
03867            << " node descriptors." << FEI_ENDL;
03868       std::abort();
03869    }
03870 
03871    int offset = 0;
03872    if (fieldDatabase_->size() == 1) {
03873      err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
03874             scatterIndices, offset);
03875      if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
03876    }
03877    else {
03878      switch (interleaveStrategy) {
03879      case 0:
03880        err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03881          scatterIndices, offset);
03882        if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
03883        break;
03884 
03885      case 1:
03886        err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03887           scatterIndices, offset);
03888        if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
03889        break;
03890 
03891      default:
03892        fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
03893        break;
03894      }
03895    }
03896 
03897    //now the element-DOF.
03898    int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
03899    std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
03900 
03901    for(int i=0; i<numElemDOF; i++) {
03902       scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
03903    }
03904 }
03905 
03906 //------------------------------------------------------------------------------
03907 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
03908             int interleaveStrategy,
03909             int* scatterIndices,
03910             int* blkScatterIndices,
03911             int* blkSizes)
03912 {
03913 //On input, scatterIndices, is assumed to be allocated by the calling code,
03914 // and be of length the number of equations per element.
03915 //
03916    BlockDescriptor& block = *(blocks_[blockIndex]);
03917    int numNodes = block.getNumNodesPerElement();
03918    int* fieldsPerNode = block.fieldsPerNodePtr();
03919    int** fieldIDs = block.fieldIDsTablePtr();
03920 
03921    work_nodePtrs_.resize(numNodes);
03922    NodeDescriptor** nodes = &work_nodePtrs_[0];
03923 
03924    int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
03925    if (err) {
03926       fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
03927            << " node descriptors." << FEI_ENDL;
03928       std::abort();
03929    }
03930 
03931    int offset = 0, blkOffset = 0;
03932    if (fieldDatabase_->size() == 1) {
03933      err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
03934          scatterIndices, offset,
03935          blkScatterIndices, blkSizes, blkOffset);
03936      if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
03937    }
03938    else {
03939      switch (interleaveStrategy) {
03940      case 0:
03941        err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03942          scatterIndices, offset,
03943          blkScatterIndices, blkSizes, blkOffset);
03944        if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
03945        break;
03946 
03947      case 1:
03948        err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
03949           scatterIndices, offset);
03950        if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
03951        break;
03952 
03953      default:
03954        fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
03955        break;
03956      }
03957    }
03958 
03959    //now the element-DOF.
03960    int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
03961    std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
03962 
03963    for(int i=0; i<numElemDOF; i++) {
03964       scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
03965       if (interleaveStrategy == 0) {
03966   blkSizes[blkOffset] = 1;
03967   blkScatterIndices[blkOffset++] = elemDOFEqns[elemIndex] + i;
03968       }
03969    }
03970 }
03971 
03972 //------------------------------------------------------------------------------
03973 int SNL_FEI_Structure::getElemNodeDescriptors(int blockIndex, int elemIndex,
03974                                              NodeDescriptor** nodes)
03975 {
03976   //Private function, called by 'getScatterIndices_index'. We can safely
03977   //assume that blockIndex and elemIndex are valid in-range indices.
03978   //
03979   //This function's task is to obtain the NodeDescriptor objects, from the
03980   //nodeDatabase, that are connected to the specified element.
03981   //
03982   int numNodes = connTables_[blockIndex]->numNodesPerElem;
03983 
03984   if (activeNodesInitialized_) {
03985     NodeDescriptor** elemNodePtrs =
03986       &((*connTables_[blockIndex]->elem_conn_ptrs)[elemIndex*numNodes]);
03987     for(int i=0; i<numNodes; ++i) nodes[i] = elemNodePtrs[i];
03988   }
03989   else {
03990     const GlobalID* elemConn = 
03991       &((*connTables_[blockIndex]->elem_conn_ids)[elemIndex*numNodes]);
03992     for(int i=0; i<numNodes; ++i) {
03993       CHK_ERR( nodeDatabase_->getNodeWithID(elemConn[i], nodes[i]));
03994     }
03995   }
03996 
03997   return(FEI_SUCCESS);
03998 }
03999 
04000 //------------------------------------------------------------------------------
04001 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
04002                int numNodes,
04003                int fieldID,
04004                int* scatterIndices,
04005                int& offset)
04006 {
04007   int fieldSize = getFieldSize(fieldID);
04008 
04009   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04010     NodeDescriptor& node = *(nodes[nodeIndex]);
04011     const int* eqnNumbers = node.getFieldEqnNumbers();
04012     int eqn = eqnNumbers[0];
04013     scatterIndices[offset++] = eqn;
04014     if (fieldSize > 1) {
04015       for(int i=1; i<fieldSize; i++) {
04016   scatterIndices[offset++] = eqn+i;
04017       }
04018     }
04019   }
04020   return(0);
04021 }
04022 
04023 //------------------------------------------------------------------------------
04024 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
04025                int numNodes,
04026                int fieldID,
04027                int* scatterIndices,
04028                int& offset,
04029                int* blkScatterIndices,
04030                int* blkSizes,
04031                int& blkOffset)
04032 {
04033   int fieldSize = getFieldSize(fieldID);
04034 
04035   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04036     NodeDescriptor& node = *(nodes[nodeIndex]);
04037     const int* eqnNumbers = node.getFieldEqnNumbers();
04038     int eqn = eqnNumbers[0];
04039     scatterIndices[offset++] = eqn;
04040     if (fieldSize > 1) {
04041       for(int i=1; i<fieldSize; i++) {
04042   scatterIndices[offset++] = eqn+i;
04043       }
04044     }
04045     blkSizes[blkOffset] = node.getNumNodalDOF();
04046     blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
04047   }
04048   return(0);
04049 }
04050 
04051 //------------------------------------------------------------------------------
04052 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04053                                           int** fieldIDs, int* fieldsPerNode,
04054                                           int* scatterIndices, int& offset)
04055 {
04056   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04057 
04058       NodeDescriptor& node = *(nodes[nodeIndex]);
04059 
04060       const int* nodeFieldIDList = node.getFieldIDList();
04061       const int* nodeEqnNums = node.getFieldEqnNumbers();
04062       int numFields = node.getNumFields();
04063 
04064       int* fieldID_ind = fieldIDs[nodeIndex];
04065 
04066       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04067          int numEqns = getFieldSize(fieldID_ind[j]);
04068          assert(numEqns >= 0);
04069 
04070          int insert = -1;
04071          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04072                                              numFields, insert);
04073 
04074          if (nind >= 0) {
04075      int eqn = nodeEqnNums[nind];
04076 
04077      if (eqn < 0) {
04078        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04079       << (int)node.getGlobalNodeID()
04080       << ", field " << nodeFieldIDList[nind]
04081       << " has equation number " << eqn << FEI_ENDL;
04082        std::abort();
04083      }
04084 
04085      for(int jj=0; jj<numEqns; jj++) {
04086        scatterIndices[offset++] = eqn+jj;
04087      }
04088          }
04089          else {
04090      if (outputLevel_ > 2) {
04091        fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
04092       << " not found for node "
04093       << (int)(node.getGlobalNodeID()) << FEI_ENDL;
04094      }
04095          }
04096       }
04097    }
04098 
04099    return(FEI_SUCCESS);
04100 }
04101 
04102 //------------------------------------------------------------------------------
04103 int SNL_FEI_Structure::getNodeBlkIndices(NodeDescriptor** nodes,
04104            int numNodes,
04105            int* scatterIndices,
04106            int& offset)
04107 {
04108   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04109     NodeDescriptor* node = nodes[nodeIndex];
04110     scatterIndices[offset++] = node->getBlkEqnNumber();
04111   }
04112   return(0);
04113 }
04114 
04115 //------------------------------------------------------------------------------
04116 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04117              int** fieldIDs, int* fieldsPerNode,
04118              int* scatterIndices, int& offset,
04119              int* blkScatterIndices,
04120              int* blkSizes,
04121              int& blkOffset)
04122 {
04123   for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04124 
04125       NodeDescriptor& node = *(nodes[nodeIndex]);
04126 
04127       const int* nodeFieldIDList = node.getFieldIDList();
04128       const int* nodeEqnNums = node.getFieldEqnNumbers();
04129       int numFields = node.getNumFields();
04130 
04131       blkSizes[blkOffset] = node.getNumNodalDOF();
04132       blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
04133 
04134       int* fieldID_ind = fieldIDs[nodeIndex];
04135 
04136       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04137          int numEqns = getFieldSize(fieldID_ind[j]);
04138          assert(numEqns >= 0);
04139 
04140          int insert = -1;
04141          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04142                                              numFields, insert);
04143 
04144          if (nind >= 0) {
04145      int eqn = nodeEqnNums[nind];
04146      if (eqn < 0) {
04147        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04148       << (int)node.getGlobalNodeID()
04149       << ", field " << nodeFieldIDList[nind]
04150       << " has equation number " << eqn << FEI_ENDL;
04151        std::abort();
04152      }
04153 
04154      for(int jj=0; jj<numEqns; jj++) {
04155        scatterIndices[offset++] = eqn+jj;
04156      }
04157          }
04158          else {
04159      if (outputLevel_ > 2) {
04160        fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
04161       << " not found for node "
04162       << (int)(node.getGlobalNodeID()) << FEI_ENDL;
04163      }
04164          }
04165       }
04166    }
04167 
04168    return(FEI_SUCCESS);
04169 }
04170 
04171 //------------------------------------------------------------------------------
04172 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
04173                                           std::vector<int>* fieldIDs,
04174                                           std::vector<int>& fieldsPerNode,
04175                                           std::vector<int>& scatterIndices)
04176 {
04177    int offset = 0;
04178    scatterIndices.resize(0);
04179 
04180    for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04181 
04182       NodeDescriptor& node = *(nodes[nodeIndex]);
04183 
04184       const int* nodeFieldIDList = node.getFieldIDList();
04185       const int* nodeEqnNums = node.getFieldEqnNumbers();
04186       int numFields = node.getNumFields();
04187 
04188       int* fieldID_ind = &(fieldIDs[nodeIndex][0]);
04189 
04190       for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
04191          int numEqns = getFieldSize(fieldID_ind[j]);
04192          assert(numEqns >= 0);
04193 
04194          int insert = -1;
04195          int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
04196                                              numFields, insert);
04197 
04198          if (nind >= 0) {
04199      int eqn = nodeEqnNums[nind];
04200 
04201      if (eqn < 0) {
04202        FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
04203       << (int)node.getGlobalNodeID()
04204       << ", field " << nodeFieldIDList[nind]
04205       << " has equation number " << eqn << FEI_ENDL;
04206        MPI_Abort(comm_, -1);
04207      }
04208 
04209      scatterIndices.resize(offset+numEqns);
04210      int* inds = &scatterIndices[0];
04211 
04212      for(int jj=0; jj<numEqns; jj++) {
04213        inds[offset+jj] = eqn+jj;
04214      }
04215      offset += numEqns;
04216          }
04217          else {
04218      if (outputLevel_ > 2) {
04219        fei::console_out() << "WARNING, field ID " << fieldID_ind[j]
04220       << " not found for node "
04221       << (int)node.getGlobalNodeID() << FEI_ENDL;
04222            }
04223          }
04224       }
04225    }
04226 
04227    return(FEI_SUCCESS);
04228 }
04229 
04230 //------------------------------------------------------------------------------
04231 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
04232                                           int** fieldIDs, int* fieldsPerNode,
04233                                           int* scatterIndices, int& offset)
04234 {
04235   //In this function we want to group equations by field, but
04236   //in what order should the fields be?
04237   //Let's just run through the fieldIDs table, and add the fields to a
04238   //flat list, in the order we encounter them, but making sure no fieldID
04239   //gets added more than once.
04240 
04241   std::vector<int> fields;
04242   for(int ii=0; ii<numNodes; ii++) {
04243     for(int j=0; j<fieldsPerNode[ii]; j++) {
04244       if (std::find(fields.begin(), fields.end(), fieldIDs[ii][j]) == fields.end()) {
04245         fields.push_back(fieldIDs[ii][j]);
04246       }
04247     }
04248   }
04249 
04250   int* fieldsPtr = fields.size()>0 ? &fields[0] : NULL;
04251 
04252   //ok, we've got our flat list of fields, so let's proceed to get the
04253   //scatter indices.
04254 
04255   for(size_t i=0; i<fields.size(); i++) {
04256     int field = fieldsPtr[i];
04257 
04258     for(int nodeIndex = 0; nodeIndex < numNodes; ++nodeIndex) {
04259       int fidx = fei::searchList(field, fieldIDs[nodeIndex],
04260             fieldsPerNode[nodeIndex]);
04261       if (fidx < 0) {
04262   continue;
04263       }
04264 
04265       NodeDescriptor* node = nodes[nodeIndex];
04266 
04267       const int* nodeFieldIDList = node->getFieldIDList();
04268       const int* nodeEqnNums = node->getFieldEqnNumbers();
04269       int numFields = node->getNumFields();
04270 
04271       int numEqns = getFieldSize(field);
04272       assert(numEqns >= 0);
04273 
04274       int insert = -1;
04275       int nind = fei::binarySearch(field, nodeFieldIDList,
04276             numFields, insert);
04277 
04278       if (nind > -1) {
04279   for(int jj=0; jj<numEqns; ++jj) {
04280     scatterIndices[offset++] = nodeEqnNums[nind]+jj;
04281   }
04282       }
04283       else {
04284   ERReturn(-1);
04285       }
04286     }
04287   }
04288 
04289   return(FEI_SUCCESS);
04290 }
04291 
04292 //------------------------------------------------------------------------------
04293 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
04294                                           std::vector<int>* fieldIDs,
04295                                           std::vector<int>& fieldsPerNode,
04296                                           std::vector<int>& scatterIndices)
04297 {
04298    //In this function we want to group equations by field, but
04299    //in what order should the fields be?
04300    //Let's just run through the fieldIDs table, and add the fields to a
04301    //flat list, in the order we encounter them, but making sure no fieldID
04302    //gets added more than once.
04303 
04304    std::vector<int> fields;
04305    for(int ii=0; ii<numNodes; ii++) {
04306       for(int j=0; j<fieldsPerNode[ii]; j++) {
04307          std::vector<int>& fldIDs = fieldIDs[ii];
04308          if (std::find(fields.begin(), fields.end(), fldIDs[j]) == fields.end()) {
04309            fields.push_back(fldIDs[j]);
04310          }
04311       }
04312    }
04313 
04314    //ok, we've got our flat list of fields, so let's proceed to get the
04315    //scatter indices.
04316 
04317    int offset = 0;
04318    scatterIndices.resize(0);
04319 
04320    for(size_t i=0; i<fields.size(); i++) {
04321       for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
04322 
04323          const int* nodeFieldIDList = nodes[nodeIndex]->getFieldIDList();
04324          const int* nodeEqnNums = nodes[nodeIndex]->getFieldEqnNumbers();
04325          int numFields = nodes[nodeIndex]->getNumFields();
04326 
04327          int numEqns = getFieldSize(fields[i]);
04328          assert(numEqns >= 0);
04329 
04330          int insert = -1;
04331          int nind = fei::binarySearch(fields[i], nodeFieldIDList,
04332                                              numFields, insert);
04333 
04334          if (nind >= 0) {
04335             for(int jj=0; jj<numEqns; jj++) {
04336                scatterIndices.push_back(nodeEqnNums[nind]+jj);
04337          offset++;
04338             }
04339          }
04340          else {
04341      if (outputLevel_ > 2) {
04342        fei::console_out() << "WARNING, field ID " << fields[i]
04343       << " not found for node "
04344       << (int)nodes[nodeIndex]->getGlobalNodeID() << FEI_ENDL;
04345      }
04346          }
04347       }
04348    }
04349 
04350    return(FEI_SUCCESS);
04351 }
04352 
04353 //------------------------------------------------------------------------------
04354 void SNL_FEI_Structure::addCR(int CRID,
04355            snl_fei::Constraint<GlobalID>*& cr,
04356          std::map<GlobalID,snl_fei::Constraint<GlobalID>* >& crDB)
04357 {
04358   std::map<GlobalID,snl_fei::Constraint<GlobalID>* >::iterator
04359     cr_iter = crDB.find(CRID);
04360 
04361   if (cr_iter == crDB.end()) {
04362     cr = new ConstraintType;
04363     crDB.insert(std::pair<GlobalID,snl_fei::Constraint<GlobalID>*>(CRID, cr));
04364   }
04365 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends