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