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