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