FEI Version of the Day
fei_VectorSpace.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include "fei_sstream.hpp"
00010 #include "fei_fstream.hpp"
00011 
00012 #include "fei_utils.hpp"
00013 
00014 #include "fei_TemplateUtils.hpp"
00015 #include "fei_chk_mpi.hpp"
00016 #include <fei_CommUtils.hpp>
00017 #include <fei_set_shared_ids.hpp>
00018 #include "snl_fei_Utils.hpp"
00019 #include "fei_Record.hpp"
00020 #include "snl_fei_RecordCollection.hpp"
00021 #include "fei_ParameterSet.hpp"
00022 #include "snl_fei_RecordMsgHandler.hpp"
00023 #include "fei_SharedIDs.hpp"
00024 #include "fei_Pattern.hpp"
00025 #include "fei_VectorSpace.hpp"
00026 #include "fei_FieldMask.hpp"
00027 #include "snl_fei_PointBlockMap.hpp"
00028 #include "fei_LogManager.hpp"
00029 
00030 #undef fei_file
00031 #define fei_file "fei_VectorSpace.cpp"
00032 #include "fei_ErrMacros.hpp"
00033 
00034 namespace fei {
00035   class RecordAttributeCounter : public Record_Operator<int> {
00036   public:
00037     RecordAttributeCounter(int proc)
00038       : numLocalDOF_(0),
00039         numLocalIDs_(0),
00040         numLocallyOwnedIDs_(0),
00041         numRemoteSharedDOF_(0),
00042         proc_(proc)
00043     {}
00044 
00045     virtual ~RecordAttributeCounter(){}
00046 
00047     void operator()(fei::Record<int>& record)
00048     {
00049       fei::FieldMask* mask = record.getFieldMask();
00050       int owner = record.getOwnerProc();
00051 
00052       if (owner != proc_) {
00053         numRemoteSharedDOF_ += mask->getNumIndices();
00054         return;
00055       }
00056       else {
00057         ++numLocallyOwnedIDs_;
00058       }
00059 
00060       ++numLocalIDs_;
00061 
00062       int numDOF = mask->getNumIndices();
00063 
00064       numLocalDOF_ += numDOF;  
00065     }
00066 
00067     int numLocalDOF_;
00068     int numLocalIDs_;
00069     int numLocallyOwnedIDs_;
00070     int numRemoteSharedDOF_;
00071 
00072   private:
00073     int proc_;
00074   };//class Record_Operator
00075 
00076   class BlkIndexAccessor : public Record_Operator<int> {
00077   public:
00078     BlkIndexAccessor(int localProc,
00079                      int lenBlkIndices,
00080                      int* globalBlkIndices,
00081                      int* blkSizes)
00082       : numBlkIndices_(0),
00083         proc_(localProc),
00084         lenBlkIndices_(lenBlkIndices),
00085         globalBlkIndices_(globalBlkIndices),
00086         blkSizes_(blkSizes)
00087     {
00088     }
00089 
00090     BlkIndexAccessor(int lenBlkIndices,
00091                      int* globalBlkIndices,
00092                      int* blkSizes)
00093       : numBlkIndices_(0),
00094         proc_(-1),
00095         lenBlkIndices_(lenBlkIndices),
00096         globalBlkIndices_(globalBlkIndices),
00097         blkSizes_(blkSizes)
00098     {
00099     }
00100 
00101     void operator()(fei::Record<int>& record)
00102     {
00103       int owner = record.getOwnerProc();
00104       if (owner != proc_ && proc_ > -1) {
00105         return;
00106       }
00107 
00108       fei::FieldMask* mask = record.getFieldMask();
00109       int blkSize = mask->getNumIndices();
00110 
00111       if (numBlkIndices_ < lenBlkIndices_) {
00112         globalBlkIndices_[numBlkIndices_] = record.getNumber();
00113         blkSizes_[numBlkIndices_] = blkSize;
00114       }
00115 
00116       ++numBlkIndices_;
00117     }
00118 
00119     int numBlkIndices_;
00120 
00121   private:
00122     int proc_;
00123     int lenBlkIndices_;
00124     int* globalBlkIndices_;
00125     int* blkSizes_;
00126   };//class BlkIndexAccessor
00127 }//namespace snl_fei
00128 
00129 //----------------------------------------------------------------------------
00130 fei::SharedPtr<fei::VectorSpace>
00131 fei::VectorSpace::Factory::createVectorSpace(MPI_Comm comm,
00132                                              const char* name)
00133 {
00134   fei::SharedPtr<fei::VectorSpace> sptr(new fei::VectorSpace(comm, name));
00135   return(sptr);
00136 }
00137 
00138 //----------------------------------------------------------------------------
00139 fei::VectorSpace::VectorSpace(MPI_Comm comm, const char* name)
00140   : fieldMasks_(),
00141     comm_(comm),
00142     idTypes_(),
00143     fieldDatabase_(),
00144     fieldDofMap_(),
00145     recordCollections_(),
00146     sharedIDTables_(),
00147     ownerPatterns_(),
00148     sharerPatterns_(),
00149     sharedRecordsSynchronized_(true),
00150     ptBlkMap_(NULL),
00151     globalOffsets_(),
00152     globalIDOffsets_(),
00153     simpleProblem_(false),
00154     firstLocalOffset_(-1),
00155     lastLocalOffset_(-1),
00156     eqnNumbers_(),
00157     newInitData_(false),
00158     name_(),
00159     dbgprefix_("VecSpc: "),
00160     checkSharedIDs_(false)
00161 {
00162   check_version();
00163 
00164 //The initializations below are redundant with the ones above in the
00165 //initializer list. For some reason, when compiled for janus with optimization,
00166 //these pointers are not being set to NULL by the above initializers.
00167 //This was causing seg faults later when other VectorSpace methods delete these
00168 //pointers if they aren't NULL.  ABW 5/19/2004
00169 
00170   ptBlkMap_ = NULL;
00171 
00172   int numProcs = fei::numProcs(comm_);
00173   globalOffsets_.assign(numProcs+1, 0);
00174   globalIDOffsets_.assign(numProcs+1, 0);
00175 
00176   setName(name);
00177 }
00178 
00179 //----------------------------------------------------------------------------
00180 fei::VectorSpace::~VectorSpace()
00181 {
00182   int i, len = fieldMasks_.size();
00183   for(i=0; i<len; ++i) delete fieldMasks_[i];
00184 
00185   len = recordCollections_.size();
00186   for(i=0; i<len; ++i) delete recordCollections_[i];
00187 
00188   std::map<int,fei::comm_map*>::iterator
00189     o_iter = ownerPatterns_.begin(), o_end = ownerPatterns_.end();
00190   for(; o_iter != o_end; ++o_iter) {
00191     delete o_iter->second;
00192   }
00193 
00194   std::map<int,fei::comm_map*>::iterator
00195     s_iter = sharerPatterns_.begin(), s_end = sharerPatterns_.end();
00196   for(; s_iter != s_end; ++s_iter) {
00197     delete s_iter->second;
00198   }
00199 
00200   delete ptBlkMap_;
00201 }
00202 
00203 //----------------------------------------------------------------------------
00204 MPI_Comm
00205 fei::VectorSpace::getCommunicator() const
00206 {
00207   return comm_;
00208 }
00209 
00210 //----------------------------------------------------------------------------
00211 void fei::VectorSpace::setParameters(const fei::ParameterSet& paramset)
00212 {
00213   const fei::Param* param = paramset.get("name");
00214   fei::Param::ParamType ptype = param != NULL ?
00215     param->getType() : fei::Param::BAD_TYPE;
00216   if (ptype == fei::Param::STRING) {
00217     setName(param->getStringValue().c_str());
00218   }
00219 
00220   param = paramset.get("FEI_OUTPUT_LEVEL");
00221   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00222   if (ptype == fei::Param::STRING) {
00223     fei::LogManager& log_manager = fei::LogManager::getLogManager();
00224     log_manager.setOutputLevel(param->getStringValue().c_str());
00225     setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
00226   }
00227 
00228   param = paramset.get("FEI_LOG_EQN");
00229   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00230   if (ptype == fei::Param::INT) {
00231     addLogEqn(param->getIntValue());
00232   }
00233 
00234   param = paramset.get("FEI_LOG_ID");
00235   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00236   if (ptype == fei::Param::INT) {
00237     addLogID(param->getIntValue());
00238   }
00239 
00240   param = paramset.get("FEI_CHECK_SHARED_IDS");
00241   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00242   if (ptype != fei::Param::BAD_TYPE) {
00243     if (ptype == fei::Param::BOOL) {
00244       checkSharedIDs_ = param->getBoolValue();
00245     }
00246     else if (ptype == fei::Param::INT) {
00247       checkSharedIDs_ = param->getIntValue() > 0 ? true : false;
00248     }
00249     else {
00250       checkSharedIDs_ = true;
00251     }
00252   }
00253   else {
00254     checkSharedIDs_ = false;
00255   }
00256 }
00257 
00258 //----------------------------------------------------------------------------
00259 void fei::VectorSpace::defineFields(int numFields,
00260                                     const int* fieldIDs,
00261                                     const int* fieldSizes,
00262                                     const int* fieldTypes)
00263 {
00264   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00265     FEI_OSTREAM& os = *output_stream_;
00266     os << dbgprefix_<<"defineFields ";
00267     for(int j=0; j<numFields; ++j) {
00268       os << "{"<<fieldIDs[j] << "," << fieldSizes[j] << "} ";
00269     }
00270     os << FEI_ENDL;
00271   }
00272 
00273   for (int i=0; i<numFields; ++i) {
00274     fieldDatabase_.insert(std::pair<int,unsigned>(fieldIDs[i], fieldSizes[i]));
00275     if (fieldIDs[i] >= 0) {
00276       if (fieldTypes != NULL) {
00277         fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i], fieldTypes[i]);
00278       }
00279       else {
00280         fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i]);
00281       }
00282     }
00283   }
00284 }
00285 
00286 //----------------------------------------------------------------------------
00287 void fei::VectorSpace::defineIDTypes(int numIDTypes,
00288                                      const int* idTypes)
00289 {
00290   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00291     FEI_OSTREAM& os = *output_stream_;
00292     os << dbgprefix_<<"defineIDTypes {";
00293     for(int j=0; j<numIDTypes; ++j) {
00294       os << idTypes[j] << " ";
00295     }
00296     os << "}"<<FEI_ENDL;
00297   }
00298 
00299   int localProc = fei::localProc(comm_);
00300   for (int i=0; i<numIDTypes; ++i) {
00301     int offset = fei::sortedListInsert(idTypes[i], idTypes_);
00302 
00303     if (offset >= 0) {
00304       recordCollections_.insert(recordCollections_.begin()+offset, new snl_fei::RecordCollection(localProc));
00305     }
00306   }
00307 }
00308 
00309 //----------------------------------------------------------------------------
00310 int fei::VectorSpace::addDOFs(int fieldID,
00311                                           int idType,
00312                                           int numIDs,
00313                                           const int* IDs)
00314 {
00315   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00316     FEI_OSTREAM& os = *output_stream_;
00317     os << dbgprefix_<<"addDOFs, fID=" << fieldID
00318        <<", idT="<<idType
00319        << " {";
00320     for(int j=0; j<numIDs; ++j) {
00321       os << IDs[j] << " ";
00322       if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
00323     }
00324     os << "}"<<FEI_ENDL;
00325   }
00326 
00327   if (numIDs <= 0) return(0);
00328 
00329   int idx = fei::binarySearch(idType, idTypes_);
00330   if (idx < 0) ERReturn(-1);
00331 
00332   unsigned fieldSize = getFieldSize(fieldID);
00333   recordCollections_[idx]->initRecords(fieldID, fieldSize,
00334                                        numIDs, IDs,
00335                                        fieldMasks_);
00336   newInitData_ = true;
00337   sharedRecordsSynchronized_ = false;
00338 
00339   return(0);
00340 }
00341 
00342 //----------------------------------------------------------------------------
00343 int fei::VectorSpace::addDOFs(int fieldID,
00344                               int idType,
00345                               int numIDs,
00346                               const int* IDs,
00347                               int* records)
00348 {
00349   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00350     FEI_OSTREAM& os = *output_stream_;
00351     os << dbgprefix_<<"addDOFs*, fID=" << fieldID
00352        <<", idT="<<idType
00353        << " {";
00354     for(int j=0; j<numIDs; ++j) {
00355       os << IDs[j] << " ";
00356       if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
00357     }
00358     os << "}"<<FEI_ENDL;
00359   }
00360 
00361   if (numIDs <= 0) return(0);
00362 
00363   int idx = fei::binarySearch(idType, idTypes_);
00364   if (idx < 0) {
00365     FEI_OSTRINGSTREAM osstr;
00366     osstr << "fei::VectorSpace::addDOFs: error, idType " << idType
00367           << " not recognized. (idTypes need to be initialized via the"
00368           << " method VectorSpace::defineIDTypes)";
00369     throw std::runtime_error(osstr.str());
00370   }
00371 
00372   unsigned fieldSize = getFieldSize(fieldID);
00373   recordCollections_[idx]->initRecords(fieldID, fieldSize,
00374                                        numIDs, IDs,
00375                                        fieldMasks_, records);
00376   newInitData_ = true;
00377   sharedRecordsSynchronized_ = false;
00378 
00379   return(0);
00380 }
00381 
00382 //----------------------------------------------------------------------------
00383 int fei::VectorSpace::addDOFs(int idType,
00384                                           int numIDs,
00385                                           const int* IDs)
00386 {
00387   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00388     FEI_OSTREAM& os = *output_stream_;
00389     os << dbgprefix_<<"addDOFs idT=" <<idType<<" {";
00390     for(int j=0; j<numIDs; ++j) {
00391       os << IDs[j] << " ";
00392       if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
00393     }
00394     os << "}"<<FEI_ENDL;
00395   }
00396 
00397   if (numIDs <= 0) return(0);
00398 
00399   int idx = fei::binarySearch(idType, idTypes_);
00400   if (idx < 0) ERReturn(-1);
00401 
00402   recordCollections_[idx]->initRecords(numIDs, IDs, fieldMasks_);
00403 
00404   newInitData_ = true;
00405   sharedRecordsSynchronized_ = false;
00406 
00407   return(0);
00408 }
00409 
00410 //----------------------------------------------------------------------------
00411 int fei::VectorSpace::addDOFs(int idType,
00412                                           int numIDs,
00413                                           const int* IDs,
00414                                           int* records)
00415 {
00416   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00417     FEI_OSTREAM& os = *output_stream_;
00418     os << dbgprefix_<<"addDOFs* idT=" <<idType<<" {";
00419     for(int j=0; j<numIDs; ++j) {
00420       os << IDs[j] << " ";
00421       if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
00422     }
00423     os << "}"<<FEI_ENDL;
00424   }
00425 
00426   if (numIDs <= 0) return(0);
00427 
00428   int idx = fei::binarySearch(idType, idTypes_);
00429   if (idx < 0) ERReturn(-1);
00430 
00431   recordCollections_[idx]->initRecords(numIDs, IDs,
00432                                        fieldMasks_, records);
00433 
00434   newInitData_ = true;
00435   sharedRecordsSynchronized_ = false;
00436 
00437   return(0);
00438 }
00439 
00440 //----------------------------------------------------------------------------
00441 int fei::VectorSpace::initSharedIDs(int numShared,
00442                                     int idType,
00443                                     const int* sharedIDs,
00444                                     const int* numSharingProcsPerID,
00445                                     const int* sharingProcs)
00446 {
00447   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00448     FEI_OSTREAM& os = *output_stream_;
00449     os << dbgprefix_<<"initSharedIDs n=" <<numShared<<", idT="<<idType<< FEI_ENDL;
00450     int offset = 0;
00451     for(int ns=0; ns<numShared; ++ns) {
00452       os << dbgprefix_<<"#sharedID="<<sharedIDs[ns] << ", nprocs=" << numSharingProcsPerID[ns] << ", procs: ";
00453       for(int sp=0; sp<numSharingProcsPerID[ns]; ++sp) {
00454         os << sharingProcs[offset++] << " ";
00455       }
00456       os << FEI_ENDL;
00457     }
00458     os << FEI_ENDL;
00459   }
00460 
00461   if (numShared == 0) return(0);
00462 
00463   int idx = fei::binarySearch(idType, idTypes_);
00464   if (idx < 0) ERReturn(-1);
00465 
00466   fei::SharedIDs<int>& shIDs = getSharedIDs_private(idType);
00467 
00468   int offset = 0;
00469   for(int i=0; i<numShared; ++i) {
00470     shIDs.addSharedID(sharedIDs[i], numSharingProcsPerID[i],
00471                                 &(sharingProcs[offset]));
00472     offset += numSharingProcsPerID[i];
00473 
00474     fei::Record<int>* rec = recordCollections_[idx]->getRecordWithID(sharedIDs[i]);
00475     if (rec == NULL) {
00476       CHK_ERR( addDOFs(idType, 1, &(sharedIDs[i])) );
00477     }
00478   }
00479 
00480   newInitData_ = true;
00481   sharedRecordsSynchronized_ = false;
00482 
00483   return(0);
00484 }
00485 
00486 //----------------------------------------------------------------------------
00487 int fei::VectorSpace::initSharedIDs(int numShared,
00488                                     int idType,
00489                                     const int* sharedIDs,
00490                                     const int* numSharingProcsPerID,
00491                                     const int* const* sharingProcs)
00492 {
00493   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00494     FEI_OSTREAM& os = *output_stream_;
00495     os << dbgprefix_<<"initSharedIDs n=" <<numShared<<", idT="<<idType<< FEI_ENDL;
00496     for(int ns=0; ns<numShared; ++ns) {
00497       os << dbgprefix_<<"#sharedID="<<sharedIDs[ns] << ", nprocs=" << numSharingProcsPerID[ns] << ", procs: ";
00498       for(int sp=0; sp<numSharingProcsPerID[ns]; ++sp) {
00499         os << sharingProcs[ns][sp] << " ";
00500       }
00501       os << FEI_ENDL;
00502     }
00503     os << FEI_ENDL;
00504   }
00505 
00506   if (numShared == 0) return(0);
00507 
00508   fei::SharedIDs<int>& shIDs = getSharedIDs_private(idType);
00509 
00510   int idx = fei::binarySearch(idType, idTypes_);
00511   if (idx < 0) ERReturn(-1);
00512 
00513   for(int i=0; i<numShared; ++i) {
00514     shIDs.addSharedID(sharedIDs[i], numSharingProcsPerID[i],
00515                                  sharingProcs[i]);
00516 
00517     fei::Record<int>* rec = recordCollections_[idx]->getRecordWithID(sharedIDs[i]);
00518     if (rec == NULL) {
00519       CHK_ERR( addDOFs(idType, 1, &(sharedIDs[i])) );
00520     }
00521   }
00522 
00523   newInitData_ = true;
00524   sharedRecordsSynchronized_ = false;
00525 
00526   return(0);
00527 }
00528 
00529 //----------------------------------------------------------------------------
00530 int fei::VectorSpace::addVectorSpace(fei::VectorSpace* inputSpace)
00531 {
00532   idTypes_ = inputSpace->idTypes_;
00533 
00534   std::map<int,unsigned>::const_iterator
00535     f_iter = inputSpace->fieldDatabase_.begin(),
00536     f_end = inputSpace->fieldDatabase_.end();
00537 
00538   for(; f_iter != f_end; ++f_iter) {
00539     const std::pair<const int,unsigned>& fpair = *f_iter;
00540     int fieldsize = fpair.second;
00541     defineFields(1, &(fpair.first), &fieldsize);
00542   }
00543 
00544   int i, len = inputSpace->fieldMasks_.size();
00545   fieldMasks_.resize(len);
00546   for(i=0; i<len; ++i) {
00547     fieldMasks_[i] = new fei::FieldMask(*(inputSpace->fieldMasks_[i]));
00548   }
00549 
00550   len = inputSpace->recordCollections_.size();
00551   recordCollections_.resize(len);
00552   for(i=0; i<len; ++i) {
00553     recordCollections_[i] =
00554       new snl_fei::RecordCollection(*(inputSpace->recordCollections_[i]));
00555   }
00556 
00557   sharedIDTables_ = inputSpace->sharedIDTables_;
00558 
00559   newInitData_ = true;
00560   sharedRecordsSynchronized_ = false;
00561 
00562   return(0);
00563 }
00564 
00565 //----------------------------------------------------------------------------
00566 void fei::VectorSpace::getSendProcs(std::vector<int>& sendProcs) const
00567 {
00568   sendProcs.clear();
00569   std::set<int> sendSet;
00570 
00571   std::map<int,fei::SharedIDs<int> >::const_iterator
00572     s_it = sharedIDTables_.begin(),
00573     s_end= sharedIDTables_.end();
00574   for(; s_it!=s_end; ++s_it) {
00575     const std::vector<int>& owners = s_it->second.getOwningProcs();
00576     for(size_t i=0; i<owners.size(); ++i) {
00577       sendSet.insert(owners[i]);
00578     }
00579   }
00580 
00581   for(std::set<int>::iterator it=sendSet.begin(); it!=sendSet.end(); ++it) {
00582     if (fei::localProc(comm_) != *it) sendProcs.push_back(*it);
00583   }
00584 }
00585 
00586 //----------------------------------------------------------------------------
00587 fei::SharedIDs<int>& fei::VectorSpace::getSharedIDs_private(int idType)
00588 {
00589   std::map<int,fei::SharedIDs<int> >::iterator
00590     iter = sharedIDTables_.find(idType);
00591 
00592   if (iter == sharedIDTables_.end()) {
00593     iter = sharedIDTables_.insert(std::make_pair(idType,fei::SharedIDs<int>())).first;
00594   }
00595 
00596   return iter->second;
00597 }
00598 
00599 //----------------------------------------------------------------------------
00600 void fei::VectorSpace::compute_shared_ids()
00601 {
00602   //For each ID-type, call the fei::set_shared_ids function, which
00603   //takes a RecordCollection as input and fills a SharedIDs object
00604   //with mappings for which processors share each ID that appears
00605   //on multiple processors.
00606   for(size_t i=0; i<idTypes_.size(); ++i) {
00607     snl_fei::RecordCollection* records = recordCollections_[i];
00608     fei::SharedIDs<int>& shIDs = getSharedIDs_private(idTypes_[i]);
00609 
00610     fei::set_shared_ids(comm_, *records, shIDs);
00611   }
00612 }
00613 
00614 //----------------------------------------------------------------------------
00615 int fei::VectorSpace::initComplete()
00616 {
00617   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00618     FEI_OSTREAM& os = *output_stream_;
00619     os <<dbgprefix_<< "initComplete" << FEI_ENDL;
00620   }
00621 
00622   simpleProblem_ = (fieldMasks_.size()==1) && (fieldMasks_[0]->getNumFields()==1);
00623 
00624   //we need to know if any processor has newInitData_.
00625   int localInitData = newInitData_ ? 1 : 0;
00626   int globalInitData = 0;
00627   CHK_ERR( fei::GlobalMax(comm_, localInitData, globalInitData) );
00628   newInitData_ = globalInitData > 0 ? true : false;
00629 
00630   if (!newInitData_) {
00631     return(0);
00632   }
00633 
00634   //if running on multiple procs and initSharedIDs(...) has not been
00635   //called, then we need to specifically figure out which IDs are
00636   //shared and which procs share them.
00637 
00638   bool need_to_compute_shared_ids = false;
00639   if (fei::numProcs(comm_) > 1 && sharedIDTables_.size() < 1) {
00640     need_to_compute_shared_ids = true;
00641   }
00642 
00643   int localNeedToCompute = need_to_compute_shared_ids ? 1 : 0;
00644   int globalNeedToCompute = 0;
00645   CHK_ERR( fei::GlobalMin(comm_, localNeedToCompute, globalNeedToCompute) );
00646   need_to_compute_shared_ids = globalNeedToCompute==1 ? true : false;
00647   if (need_to_compute_shared_ids) {
00648     compute_shared_ids();
00649   }
00650 
00651   //setOwners_lowestSharing is a local operation (no communication), which
00652   //assumes that each processor holds CORRECT (globally symmetric)
00653   //shared-id/sharing-proc tables. No correctness-checking is performed here.
00654   setOwners_lowestSharing();
00655 
00656   //synchronizeSharedRecords ensures that each sharing processor has the same
00657   //view of the shared records, with respect to the layout of fields, which
00658   //determines how many DOFs and equation-numbers reside at each ID.
00659   //This involves inter-processor communication.
00660   if ( synchronizeSharedRecords() != 0) {
00661     return(-1);
00662   }
00663 
00664   //calculateGlobalIndices is also a global operation.
00665   if ( calculateGlobalIndices() != 0) {
00666     return(-1);
00667   }
00668 
00669   //finally we need to exchange global indices for shared records. i.e., the
00670   //processors that own shared records, need to send the global indices for
00671   //those records to the sharing-but-not-owning processors.
00672   if (fei::numProcs(comm_) > 1) {
00673     if ( exchangeGlobalIndices() != 0) {
00674       return(-1);
00675     }
00676   }
00677 
00678   newInitData_ = false;
00679 
00680   return(0);
00681 }
00682 
00683 //----------------------------------------------------------------------------
00684 int fei::VectorSpace::getGlobalIndex(int idType,
00685                                      int ID,
00686                                      int fieldID,
00687                                      int fieldOffset,
00688                                      int whichComponentOfField,
00689                                      int& globalIndex)
00690 {
00691   int idindex = fei::binarySearch(idType, idTypes_);
00692   if (idindex < 0) return(-1);
00693 
00694   unsigned fieldSize = 0;
00695   if (fieldOffset > 0) {
00696     fieldSize = getFieldSize(fieldID);
00697   }
00698 
00699   try {
00700     globalIndex = recordCollections_[idindex]->getGlobalIndex(ID,
00701                                                          fieldID,
00702                                                          fieldSize,
00703                                                          fieldOffset,
00704                                                          whichComponentOfField,
00705                                                          &eqnNumbers_[0]);
00706   }
00707   catch (std::runtime_error& exc) {
00708     FEI_OSTRINGSTREAM osstr;
00709     osstr << "VectorSpace::getGlobalIndex caught exception: " << exc.what();
00710     fei::console_out() << osstr.str()<<FEI_ENDL;
00711     ERReturn(-1);
00712   }
00713 
00714   return(0);
00715 }
00716 
00717 //----------------------------------------------------------------------------
00718 int fei::VectorSpace::getGlobalIndex(int idType,
00719                                      int ID,
00720                                      int fieldID,
00721                                      int& globalIndex)
00722 {
00723   return( getGlobalIndex(idType, ID, fieldID, 0, 0, globalIndex) );
00724 }
00725 
00726 //----------------------------------------------------------------------------
00727 int fei::VectorSpace::getGlobalBlkIndex(int idType,
00728                                         int ID,
00729                                         int& globalBlkIndex)
00730 {
00731   int idindex = fei::binarySearch(idType, idTypes_);
00732   if (idindex < 0) return(-1);
00733 
00734   CHK_ERR(recordCollections_[idindex]->getGlobalBlkIndex(ID, globalBlkIndex));
00735 
00736   return(0);
00737 }
00738 
00739 //----------------------------------------------------------------------------
00740 int fei::VectorSpace::getGlobalIndices(int numIDs,
00741                                        const int* IDs,
00742                                        int idType,
00743                                        int fieldID,
00744                                        int* globalIndices)
00745 {
00746   int idindex = fei::binarySearch(idType, idTypes_);
00747   if (idindex < 0) return(-1);
00748 
00749   unsigned fieldSize = getFieldSize(fieldID);
00750   int offset = 0;
00751 
00752   for(int i=0; i<numIDs; ++i) {
00753     try {
00754       globalIndices[offset] =
00755         recordCollections_[idindex]->getGlobalIndex(IDs[i],
00756                                                     fieldID,
00757                                                     fieldSize,
00758                                                     0, 0,
00759                                                     &eqnNumbers_[0]);
00760       if (fieldSize>1) {
00761         int eqn = globalIndices[offset];
00762         for(unsigned j=1; j<fieldSize; ++j) {
00763           globalIndices[offset+j] = eqn+j;
00764         }
00765       }
00766     }
00767     catch (...) {
00768       for(unsigned j=0; j<fieldSize; ++j) {
00769         globalIndices[offset+j] = -1;
00770       }
00771     }
00772 
00773     offset += fieldSize;
00774   }
00775 
00776   return(0);
00777 }
00778 
00779 //----------------------------------------------------------------------------
00780 int fei::VectorSpace::getGlobalBlkIndices(int numIDs,
00781                                           const int* IDs,
00782                                           int idType,
00783                                           int* globalBlkIndices)
00784 {
00785   int err;
00786   int idindex = fei::binarySearch(idType, idTypes_);
00787   if (idindex < 0) return(-1);
00788 
00789   int offset = 0;
00790 
00791   for(int i=0; i<numIDs; ++i) {
00792     err = recordCollections_[idindex]->getGlobalBlkIndex(IDs[i],
00793                                                      globalBlkIndices[offset]);
00794     if (err != 0) {
00795       globalBlkIndices[offset] = -1;
00796     }
00797     ++offset;
00798   }
00799 
00800   return(0);
00801 }
00802 
00803 //----------------------------------------------------------------------------
00804 int fei::VectorSpace::getGlobalIndices(int numIDs,
00805                                        const int* IDs,
00806                                        const int* idTypes,
00807                                        const int* fieldIDs,
00808                                        int* globalIndices)
00809 {
00810   int err;
00811   int offset = 0;
00812   for(int i=0; i<numIDs; ++i) {
00813     unsigned fieldSize = getFieldSize(fieldIDs[i]);
00814     err = getGlobalIndex(idTypes[i], IDs[i], fieldIDs[i], 0, 0,
00815                          globalIndices[offset]);
00816     if (err) {
00817       for(unsigned j=1; j<fieldSize; ++j) {
00818         globalIndices[offset+j] = -1;
00819       }
00820     }
00821     else {
00822       if (fieldSize>1) {
00823         int eqn = globalIndices[offset];
00824         for(unsigned j=1; j<fieldSize; ++j) {
00825           globalIndices[offset+j] = eqn+j;
00826         }
00827       }
00828     }
00829     offset += fieldSize;
00830   }
00831 
00832   return(0);
00833 }
00834 
00835 //----------------------------------------------------------------------------
00836 void fei::VectorSpace::getGlobalBlkIndices(const fei::Pattern* pattern,
00837                                            const fei::Record<int>*const* records,
00838                                            std::vector<int>& indices)
00839 {
00840   int numRecords = pattern->getNumIDs();
00841   indices.resize(numRecords);
00842   int numIndices;
00843   getGlobalBlkIndices(numRecords, records, numRecords, &indices[0],
00844                       numIndices);
00845 }
00846 
00847 //----------------------------------------------------------------------------
00848 void fei::VectorSpace::getGlobalIndices(const fei::Pattern* pattern,
00849                                         const fei::Record<int>*const* records,
00850                                         std::vector<int>& indices)
00851 {
00852   int numRecords = pattern->getNumIDs();
00853   int numIndices = pattern->getNumIndices();
00854   indices.resize(numIndices);
00855   int* indices_ptr = &indices[0];
00856 
00857   fei::Pattern::PatternType pType = pattern->getPatternType();
00858 
00859   if (pType == fei::Pattern::GENERAL ||
00860       pType == fei::Pattern::SINGLE_IDTYPE) {
00861     const int* numFieldsPerID = pattern->getNumFieldsPerID();
00862     const int* fieldIDs = pattern->getFieldIDs();
00863     int totalNumFields = pattern->getTotalNumFields();
00864 
00865     std::vector<int> fieldSizes(totalNumFields);
00866 
00867     for(int j=0; j<totalNumFields; ++j) {
00868       fieldSizes[j] = getFieldSize(fieldIDs[j]);
00869     }
00870 
00871     getGlobalIndices(numRecords, records, numFieldsPerID,
00872                      fieldIDs, &(fieldSizes[0]),
00873                      numIndices, indices_ptr, numIndices);
00874   }
00875   else if (pType == fei::Pattern::SIMPLE) {
00876     const int* fieldIDs = pattern->getFieldIDs();
00877 
00878     int fieldID = fieldIDs[0];
00879     unsigned fieldSize = getFieldSize(fieldID);
00880 
00881     getGlobalIndices(numRecords, records,
00882                      fieldID, fieldSize,
00883                      numIndices, indices_ptr, numIndices);
00884   }
00885   else if (pType == fei::Pattern::NO_FIELD) {
00886     getGlobalBlkIndices(numRecords, records, numIndices, indices_ptr, numIndices);
00887   }
00888 }
00889 
00890 //----------------------------------------------------------------------------
00891 void fei::VectorSpace::getGlobalIndicesL(const fei::Pattern* pattern,
00892                                         const int* records,
00893                                         std::vector<int>& indices)
00894 {
00895   int numRecords = pattern->getNumIDs();
00896   int numIndices = pattern->getNumIndices();
00897   const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
00898   indices.resize(numIndices);
00899   int* indices_ptr = &indices[0];
00900 
00901   fei::Pattern::PatternType pType = pattern->getPatternType();
00902 
00903   if (pType == fei::Pattern::GENERAL ||
00904       pType == fei::Pattern::SINGLE_IDTYPE) {
00905     const int* numFieldsPerID = pattern->getNumFieldsPerID();
00906     const int* fieldIDs = pattern->getFieldIDs();
00907     int totalNumFields = pattern->getTotalNumFields();
00908 
00909     std::vector<int> fieldSizes(totalNumFields);
00910 
00911     for(int j=0; j<totalNumFields; ++j) {
00912       fieldSizes[j] = getFieldSize(fieldIDs[j]);
00913     }
00914 
00915     getGlobalIndicesL(numRecords, recordCollections, records, numFieldsPerID,
00916                      fieldIDs, &(fieldSizes[0]),
00917                      numIndices, indices_ptr, numIndices);
00918   }
00919   else if (pType == fei::Pattern::SIMPLE) {
00920     const int* fieldIDs = pattern->getFieldIDs();
00921     const int* idTypes = pattern->getIDTypes();
00922 
00923     int fieldID = fieldIDs[0];
00924     unsigned fieldSize = getFieldSize(fieldID);
00925 
00926     getGlobalIndicesL(numRecords, recordCollections, records,
00927                      fieldID, fieldSize,
00928                      numIndices, indices_ptr, numIndices);
00929   }
00930   else if (pType == fei::Pattern::NO_FIELD) {
00931     const int* idTypes = pattern->getIDTypes();
00932     getGlobalBlkIndicesL(numRecords, recordCollections, records, numIndices, indices_ptr, numIndices);
00933   }
00934 }
00935 
00936 //----------------------------------------------------------------------------
00937 void fei::VectorSpace::getGlobalIndices(int numRecords,
00938                                         const fei::Record<int>*const* records,
00939                                         int fieldID,
00940                                         int fieldSize,
00941                                         int indicesAllocLen,
00942                                         int* indices,
00943                                         int& numIndices)
00944 {
00945   numIndices = 0;
00946   int* eqnPtr = &eqnNumbers_[0];
00947 
00948   int len = numRecords;
00949   if (len*fieldSize >= indicesAllocLen) {
00950     len = indicesAllocLen/fieldSize;
00951   }
00952 
00953   if (fieldSize == 1 && simpleProblem_) {
00954     for(int i=0; i<len; ++i) {
00955       const fei::Record<int>* record = records[i];
00956       indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
00957     }    
00958     return;
00959   }
00960 
00961   if (fieldSize == 1) {
00962     for(int i=0; i<len; ++i) {
00963       const fei::Record<int>* record = records[i];
00964 
00965       int eqnOffset = 0;
00966       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
00967 
00968       const fei::FieldMask* fieldMask = record->getFieldMask();
00969       fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
00970       indices[numIndices++] = eqnNumbers[eqnOffset];
00971     }
00972   }
00973   else {
00974     for(int i=0; i<len; ++i) {
00975       const fei::Record<int>* record = records[i];
00976 
00977       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
00978 
00979       int eqnOffset = 0;
00980       if (!simpleProblem_) {
00981         const fei::FieldMask* fieldMask = record->getFieldMask();
00982         fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
00983       }
00984 
00985       for(int fs=0; fs<fieldSize; ++fs) {
00986         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
00987       }
00988     }
00989   }
00990 }
00991 
00992 //----------------------------------------------------------------------------
00993 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
00994                              const snl_fei::RecordCollection*const* recordCollections,
00995                                         const int* records,
00996                                         int fieldID,
00997                                         int fieldSize,
00998                                         int indicesAllocLen,
00999                                         int* indices,
01000                                         int& numIndices)
01001 {
01002   numIndices = 0;
01003   int* eqnPtr = &eqnNumbers_[0];
01004 
01005   int len = numRecords;
01006   if (len*fieldSize >= indicesAllocLen) {
01007     len = indicesAllocLen/fieldSize;
01008   }
01009 
01010   if (fieldSize == 1 && simpleProblem_) {
01011     for(int i=0; i<len; ++i) {
01012       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01013       indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
01014     }    
01015     return;
01016   }
01017 
01018   if (fieldSize == 1) {
01019     for(int i=0; i<len; ++i) {
01020       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01021 
01022       int eqnOffset = 0;
01023       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01024 
01025       const fei::FieldMask* fieldMask = record->getFieldMask();
01026       fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01027       indices[numIndices++] = eqnNumbers[eqnOffset];
01028     }
01029   }
01030   else {
01031     for(int i=0; i<len; ++i) {
01032       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01033 
01034       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01035 
01036       int eqnOffset = 0;
01037       if (!simpleProblem_) {
01038         const fei::FieldMask* fieldMask = record->getFieldMask();
01039         fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01040       }
01041 
01042       for(int fs=0; fs<fieldSize; ++fs) {
01043         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01044       }
01045     }
01046   }
01047 }
01048 
01049 //----------------------------------------------------------------------------
01050 void fei::VectorSpace::getGlobalIndices(int numRecords,
01051                                         const fei::Record<int>*const* records,
01052                                         const int* numFieldsPerID,
01053                                         const int* fieldIDs,
01054                                         const int* fieldSizes,
01055                                         int indicesAllocLen,
01056                                         int* indices,
01057                                         int& numIndices)
01058 {
01059   numIndices = 0;
01060   int fld_offset = 0;
01061   int* eqnPtr = &eqnNumbers_[0];
01062 
01063   for(int i=0; i<numRecords; ++i) {
01064     const fei::Record<int>* record = records[i];
01065 
01066     const fei::FieldMask* fieldMask = record->getFieldMask();
01067     int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
01068 
01069     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
01070       int eqnOffset = 0;
01071       if (!simpleProblem_) {
01072         fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
01073       }
01074 
01075       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01076         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01077       }
01078 
01079       ++fld_offset;
01080     }
01081   }
01082 }
01083 
01084 //----------------------------------------------------------------------------
01085 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
01086                                         const snl_fei::RecordCollection*const* recordCollections,
01087                                         const int* records,
01088                                         const int* numFieldsPerID,
01089                                         const int* fieldIDs,
01090                                         const int* fieldSizes,
01091                                         int indicesAllocLen,
01092                                         int* indices,
01093                                         int& numIndices)
01094 {
01095   numIndices = 0;
01096   int fld_offset = 0;
01097   int* eqnPtr = &eqnNumbers_[0];
01098 
01099   for(int i=0; i<numRecords; ++i) {
01100     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01101 
01102     const fei::FieldMask* fieldMask = record->getFieldMask();
01103     int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
01104 
01105     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
01106       int eqnOffset = 0;
01107       if (!simpleProblem_) {
01108         fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
01109       }
01110 
01111       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01112         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01113       }
01114 
01115       ++fld_offset;
01116     }
01117   }
01118 }
01119 
01120 //----------------------------------------------------------------------------
01121 void fei::VectorSpace::getGlobalBlkIndices(int numRecords,
01122                                            const fei::Record<int>*const* records,
01123                                            int indicesAllocLen,
01124                                            int* indices,
01125                                            int& numIndices)
01126 {
01127   numIndices = 0;
01128   for(int i=0; i<numRecords; ++i) {
01129     if (numIndices < indicesAllocLen) {
01130       indices[numIndices++] = records[i]->getNumber();
01131     }
01132     else ++numIndices;
01133   }
01134 }
01135 
01136 //----------------------------------------------------------------------------
01137 void fei::VectorSpace::getGlobalBlkIndicesL(int numRecords,
01138                              const snl_fei::RecordCollection*const* recordCollections,
01139                                            const int* records,
01140                                            int indicesAllocLen,
01141                                            int* indices,
01142                                            int& numIndices)
01143 {
01144   numIndices = 0;
01145   for(int i=0; i<numRecords; ++i) {
01146     if (numIndices < indicesAllocLen) {
01147       indices[numIndices++] = recordCollections[i]->getRecordWithLocalID(records[i])->getNumber();
01148     }
01149     else ++numIndices;
01150   }
01151 }
01152 
01153 //----------------------------------------------------------------------------
01154 int fei::VectorSpace::getGlobalIndex(int idType,
01155                                      int ID,
01156                                      int& globalIndex)
01157 {
01158   int idindex = fei::binarySearch(idType, idTypes_);
01159   if (idindex < 0) return(-1);
01160 
01161   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01162   if (record == NULL) {
01163     ERReturn(-1);
01164   }
01165 
01166   const int* eqnNums = &eqnNumbers_[0]
01167                      + record->getOffsetIntoEqnNumbers();
01168 
01169   if (eqnNums != NULL) { globalIndex = eqnNums[0]; return(0); }
01170   return(-1);
01171 }
01172 
01173 //----------------------------------------------------------------------------
01174 int fei::VectorSpace::getNumDegreesOfFreedom(int idType,
01175                                              int ID)
01176 {
01177   int idindex = fei::binarySearch(idType, idTypes_);
01178   if (idindex < 0) return(0);
01179 
01180   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01181   if (record == NULL) {
01182     ERReturn(-1);
01183   }
01184 
01185   return( record->getFieldMask()->getNumIndices() );
01186 }
01187 
01188 //----------------------------------------------------------------------------
01189 int fei::VectorSpace::getNumFields()
01190 {
01191   return(fieldDatabase_.size());
01192 }
01193 
01194 //----------------------------------------------------------------------------
01195 void fei::VectorSpace::getFields(std::vector<int>& fieldIDs)
01196 {
01197   unsigned numFields = fieldDatabase_.size();
01198 
01199   fieldIDs.resize(numFields);
01200 
01201   fei::copyKeysToArray<std::map<int,unsigned> >(fieldDatabase_, numFields,
01202                                                     &fieldIDs[0]);
01203 }
01204 
01205 //----------------------------------------------------------------------------
01206 size_t fei::VectorSpace::getNumIDTypes()
01207 {
01208   return(idTypes_.size());
01209 }
01210 
01211 //----------------------------------------------------------------------------
01212 void fei::VectorSpace::getIDTypes(std::vector<int>& idTypes) const
01213 {
01214   size_t numIDTypes = idTypes_.size();
01215   idTypes.resize(numIDTypes);
01216   for(size_t i=0; i<numIDTypes; ++i) idTypes[i] = idTypes_[i];
01217 }
01218 
01219 //----------------------------------------------------------------------------
01220 int fei::VectorSpace::getNumFields(int idType,
01221                                    int ID)
01222 {
01223   int idindex = fei::binarySearch(idType, idTypes_);
01224   if (idindex < 0) return(0);
01225 
01226   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01227   if (record == NULL) {
01228     ERReturn(-1);
01229   }
01230 
01231   return( record->getFieldMask()->getNumFields() );
01232 }
01233 
01234 //----------------------------------------------------------------------------
01235 bool fei::VectorSpace::isLocal(int idType, int ID)
01236 {
01237   int idindex = fei::binarySearch(idType, idTypes_);
01238   if (idindex < 0) return(false);
01239 
01240   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01241   if (record == NULL) {
01242     return(false);
01243   }
01244 
01245   return(true);
01246 }
01247 
01248 //----------------------------------------------------------------------------
01249 bool fei::VectorSpace::isLocallyOwned(int idType, int ID)
01250 {
01251   int idindex = fei::binarySearch(idType, idTypes_);
01252   if (idindex < 0) return(false);
01253 
01254   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01255   if (record == NULL) {
01256     return(false);
01257   }
01258   if (record->getOwnerProc() == fei::localProc(comm_)) {
01259     return(true);
01260   }
01261 
01262   return(false);
01263 }
01264 
01265 //----------------------------------------------------------------------------
01266 unsigned fei::VectorSpace::getFieldSize(int fieldID)
01267 {
01268   std::map<int,unsigned>::const_iterator
01269     f_iter = fieldDatabase_.find(fieldID);
01270 
01271   if (f_iter == fieldDatabase_.end()) {
01272     FEI_OSTRINGSTREAM osstr;
01273     osstr << "fei::VectorSpace";
01274     if (name_.length() > 0) {
01275       osstr << "(name: "<<name_<<")";
01276     }
01277     osstr << "::getFieldSize: fieldID " << fieldID << " not found.";
01278     throw std::runtime_error(osstr.str());
01279   }
01280 
01281   return((*f_iter).second);
01282 }
01283 
01284 //----------------------------------------------------------------------------
01285 void fei::VectorSpace::getFields(int idType, int ID, std::vector<int>& fieldIDs)
01286 {
01287   int idindex = fei::binarySearch(idType, idTypes_);
01288   if (idindex < 0) {
01289     fieldIDs.resize(0);
01290     return;
01291   }
01292 
01293   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01294   if (record == NULL) {
01295     voidERReturn;
01296   }
01297 
01298   fei::FieldMask* fieldMask = record->getFieldMask();
01299   std::vector<int>& maskFieldIDs = fieldMask->getFieldIDs();
01300   int numFields = maskFieldIDs.size();
01301   fieldIDs.resize(numFields);
01302   for(int i=0; i<numFields; ++i) {
01303     fieldIDs[i] = maskFieldIDs[i];
01304   }
01305 }
01306 
01307 //----------------------------------------------------------------------------
01308 void fei::VectorSpace::getGlobalIndexOffsets(std::vector<int>& globalOffsets) const
01309 {
01310   globalOffsets = globalOffsets_;
01311 }
01312 
01313 //----------------------------------------------------------------------------
01314 void fei::VectorSpace::getGlobalBlkIndexOffsets(std::vector<int>& globalBlkOffsets) const
01315 {
01316   globalBlkOffsets = globalIDOffsets_;
01317 }
01318 
01319 //----------------------------------------------------------------------------
01320 int fei::VectorSpace::getOwnerProcPtIndex(int globalIndex)
01321 {
01322   if (globalIndex < 0) return(-1);
01323 
01324   unsigned len = globalOffsets_.size();
01325   for(int i=0; i<(int)(len-1); ++i) {
01326     if (globalIndex < globalOffsets_[i+1]) {
01327       return(i);
01328     }
01329   }
01330 
01331   return(-1);
01332 }
01333 
01334 //----------------------------------------------------------------------------
01335 int fei::VectorSpace::getOwnerProcBlkIndex(int globalIndex)
01336 {
01337   if (globalIndex < 0) return(-1);
01338 
01339   unsigned len = globalOffsets_.size();
01340   for(int i=0; i<(int)(len-1); ++i) {
01341     if (globalIndex < globalIDOffsets_[i+1]) {
01342       return(i);
01343     }
01344   }
01345 
01346   return(-1);
01347 }
01348 
01349 //----------------------------------------------------------------------------
01350 int fei::VectorSpace::getNumOwnedAndSharedIDs(int idType)
01351 {
01352   int idx = fei::binarySearch(idType, idTypes_);
01353   if (idx < 0) return(0);
01354 
01355   return( recordCollections_[idx]->getNumRecords() );
01356 }
01357 
01358 //----------------------------------------------------------------------------
01359 int fei::VectorSpace::getNumOwnedIDs(int idType)
01360 {
01361   int idx = fei::binarySearch(idType, idTypes_);
01362   if (idx < 0) return(0);
01363 
01364   fei::RecordAttributeCounter attrCounter(fei::localProc(comm_));
01365   runRecords(attrCounter);
01366 
01367   return( attrCounter.numLocallyOwnedIDs_ );
01368 }
01369 
01370 //----------------------------------------------------------------------------
01371 int fei::VectorSpace::getNumSharedIDs(int idType, int& numShared)
01372 {
01373   numShared = sharedIDTables_[idType].getSharedIDs().size();
01374 
01375   return(0);
01376 }
01377 
01378 //----------------------------------------------------------------------------
01379 int fei::VectorSpace::getOwnedAndSharedIDs(int idType,
01380                                       int lenList,
01381                                       int* IDs,
01382                                       int& numLocalIDs)
01383 {
01384   int idx = fei::binarySearch(idType, idTypes_);
01385   if (idx < 0) return(-1);
01386 
01387   snl_fei::RecordCollection* records = recordCollections_[idx];
01388   std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
01389   std::map<int,int>::iterator
01390     it = global_to_local.begin(),
01391     end = global_to_local.end();
01392   numLocalIDs = records->getNumRecords();
01393   int len = numLocalIDs;
01394   if (lenList < len) len = lenList;
01395   int i=0;
01396   for(; it!=end; ++it) {
01397     int local_id = it->second;
01398     IDs[i++] = records->getRecordWithLocalID(local_id)->getID();
01399     if (i >= len) break;
01400   }
01401 
01402   return(0);
01403 }
01404 
01405 //----------------------------------------------------------------------------
01406 int fei::VectorSpace::getOwnedIDs(int idType,
01407                                          int lenList,
01408                                          int* IDs,
01409                                          int& numLocalIDs)
01410 {
01411   int idx = fei::binarySearch(idType, idTypes_);
01412   if (idx < 0) return(-1);
01413 
01414   snl_fei::RecordCollection* records = recordCollections_[idx];
01415 
01416   std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
01417   std::map<int,int>::iterator
01418     it = global_to_local.begin(),
01419     end = global_to_local.end();
01420 
01421   numLocalIDs = 0;
01422 
01423   for(; it!=end; ++it) {
01424     fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
01425 
01426     if (thisrecord.getOwnerProc() == fei::localProc(comm_)) {
01427       if (numLocalIDs < lenList) {
01428         IDs[numLocalIDs] = thisrecord.getID();
01429       }
01430       ++numLocalIDs;
01431     }
01432   }
01433 
01434   return(0);
01435 }
01436 
01437 //----------------------------------------------------------------------------
01438 int fei::VectorSpace::getNumIndices_SharedAndOwned() const
01439 {
01440   return(eqnNumbers_.size());
01441 }
01442 
01443 //----------------------------------------------------------------------------
01444 int fei::VectorSpace::getIndices_SharedAndOwned(std::vector<int>& globalIndices) const
01445 {
01446   if (eqnNumbers_.size() == 0) {
01447     globalIndices.resize(0);
01448     return(0);
01449   }
01450 
01451   size_t numIndices = eqnNumbers_.size();
01452   const int* indicesPtr = &eqnNumbers_[0];
01453 
01454   globalIndices.resize(numIndices);
01455   for(size_t i=0; i<numIndices; ++i) {
01456     globalIndices[i] = indicesPtr[i];
01457   }
01458 
01459   return(0);
01460 }
01461 
01462 //----------------------------------------------------------------------------
01463 int fei::VectorSpace::getNumBlkIndices_SharedAndOwned(int& numBlkIndices) const
01464 {
01465   numBlkIndices = 0;
01466   for(size_t i=0; i<recordCollections_.size(); ++i) {
01467     numBlkIndices += recordCollections_[i]->getNumRecords();
01468   }
01469 
01470   return(0);
01471 }
01472 
01473 //----------------------------------------------------------------------------
01474 int fei::VectorSpace::getBlkIndices_SharedAndOwned(int lenBlkIndices,
01475                                                    int* globalBlkIndices,
01476                                                    int* blkSizes,
01477                                                    int& numBlkIndices)
01478 {
01479   if (!sharedRecordsSynchronized_) {
01480     numBlkIndices = 0;
01481     return(-1);
01482   }
01483 
01484   fei::BlkIndexAccessor blkIndAccessor(lenBlkIndices,
01485                                   globalBlkIndices, blkSizes);
01486   runRecords(blkIndAccessor);
01487 
01488   numBlkIndices = blkIndAccessor.numBlkIndices_;
01489 
01490   return(0);
01491 }
01492 
01493 //----------------------------------------------------------------------------
01494 int fei::VectorSpace::getGlobalNumIndices() const
01495 {
01496   if (globalOffsets_.size() < 1) return(0);
01497   return(globalOffsets_[globalOffsets_.size()-1]);
01498 }
01499 
01500 //----------------------------------------------------------------------------
01501 int fei::VectorSpace::getNumIndices_Owned() const
01502 {
01503   if (!sharedRecordsSynchronized_) {
01504     return(-1);
01505   }
01506 
01507   int localProc = fei::localProc(comm_);
01508   int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01509 
01510   return(numIndices);
01511 }
01512 
01513 //----------------------------------------------------------------------------
01514 int fei::VectorSpace::getIndices_Owned(std::vector<int>& globalIndices) const
01515 {
01516   if (!sharedRecordsSynchronized_) {
01517     globalIndices.clear();
01518     return(-1);
01519   }
01520 
01521   int localProc = fei::localProc(comm_);
01522   int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01523 
01524   globalIndices.resize(numIndices);
01525 
01526   int firstOffset = globalOffsets_[localProc];
01527   for(int i=0; i<numIndices; ++i) {
01528     globalIndices[i] = firstOffset+i;
01529   }
01530 
01531   return(0);
01532 }
01533 
01534 //----------------------------------------------------------------------------
01535 int fei::VectorSpace::getIndices_Owned(int lenIndices, int* globalIndices, int& numIndices) const
01536 {
01537   if (!sharedRecordsSynchronized_) {
01538     numIndices = 0;
01539     return(-1);
01540   }
01541 
01542   int localProc = fei::localProc(comm_);
01543   numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01544 
01545   int len = lenIndices >= numIndices ? numIndices : lenIndices;
01546 
01547   int firstOffset = globalOffsets_[localProc];
01548   for(int i=0; i<len; ++i) {
01549     globalIndices[i] = firstOffset+i;
01550   }
01551 
01552   return(0);
01553 }
01554 
01555 //----------------------------------------------------------------------------
01556 int fei::VectorSpace::getNumBlkIndices_Owned() const
01557 {
01558   if (!sharedRecordsSynchronized_) {
01559     return(-1);
01560   }
01561 
01562   int localProc = fei::localProc(comm_);
01563   int numBlkIndices = globalIDOffsets_[localProc+1]-globalIDOffsets_[localProc];
01564 
01565   return(numBlkIndices);
01566 }
01567 
01568 //----------------------------------------------------------------------------
01569 int fei::VectorSpace::getGlobalNumBlkIndices() const
01570 {
01571   int numBlkIndices = 0;
01572   unsigned len = globalIDOffsets_.size();
01573   if (len > 0) {
01574     numBlkIndices = globalIDOffsets_[len-1];
01575   }
01576 
01577   return(numBlkIndices);
01578 }
01579 
01580 //----------------------------------------------------------------------------
01581 int fei::VectorSpace::getBlkIndices_Owned(int lenBlkIndices,
01582                                           int* globalBlkIndices,
01583                                           int* blkSizes,
01584                                           int& numBlkIndices)
01585 {
01586   if (!sharedRecordsSynchronized_) {
01587     numBlkIndices = 0;
01588     return(-1);
01589   }
01590 
01591   int localProc = fei::localProc(comm_);
01592   fei::BlkIndexAccessor blkIndAccessor(localProc, lenBlkIndices,
01593                                            globalBlkIndices, blkSizes);
01594   runRecords(blkIndAccessor);
01595 
01596   numBlkIndices = blkIndAccessor.numBlkIndices_;
01597 
01598   return(0);
01599 }
01600 
01601 //----------------------------------------------------------------------------
01602 int fei::VectorSpace::getRecordCollection(int idType,
01603                                           snl_fei::RecordCollection*& records)
01604 {
01605   int idx = fei::binarySearch(idType, idTypes_);
01606   if (idx < 0) return(-1);
01607 
01608   records = recordCollections_[idx];
01609   return(0);
01610 }
01611 
01612 //----------------------------------------------------------------------------
01613 int fei::VectorSpace::getRecordCollection(int idType,
01614                                           const snl_fei::RecordCollection*& records) const
01615 {
01616   int idx = fei::binarySearch(idType, idTypes_);
01617   if (idx < 0) return(-1);
01618 
01619   records = recordCollections_[idx];
01620   return(0);
01621 }
01622 
01623 //----------------------------------------------------------------------------
01624 void fei::VectorSpace::setOwners_lowestSharing()
01625 {
01626   //first, add localProc to each of the sharing-proc lists, in case it wasn't
01627   //included via initSharedIDs().
01628   int localProc = fei::localProc(comm_);
01629 
01630   std::map<int, fei::SharedIDs<int> >::iterator
01631     s_iter = sharedIDTables_.begin(), s_end = sharedIDTables_.end();
01632 
01633   for(; s_iter != s_end; ++s_iter) {
01634     fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
01635 
01636     fei::SharedIDs<int>::map_type::iterator
01637       t_iter = shid_table.begin(),
01638       t_end = shid_table.end();
01639 
01640     for(; t_iter != t_end; ++t_iter) {
01641       std::set<int>& shProcs = t_iter->second;
01642       shProcs.insert(localProc);
01643     }
01644   }
01645 
01646   //now set the owningProcs for the SharedIDs records, and the owning procs on
01647   //the appropriate records in the recordCollections. Set the owner to be the
01648   //lowest-numbered sharing proc in all cases.
01649   s_iter = sharedIDTables_.begin();
01650   for(; s_iter != s_end; ++s_iter) {
01651     fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
01652 
01653     std::vector<int>& owningProcs = s_iter->second.getOwningProcs();
01654 
01655     int len = shid_table.size();
01656     owningProcs.resize(len);
01657     int j = 0;
01658 
01659     fei::SharedIDs<int>::map_type::iterator
01660       t_iter = shid_table.begin(),
01661       t_end = shid_table.end();
01662 
01663     for(; t_iter != t_end; ++t_iter, ++j) {
01664       std::set<int>& shProcs = (*t_iter).second;
01665       int lowest = *(shProcs.begin());
01666       owningProcs[j] = lowest;
01667     }
01668     
01669     int idx = fei::binarySearch(s_iter->first, idTypes_);
01670     if (idx < 0) {
01671       throw std::runtime_error("internal error in fei::VectorSapce::setOwners_lowestSharing");
01672     }
01673 
01674     if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01675       recordCollections_[idx]->setDebugOutput(output_stream_);
01676     }
01677 
01678     recordCollections_[idx]->setOwners_lowestSharing(s_iter->second);
01679   }
01680 }
01681 
01682 //----------------------------------------------------------------------------
01683 int fei::VectorSpace::calculateGlobalIndices()
01684 {
01685   int localProc = fei::localProc(comm_);
01686   int numProcs = fei::numProcs(comm_);
01687   std::vector<int> localOffsets(numProcs +1, 0);
01688   std::vector<int> localIDOffsets(numProcs+1, 0);
01689   globalOffsets_.resize(numProcs + 1, 0);
01690 
01691   globalIDOffsets_.assign(globalIDOffsets_.size(), 0);
01692 
01693   //first we'll calculate the number of local degrees of freedom, and the
01694   //number of local identifiers.
01695   fei::RecordAttributeCounter counter(localProc);
01696   runRecords(counter);
01697 
01698   int numLocalDOF = counter.numLocalDOF_;
01699   int numLocalIDs = counter.numLocalIDs_;
01700   int numRemoteSharedDOF = counter.numRemoteSharedDOF_;
01701 
01702   eqnNumbers_.resize(numLocalDOF+numRemoteSharedDOF);
01703 
01704   localOffsets[localProc] = numLocalDOF;
01705   CHK_MPI( fei::GlobalMax(comm_, localOffsets, globalOffsets_) );
01706 
01707   localIDOffsets[localProc] = numLocalIDs;
01708   CHK_MPI( fei::GlobalMax(comm_, localIDOffsets, globalIDOffsets_) );
01709 
01710   //Now globalOffsets_ contains numLocalDOF for proc i, in the i-th position.
01711   //(and similarly for globalIDOffsets_)
01712   //So all we need to do is turn that data into global-offsets (i.e., the
01713   //starting global-offset for each processor).
01714   int localOffset = 0;
01715   int localIDOffset = 0;
01716   for(int p=0; p<numProcs; ++p) {
01717     numLocalDOF = globalOffsets_[p];
01718     globalOffsets_[p] = localOffset;
01719     localOffset += numLocalDOF;
01720     numLocalIDs = globalIDOffsets_[p];
01721     globalIDOffsets_[p] = localIDOffset;
01722     localIDOffset += numLocalIDs;
01723   }
01724   globalOffsets_[numProcs] = localOffset;
01725   globalIDOffsets_[numProcs] = localIDOffset;
01726 
01727   firstLocalOffset_ = globalOffsets_[localProc];
01728   lastLocalOffset_ = globalOffsets_[localProc+1] - 1;
01729 
01730   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01731     FEI_OSTREAM& os = *output_stream_;
01732     os <<dbgprefix_<< "  firstLocalOffset_: " << firstLocalOffset_ << ", lastLocalOffset_: "
01733       << lastLocalOffset_ << FEI_ENDL;
01734   }
01735 
01736   //Now we're ready to set the equation-numbers on all local records.
01737   CHK_ERR( setLocalEqnNumbers() );
01738 
01739   return(0);
01740 }
01741 
01742 //----------------------------------------------------------------------------
01743 int fei::VectorSpace::synchronizeSharedRecords()
01744 {
01745   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01746     FEI_OSTREAM& os = *output_stream_;
01747     os<<dbgprefix_<<"#synchronizeSharedRecords num-field-masks: "<<fieldMasks_.size()<<FEI_ENDL;
01748     for(unsigned fm=0; fm<fieldMasks_.size(); ++fm) {
01749       os << dbgprefix_<<"#     maskID["<<fm<<"]: " << fieldMasks_[fm]->getMaskID() << FEI_ENDL;
01750     }
01751   }
01752 
01753   bool safetyCheck = checkSharedIDs_;
01754 
01755   int numProcs = fei::numProcs(comm_);
01756   int localProc = fei::localProc(comm_);
01757   int numShTables = sharedIDTables_.size();
01758 
01759   if (numProcs < 2) {
01760     sharedRecordsSynchronized_ = true;
01761     return(0);
01762   }
01763 
01764   if (safetyCheck) {
01765     if (localProc == 0 && output_level_ >= fei::BRIEF_LOGS) {
01766       FEI_COUT << "fei::VectorSpace: global consistency-check of shared ID"
01767            << " data (involves all-to-all communication). This is done "
01768            << "only if requested by parameter 'FEI_CHECK_SHARED_IDS true'."
01769            << FEI_ENDL;
01770     }
01771 
01772     int totalNumShTables = 0;
01773     CHK_ERR( fei::GlobalSum(comm_, numShTables, totalNumShTables) );
01774     if (totalNumShTables != numShTables * numProcs) {
01775       //not all processors have the same number of shared-id tables. This means
01776       //that one or more processors is 'empty', which may not be an error. But
01777       //it does mean that the safety check can't be performed because the safety
01778       //check involves all-to-all communication and if one or more processors
01779       //don't participate then we'll hang...
01780       safetyCheck = false;
01781     }
01782   }
01783 
01784   //Create a list of fei::comm_maps which will be the communication-patterns
01785   //for each of the sharedIDTables. A sharedIDTable maps lists of processors
01786   //to each shared ID. The communication-pattern will be a transpose of
01787   //that, mapping lists of IDs to owning or sharing processors.
01788   //
01789 
01790   int local_err = 0;
01791 
01792   for(size_t i=0; i<idTypes_.size(); ++i) {
01793 
01794     fei::SharedIDs<int>& shared = getSharedIDs_private(idTypes_[i]);
01795     
01796     //now create/initialize ownerPatterns which map owning processors to lists
01797     //of ids that are shared locally, and sharerPatterns which map sharing
01798     //processors to lists of ids that are owned locally.
01799     fei::comm_map* ownerPattern = new fei::comm_map(1, numProcs);
01800 
01801     std::map<int,fei::comm_map*>::iterator iter = ownerPatterns_.find(idTypes_[i]);
01802     if (iter == ownerPatterns_.end()) {
01803       ownerPatterns_.insert(std::make_pair(idTypes_[i], ownerPattern));
01804     }
01805     else {
01806       delete iter->second;
01807       iter->second = ownerPattern;
01808     }
01809 
01810     fei::comm_map* sharerPattern = new fei::comm_map(1, numProcs);
01811 
01812     iter = sharerPatterns_.find(idTypes_[i]);
01813     if (iter == sharerPatterns_.end()) {
01814       sharerPatterns_.insert(std::make_pair(idTypes_[i], sharerPattern));
01815     }
01816     else {
01817       delete iter->second;
01818       iter->second = sharerPattern;
01819     }
01820 
01821     std::vector<int>& owningProcs = shared.getOwningProcs();
01822 
01823     fei::SharedIDs<int>::map_type& shtable = shared.getSharedIDs();
01824 
01825     fei::SharedIDs<int>::map_type::iterator
01826       s_iter = shtable.begin(),
01827       s_end = shtable.end();
01828     
01829     int j = 0;
01830     for(; s_iter != s_end; ++s_iter, ++j) {
01831       int ID = (*s_iter).first;
01832       std::set<int>& shProcs = (*s_iter).second;
01833 
01834       int owner = owningProcs[j];
01835 
01836       if (owner == localProc) {
01837         std::set<int>::const_iterator
01838           p_iter = shProcs.begin(),
01839           p_end = shProcs.end();
01840         for(; p_iter != p_end; ++p_iter) {
01841           if (*p_iter != localProc) {
01842             sharerPattern->addIndices(*p_iter, 1, &ID);
01843           }
01844         }
01845       }
01846       else {
01847         ownerPattern->addIndices(owner, 1, &ID);
01848       }
01849     }
01850 
01851     if (safetyCheck) {
01852       fei::comm_map* checkPattern = NULL;
01853       CHK_ERR( fei::mirrorCommPattern(comm_, ownerPattern, checkPattern) );
01854       fei::Barrier(comm_);
01855       if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01856         FEI_OSTREAM& os = *output_stream_;
01857         fei::comm_map::map_type& owner_map = ownerPattern->getMap();
01858         int numKeys = owner_map.size();
01859         fei::comm_map::map_type::const_iterator
01860           omap_iter = owner_map.begin(),
01861           omap_end = owner_map.end();
01862 
01863         os << dbgprefix_<<"#  synchronizeSharedRecords" << FEI_ENDL
01864            << dbgprefix_<<"#  ownerPattern, num-procs-to-send-to: " << numKeys << FEI_ENDL;
01865         for(int sk=0; omap_iter != omap_end; ++sk, ++omap_iter) {
01866           os << dbgprefix_<<"#    sendProc["<<sk<<"]: " << omap_iter->first << " IDs: ";
01867           fei::comm_map::row_type::const_iterator
01868             val_iter = omap_iter->second->begin(),
01869             val_end =  omap_iter->second->end();
01870           for(; val_iter != val_end; ++val_iter) {
01871             os << *val_iter << " ";
01872           }
01873           os << FEI_ENDL;
01874         }
01875 
01876         fei::comm_map::map_type& check_map = checkPattern->getMap();
01877         int numCheckKeys = check_map.size();
01878         fei::comm_map::map_type::const_iterator
01879           cmap_iter = check_map.begin(),
01880           cmap_end = check_map.end();
01881 
01882         os <<dbgprefix_<< "#  synchronizeSharedRecords" << FEI_ENDL
01883            <<dbgprefix_<< "#  checkPattern (send mirror), num-procs: "
01884            << numCheckKeys << FEI_ENDL;
01885         for(int sk=0; cmap_iter != cmap_end; ++sk, ++cmap_iter) {
01886           os <<dbgprefix_<< "#    proc["<<sk<<"]: " << cmap_iter->first << " IDs: ";
01887           fei::comm_map::row_type::const_iterator
01888           val_iter = cmap_iter->second->begin(),
01889             val_end = cmap_iter->second->end();
01890           for(; val_iter != val_end; ++val_iter) {
01891             os << *val_iter << " ";
01892           }
01893           os << FEI_ENDL;
01894         }
01895       }
01896 
01897       int err = 0;
01898       bool quiet = false;
01899       if (!checkPattern->equal(*sharerPattern, quiet)) {
01900         //fei::console_out() << "shared-ID safety-check failed..."
01901         //     << FEI_ENDL;
01902         err = -1;
01903       }
01904       int globalErr = 0;
01905       CHK_ERR( fei::GlobalSum(comm_, err, globalErr) );
01906 
01907       delete checkPattern;
01908 
01909       if (globalErr != 0) {
01910         return(globalErr);
01911       }
01912     }
01913 
01914     local_err += exchangeFieldInfo(ownerPattern, sharerPattern,
01915                                    recordCollections_[i], fieldMasks_);
01916   }
01917 
01918   int global_err = 0;
01919   CHK_ERR( fei::GlobalSum(comm_, local_err, global_err) );
01920 
01921   if (global_err != 0) {
01922     ERReturn(-1);
01923   }
01924 
01925   sharedRecordsSynchronized_ = true;
01926 
01927   return(0);
01928 }
01929 
01930 //----------------------------------------------------------------------------
01931 int fei::VectorSpace::exchangeGlobalIndices()
01932 {
01933   for(size_t i=0; i<idTypes_.size(); ++i) {
01934     snl_fei::RecordMsgHandler recmsghndlr(fei::localProc(comm_),
01935                                  recordCollections_[i], *ptBlkMap_,
01936                                           fieldMasks_, eqnNumbers_);
01937     recmsghndlr.setTask(snl_fei::RecordMsgHandler::_EqnNumbers_);
01938 
01939     recmsghndlr.setSendPattern(sharerPatterns_[idTypes_[i]]);
01940     recmsghndlr.setRecvPattern(ownerPatterns_[idTypes_[i]]);
01941     CHK_ERR( fei::exchange(comm_, &recmsghndlr) );
01942   }
01943 
01944   return(0);
01945 }
01946 
01947 //----------------------------------------------------------------------------
01948 void fei::VectorSpace::runRecords(fei::Record_Operator<int>& record_op)
01949 {
01950   for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
01951     snl_fei::RecordCollection* records = recordCollections_[rec];
01952     std::map<int,int>& g2l = records->getGlobalToLocalMap();
01953     std::map<int,int>::iterator
01954       it = g2l.begin(),
01955       end= g2l.end();
01956 
01957     for(; it!=end; ++it) {
01958       fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
01959 
01960       record_op(thisrecord);
01961     }
01962   }
01963 }
01964 
01965 //----------------------------------------------------------------------------
01966 int fei::VectorSpace::setLocalEqnNumbers()
01967 {
01968   int proc = fei::localProc(comm_);
01969   int eqnNumber = firstLocalOffset_;
01970   int idNumber = globalIDOffsets_[proc];
01971 
01972   int numProcs = fei::numProcs(comm_);
01973   int localProc = fei::localProc(comm_);
01974 
01975   if (ptBlkMap_ != NULL) {
01976     delete ptBlkMap_;
01977   }
01978   ptBlkMap_ = new snl_fei::PointBlockMap;
01979 
01980   int maxNumIndices = 0;
01981   for(unsigned i=0; i<fieldMasks_.size(); ++i) {
01982     if (fieldMasks_[i]->getNumIndices() > maxNumIndices) {
01983       maxNumIndices = fieldMasks_[i]->getNumIndices();
01984     }
01985   }
01986 
01987   if (maxNumIndices == 1) {
01988     ptBlkMap_->setPtEqualBlk();
01989   }
01990 
01991   FEI_OSTREAM* id2eqnStream = NULL;
01992   if (output_level_ >= fei::BRIEF_LOGS) {
01993     std::string path = fei::LogManager::getLogManager().getOutputPath();
01994     if (path == "") path = ".";
01995     FEI_OSTRINGSTREAM osstr;
01996     osstr << path << "/fei_ID2Eqn";
01997     if (name_.size() > 0) osstr << "_" << name_;
01998     osstr << "." <<numProcs<< "." << localProc;
01999 
02000     id2eqnStream = new FEI_OFSTREAM(osstr.str().c_str(), IOS_OUT);
02001     FEI_OSTREAM& os = *id2eqnStream;
02002     os << "# Each line contains:\n#   ID   blk-eqn   eqn" << FEI_ENDL;
02003   }
02004 
02005   int eqnNumberOffset = 0;
02006   int maxNumDOF = 0;
02007 
02008   for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
02009     snl_fei::RecordCollection* records = recordCollections_[rec];
02010 
02011     const std::map<int,int>& rmap = records->getGlobalToLocalMap();
02012     std::map<int,int>::const_iterator
02013       it = rmap.begin(), it_end = rmap.end();
02014 
02015     int* eqnNumPtr = &eqnNumbers_[0];
02016 
02017     for(; it!=it_end; ++it) {
02018       fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
02019 
02020       fei::FieldMask* mask = thisrecord.getFieldMask();
02021 
02022       thisrecord.setOffsetIntoEqnNumbers(eqnNumberOffset);
02023 
02024       int owner = thisrecord.getOwnerProc();
02025       if (owner == proc) {
02026         thisrecord.setNumber(idNumber++);
02027       }
02028 
02029       int* eqnNumbers = eqnNumPtr
02030                       + thisrecord.getOffsetIntoEqnNumbers();
02031 
02032       int numDOF = mask->getNumIndices();
02033       eqnNumberOffset += numDOF;
02034 
02035       if (output_level_ >= fei::BRIEF_LOGS) {
02036         for(int ii=0; ii<numDOF; ++ii) {
02037           if (isLogEqn(eqnNumber+ii) && output_stream_ != NULL) {
02038             FEI_OSTREAM& os = *output_stream_;
02039             os <<dbgprefix_<< "setLocalEqnNumbers: ID "<<thisrecord.getID()
02040                << " <--> eqn " << eqnNumber+ii<<FEI_ENDL;
02041           }
02042         }
02043       }
02044 
02045       if (owner == proc) {
02046         int offset = 0;
02047         for(int n=0; n<numDOF; ++n) {
02048           eqnNumbers[offset++] = eqnNumber++;
02049         }
02050       }
02051 
02052       if (numDOF > maxNumDOF) maxNumDOF = numDOF;
02053 
02054       if (owner == proc) {
02055         int thiseqn = eqnNumber-numDOF;
02056         int thisrecordnumber = thisrecord.getNumber();
02057         if (maxNumIndices > 1) {
02058           CHK_ERR( ptBlkMap_->setEqn(thiseqn, thisrecordnumber, numDOF) );
02059           if (numDOF > 1) {
02060             for(int i=1; i<numDOF; ++i) {
02061               CHK_ERR( ptBlkMap_->setEqn(thiseqn+i, thisrecordnumber, numDOF) );
02062             }
02063           }
02064         }
02065       }
02066 
02067       if (id2eqnStream != NULL) {
02068         if (owner == proc) {
02069           FEI_OSTREAM& os = *id2eqnStream;
02070           for(int n=0; n<numDOF; ++n) {
02071             os << thisrecord.getID() << " " << thisrecord.getNumber() << " "
02072                << eqnNumber-numDOF+n<<FEI_ENDL;
02073           }
02074         }
02075       }
02076     }
02077 
02078   }
02079 
02080   ptBlkMap_->setMaxBlkEqnSize(maxNumDOF);
02081 
02082   int globalMaxNumDOF;
02083   CHK_ERR( fei::GlobalMax(comm_, maxNumDOF, globalMaxNumDOF) );
02084 
02085   if (globalMaxNumDOF == 1) {
02086     ptBlkMap_->setPtEqualBlk();
02087   }
02088 
02089   delete id2eqnStream;
02090 
02091   return(0);
02092 }
02093 
02094 //----------------------------------------------------------------------------
02095 int fei::VectorSpace::exchangeFieldInfo(fei::comm_map* ownerPattern,
02096                                         fei::comm_map* sharerPattern,
02097                                   snl_fei::RecordCollection* recordCollection,
02098                                       std::vector<fei::FieldMask*>& fieldMasks)
02099 {
02100   //ownerPattern maps owning processors to lists of IDs that we (the local
02101   //processor) share.
02102   //
02103   //sharerPattern maps sharing processors to lists of IDs that we own.
02104   //
02105   //In this function we need to perform these tasks:
02106   //1. processors exchange and combine their sets of field-masks so that all
02107   //   processors will have the super-set of field-masks that they need.
02108   //2. sharing processors send maskIDs for their shared IDs to the owning
02109   //   processors. The owning processors then combine masks as necessary to
02110   //   make sure each shared ID has the union of field-masks that are held by
02111   //   each of the sharing processors. all owning processors now send their
02112   //   maskIDs out to the sharing processors. At this point all processors
02113   //   should have the same view of the masks that belong on shared IDs.
02114   //3. exchange info describing inactive DOF offsets for shared records.
02115   //
02116 
02117   int numProcs = fei::numProcs(comm_);
02118   if (numProcs < 2) return(0);
02119 
02120   if (ptBlkMap_ == 0) ptBlkMap_ = new snl_fei::PointBlockMap;
02121 
02122   snl_fei::RecordMsgHandler recMsgHandler(fei::localProc(comm_),
02123                                           recordCollection,
02124                                           *ptBlkMap_,
02125                                           fieldMasks, eqnNumbers_);
02126 
02127   //Step 1a.
02128   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
02129 
02130   recMsgHandler.setSendPattern(ownerPattern);
02131   recMsgHandler.setRecvPattern(sharerPattern);
02132   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02133 
02134   //Step 2a.
02135   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
02136 
02137   recMsgHandler.setSendPattern(ownerPattern);
02138   recMsgHandler.setRecvPattern(sharerPattern);
02139   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02140 
02141   //Step 1b.
02142   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
02143 
02144   recMsgHandler.setSendPattern(sharerPattern);
02145   recMsgHandler.setRecvPattern(ownerPattern);
02146   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02147 
02148   //Step 2b.
02149   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
02150 
02151   recMsgHandler.setSendPattern(sharerPattern);
02152   recMsgHandler.setRecvPattern(ownerPattern);
02153   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02154 
02155   return(0);
02156 }
02157 
02158 //----------------------------------------------------------------------------
02159 fei::FieldDofMap<int>&
02160 fei::VectorSpace::getFieldDofMap()
02161 {
02162   return fieldDofMap_;
02163 }
02164 
02165 //----------------------------------------------------------------------------
02166 void fei::VectorSpace::setName(const char* name)
02167 {
02168   if (name == NULL) return;
02169 
02170   if (name_ == name) return;
02171 
02172   name_ = name;
02173   dbgprefix_ = "VecSpc_"+name_+": ";
02174 }
02175 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends