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