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

Generated on Wed May 12 21:30:42 2010 for FEI by  doxygen 1.4.7