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   globalIndex = recordCollections_[idindex]->getGlobalIndex(ID,
00749                                                          fieldID,
00750                                                          fieldSize,
00751                                                          fieldOffset,
00752                                                          whichComponentOfField,
00753                                                          &eqnNumbers_[0]);
00754   int return_value = globalIndex >= 0 ? 0 : -1;
00755   return return_value;
00756 }
00757 
00758 //----------------------------------------------------------------------------
00759 int fei::VectorSpace::getGlobalIndex(int idType,
00760                                      int ID,
00761                                      int fieldID,
00762                                      int& globalIndex)
00763 {
00764   return( getGlobalIndex(idType, ID, fieldID, 0, 0, globalIndex) );
00765 }
00766 
00767 //----------------------------------------------------------------------------
00768 int fei::VectorSpace::getGlobalBlkIndex(int idType,
00769                                         int ID,
00770                                         int& globalBlkIndex)
00771 {
00772   int idindex = fei::binarySearch(idType, idTypes_);
00773   if (idindex < 0) return(-1);
00774 
00775   CHK_ERR(recordCollections_[idindex]->getGlobalBlkIndex(ID, globalBlkIndex));
00776 
00777   return(0);
00778 }
00779 
00780 //----------------------------------------------------------------------------
00781 int fei::VectorSpace::getGlobalIndices(int numIDs,
00782                                        const int* IDs,
00783                                        int idType,
00784                                        int fieldID,
00785                                        int* globalIndices)
00786 {
00787   int idindex = fei::binarySearch(idType, idTypes_);
00788   if (idindex < 0) return(-1);
00789 
00790   unsigned fieldSize = getFieldSize(fieldID);
00791   int offset = 0;
00792 
00793   for(int i=0; i<numIDs; ++i) {
00794     globalIndices[offset] = recordCollections_[idindex]->getGlobalIndex(IDs[i],
00795                                                     fieldID, fieldSize,
00796                                                     0, 0, &eqnNumbers_[0]);
00797     if (fieldSize>1) {
00798       if (globalIndices[offset] >= 0) {
00799         int eqn = globalIndices[offset];
00800         for(unsigned j=1; j<fieldSize; ++j) {
00801           globalIndices[offset+j] = eqn+j;
00802         }
00803       }
00804       else {
00805         for(unsigned j=0; j<fieldSize; ++j) {
00806           globalIndices[offset+j] = -1;
00807         }
00808       }
00809     }
00810 
00811     offset += fieldSize;
00812   }
00813 
00814   return(0);
00815 }
00816 
00817 //----------------------------------------------------------------------------
00818 int fei::VectorSpace::getGlobalIndicesLocalIDs(int numIDs,
00819                                        const int* localIDs,
00820                                        int idType,
00821                                        int fieldID,
00822                                        int* globalIndices)
00823 {
00824   int idindex = fei::binarySearch(idType, idTypes_);
00825   if (idindex < 0) return(-1);
00826 
00827   unsigned fieldSize = getFieldSize(fieldID);
00828   int offset = 0;
00829 
00830   for(int i=0; i<numIDs; ++i) {
00831     globalIndices[offset] = recordCollections_[idindex]->getGlobalIndexLocalID(localIDs[i],
00832                                                     fieldID, fieldSize,
00833                                                     0, 0, &eqnNumbers_[0]);
00834     if (globalIndices[offset] >= 0) {
00835       if (fieldSize>1) {
00836         int eqn = globalIndices[offset];
00837         for(unsigned j=1; j<fieldSize; ++j) {
00838           globalIndices[offset+j] = eqn+j;
00839         }
00840       }
00841     }
00842     else {
00843       for(unsigned j=0; j<fieldSize; ++j) {
00844         globalIndices[offset+j] = -1;
00845       }
00846     }
00847 
00848     offset += fieldSize;
00849   }
00850 
00851   return(0);
00852 }
00853 
00854 //----------------------------------------------------------------------------
00855 int fei::VectorSpace::getGlobalBlkIndices(int numIDs,
00856                                           const int* IDs,
00857                                           int idType,
00858                                           int* globalBlkIndices)
00859 {
00860   int err;
00861   int idindex = fei::binarySearch(idType, idTypes_);
00862   if (idindex < 0) return(-1);
00863 
00864   int offset = 0;
00865 
00866   for(int i=0; i<numIDs; ++i) {
00867     err = recordCollections_[idindex]->getGlobalBlkIndex(IDs[i],
00868                                                      globalBlkIndices[offset]);
00869     if (err != 0) {
00870       globalBlkIndices[offset] = -1;
00871     }
00872     ++offset;
00873   }
00874 
00875   return(0);
00876 }
00877 
00878 //----------------------------------------------------------------------------
00879 int fei::VectorSpace::getGlobalIndices(int numIDs,
00880                                        const int* IDs,
00881                                        const int* idTypes,
00882                                        const int* fieldIDs,
00883                                        int* globalIndices)
00884 {
00885   int err;
00886   int offset = 0;
00887   for(int i=0; i<numIDs; ++i) {
00888     unsigned fieldSize = getFieldSize(fieldIDs[i]);
00889     err = getGlobalIndex(idTypes[i], IDs[i], fieldIDs[i], 0, 0,
00890                          globalIndices[offset]);
00891     if (err) {
00892       for(unsigned j=1; j<fieldSize; ++j) {
00893         globalIndices[offset+j] = -1;
00894       }
00895     }
00896     else {
00897       if (fieldSize>1) {
00898         int eqn = globalIndices[offset];
00899         for(unsigned j=1; j<fieldSize; ++j) {
00900           globalIndices[offset+j] = eqn+j;
00901         }
00902       }
00903     }
00904     offset += fieldSize;
00905   }
00906 
00907   return(0);
00908 }
00909 
00910 //----------------------------------------------------------------------------
00911 void fei::VectorSpace::getGlobalBlkIndices(const fei::Pattern* pattern,
00912                                            const fei::Record<int>*const* records,
00913                                            std::vector<int>& indices)
00914 {
00915   int numRecords = pattern->getNumIDs();
00916   indices.resize(numRecords);
00917   int numIndices;
00918   getGlobalBlkIndices(numRecords, records, numRecords, &indices[0],
00919                       numIndices);
00920 }
00921 
00922 //----------------------------------------------------------------------------
00923 void fei::VectorSpace::getGlobalIndices(const fei::Pattern* pattern,
00924                                         const fei::Record<int>*const* records,
00925                                         std::vector<int>& indices)
00926 {
00927   int numRecords = pattern->getNumIDs();
00928   int numIndices = pattern->getNumIndices();
00929   indices.resize(numIndices);
00930   int* indices_ptr = &indices[0];
00931 
00932   fei::Pattern::PatternType pType = pattern->getPatternType();
00933 
00934   if (pType == fei::Pattern::GENERAL ||
00935       pType == fei::Pattern::SINGLE_IDTYPE) {
00936     const int* numFieldsPerID = pattern->getNumFieldsPerID();
00937     const int* fieldIDs = pattern->getFieldIDs();
00938     int totalNumFields = pattern->getTotalNumFields();
00939 
00940     std::vector<int> fieldSizes(totalNumFields);
00941 
00942     for(int j=0; j<totalNumFields; ++j) {
00943       fieldSizes[j] = getFieldSize(fieldIDs[j]);
00944     }
00945 
00946     getGlobalIndices(numRecords, records, numFieldsPerID,
00947                      fieldIDs, &(fieldSizes[0]),
00948                      numIndices, indices_ptr, numIndices);
00949   }
00950   else if (pType == fei::Pattern::SIMPLE) {
00951     const int* fieldIDs = pattern->getFieldIDs();
00952 
00953     int fieldID = fieldIDs[0];
00954     unsigned fieldSize = getFieldSize(fieldID);
00955 
00956     getGlobalIndices(numRecords, records,
00957                      fieldID, fieldSize,
00958                      numIndices, indices_ptr, numIndices);
00959   }
00960   else if (pType == fei::Pattern::NO_FIELD) {
00961     getGlobalBlkIndices(numRecords, records, numIndices, indices_ptr, numIndices);
00962   }
00963 }
00964 
00965 //----------------------------------------------------------------------------
00966 void fei::VectorSpace::getGlobalIndicesL(const fei::Pattern* pattern,
00967                                         const int* records,
00968                                         std::vector<int>& indices)
00969 {
00970   int numRecords = pattern->getNumIDs();
00971   int numIndices = pattern->getNumIndices();
00972   const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
00973   indices.resize(numIndices);
00974   int* indices_ptr = &indices[0];
00975 
00976   fei::Pattern::PatternType pType = pattern->getPatternType();
00977 
00978   if (pType == fei::Pattern::GENERAL ||
00979       pType == fei::Pattern::SINGLE_IDTYPE) {
00980     const int* numFieldsPerID = pattern->getNumFieldsPerID();
00981     const int* fieldIDs = pattern->getFieldIDs();
00982     int totalNumFields = pattern->getTotalNumFields();
00983 
00984     std::vector<int> fieldSizes(totalNumFields);
00985 
00986     for(int j=0; j<totalNumFields; ++j) {
00987       fieldSizes[j] = getFieldSize(fieldIDs[j]);
00988     }
00989 
00990     getGlobalIndicesL(numRecords, recordCollections, records, numFieldsPerID,
00991                      fieldIDs, &(fieldSizes[0]),
00992                      numIndices, indices_ptr, numIndices);
00993   }
00994   else if (pType == fei::Pattern::SIMPLE) {
00995     const int* fieldIDs = pattern->getFieldIDs();
00996 
00997     int fieldID = fieldIDs[0];
00998     unsigned fieldSize = getFieldSize(fieldID);
00999 
01000     getGlobalIndicesL(numRecords, recordCollections, records,
01001                      fieldID, fieldSize,
01002                      numIndices, indices_ptr, numIndices);
01003   }
01004   else if (pType == fei::Pattern::NO_FIELD) {
01005     getGlobalBlkIndicesL(numRecords, recordCollections, records, numIndices, indices_ptr, numIndices);
01006   }
01007 }
01008 
01009 //----------------------------------------------------------------------------
01010 void fei::VectorSpace::getGlobalIndices(int numRecords,
01011                                         const fei::Record<int>*const* records,
01012                                         int fieldID,
01013                                         int fieldSize,
01014                                         int indicesAllocLen,
01015                                         int* indices,
01016                                         int& numIndices)
01017 {
01018   numIndices = 0;
01019   int* eqnPtr = &eqnNumbers_[0];
01020 
01021   int len = numRecords;
01022   if (len*fieldSize >= indicesAllocLen) {
01023     len = indicesAllocLen/fieldSize;
01024   }
01025 
01026   if (fieldSize == 1 && simpleProblem_) {
01027     for(int i=0; i<len; ++i) {
01028       const fei::Record<int>* record = records[i];
01029       indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
01030     }    
01031     return;
01032   }
01033 
01034   if (fieldSize == 1) {
01035     for(int i=0; i<len; ++i) {
01036       const fei::Record<int>* record = records[i];
01037 
01038       int eqnOffset = 0;
01039       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01040 
01041       const fei::FieldMask* fieldMask = record->getFieldMask();
01042       int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01043       indices[numIndices++] = err == 0 ? eqnNumbers[eqnOffset] : -1;
01044     }
01045   }
01046   else {
01047     for(int i=0; i<len; ++i) {
01048       const fei::Record<int>* record = records[i];
01049 
01050       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01051 
01052       int eqnOffset = 0;
01053       if (!simpleProblem_) {
01054         const fei::FieldMask* fieldMask = record->getFieldMask();
01055         int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01056         if (err != 0) {
01057           for(int fs=0; fs<fieldSize; ++fs) {
01058             indices[numIndices++] = -1;
01059           }
01060           continue;
01061         }
01062       }
01063 
01064       for(int fs=0; fs<fieldSize; ++fs) {
01065         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01066       }
01067     }
01068   }
01069 }
01070 
01071 //----------------------------------------------------------------------------
01072 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
01073                              const snl_fei::RecordCollection*const* recordCollections,
01074                                         const int* records,
01075                                         int fieldID,
01076                                         int fieldSize,
01077                                         int indicesAllocLen,
01078                                         int* indices,
01079                                         int& numIndices)
01080 {
01081   numIndices = 0;
01082   int* eqnPtr = &eqnNumbers_[0];
01083 
01084   int len = numRecords;
01085   if (len*fieldSize >= indicesAllocLen) {
01086     len = indicesAllocLen/fieldSize;
01087   }
01088 
01089   if (fieldSize == 1 && simpleProblem_) {
01090     for(int i=0; i<len; ++i) {
01091       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01092       indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
01093     }    
01094     return;
01095   }
01096 
01097   if (fieldSize == 1) {
01098     for(int i=0; i<len; ++i) {
01099       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01100 
01101       int eqnOffset = 0;
01102       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01103 
01104       const fei::FieldMask* fieldMask = record->getFieldMask();
01105       int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01106       indices[numIndices++] = err == 0 ? eqnNumbers[eqnOffset] : -1;
01107     }
01108   }
01109   else {
01110     for(int i=0; i<len; ++i) {
01111       const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01112 
01113       int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
01114 
01115       int eqnOffset = 0;
01116       if (!simpleProblem_) {
01117         const fei::FieldMask* fieldMask = record->getFieldMask();
01118         int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
01119         if (err != 0) {
01120           for(int fs=0; fs<fieldSize; ++fs) {
01121             indices[numIndices++] = -1;
01122           }
01123           continue;
01124         }
01125       }
01126 
01127       for(int fs=0; fs<fieldSize; ++fs) {
01128         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01129       }
01130     }
01131   }
01132 }
01133 
01134 //----------------------------------------------------------------------------
01135 void fei::VectorSpace::getGlobalIndices(int numRecords,
01136                                         const fei::Record<int>*const* records,
01137                                         const int* numFieldsPerID,
01138                                         const int* fieldIDs,
01139                                         const int* fieldSizes,
01140                                         int indicesAllocLen,
01141                                         int* indices,
01142                                         int& numIndices)
01143 {
01144   numIndices = 0;
01145   int fld_offset = 0;
01146   int* eqnPtr = &eqnNumbers_[0];
01147 
01148   for(int i=0; i<numRecords; ++i) {
01149     const fei::Record<int>* record = records[i];
01150 
01151     const fei::FieldMask* fieldMask = record->getFieldMask();
01152     int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
01153 
01154     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
01155       int eqnOffset = 0;
01156       if (!simpleProblem_) {
01157         int err = fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
01158         if (err != 0) {
01159           for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01160             indices[numIndices++] = -1;
01161           }
01162           continue;
01163         }
01164       }
01165 
01166       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01167         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01168       }
01169 
01170       ++fld_offset;
01171     }
01172   }
01173 }
01174 
01175 //----------------------------------------------------------------------------
01176 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
01177                                         const snl_fei::RecordCollection*const* recordCollections,
01178                                         const int* records,
01179                                         const int* numFieldsPerID,
01180                                         const int* fieldIDs,
01181                                         const int* fieldSizes,
01182                                         int indicesAllocLen,
01183                                         int* indices,
01184                                         int& numIndices)
01185 {
01186   numIndices = 0;
01187   int fld_offset = 0;
01188   int* eqnPtr = &eqnNumbers_[0];
01189 
01190   for(int i=0; i<numRecords; ++i) {
01191     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
01192 
01193     const fei::FieldMask* fieldMask = record->getFieldMask();
01194     int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
01195 
01196     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
01197       int eqnOffset = 0;
01198       if (!simpleProblem_) {
01199         int err = fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
01200         if (err != 0) {
01201           for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01202             indices[numIndices++] = -1;
01203           }
01204           continue;
01205         }
01206       }
01207 
01208       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
01209         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
01210       }
01211 
01212       ++fld_offset;
01213     }
01214   }
01215 }
01216 
01217 //----------------------------------------------------------------------------
01218 void fei::VectorSpace::getGlobalBlkIndices(int numRecords,
01219                                            const fei::Record<int>*const* records,
01220                                            int indicesAllocLen,
01221                                            int* indices,
01222                                            int& numIndices)
01223 {
01224   numIndices = 0;
01225   for(int i=0; i<numRecords; ++i) {
01226     if (numIndices < indicesAllocLen) {
01227       indices[numIndices++] = records[i]->getNumber();
01228     }
01229     else ++numIndices;
01230   }
01231 }
01232 
01233 //----------------------------------------------------------------------------
01234 void fei::VectorSpace::getGlobalBlkIndicesL(int numRecords,
01235                              const snl_fei::RecordCollection*const* recordCollections,
01236                                            const int* records,
01237                                            int indicesAllocLen,
01238                                            int* indices,
01239                                            int& numIndices)
01240 {
01241   numIndices = 0;
01242   for(int i=0; i<numRecords; ++i) {
01243     if (numIndices < indicesAllocLen) {
01244       indices[numIndices++] = recordCollections[i]->getRecordWithLocalID(records[i])->getNumber();
01245     }
01246     else ++numIndices;
01247   }
01248 }
01249 
01250 //----------------------------------------------------------------------------
01251 int fei::VectorSpace::getGlobalIndex(int idType,
01252                                      int ID,
01253                                      int& globalIndex)
01254 {
01255   int idindex = fei::binarySearch(idType, idTypes_);
01256   if (idindex < 0) return(-1);
01257 
01258   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01259   if (record == NULL) {
01260     ERReturn(-1);
01261   }
01262 
01263   const int* eqnNums = &eqnNumbers_[0]
01264                      + record->getOffsetIntoEqnNumbers();
01265 
01266   if (eqnNums != NULL) { globalIndex = eqnNums[0]; return(0); }
01267   return(-1);
01268 }
01269 
01270 //----------------------------------------------------------------------------
01271 int fei::VectorSpace::getNumDegreesOfFreedom(int idType,
01272                                              int ID)
01273 {
01274   int idindex = fei::binarySearch(idType, idTypes_);
01275   if (idindex < 0) return(0);
01276 
01277   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01278   if (record == NULL) {
01279     ERReturn(-1);
01280   }
01281 
01282   return( record->getFieldMask()->getNumIndices() );
01283 }
01284 
01285 //----------------------------------------------------------------------------
01286 int fei::VectorSpace::getNumFields()
01287 {
01288   return(fieldDatabase_.size());
01289 }
01290 
01291 //----------------------------------------------------------------------------
01292 void fei::VectorSpace::getFields(std::vector<int>& fieldIDs)
01293 {
01294   unsigned numFields = fieldDatabase_.size();
01295 
01296   fieldIDs.resize(numFields);
01297 
01298   fei::copyKeysToArray<std::map<int,unsigned> >(fieldDatabase_, numFields,
01299                                                     &fieldIDs[0]);
01300 }
01301 
01302 //----------------------------------------------------------------------------
01303 size_t fei::VectorSpace::getNumIDTypes()
01304 {
01305   return(idTypes_.size());
01306 }
01307 
01308 //----------------------------------------------------------------------------
01309 void fei::VectorSpace::getIDTypes(std::vector<int>& idTypes) const
01310 {
01311   size_t numIDTypes = idTypes_.size();
01312   idTypes.resize(numIDTypes);
01313   for(size_t i=0; i<numIDTypes; ++i) idTypes[i] = idTypes_[i];
01314 }
01315 
01316 //----------------------------------------------------------------------------
01317 int fei::VectorSpace::getNumFields(int idType,
01318                                    int ID)
01319 {
01320   int idindex = fei::binarySearch(idType, idTypes_);
01321   if (idindex < 0) return(0);
01322 
01323   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01324   if (record == NULL) {
01325     ERReturn(-1);
01326   }
01327 
01328   return( record->getFieldMask()->getNumFields() );
01329 }
01330 
01331 //----------------------------------------------------------------------------
01332 bool fei::VectorSpace::isLocal(int idType, int ID)
01333 {
01334   int idindex = fei::binarySearch(idType, idTypes_);
01335   if (idindex < 0) return(false);
01336 
01337   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01338   if (record == NULL) {
01339     return(false);
01340   }
01341 
01342   return(true);
01343 }
01344 
01345 //----------------------------------------------------------------------------
01346 bool fei::VectorSpace::isLocallyOwned(int idType, int ID)
01347 {
01348   int idindex = fei::binarySearch(idType, idTypes_);
01349   if (idindex < 0) return(false);
01350 
01351   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01352   if (record == NULL) {
01353     return(false);
01354   }
01355   if (record->getOwnerProc() == fei::localProc(comm_)) {
01356     return(true);
01357   }
01358 
01359   return(false);
01360 }
01361 
01362 //----------------------------------------------------------------------------
01363 unsigned fei::VectorSpace::getFieldSize(int fieldID)
01364 {
01365   std::map<int,unsigned>::const_iterator
01366     f_iter = fieldDatabase_.find(fieldID);
01367 
01368   if (f_iter == fieldDatabase_.end()) {
01369     FEI_OSTRINGSTREAM osstr;
01370     osstr << "fei::VectorSpace";
01371     if (name_.length() > 0) {
01372       osstr << "(name: "<<name_<<")";
01373     }
01374     osstr << "::getFieldSize: fieldID " << fieldID << " not found.";
01375     throw std::runtime_error(osstr.str());
01376   }
01377 
01378   return((*f_iter).second);
01379 }
01380 
01381 //----------------------------------------------------------------------------
01382 void fei::VectorSpace::getFields(int idType, int ID, std::vector<int>& fieldIDs)
01383 {
01384   int idindex = fei::binarySearch(idType, idTypes_);
01385   if (idindex < 0) {
01386     fieldIDs.resize(0);
01387     return;
01388   }
01389 
01390   fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
01391   if (record == NULL) {
01392     voidERReturn;
01393   }
01394 
01395   fei::FieldMask* fieldMask = record->getFieldMask();
01396   std::vector<int>& maskFieldIDs = fieldMask->getFieldIDs();
01397   int numFields = maskFieldIDs.size();
01398   fieldIDs.resize(numFields);
01399   for(int i=0; i<numFields; ++i) {
01400     fieldIDs[i] = maskFieldIDs[i];
01401   }
01402 }
01403 
01404 //----------------------------------------------------------------------------
01405 void fei::VectorSpace::getGlobalIndexOffsets(std::vector<int>& globalOffsets) const
01406 {
01407   globalOffsets = globalOffsets_;
01408 }
01409 
01410 //----------------------------------------------------------------------------
01411 void fei::VectorSpace::getGlobalBlkIndexOffsets(std::vector<int>& globalBlkOffsets) const
01412 {
01413   globalBlkOffsets = globalIDOffsets_;
01414 }
01415 
01416 //----------------------------------------------------------------------------
01417 int fei::VectorSpace::getOwnerProcPtIndex(int globalIndex)
01418 {
01419   if (globalIndex < 0) return(-1);
01420 
01421   unsigned len = globalOffsets_.size();
01422   for(int i=0; i<(int)(len-1); ++i) {
01423     if (globalIndex < globalOffsets_[i+1]) {
01424       return(i);
01425     }
01426   }
01427 
01428   return(-1);
01429 }
01430 
01431 //----------------------------------------------------------------------------
01432 int fei::VectorSpace::getOwnerProcBlkIndex(int globalIndex)
01433 {
01434   if (globalIndex < 0) return(-1);
01435 
01436   unsigned len = globalOffsets_.size();
01437   for(int i=0; i<(int)(len-1); ++i) {
01438     if (globalIndex < globalIDOffsets_[i+1]) {
01439       return(i);
01440     }
01441   }
01442 
01443   return(-1);
01444 }
01445 
01446 //----------------------------------------------------------------------------
01447 int fei::VectorSpace::getNumOwnedAndSharedIDs(int idType)
01448 {
01449   int idx = fei::binarySearch(idType, idTypes_);
01450   if (idx < 0) return(0);
01451 
01452   return( recordCollections_[idx]->getNumRecords() );
01453 }
01454 
01455 //----------------------------------------------------------------------------
01456 int fei::VectorSpace::getNumOwnedIDs(int idType)
01457 {
01458   int idx = fei::binarySearch(idType, idTypes_);
01459   if (idx < 0) return(0);
01460 
01461   fei::RecordAttributeCounter attrCounter(fei::localProc(comm_));
01462   runRecords(attrCounter);
01463 
01464   return( attrCounter.numLocallyOwnedIDs_ );
01465 }
01466 
01467 //----------------------------------------------------------------------------
01468 int fei::VectorSpace::getNumSharedIDs(int idType, int& numShared)
01469 {
01470   numShared = sharedIDTables_[idType].getSharedIDs().size();
01471 
01472   return(0);
01473 }
01474 
01475 //----------------------------------------------------------------------------
01476 int fei::VectorSpace::getOwnedAndSharedIDs(int idType,
01477                                       int lenList,
01478                                       int* IDs,
01479                                       int& numLocalIDs)
01480 {
01481   int idx = fei::binarySearch(idType, idTypes_);
01482   if (idx < 0) return(-1);
01483 
01484   snl_fei::RecordCollection* records = recordCollections_[idx];
01485   std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
01486   std::map<int,int>::iterator
01487     it = global_to_local.begin(),
01488     end = global_to_local.end();
01489   numLocalIDs = records->getNumRecords();
01490   int len = numLocalIDs;
01491   if (lenList < len) len = lenList;
01492   int i=0;
01493   for(; it!=end; ++it) {
01494     int local_id = it->second;
01495     IDs[i++] = records->getRecordWithLocalID(local_id)->getID();
01496     if (i >= len) break;
01497   }
01498 
01499   return(0);
01500 }
01501 
01502 //----------------------------------------------------------------------------
01503 int fei::VectorSpace::getOwnedIDs(int idType,
01504                                          int lenList,
01505                                          int* IDs,
01506                                          int& numLocalIDs)
01507 {
01508   int idx = fei::binarySearch(idType, idTypes_);
01509   if (idx < 0) return(-1);
01510 
01511   snl_fei::RecordCollection* records = recordCollections_[idx];
01512 
01513   std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
01514   std::map<int,int>::iterator
01515     it = global_to_local.begin(),
01516     end = global_to_local.end();
01517 
01518   numLocalIDs = 0;
01519 
01520   for(; it!=end; ++it) {
01521     fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
01522 
01523     if (thisrecord.getOwnerProc() == fei::localProc(comm_)) {
01524       if (numLocalIDs < lenList) {
01525         IDs[numLocalIDs] = thisrecord.getID();
01526       }
01527       ++numLocalIDs;
01528     }
01529   }
01530 
01531   return(0);
01532 }
01533 
01534 //----------------------------------------------------------------------------
01535 int fei::VectorSpace::getNumIndices_SharedAndOwned() const
01536 {
01537   return(eqnNumbers_.size());
01538 }
01539 
01540 //----------------------------------------------------------------------------
01541 int fei::VectorSpace::getIndices_SharedAndOwned(std::vector<int>& globalIndices) const
01542 {
01543   if (eqnNumbers_.size() == 0) {
01544     globalIndices.resize(0);
01545     return(0);
01546   }
01547 
01548   size_t numIndices = eqnNumbers_.size();
01549   const int* indicesPtr = &eqnNumbers_[0];
01550 
01551   globalIndices.resize(numIndices);
01552   for(size_t i=0; i<numIndices; ++i) {
01553     globalIndices[i] = indicesPtr[i];
01554   }
01555 
01556   return(0);
01557 }
01558 
01559 //----------------------------------------------------------------------------
01560 int fei::VectorSpace::getNumBlkIndices_SharedAndOwned(int& numBlkIndices) const
01561 {
01562   numBlkIndices = 0;
01563   for(size_t i=0; i<recordCollections_.size(); ++i) {
01564     numBlkIndices += recordCollections_[i]->getNumRecords();
01565   }
01566 
01567   return(0);
01568 }
01569 
01570 //----------------------------------------------------------------------------
01571 int fei::VectorSpace::getBlkIndices_SharedAndOwned(int lenBlkIndices,
01572                                                    int* globalBlkIndices,
01573                                                    int* blkSizes,
01574                                                    int& numBlkIndices)
01575 {
01576   if (!sharedRecordsSynchronized_) {
01577     numBlkIndices = 0;
01578     return(-1);
01579   }
01580 
01581   fei::BlkIndexAccessor blkIndAccessor(lenBlkIndices,
01582                                   globalBlkIndices, blkSizes);
01583   runRecords(blkIndAccessor);
01584 
01585   numBlkIndices = blkIndAccessor.numBlkIndices_;
01586 
01587   return(0);
01588 }
01589 
01590 //----------------------------------------------------------------------------
01591 int fei::VectorSpace::getGlobalNumIndices() const
01592 {
01593   if (globalOffsets_.size() < 1) return(0);
01594   return(globalOffsets_[globalOffsets_.size()-1]);
01595 }
01596 
01597 //----------------------------------------------------------------------------
01598 int fei::VectorSpace::getNumIndices_Owned() const
01599 {
01600   if (!sharedRecordsSynchronized_) {
01601     return(-1);
01602   }
01603 
01604   int localProc = fei::localProc(comm_);
01605   int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01606 
01607   return(numIndices);
01608 }
01609 
01610 //----------------------------------------------------------------------------
01611 int fei::VectorSpace::getIndices_Owned(std::vector<int>& globalIndices) const
01612 {
01613   if (!sharedRecordsSynchronized_) {
01614     globalIndices.clear();
01615     return(-1);
01616   }
01617 
01618   int localProc = fei::localProc(comm_);
01619   int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01620 
01621   globalIndices.resize(numIndices);
01622 
01623   int firstOffset = globalOffsets_[localProc];
01624   for(int i=0; i<numIndices; ++i) {
01625     globalIndices[i] = firstOffset+i;
01626   }
01627 
01628   return(0);
01629 }
01630 
01631 //----------------------------------------------------------------------------
01632 int fei::VectorSpace::getIndices_Owned(int lenIndices, int* globalIndices, int& numIndices) const
01633 {
01634   if (!sharedRecordsSynchronized_) {
01635     numIndices = 0;
01636     return(-1);
01637   }
01638 
01639   int localProc = fei::localProc(comm_);
01640   numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
01641 
01642   int len = lenIndices >= numIndices ? numIndices : lenIndices;
01643 
01644   int firstOffset = globalOffsets_[localProc];
01645   for(int i=0; i<len; ++i) {
01646     globalIndices[i] = firstOffset+i;
01647   }
01648 
01649   return(0);
01650 }
01651 
01652 //----------------------------------------------------------------------------
01653 int fei::VectorSpace::getNumBlkIndices_Owned() const
01654 {
01655   if (!sharedRecordsSynchronized_) {
01656     return(-1);
01657   }
01658 
01659   int localProc = fei::localProc(comm_);
01660   int numBlkIndices = globalIDOffsets_[localProc+1]-globalIDOffsets_[localProc];
01661 
01662   return(numBlkIndices);
01663 }
01664 
01665 //----------------------------------------------------------------------------
01666 int fei::VectorSpace::getGlobalNumBlkIndices() const
01667 {
01668   int numBlkIndices = 0;
01669   unsigned len = globalIDOffsets_.size();
01670   if (len > 0) {
01671     numBlkIndices = globalIDOffsets_[len-1];
01672   }
01673 
01674   return(numBlkIndices);
01675 }
01676 
01677 //----------------------------------------------------------------------------
01678 int fei::VectorSpace::getBlkIndices_Owned(int lenBlkIndices,
01679                                           int* globalBlkIndices,
01680                                           int* blkSizes,
01681                                           int& numBlkIndices)
01682 {
01683   if (!sharedRecordsSynchronized_) {
01684     numBlkIndices = 0;
01685     return(-1);
01686   }
01687 
01688   int localProc = fei::localProc(comm_);
01689   fei::BlkIndexAccessor blkIndAccessor(localProc, lenBlkIndices,
01690                                            globalBlkIndices, blkSizes);
01691   runRecords(blkIndAccessor);
01692 
01693   numBlkIndices = blkIndAccessor.numBlkIndices_;
01694 
01695   return(0);
01696 }
01697 
01698 //----------------------------------------------------------------------------
01699 int fei::VectorSpace::getRecordCollection(int idType,
01700                                           snl_fei::RecordCollection*& records)
01701 {
01702   int idx = fei::binarySearch(idType, idTypes_);
01703   if (idx < 0) return(-1);
01704 
01705   records = recordCollections_[idx];
01706   return(0);
01707 }
01708 
01709 //----------------------------------------------------------------------------
01710 int fei::VectorSpace::getRecordCollection(int idType,
01711                                           const snl_fei::RecordCollection*& records) const
01712 {
01713   int idx = fei::binarySearch(idType, idTypes_);
01714   if (idx < 0) return(-1);
01715 
01716   records = recordCollections_[idx];
01717   return(0);
01718 }
01719 
01720 //----------------------------------------------------------------------------
01721 void fei::VectorSpace::setOwners_shared()
01722 {
01723   //first, add localProc to each of the sharing-proc lists, in case it wasn't
01724   //included via initSharedIDs().
01725   int localProc = fei::localProc(comm_);
01726 
01727   std::map<int, fei::SharedIDs<int> >::iterator
01728     s_iter = sharedIDTables_.begin(), s_end = sharedIDTables_.end();
01729 
01730   for(; s_iter != s_end; ++s_iter) {
01731     fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
01732 
01733     fei::SharedIDs<int>::map_type::iterator
01734       t_iter = shid_table.begin(),
01735       t_end = shid_table.end();
01736 
01737     for(; t_iter != t_end; ++t_iter) {
01738       std::set<int>& shProcs = t_iter->second;
01739       shProcs.insert(localProc);
01740     }
01741   }
01742 
01743   //now set the owningProcs for the SharedIDs records, and the owning procs on
01744   //the appropriate records in the recordCollections. Set the owner to be the
01745   //lowest-numbered sharing proc in all cases.
01746   s_iter = sharedIDTables_.begin();
01747   for(; s_iter != s_end; ++s_iter) {
01748     fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
01749 
01750     std::vector<int>& owningProcs = s_iter->second.getOwningProcs();
01751 
01752     int len = shid_table.size();
01753     owningProcs.resize(len);
01754     int j = 0;
01755 
01756     fei::SharedIDs<int>::map_type::iterator
01757       t_iter = shid_table.begin(),
01758       t_end = shid_table.end();
01759 
01760     for(; t_iter != t_end; ++t_iter, ++j) {
01761       std::set<int>& shProcs = (*t_iter).second;
01762       int lowest = *(shProcs.begin());
01763       owningProcs[j] = lowest;
01764     }
01765     
01766     int idx = fei::binarySearch(s_iter->first, idTypes_);
01767     if (idx < 0) {
01768       throw std::runtime_error("internal error in fei::VectorSapce::setOwners_lowestSharing");
01769     }
01770 
01771     if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01772       recordCollections_[idx]->setDebugOutput(output_stream_);
01773     }
01774 
01775     recordCollections_[idx]->setOwners_lowestSharing(s_iter->second);
01776   }
01777 }
01778 
01779 //----------------------------------------------------------------------------
01780 int fei::VectorSpace::calculateGlobalIndices()
01781 {
01782   int localProc = fei::localProc(comm_);
01783   int numProcs = fei::numProcs(comm_);
01784   std::vector<int> localOffsets(numProcs +1, 0);
01785   std::vector<int> localIDOffsets(numProcs+1, 0);
01786   globalOffsets_.resize(numProcs + 1, 0);
01787 
01788   globalIDOffsets_.assign(globalIDOffsets_.size(), 0);
01789 
01790   //first we'll calculate the number of local degrees of freedom, and the
01791   //number of local identifiers.
01792   fei::RecordAttributeCounter counter(localProc);
01793   runRecords(counter);
01794 
01795   int numLocalDOF = counter.numLocalDOF_;
01796   int numLocalIDs = counter.numLocalIDs_;
01797   int numRemoteSharedDOF = counter.numRemoteSharedDOF_;
01798 
01799   eqnNumbers_.resize(numLocalDOF+numRemoteSharedDOF);
01800 
01801   localOffsets[localProc] = numLocalDOF;
01802   CHK_MPI( fei::GlobalMax(comm_, localOffsets, globalOffsets_) );
01803 
01804   localIDOffsets[localProc] = numLocalIDs;
01805   CHK_MPI( fei::GlobalMax(comm_, localIDOffsets, globalIDOffsets_) );
01806 
01807   //Now globalOffsets_ contains numLocalDOF for proc i, in the i-th position.
01808   //(and similarly for globalIDOffsets_)
01809   //So all we need to do is turn that data into global-offsets (i.e., the
01810   //starting global-offset for each processor).
01811   int localOffset = 0;
01812   int localIDOffset = 0;
01813   for(int p=0; p<numProcs; ++p) {
01814     numLocalDOF = globalOffsets_[p];
01815     globalOffsets_[p] = localOffset;
01816     localOffset += numLocalDOF;
01817     numLocalIDs = globalIDOffsets_[p];
01818     globalIDOffsets_[p] = localIDOffset;
01819     localIDOffset += numLocalIDs;
01820   }
01821   globalOffsets_[numProcs] = localOffset;
01822   globalIDOffsets_[numProcs] = localIDOffset;
01823 
01824   firstLocalOffset_ = globalOffsets_[localProc];
01825   lastLocalOffset_ = globalOffsets_[localProc+1] - 1;
01826 
01827   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01828     FEI_OSTREAM& os = *output_stream_;
01829     os <<dbgprefix_<< "  firstLocalOffset_: " << firstLocalOffset_ << ", lastLocalOffset_: "
01830       << lastLocalOffset_ << FEI_ENDL;
01831   }
01832 
01833   //Now we're ready to set the equation-numbers on all local records.
01834   CHK_ERR( setLocalEqnNumbers() );
01835 
01836   return(0);
01837 }
01838 
01839 //----------------------------------------------------------------------------
01840 int fei::VectorSpace::synchronizeSharedRecords()
01841 {
01842   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01843     FEI_OSTREAM& os = *output_stream_;
01844     os<<dbgprefix_<<"#synchronizeSharedRecords num-field-masks: "<<fieldMasks_.size()<<FEI_ENDL;
01845     for(unsigned fm=0; fm<fieldMasks_.size(); ++fm) {
01846       os << dbgprefix_<<"#     maskID["<<fm<<"]: " << fieldMasks_[fm]->getMaskID() << FEI_ENDL;
01847     }
01848   }
01849 
01850   bool safetyCheck = checkSharedIDs_;
01851 
01852   int numProcs = fei::numProcs(comm_);
01853   int localProc = fei::localProc(comm_);
01854   int numShTables = sharedIDTables_.size();
01855 
01856   if (numProcs < 2) {
01857     sharedRecordsSynchronized_ = true;
01858     return(0);
01859   }
01860 
01861   if (safetyCheck) {
01862     if (localProc == 0 && output_level_ >= fei::BRIEF_LOGS) {
01863       FEI_COUT << "fei::VectorSpace: global consistency-check of shared ID"
01864            << " data (involves all-to-all communication). This is done "
01865            << "only if requested by parameter 'FEI_CHECK_SHARED_IDS true'."
01866            << FEI_ENDL;
01867     }
01868 
01869     int totalNumShTables = 0;
01870     CHK_ERR( fei::GlobalSum(comm_, numShTables, totalNumShTables) );
01871     if (totalNumShTables != numShTables * numProcs) {
01872       //not all processors have the same number of shared-id tables. This means
01873       //that one or more processors is 'empty', which may not be an error. But
01874       //it does mean that the safety check can't be performed because the safety
01875       //check involves all-to-all communication and if one or more processors
01876       //don't participate then we'll hang...
01877       safetyCheck = false;
01878     }
01879   }
01880 
01881   //Create a list of fei::comm_maps which will be the communication-patterns
01882   //for each of the sharedIDTables. A sharedIDTable maps lists of processors
01883   //to each shared ID. The communication-pattern will be a transpose of
01884   //that, mapping lists of IDs to owning or sharing processors.
01885   //
01886 
01887   int local_err = 0;
01888 
01889   for(size_t i=0; i<idTypes_.size(); ++i) {
01890 
01891     fei::SharedIDs<int>& shared = getSharedIDs(idTypes_[i]);
01892     
01893     //now create/initialize ownerPatterns which map owning processors to lists
01894     //of ids that are shared locally, and sharerPatterns which map sharing
01895     //processors to lists of ids that are owned locally.
01896     fei::comm_map* ownerPattern = new fei::comm_map(1, numProcs);
01897 
01898     std::map<int,fei::comm_map*>::iterator iter = ownerPatterns_.find(idTypes_[i]);
01899     if (iter == ownerPatterns_.end()) {
01900       ownerPatterns_.insert(std::make_pair(idTypes_[i], ownerPattern));
01901     }
01902     else {
01903       delete iter->second;
01904       iter->second = ownerPattern;
01905     }
01906 
01907     fei::comm_map* sharerPattern = new fei::comm_map(1, numProcs);
01908 
01909     iter = sharerPatterns_.find(idTypes_[i]);
01910     if (iter == sharerPatterns_.end()) {
01911       sharerPatterns_.insert(std::make_pair(idTypes_[i], sharerPattern));
01912     }
01913     else {
01914       delete iter->second;
01915       iter->second = sharerPattern;
01916     }
01917 
01918     std::vector<int>& owningProcs = shared.getOwningProcs();
01919 
01920     fei::SharedIDs<int>::map_type& shtable = shared.getSharedIDs();
01921 
01922     fei::SharedIDs<int>::map_type::iterator
01923       s_iter = shtable.begin(),
01924       s_end = shtable.end();
01925     
01926     int j = 0;
01927     for(; s_iter != s_end; ++s_iter, ++j) {
01928       int ID = (*s_iter).first;
01929       std::set<int>& shProcs = (*s_iter).second;
01930 
01931       int owner = owningProcs[j];
01932 
01933       if (owner == localProc) {
01934         std::set<int>::const_iterator
01935           p_iter = shProcs.begin(),
01936           p_end = shProcs.end();
01937         for(; p_iter != p_end; ++p_iter) {
01938           if (*p_iter != localProc) {
01939             sharerPattern->addIndices(*p_iter, 1, &ID);
01940           }
01941         }
01942       }
01943       else {
01944         ownerPattern->addIndices(owner, 1, &ID);
01945       }
01946     }
01947 
01948     if (safetyCheck) {
01949       fei::comm_map* checkPattern = NULL;
01950       CHK_ERR( fei::mirrorCommPattern(comm_, ownerPattern, checkPattern) );
01951       fei::Barrier(comm_);
01952       if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01953         FEI_OSTREAM& os = *output_stream_;
01954         fei::comm_map::map_type& owner_map = ownerPattern->getMap();
01955         int numKeys = owner_map.size();
01956         fei::comm_map::map_type::const_iterator
01957           omap_iter = owner_map.begin(),
01958           omap_end = owner_map.end();
01959 
01960         os << dbgprefix_<<"#  synchronizeSharedRecords" << FEI_ENDL
01961            << dbgprefix_<<"#  ownerPattern, num-procs-to-send-to: " << numKeys << FEI_ENDL;
01962         for(int sk=0; omap_iter != omap_end; ++sk, ++omap_iter) {
01963           os << dbgprefix_<<"#    sendProc["<<sk<<"]: " << omap_iter->first << " IDs: ";
01964           fei::comm_map::row_type::const_iterator
01965             val_iter = omap_iter->second->begin(),
01966             val_end =  omap_iter->second->end();
01967           for(; val_iter != val_end; ++val_iter) {
01968             os << *val_iter << " ";
01969           }
01970           os << FEI_ENDL;
01971         }
01972 
01973         fei::comm_map::map_type& check_map = checkPattern->getMap();
01974         int numCheckKeys = check_map.size();
01975         fei::comm_map::map_type::const_iterator
01976           cmap_iter = check_map.begin(),
01977           cmap_end = check_map.end();
01978 
01979         os <<dbgprefix_<< "#  synchronizeSharedRecords" << FEI_ENDL
01980            <<dbgprefix_<< "#  checkPattern (send mirror), num-procs: "
01981            << numCheckKeys << FEI_ENDL;
01982         for(int sk=0; cmap_iter != cmap_end; ++sk, ++cmap_iter) {
01983           os <<dbgprefix_<< "#    proc["<<sk<<"]: " << cmap_iter->first << " IDs: ";
01984           fei::comm_map::row_type::const_iterator
01985           val_iter = cmap_iter->second->begin(),
01986             val_end = cmap_iter->second->end();
01987           for(; val_iter != val_end; ++val_iter) {
01988             os << *val_iter << " ";
01989           }
01990           os << FEI_ENDL;
01991         }
01992       }
01993 
01994       int err = 0;
01995       bool quiet = false;
01996       if (!checkPattern->equal(*sharerPattern, quiet)) {
01997         //fei::console_out() << "shared-ID safety-check failed..."
01998         //     << FEI_ENDL;
01999         err = -1;
02000       }
02001       int globalErr = 0;
02002       CHK_ERR( fei::GlobalSum(comm_, err, globalErr) );
02003 
02004       delete checkPattern;
02005 
02006       if (globalErr != 0) {
02007         return(globalErr);
02008       }
02009     }
02010 
02011     local_err += exchangeFieldInfo(ownerPattern, sharerPattern,
02012                                    recordCollections_[i], fieldMasks_);
02013   }
02014 
02015   int global_err = 0;
02016   CHK_ERR( fei::GlobalSum(comm_, local_err, global_err) );
02017 
02018   if (global_err != 0) {
02019     ERReturn(-1);
02020   }
02021 
02022   sharedRecordsSynchronized_ = true;
02023 
02024   return(0);
02025 }
02026 
02027 //----------------------------------------------------------------------------
02028 int fei::VectorSpace::exchangeGlobalIndices()
02029 {
02030   for(size_t i=0; i<idTypes_.size(); ++i) {
02031     snl_fei::RecordMsgHandler recmsghndlr(fei::localProc(comm_),
02032                                  recordCollections_[i], *ptBlkMap_,
02033                                           fieldMasks_, eqnNumbers_);
02034     recmsghndlr.setTask(snl_fei::RecordMsgHandler::_EqnNumbers_);
02035 
02036     recmsghndlr.setSendPattern(sharerPatterns_[idTypes_[i]]);
02037     recmsghndlr.setRecvPattern(ownerPatterns_[idTypes_[i]]);
02038     CHK_ERR( fei::exchange(comm_, &recmsghndlr) );
02039   }
02040 
02041   return(0);
02042 }
02043 
02044 //----------------------------------------------------------------------------
02045 void fei::VectorSpace::runRecords(fei::Record_Operator<int>& record_op)
02046 {
02047   for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
02048     snl_fei::RecordCollection* records = recordCollections_[rec];
02049     std::map<int,int>& g2l = records->getGlobalToLocalMap();
02050     std::map<int,int>::iterator
02051       it = g2l.begin(),
02052       end= g2l.end();
02053 
02054     for(; it!=end; ++it) {
02055       fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
02056 
02057       record_op(thisrecord);
02058     }
02059   }
02060 }
02061 
02062 //----------------------------------------------------------------------------
02063 int fei::VectorSpace::setLocalEqnNumbers()
02064 {
02065   int proc = fei::localProc(comm_);
02066   int eqnNumber = firstLocalOffset_;
02067   int idNumber = globalIDOffsets_[proc];
02068 
02069   int numProcs = fei::numProcs(comm_);
02070   int localProc = fei::localProc(comm_);
02071 
02072   if (ptBlkMap_ != NULL) {
02073     delete ptBlkMap_;
02074   }
02075   ptBlkMap_ = new snl_fei::PointBlockMap;
02076 
02077   int maxNumIndices = 0;
02078   for(unsigned i=0; i<fieldMasks_.size(); ++i) {
02079     if (fieldMasks_[i]->getNumIndices() > maxNumIndices) {
02080       maxNumIndices = fieldMasks_[i]->getNumIndices();
02081     }
02082   }
02083 
02084   if (maxNumIndices == 1) {
02085     ptBlkMap_->setPtEqualBlk();
02086   }
02087 
02088   FEI_OSTREAM* id2eqnStream = NULL;
02089   if (output_level_ >= fei::BRIEF_LOGS) {
02090     std::string path = fei::LogManager::getLogManager().getOutputPath();
02091     if (path == "") path = ".";
02092     FEI_OSTRINGSTREAM osstr;
02093     osstr << path << "/fei_ID2Eqn";
02094     if (name_.size() > 0) osstr << "_" << name_;
02095     osstr << "." <<numProcs<< "." << localProc;
02096 
02097     id2eqnStream = new FEI_OFSTREAM(osstr.str().c_str(), IOS_OUT);
02098     FEI_OSTREAM& os = *id2eqnStream;
02099     os << "# Each line contains:\n#   ID   blk-eqn   eqn" << FEI_ENDL;
02100   }
02101 
02102   int eqnNumberOffset = 0;
02103   int maxNumDOF = 0;
02104 
02105   for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
02106     snl_fei::RecordCollection* records = recordCollections_[rec];
02107 
02108     const std::map<int,int>& rmap = records->getGlobalToLocalMap();
02109     std::map<int,int>::const_iterator
02110       it = rmap.begin(), it_end = rmap.end();
02111 
02112     int* eqnNumPtr = eqnNumbers_.empty() ? NULL : &eqnNumbers_[0];
02113 
02114     for(; it!=it_end; ++it) {
02115       fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
02116 
02117       fei::FieldMask* mask = thisrecord.getFieldMask();
02118 
02119       thisrecord.setOffsetIntoEqnNumbers(eqnNumberOffset);
02120 
02121       int owner = thisrecord.getOwnerProc();
02122       if (owner == proc) {
02123         thisrecord.setNumber(idNumber++);
02124       }
02125 
02126       int* eqnNumbers = eqnNumPtr
02127                       + thisrecord.getOffsetIntoEqnNumbers();
02128 
02129       int numDOF = mask->getNumIndices();
02130       eqnNumberOffset += numDOF;
02131 
02132       if (output_level_ >= fei::BRIEF_LOGS) {
02133         for(int ii=0; ii<numDOF; ++ii) {
02134           if (isLogEqn(eqnNumber+ii) && output_stream_ != NULL) {
02135             FEI_OSTREAM& os = *output_stream_;
02136             os <<dbgprefix_<< "setLocalEqnNumbers: ID "<<thisrecord.getID()
02137                << " <--> eqn " << eqnNumber+ii<<FEI_ENDL;
02138           }
02139         }
02140       }
02141 
02142       if (owner == proc) {
02143         int offset = 0;
02144         for(int n=0; n<numDOF; ++n) {
02145           eqnNumbers[offset++] = eqnNumber++;
02146         }
02147       }
02148 
02149       if (numDOF > maxNumDOF) maxNumDOF = numDOF;
02150 
02151       if (owner == proc) {
02152         int thiseqn = eqnNumber-numDOF;
02153         int thisrecordnumber = thisrecord.getNumber();
02154         if (maxNumIndices > 1) {
02155           CHK_ERR( ptBlkMap_->setEqn(thiseqn, thisrecordnumber, numDOF) );
02156           if (numDOF > 1) {
02157             for(int i=1; i<numDOF; ++i) {
02158               CHK_ERR( ptBlkMap_->setEqn(thiseqn+i, thisrecordnumber, numDOF) );
02159             }
02160           }
02161         }
02162       }
02163 
02164       if (id2eqnStream != NULL) {
02165         if (owner == proc) {
02166           FEI_OSTREAM& os = *id2eqnStream;
02167           for(int n=0; n<numDOF; ++n) {
02168             os << thisrecord.getID() << " " << thisrecord.getNumber() << " "
02169                << eqnNumber-numDOF+n<<FEI_ENDL;
02170           }
02171         }
02172       }
02173     }
02174 
02175   }
02176 
02177   ptBlkMap_->setMaxBlkEqnSize(maxNumDOF);
02178 
02179   int globalMaxNumDOF;
02180   CHK_ERR( fei::GlobalMax(comm_, maxNumDOF, globalMaxNumDOF) );
02181 
02182   if (globalMaxNumDOF == 1) {
02183     ptBlkMap_->setPtEqualBlk();
02184   }
02185 
02186   delete id2eqnStream;
02187 
02188   return(0);
02189 }
02190 
02191 //----------------------------------------------------------------------------
02192 int fei::VectorSpace::exchangeFieldInfo(fei::comm_map* ownerPattern,
02193                                         fei::comm_map* sharerPattern,
02194                                   snl_fei::RecordCollection* recordCollection,
02195                                       std::vector<fei::FieldMask*>& fieldMasks)
02196 {
02197   //ownerPattern maps owning processors to lists of IDs that we (the local
02198   //processor) share.
02199   //
02200   //sharerPattern maps sharing processors to lists of IDs that we own.
02201   //
02202   //In this function we need to perform these tasks:
02203   //1. processors exchange and combine their sets of field-masks so that all
02204   //   processors will have the super-set of field-masks that they need.
02205   //2. sharing processors send maskIDs for their shared IDs to the owning
02206   //   processors. The owning processors then combine masks as necessary to
02207   //   make sure each shared ID has the union of field-masks that are held by
02208   //   each of the sharing processors. all owning processors now send their
02209   //   maskIDs out to the sharing processors. At this point all processors
02210   //   should have the same view of the masks that belong on shared IDs.
02211   //3. exchange info describing inactive DOF offsets for shared records.
02212   //
02213 
02214   int numProcs = fei::numProcs(comm_);
02215   if (numProcs < 2) return(0);
02216 
02217   if (ptBlkMap_ == 0) ptBlkMap_ = new snl_fei::PointBlockMap;
02218 
02219   snl_fei::RecordMsgHandler recMsgHandler(fei::localProc(comm_),
02220                                           recordCollection,
02221                                           *ptBlkMap_,
02222                                           fieldMasks, eqnNumbers_);
02223 
02224   //Step 1a.
02225   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
02226 
02227   recMsgHandler.setSendPattern(ownerPattern);
02228   recMsgHandler.setRecvPattern(sharerPattern);
02229   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02230 
02231   //Step 2a.
02232   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
02233 
02234   recMsgHandler.setSendPattern(ownerPattern);
02235   recMsgHandler.setRecvPattern(sharerPattern);
02236   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02237 
02238   //Step 1b.
02239   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
02240 
02241   recMsgHandler.setSendPattern(sharerPattern);
02242   recMsgHandler.setRecvPattern(ownerPattern);
02243   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02244 
02245   //Step 2b.
02246   recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
02247 
02248   recMsgHandler.setSendPattern(sharerPattern);
02249   recMsgHandler.setRecvPattern(ownerPattern);
02250   CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
02251 
02252   return(0);
02253 }
02254 
02255 //----------------------------------------------------------------------------
02256 fei::FieldDofMap<int>&
02257 fei::VectorSpace::getFieldDofMap()
02258 {
02259   return fieldDofMap_;
02260 }
02261 
02262 //----------------------------------------------------------------------------
02263 void fei::VectorSpace::setName(const char* name)
02264 {
02265   if (name == NULL) return;
02266 
02267   if (name_ == name) return;
02268 
02269   name_ = name;
02270   dbgprefix_ = "VecSpc_"+name_+": ";
02271 }
02272 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends