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