FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fei_VectorSpace.cpp
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2005 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 #include "fei_sstream.hpp"
10 #include "fei_fstream.hpp"
11 
12 #include "fei_utils.hpp"
13 
14 #include "fei_TemplateUtils.hpp"
15 #include "fei_chk_mpi.hpp"
16 #include <fei_CommUtils.hpp>
17 #include <fei_set_shared_ids.hpp>
18 #include "snl_fei_Utils.hpp"
19 #include "fei_Record.hpp"
20 #include "snl_fei_RecordCollection.hpp"
21 #include "fei_ParameterSet.hpp"
22 #include "snl_fei_RecordMsgHandler.hpp"
23 #include "fei_SharedIDs.hpp"
24 #include "fei_Pattern.hpp"
25 #include "fei_VectorSpace.hpp"
26 #include "fei_FieldMask.hpp"
27 #include "snl_fei_PointBlockMap.hpp"
28 #include "fei_LogManager.hpp"
29 
30 #undef fei_file
31 #define fei_file "fei_VectorSpace.cpp"
32 #include "fei_ErrMacros.hpp"
33 
34 namespace fei {
35  class RecordAttributeCounter : public Record_Operator<int> {
36  public:
37  RecordAttributeCounter(int proc)
38  : numLocalDOF_(0),
39  numLocalIDs_(0),
40  numLocallyOwnedIDs_(0),
41  numRemoteSharedDOF_(0),
42  proc_(proc)
43  {}
44 
45  virtual ~RecordAttributeCounter(){}
46 
47  void operator()(fei::Record<int>& record)
48  {
49  fei::FieldMask* mask = record.getFieldMask();
50  int owner = record.getOwnerProc();
51 
52  if (owner != proc_) {
53  numRemoteSharedDOF_ += mask->getNumIndices();
54  return;
55  }
56  else {
57  ++numLocallyOwnedIDs_;
58  }
59 
60  ++numLocalIDs_;
61 
62  int numDOF = mask->getNumIndices();
63 
64  numLocalDOF_ += numDOF;
65  }
66 
67  int numLocalDOF_;
68  int numLocalIDs_;
69  int numLocallyOwnedIDs_;
70  int numRemoteSharedDOF_;
71 
72  private:
73  int proc_;
74  };//class Record_Operator
75 
76  class BlkIndexAccessor : public Record_Operator<int> {
77  public:
78  BlkIndexAccessor(int localProc,
79  int lenBlkIndices,
80  int* globalBlkIndices,
81  int* blkSizes)
82  : numBlkIndices_(0),
83  proc_(localProc),
84  lenBlkIndices_(lenBlkIndices),
85  globalBlkIndices_(globalBlkIndices),
86  blkSizes_(blkSizes)
87  {
88  }
89 
90  BlkIndexAccessor(int lenBlkIndices,
91  int* globalBlkIndices,
92  int* blkSizes)
93  : numBlkIndices_(0),
94  proc_(-1),
95  lenBlkIndices_(lenBlkIndices),
96  globalBlkIndices_(globalBlkIndices),
97  blkSizes_(blkSizes)
98  {
99  }
100 
101  void operator()(fei::Record<int>& record)
102  {
103  int owner = record.getOwnerProc();
104  if (owner != proc_ && proc_ > -1) {
105  return;
106  }
107 
108  fei::FieldMask* mask = record.getFieldMask();
109  int blkSize = mask->getNumIndices();
110 
111  if (numBlkIndices_ < lenBlkIndices_) {
112  globalBlkIndices_[numBlkIndices_] = record.getNumber();
113  blkSizes_[numBlkIndices_] = blkSize;
114  }
115 
116  ++numBlkIndices_;
117  }
118 
119  int numBlkIndices_;
120 
121  private:
122  int proc_;
123  int lenBlkIndices_;
124  int* globalBlkIndices_;
125  int* blkSizes_;
126  };//class BlkIndexAccessor
127 }//namespace snl_fei
128 
129 //----------------------------------------------------------------------------
132  const char* name)
133 {
135  return(sptr);
136 }
137 
138 //----------------------------------------------------------------------------
139 fei::VectorSpace::VectorSpace(MPI_Comm comm, const char* name)
140  : fieldMasks_(),
141  comm_(comm),
142  idTypes_(),
143  fieldDatabase_(),
144  fieldDofMap_(),
145  recordCollections_(),
146  sharedIDTables_(),
147  ownerPatterns_(),
148  sharerPatterns_(),
149  sharedRecordsSynchronized_(true),
150  ptBlkMap_(NULL),
151  globalOffsets_(),
152  globalIDOffsets_(),
153  simpleProblem_(false),
154  firstLocalOffset_(-1),
155  lastLocalOffset_(-1),
156  eqnNumbers_(),
157  newInitData_(false),
158  initCompleteAlreadyCalled_(false),
159  name_(),
160  dbgprefix_("VecSpc: "),
161  checkSharedIDs_(false)
162 {
163  check_version();
164 
165 //The initializations below are redundant with the ones above in the
166 //initializer list. For some reason, when compiled for janus with optimization,
167 //these pointers are not being set to NULL by the above initializers.
168 //This was causing seg faults later when other VectorSpace methods delete these
169 //pointers if they aren't NULL. ABW 5/19/2004
170 
171  ptBlkMap_ = NULL;
172 
173  int numProcs = fei::numProcs(comm_);
174  globalOffsets_.assign(numProcs+1, 0);
175  globalIDOffsets_.assign(numProcs+1, 0);
176 
177  setName(name);
178 }
179 
180 //----------------------------------------------------------------------------
182 {
183  int i, len = fieldMasks_.size();
184  for(i=0; i<len; ++i) delete fieldMasks_[i];
185 
186  len = recordCollections_.size();
187  for(i=0; i<len; ++i) delete recordCollections_[i];
188 
189  std::map<int,fei::comm_map*>::iterator
190  o_iter = ownerPatterns_.begin(), o_end = ownerPatterns_.end();
191  for(; o_iter != o_end; ++o_iter) {
192  delete o_iter->second;
193  }
194 
195  std::map<int,fei::comm_map*>::iterator
196  s_iter = sharerPatterns_.begin(), s_end = sharerPatterns_.end();
197  for(; s_iter != s_end; ++s_iter) {
198  delete s_iter->second;
199  }
200 
201  delete ptBlkMap_;
202 }
203 
204 //----------------------------------------------------------------------------
205 MPI_Comm
207 {
208  return comm_;
209 }
210 
211 //----------------------------------------------------------------------------
213 {
214  const fei::Param* param = paramset.get("name");
215  fei::Param::ParamType ptype = param != NULL ?
216  param->getType() : fei::Param::BAD_TYPE;
217  if (ptype == fei::Param::STRING) {
218  setName(param->getStringValue().c_str());
219  }
220 
221  param = paramset.get("FEI_OUTPUT_LEVEL");
222  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
223  if (ptype == fei::Param::STRING) {
225  log_manager.setOutputLevel(param->getStringValue().c_str());
226  setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
227  }
228 
229  param = paramset.get("FEI_LOG_EQN");
230  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
231  if (ptype == fei::Param::INT) {
232  addLogEqn(param->getIntValue());
233  }
234 
235  param = paramset.get("FEI_LOG_ID");
236  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
237  if (ptype == fei::Param::INT) {
238  addLogID(param->getIntValue());
239  }
240 
241  param = paramset.get("FEI_CHECK_SHARED_IDS");
242  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
243  if (ptype != fei::Param::BAD_TYPE) {
244  if (ptype == fei::Param::BOOL) {
245  checkSharedIDs_ = param->getBoolValue();
246  }
247  else if (ptype == fei::Param::INT) {
248  checkSharedIDs_ = param->getIntValue() > 0 ? true : false;
249  }
250  else {
251  checkSharedIDs_ = true;
252  }
253  }
254  else {
255  checkSharedIDs_ = false;
256  }
257 }
258 
259 //----------------------------------------------------------------------------
261  const int* fieldIDs,
262  const int* fieldSizes,
263  const int* fieldTypes)
264 {
265  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
266  FEI_OSTREAM& os = *output_stream_;
267  os << dbgprefix_<<"defineFields ";
268  for(int j=0; j<numFields; ++j) {
269  os << "{"<<fieldIDs[j] << "," << fieldSizes[j] << "} ";
270  }
271  os << FEI_ENDL;
272  }
273 
274  for (int i=0; i<numFields; ++i) {
275  fieldDatabase_.insert(std::pair<int,unsigned>(fieldIDs[i], fieldSizes[i]));
276  if (fieldIDs[i] >= 0) {
277  if (fieldTypes != NULL) {
278  fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i], fieldTypes[i]);
279  }
280  else {
281  fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i]);
282  }
283  }
284  }
285 }
286 
287 //----------------------------------------------------------------------------
289  const int* idTypes)
290 {
291  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
292  FEI_OSTREAM& os = *output_stream_;
293  os << dbgprefix_<<"defineIDTypes {";
294  for(int j=0; j<numIDTypes; ++j) {
295  os << idTypes[j] << " ";
296  }
297  os << "}"<<FEI_ENDL;
298  }
299 
300  int localProc = fei::localProc(comm_);
301  for (int i=0; i<numIDTypes; ++i) {
302  int offset = fei::sortedListInsert(idTypes[i], idTypes_);
303 
304  if (offset >= 0) {
305  recordCollections_.insert(recordCollections_.begin()+offset, new snl_fei::RecordCollection(localProc));
306  }
307  }
308 }
309 
310 //----------------------------------------------------------------------------
311 void fei::VectorSpace::setIDMap(int idType,
312  const int* localIDs_begin, const int* localIDs_end,
313  const int* globalIDs_begin, const int* globalIDs_end)
314 {
315  int idx = fei::binarySearch(idType, idTypes_);
316  if (idx < 0) {
317  throw std::runtime_error("fei::VectorSpace::setIDMap ERROR, idType not found.");
318  }
319 
320  recordCollections_[idx]->setIDMap(localIDs_begin, localIDs_end,
321  globalIDs_begin, globalIDs_end);
322 }
323 
324 //----------------------------------------------------------------------------
326  int idType,
327  int numIDs,
328  const int* IDs)
329 {
330  if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
331  FEI_OSTREAM& os = *output_stream_;
332  os << dbgprefix_<<"addDOFs, fID=" << fieldID
333  <<", idT="<<idType
334  << " {";
335  for(int j=0; j<numIDs; ++j) {
336  os << IDs[j] << " ";
337  if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
338  }
339  os << "}"<<FEI_ENDL;
340  }
341 
342  if (numIDs <= 0) return(0);
343 
344  int idx = fei::binarySearch(idType, idTypes_);
345  if (idx < 0) {
346  ERReturn(-1);
347  }
348 
349  unsigned fieldSize = getFieldSize(fieldID);
350  recordCollections_[idx]->initRecords(fieldID, fieldSize,
351  numIDs, IDs,
352  fieldMasks_);
353  newInitData_ = true;
354  sharedRecordsSynchronized_ = false;
355 
356  return(0);
357 }
358 
359 //----------------------------------------------------------------------------
360 int fei::VectorSpace::addDOFs(int fieldID,
361  int idType,
362  int numIDs,
363  const int* IDs,
364  int* records)
365 {
366  if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
367  FEI_OSTREAM& os = *output_stream_;
368  os << dbgprefix_<<"addDOFs*, fID=" << fieldID
369  <<", idT="<<idType
370  << " {";
371  for(int j=0; j<numIDs; ++j) {
372  os << IDs[j] << " ";
373  if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
374  }
375  os << "}"<<FEI_ENDL;
376  }
377 
378  if (numIDs <= 0) return(0);
379 
380  int idx = fei::binarySearch(idType, idTypes_);
381  if (idx < 0) {
382  FEI_OSTRINGSTREAM osstr;
383  osstr << "fei::VectorSpace::addDOFs: error, idType " << idType
384  << " not recognized. (idTypes need to be initialized via the"
385  << " method VectorSpace::defineIDTypes)";
386  throw std::runtime_error(osstr.str());
387  }
388 
389  unsigned fieldSize = getFieldSize(fieldID);
390  recordCollections_[idx]->initRecords(fieldID, fieldSize,
391  numIDs, IDs,
392  fieldMasks_, records);
393  newInitData_ = true;
394  sharedRecordsSynchronized_ = false;
395 
396  return(0);
397 }
398 
399 //----------------------------------------------------------------------------
401  int numIDs,
402  const int* IDs)
403 {
404  if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
405  FEI_OSTREAM& os = *output_stream_;
406  os << dbgprefix_<<"addDOFs idT=" <<idType<<" {";
407  for(int j=0; j<numIDs; ++j) {
408  os << IDs[j] << " ";
409  if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
410  }
411  os << "}"<<FEI_ENDL;
412  }
413 
414  if (numIDs <= 0) return(0);
415 
416  int idx = fei::binarySearch(idType, idTypes_);
417  if (idx < 0) ERReturn(-1);
418 
419  recordCollections_[idx]->initRecords(numIDs, IDs, fieldMasks_);
420 
421  newInitData_ = true;
422  sharedRecordsSynchronized_ = false;
423 
424  return(0);
425 }
426 
427 //----------------------------------------------------------------------------
428 int fei::VectorSpace::addDOFs(int idType,
429  int numIDs,
430  const int* IDs,
431  int* records)
432 {
433  if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
434  FEI_OSTREAM& os = *output_stream_;
435  os << dbgprefix_<<"addDOFs* idT=" <<idType<<" {";
436  for(int j=0; j<numIDs; ++j) {
437  os << IDs[j] << " ";
438  if (j>0 && j%20==0) os << FEI_ENDL << dbgprefix_;
439  }
440  os << "}"<<FEI_ENDL;
441  }
442 
443  if (numIDs <= 0) return(0);
444 
445  int idx = fei::binarySearch(idType, idTypes_);
446  if (idx < 0) ERReturn(-1);
447 
448  recordCollections_[idx]->initRecords(numIDs, IDs,
449  fieldMasks_, records);
450 
451  newInitData_ = true;
452  sharedRecordsSynchronized_ = false;
453 
454  return(0);
455 }
456 
457 //----------------------------------------------------------------------------
459  int idType,
460  const int* sharedIDs,
461  const int* numSharingProcsPerID,
462  const int* sharingProcs)
463 {
464  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
465  FEI_OSTREAM& os = *output_stream_;
466  os << dbgprefix_<<"initSharedIDs n=" <<numShared<<", idT="<<idType<< FEI_ENDL;
467  int offset = 0;
468  for(int ns=0; ns<numShared; ++ns) {
469  os << dbgprefix_<<"#sharedID="<<sharedIDs[ns] << ", nprocs=" << numSharingProcsPerID[ns] << ", procs: ";
470  for(int sp=0; sp<numSharingProcsPerID[ns]; ++sp) {
471  os << sharingProcs[offset++] << " ";
472  }
473  os << FEI_ENDL;
474  }
475  os << FEI_ENDL;
476  }
477 
478  if (numShared == 0) return(0);
479 
480  int idx = fei::binarySearch(idType, idTypes_);
481  if (idx < 0) ERReturn(-1);
482 
483  fei::SharedIDs<int>& shIDs = getSharedIDs(idType);
484 
485  int offset = 0;
486  for(int i=0; i<numShared; ++i) {
487  shIDs.addSharedID(sharedIDs[i], numSharingProcsPerID[i],
488  &(sharingProcs[offset]));
489  offset += numSharingProcsPerID[i];
490 
491  fei::Record<int>* rec = recordCollections_[idx]->getRecordWithID(sharedIDs[i]);
492  if (rec == NULL) {
493  CHK_ERR( addDOFs(idType, 1, &(sharedIDs[i])) );
494  }
495  }
496 
497  newInitData_ = true;
498  sharedRecordsSynchronized_ = false;
499 
500  return(0);
501 }
502 
503 //----------------------------------------------------------------------------
505  int idType,
506  const int* sharedIDs,
507  const int* numSharingProcsPerID,
508  const int* const* sharingProcs)
509 {
510  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
511  FEI_OSTREAM& os = *output_stream_;
512  os << dbgprefix_<<"initSharedIDs n=" <<numShared<<", idT="<<idType<< FEI_ENDL;
513  for(int ns=0; ns<numShared; ++ns) {
514  os << dbgprefix_<<"#sharedID="<<sharedIDs[ns] << ", nprocs=" << numSharingProcsPerID[ns] << ", procs: ";
515  for(int sp=0; sp<numSharingProcsPerID[ns]; ++sp) {
516  os << sharingProcs[ns][sp] << " ";
517  }
518  os << FEI_ENDL;
519  }
520  os << FEI_ENDL;
521  }
522 
523  if (numShared == 0) return(0);
524 
525  fei::SharedIDs<int>& shIDs = getSharedIDs(idType);
526 
527  int idx = fei::binarySearch(idType, idTypes_);
528  if (idx < 0) ERReturn(-1);
529 
530  for(int i=0; i<numShared; ++i) {
531  shIDs.addSharedID(sharedIDs[i], numSharingProcsPerID[i],
532  sharingProcs[i]);
533 
534  fei::Record<int>* rec = recordCollections_[idx]->getRecordWithID(sharedIDs[i]);
535  if (rec == NULL) {
536  CHK_ERR( addDOFs(idType, 1, &(sharedIDs[i])) );
537  }
538  }
539 
540  newInitData_ = true;
541  sharedRecordsSynchronized_ = false;
542 
543  return(0);
544 }
545 
546 int fei::VectorSpace::setOwners(int numShared, int idType, const int* sharedIDs, const int* owners)
547 {
548  int idx = fei::binarySearch(idType, idTypes_);
549  if (idx < 0) ERReturn(-1);
550  for(int i=0; i<numShared; ++i) {
551  fei::Record<int>* rec = recordCollections_[idx]->getRecordWithID(sharedIDs[i]);
552  if (rec == NULL) {
553  throw std::runtime_error("shared id not found in setOwners");
554  }
555  rec->setOwnerProc(owners[i]);
556  }
557 
558  return 0;
559 }
560 
561 //----------------------------------------------------------------------------
563 {
564  idTypes_ = inputSpace->idTypes_;
565 
566  std::map<int,unsigned>::const_iterator
567  f_iter = inputSpace->fieldDatabase_.begin(),
568  f_end = inputSpace->fieldDatabase_.end();
569 
570  for(; f_iter != f_end; ++f_iter) {
571  const std::pair<const int,unsigned>& fpair = *f_iter;
572  int fieldsize = fpair.second;
573  defineFields(1, &(fpair.first), &fieldsize);
574  }
575 
576  int i, len = inputSpace->fieldMasks_.size();
577  fieldMasks_.resize(len);
578  for(i=0; i<len; ++i) {
579  fieldMasks_[i] = new fei::FieldMask(*(inputSpace->fieldMasks_[i]));
580  }
581 
582  len = inputSpace->recordCollections_.size();
583  recordCollections_.resize(len);
584  for(i=0; i<len; ++i) {
585  recordCollections_[i] =
586  new snl_fei::RecordCollection(*(inputSpace->recordCollections_[i]));
587  }
588 
589  sharedIDTables_ = inputSpace->sharedIDTables_;
590 
591  newInitData_ = true;
592  sharedRecordsSynchronized_ = false;
593 
594  return(0);
595 }
596 
597 //----------------------------------------------------------------------------
598 void fei::VectorSpace::getSendProcs(std::vector<int>& sendProcs) const
599 {
600  sendProcs.clear();
601  std::set<int> sendSet;
602 
603  std::map<int,fei::SharedIDs<int> >::const_iterator
604  s_it = sharedIDTables_.begin(),
605  s_end= sharedIDTables_.end();
606  for(; s_it!=s_end; ++s_it) {
607  const std::vector<int>& owners = s_it->second.getOwningProcs();
608  for(size_t i=0; i<owners.size(); ++i) {
609  sendSet.insert(owners[i]);
610  }
611  }
612 
613  for(std::set<int>::iterator it=sendSet.begin(); it!=sendSet.end(); ++it) {
614  if (fei::localProc(comm_) != *it) sendProcs.push_back(*it);
615  }
616 }
617 
618 //----------------------------------------------------------------------------
619 fei::SharedIDs<int>& fei::VectorSpace::getSharedIDs(int idType)
620 {
621  std::map<int,fei::SharedIDs<int> >::iterator
622  iter = sharedIDTables_.find(idType);
623 
624  if (iter == sharedIDTables_.end()) {
625  iter = sharedIDTables_.insert(std::make_pair(idType,fei::SharedIDs<int>())).first;
626  }
627 
628  return iter->second;
629 }
630 
631 //----------------------------------------------------------------------------
632 void fei::VectorSpace::compute_shared_ids(const std::vector<int>& global_min, const std::vector<int>& global_max)
633 {
634  //For each ID-type, call the fei::set_shared_ids function, which
635  //takes a RecordCollection as input and fills a SharedIDs object
636  //with mappings for which processors share each ID that appears
637  //on multiple processors.
638  for(size_t i=0; i<idTypes_.size(); ++i) {
639  snl_fei::RecordCollection* records = recordCollections_[i];
640  fei::SharedIDs<int>& shIDs = getSharedIDs(idTypes_[i]);
641 
642  fei::set_shared_ids(comm_, *records, shIDs, global_min[i], global_max[i]);
643  }
644 }
645 
646 //----------------------------------------------------------------------------
648 {
649  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
650  FEI_OSTREAM& os = *output_stream_;
651  os <<dbgprefix_<< "initComplete" << FEI_ENDL;
652  }
653 
654  simpleProblem_ = (fieldMasks_.size()==1) && (fieldMasks_[0]->getNumFields()==1);
655 
656  std::vector<int> local_data(recordCollections_.size()+1);
657  for(size_t i=0; i<recordCollections_.size(); ++i) {
658  local_data[i] = recordCollections_[i]->getMaxID();
659  }
660 
661  std::vector<int> global_max_data(recordCollections_.size()+1, 0);
662 
663  //we need to know if any processor has newInitData_.
664  int localInitData = newInitData_ ? 1 : 0;
665  local_data[local_data.size()-1] = localInitData;
666 
667  CHK_ERR( fei::GlobalMax(comm_, local_data, global_max_data) );
668  int globalInitData = global_max_data[global_max_data.size()-1];
669  newInitData_ = globalInitData > 0 ? true : false;
670 
671  if (!newInitData_) {
672  return(0);
673  }
674 
675  //if running on multiple procs and initSharedIDs(...) has not been
676  //called, then we need to specifically figure out which IDs are
677  //shared and which procs share them.
678 
679  bool need_to_compute_shared_ids = false;
680  if (fei::numProcs(comm_) > 1 && sharedIDTables_.size() < 1) {
681  need_to_compute_shared_ids = true;
682  }
683 
684  for(size_t i=0; i<recordCollections_.size(); ++i) {
685  local_data[i] = recordCollections_[i]->getMinID();
686  }
687  int localNeedToCompute = need_to_compute_shared_ids ? 1 : 0;
688  local_data[local_data.size()-1] = localNeedToCompute;
689 
690  std::vector<int> global_min_data(recordCollections_.size()+1);
691  CHK_ERR( fei::GlobalMin(comm_, local_data, global_min_data) );
692  int globalNeedToCompute = global_min_data[global_min_data.size()-1];
693  need_to_compute_shared_ids = globalNeedToCompute==1 ? true : false;
694  if (need_to_compute_shared_ids) {
695  compute_shared_ids(global_min_data, global_max_data);
696  }
697 
698  //setOwners_shared is a local operation (no communication), which
699  //assumes that each processor holds CORRECT (globally symmetric)
700  //shared-id/sharing-proc tables. No correctness-checking is performed here.
701  setOwners_shared();
702 
703  for(size_t i=0; i<recordCollections_.size(); ++i) {
704  recordCollections_[i]->setOwners_local();
705  }
706 
707  //synchronizeSharedRecords ensures that each sharing processor has the same
708  //view of the shared records, with respect to the layout of fields, which
709  //determines how many DOFs and equation-numbers reside at each ID.
710  //This involves inter-processor communication.
711  if ( synchronizeSharedRecords() != 0) {
712  return(-1);
713  }
714 
715  //calculateGlobalIndices is also a global operation.
716  if ( calculateGlobalIndices() != 0) {
717  return(-1);
718  }
719 
720  //finally we need to exchange global indices for shared records. i.e., the
721  //processors that own shared records, need to send the global indices for
722  //those records to the sharing-but-not-owning processors.
723  if (fei::numProcs(comm_) > 1) {
724  if ( exchangeGlobalIndices() != 0) {
725  return(-1);
726  }
727  }
728 
729  newInitData_ = false;
730  initCompleteAlreadyCalled_ = true;
731  return(0);
732 }
733 
734 //----------------------------------------------------------------------------
736  int ID,
737  int fieldID,
738  int fieldOffset,
739  int whichComponentOfField,
740  int& globalIndex)
741 {
742  int idindex = fei::binarySearch(idType, idTypes_);
743  if (idindex < 0) return(-1);
744 
745  unsigned fieldSize = 0;
746  if (fieldOffset > 0) {
747  fieldSize = getFieldSize(fieldID);
748  }
749 
750  globalIndex = recordCollections_[idindex]->getGlobalIndex(ID,
751  fieldID,
752  fieldSize,
753  fieldOffset,
754  whichComponentOfField,
755  &eqnNumbers_[0]);
756  int return_value = globalIndex >= 0 ? 0 : -1;
757  return return_value;
758 }
759 
760 //----------------------------------------------------------------------------
762  int ID,
763  int fieldID,
764  int& globalIndex)
765 {
766  return( getGlobalIndex(idType, ID, fieldID, 0, 0, globalIndex) );
767 }
768 
769 //----------------------------------------------------------------------------
771  int ID,
772  int& globalBlkIndex)
773 {
774  int idindex = fei::binarySearch(idType, idTypes_);
775  if (idindex < 0) return(-1);
776 
777  CHK_ERR(recordCollections_[idindex]->getGlobalBlkIndex(ID, globalBlkIndex));
778 
779  return(0);
780 }
781 
782 //----------------------------------------------------------------------------
784  const int* IDs,
785  int idType,
786  int fieldID,
787  int* globalIndices)
788 {
789  int idindex = fei::binarySearch(idType, idTypes_);
790  if (idindex < 0) return(-1);
791 
792  unsigned fieldSize = getFieldSize(fieldID);
793  int offset = 0;
794 
795  for(int i=0; i<numIDs; ++i) {
796  globalIndices[offset] = recordCollections_[idindex]->getGlobalIndex(IDs[i],
797  fieldID, fieldSize,
798  0, 0, &eqnNumbers_[0]);
799  if (fieldSize>1) {
800  if (globalIndices[offset] >= 0) {
801  int eqn = globalIndices[offset];
802  for(unsigned j=1; j<fieldSize; ++j) {
803  globalIndices[offset+j] = eqn+j;
804  }
805  }
806  else {
807  for(unsigned j=0; j<fieldSize; ++j) {
808  globalIndices[offset+j] = -1;
809  }
810  }
811  }
812 
813  offset += fieldSize;
814  }
815 
816  return(0);
817 }
818 
819 //----------------------------------------------------------------------------
821  const int* localIDs,
822  int idType,
823  int fieldID,
824  int* globalIndices)
825 {
826  int idindex = fei::binarySearch(idType, idTypes_);
827  if (idindex < 0) return(-1);
828 
829  unsigned fieldSize = getFieldSize(fieldID);
830  int offset = 0;
831 
832  for(int i=0; i<numIDs; ++i) {
833  globalIndices[offset] = recordCollections_[idindex]->getGlobalIndexLocalID(localIDs[i],
834  fieldID, fieldSize,
835  0, 0, &eqnNumbers_[0]);
836  if (globalIndices[offset] >= 0) {
837  if (fieldSize>1) {
838  int eqn = globalIndices[offset];
839  for(unsigned j=1; j<fieldSize; ++j) {
840  globalIndices[offset+j] = eqn+j;
841  }
842  }
843  }
844  else {
845  for(unsigned j=0; j<fieldSize; ++j) {
846  globalIndices[offset+j] = -1;
847  }
848  }
849 
850  offset += fieldSize;
851  }
852 
853  return(0);
854 }
855 
856 //----------------------------------------------------------------------------
858  const int* IDs,
859  int idType,
860  int* globalBlkIndices)
861 {
862  int err;
863  int idindex = fei::binarySearch(idType, idTypes_);
864  if (idindex < 0) return(-1);
865 
866  int offset = 0;
867 
868  for(int i=0; i<numIDs; ++i) {
869  err = recordCollections_[idindex]->getGlobalBlkIndex(IDs[i],
870  globalBlkIndices[offset]);
871  if (err != 0) {
872  globalBlkIndices[offset] = -1;
873  }
874  ++offset;
875  }
876 
877  return(0);
878 }
879 
880 //----------------------------------------------------------------------------
882  const int* IDs,
883  const int* idTypes,
884  const int* fieldIDs,
885  int* globalIndices)
886 {
887  int err;
888  int offset = 0;
889  for(int i=0; i<numIDs; ++i) {
890  unsigned fieldSize = getFieldSize(fieldIDs[i]);
891  err = getGlobalIndex(idTypes[i], IDs[i], fieldIDs[i], 0, 0,
892  globalIndices[offset]);
893  if (err) {
894  for(unsigned j=1; j<fieldSize; ++j) {
895  globalIndices[offset+j] = -1;
896  }
897  }
898  else {
899  if (fieldSize>1) {
900  int eqn = globalIndices[offset];
901  for(unsigned j=1; j<fieldSize; ++j) {
902  globalIndices[offset+j] = eqn+j;
903  }
904  }
905  }
906  offset += fieldSize;
907  }
908 
909  return(0);
910 }
911 
912 //----------------------------------------------------------------------------
914  const fei::Record<int>*const* records,
915  std::vector<int>& indices)
916 {
917  int numRecords = pattern->getNumIDs();
918  indices.resize(numRecords);
919  int numIndices;
920  getGlobalBlkIndices(numRecords, records, numRecords, &indices[0],
921  numIndices);
922 }
923 
924 //----------------------------------------------------------------------------
926  const fei::Record<int>*const* records,
927  std::vector<int>& indices)
928 {
929  int numRecords = pattern->getNumIDs();
930  int numIndices = pattern->getNumIndices();
931  indices.resize(numIndices);
932  int* indices_ptr = &indices[0];
933 
934  fei::Pattern::PatternType pType = pattern->getPatternType();
935 
936  if (pType == fei::Pattern::GENERAL ||
937  pType == fei::Pattern::SINGLE_IDTYPE) {
938  const int* numFieldsPerID = pattern->getNumFieldsPerID();
939  const int* fieldIDs = pattern->getFieldIDs();
940  int totalNumFields = pattern->getTotalNumFields();
941 
942  std::vector<int> fieldSizes(totalNumFields);
943 
944  for(int j=0; j<totalNumFields; ++j) {
945  fieldSizes[j] = getFieldSize(fieldIDs[j]);
946  }
947 
948  getGlobalIndices(numRecords, records, numFieldsPerID,
949  fieldIDs, &(fieldSizes[0]),
950  numIndices, indices_ptr, numIndices);
951  }
952  else if (pType == fei::Pattern::SIMPLE) {
953  const int* fieldIDs = pattern->getFieldIDs();
954 
955  int fieldID = fieldIDs[0];
956  unsigned fieldSize = getFieldSize(fieldID);
957 
958  getGlobalIndices(numRecords, records,
959  fieldID, fieldSize,
960  numIndices, indices_ptr, numIndices);
961  }
962  else if (pType == fei::Pattern::NO_FIELD) {
963  getGlobalBlkIndices(numRecords, records, numIndices, indices_ptr, numIndices);
964  }
965 }
966 
967 //----------------------------------------------------------------------------
968 void fei::VectorSpace::getGlobalIndicesL(const fei::Pattern* pattern,
969  const int* records,
970  std::vector<int>& indices)
971 {
972  int numRecords = pattern->getNumIDs();
973  int numIndices = pattern->getNumIndices();
974  const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
975  indices.resize(numIndices);
976  int* indices_ptr = &indices[0];
977 
978  fei::Pattern::PatternType pType = pattern->getPatternType();
979 
980  if (pType == fei::Pattern::GENERAL ||
981  pType == fei::Pattern::SINGLE_IDTYPE) {
982  const int* numFieldsPerID = pattern->getNumFieldsPerID();
983  const int* fieldIDs = pattern->getFieldIDs();
984  int totalNumFields = pattern->getTotalNumFields();
985 
986  std::vector<int> fieldSizes(totalNumFields);
987 
988  for(int j=0; j<totalNumFields; ++j) {
989  fieldSizes[j] = getFieldSize(fieldIDs[j]);
990  }
991 
992  getGlobalIndicesL(numRecords, recordCollections, records, numFieldsPerID,
993  fieldIDs, &(fieldSizes[0]),
994  numIndices, indices_ptr, numIndices);
995  }
996  else if (pType == fei::Pattern::SIMPLE) {
997  const int* fieldIDs = pattern->getFieldIDs();
998 
999  int fieldID = fieldIDs[0];
1000  unsigned fieldSize = getFieldSize(fieldID);
1001 
1002  getGlobalIndicesL(numRecords, recordCollections, records,
1003  fieldID, fieldSize,
1004  numIndices, indices_ptr, numIndices);
1005  }
1006  else if (pType == fei::Pattern::NO_FIELD) {
1007  getGlobalBlkIndicesL(numRecords, recordCollections, records, numIndices, indices_ptr, numIndices);
1008  }
1009 }
1010 
1011 //----------------------------------------------------------------------------
1012 void fei::VectorSpace::getGlobalIndices(int numRecords,
1013  const fei::Record<int>*const* records,
1014  int fieldID,
1015  int fieldSize,
1016  int indicesAllocLen,
1017  int* indices,
1018  int& numIndices)
1019 {
1020  numIndices = 0;
1021  int* eqnPtr = &eqnNumbers_[0];
1022 
1023  int len = numRecords;
1024  if (len*fieldSize >= indicesAllocLen) {
1025  len = indicesAllocLen/fieldSize;
1026  }
1027 
1028  if (fieldSize == 1 && simpleProblem_) {
1029  for(int i=0; i<len; ++i) {
1030  const fei::Record<int>* record = records[i];
1031  indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
1032  }
1033  return;
1034  }
1035 
1036  if (fieldSize == 1) {
1037  for(int i=0; i<len; ++i) {
1038  const fei::Record<int>* record = records[i];
1039 
1040  int eqnOffset = 0;
1041  int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
1042 
1043  const fei::FieldMask* fieldMask = record->getFieldMask();
1044  int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
1045  indices[numIndices++] = err == 0 ? eqnNumbers[eqnOffset] : -1;
1046  }
1047  }
1048  else {
1049  for(int i=0; i<len; ++i) {
1050  const fei::Record<int>* record = records[i];
1051 
1052  int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
1053 
1054  int eqnOffset = 0;
1055  if (!simpleProblem_) {
1056  const fei::FieldMask* fieldMask = record->getFieldMask();
1057  int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
1058  if (err != 0) {
1059  for(int fs=0; fs<fieldSize; ++fs) {
1060  indices[numIndices++] = -1;
1061  }
1062  continue;
1063  }
1064  }
1065 
1066  for(int fs=0; fs<fieldSize; ++fs) {
1067  indices[numIndices++] = eqnNumbers[eqnOffset+fs];
1068  }
1069  }
1070  }
1071 }
1072 
1073 //----------------------------------------------------------------------------
1074 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
1075  const snl_fei::RecordCollection*const* recordCollections,
1076  const int* records,
1077  int fieldID,
1078  int fieldSize,
1079  int indicesAllocLen,
1080  int* indices,
1081  int& numIndices)
1082 {
1083  numIndices = 0;
1084  int* eqnPtr = &eqnNumbers_[0];
1085 
1086  int len = numRecords;
1087  if (len*fieldSize >= indicesAllocLen) {
1088  len = indicesAllocLen/fieldSize;
1089  }
1090 
1091  if (fieldSize == 1 && simpleProblem_) {
1092  for(int i=0; i<len; ++i) {
1093  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
1094  indices[numIndices++] = *(eqnPtr+record->getOffsetIntoEqnNumbers());
1095  }
1096  return;
1097  }
1098 
1099  if (fieldSize == 1) {
1100  for(int i=0; i<len; ++i) {
1101  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
1102 
1103  int eqnOffset = 0;
1104  int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
1105 
1106  const fei::FieldMask* fieldMask = record->getFieldMask();
1107  int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
1108  indices[numIndices++] = err == 0 ? eqnNumbers[eqnOffset] : -1;
1109  }
1110  }
1111  else {
1112  for(int i=0; i<len; ++i) {
1113  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
1114 
1115  int* eqnNumbers = eqnPtr+record->getOffsetIntoEqnNumbers();
1116 
1117  int eqnOffset = 0;
1118  if (!simpleProblem_) {
1119  const fei::FieldMask* fieldMask = record->getFieldMask();
1120  int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
1121  if (err != 0) {
1122  for(int fs=0; fs<fieldSize; ++fs) {
1123  indices[numIndices++] = -1;
1124  }
1125  continue;
1126  }
1127  }
1128 
1129  for(int fs=0; fs<fieldSize; ++fs) {
1130  indices[numIndices++] = eqnNumbers[eqnOffset+fs];
1131  }
1132  }
1133  }
1134 }
1135 
1136 //----------------------------------------------------------------------------
1137 void fei::VectorSpace::getGlobalIndices(int numRecords,
1138  const fei::Record<int>*const* records,
1139  const int* numFieldsPerID,
1140  const int* fieldIDs,
1141  const int* fieldSizes,
1142  int indicesAllocLen,
1143  int* indices,
1144  int& numIndices)
1145 {
1146  numIndices = 0;
1147  int fld_offset = 0;
1148  int* eqnPtr = &eqnNumbers_[0];
1149 
1150  for(int i=0; i<numRecords; ++i) {
1151  const fei::Record<int>* record = records[i];
1152 
1153  const fei::FieldMask* fieldMask = record->getFieldMask();
1154  int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
1155 
1156  for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
1157  int eqnOffset = 0;
1158  if (!simpleProblem_) {
1159  int err = fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
1160  if (err != 0) {
1161  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
1162  indices[numIndices++] = -1;
1163  }
1164  continue;
1165  }
1166  }
1167 
1168  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
1169  indices[numIndices++] = eqnNumbers[eqnOffset+fs];
1170  }
1171 
1172  ++fld_offset;
1173  }
1174  }
1175 }
1176 
1177 //----------------------------------------------------------------------------
1178 void fei::VectorSpace::getGlobalIndicesL(int numRecords,
1179  const snl_fei::RecordCollection*const* recordCollections,
1180  const int* records,
1181  const int* numFieldsPerID,
1182  const int* fieldIDs,
1183  const int* fieldSizes,
1184  int indicesAllocLen,
1185  int* indices,
1186  int& numIndices)
1187 {
1188  numIndices = 0;
1189  int fld_offset = 0;
1190  int* eqnPtr = &eqnNumbers_[0];
1191 
1192  for(int i=0; i<numRecords; ++i) {
1193  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
1194 
1195  const fei::FieldMask* fieldMask = record->getFieldMask();
1196  int* eqnNumbers = eqnPtr + record->getOffsetIntoEqnNumbers();
1197 
1198  for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
1199  int eqnOffset = 0;
1200  if (!simpleProblem_) {
1201  int err = fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
1202  if (err != 0) {
1203  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
1204  indices[numIndices++] = -1;
1205  }
1206  continue;
1207  }
1208  }
1209 
1210  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
1211  indices[numIndices++] = eqnNumbers[eqnOffset+fs];
1212  }
1213 
1214  ++fld_offset;
1215  }
1216  }
1217 }
1218 
1219 //----------------------------------------------------------------------------
1220 void fei::VectorSpace::getGlobalBlkIndices(int numRecords,
1221  const fei::Record<int>*const* records,
1222  int indicesAllocLen,
1223  int* indices,
1224  int& numIndices)
1225 {
1226  numIndices = 0;
1227  for(int i=0; i<numRecords; ++i) {
1228  if (numIndices < indicesAllocLen) {
1229  indices[numIndices++] = records[i]->getNumber();
1230  }
1231  else ++numIndices;
1232  }
1233 }
1234 
1235 //----------------------------------------------------------------------------
1236 void fei::VectorSpace::getGlobalBlkIndicesL(int numRecords,
1237  const snl_fei::RecordCollection*const* recordCollections,
1238  const int* records,
1239  int indicesAllocLen,
1240  int* indices,
1241  int& numIndices)
1242 {
1243  numIndices = 0;
1244  for(int i=0; i<numRecords; ++i) {
1245  if (numIndices < indicesAllocLen) {
1246  indices[numIndices++] = recordCollections[i]->getRecordWithLocalID(records[i])->getNumber();
1247  }
1248  else ++numIndices;
1249  }
1250 }
1251 
1252 //----------------------------------------------------------------------------
1254  int ID,
1255  int& globalIndex)
1256 {
1257  int idindex = fei::binarySearch(idType, idTypes_);
1258  if (idindex < 0) return(-1);
1259 
1260  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1261  if (record == NULL) {
1262  ERReturn(-1);
1263  }
1264 
1265  const int* eqnNums = &eqnNumbers_[0]
1266  + record->getOffsetIntoEqnNumbers();
1267 
1268  if (eqnNums != NULL) { globalIndex = eqnNums[0]; return(0); }
1269  return(-1);
1270 }
1271 
1272 //----------------------------------------------------------------------------
1274  int ID)
1275 {
1276  int idindex = fei::binarySearch(idType, idTypes_);
1277  if (idindex < 0) return(0);
1278 
1279  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1280  if (record == NULL) {
1281  ERReturn(-1);
1282  }
1283 
1284  return( record->getFieldMask()->getNumIndices() );
1285 }
1286 
1287 //----------------------------------------------------------------------------
1289 {
1290  return(fieldDatabase_.size());
1291 }
1292 
1293 //----------------------------------------------------------------------------
1294 void fei::VectorSpace::getFields(std::vector<int>& fieldIDs)
1295 {
1296  unsigned numFields = fieldDatabase_.size();
1297 
1298  fieldIDs.resize(numFields);
1299 
1300  fei::copyKeysToArray<std::map<int,unsigned> >(fieldDatabase_, numFields,
1301  &fieldIDs[0]);
1302 }
1303 
1304 //----------------------------------------------------------------------------
1306 {
1307  return(idTypes_.size());
1308 }
1309 
1310 //----------------------------------------------------------------------------
1311 void fei::VectorSpace::getIDTypes(std::vector<int>& idTypes) const
1312 {
1313  size_t numIDTypes = idTypes_.size();
1314  idTypes.resize(numIDTypes);
1315  for(size_t i=0; i<numIDTypes; ++i) idTypes[i] = idTypes_[i];
1316 }
1317 
1318 //----------------------------------------------------------------------------
1320  int ID)
1321 {
1322  int idindex = fei::binarySearch(idType, idTypes_);
1323  if (idindex < 0) return(0);
1324 
1325  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1326  if (record == NULL) {
1327  ERReturn(-1);
1328  }
1329 
1330  return( record->getFieldMask()->getNumFields() );
1331 }
1332 
1333 //----------------------------------------------------------------------------
1334 bool fei::VectorSpace::isLocal(int idType, int ID)
1335 {
1336  int idindex = fei::binarySearch(idType, idTypes_);
1337  if (idindex < 0) return(false);
1338 
1339  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1340  if (record == NULL) {
1341  return(false);
1342  }
1343 
1344  return(true);
1345 }
1346 
1347 //----------------------------------------------------------------------------
1348 bool fei::VectorSpace::isLocallyOwned(int idType, int ID)
1349 {
1350  int idindex = fei::binarySearch(idType, idTypes_);
1351  if (idindex < 0) return(false);
1352 
1353  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1354  if (record == NULL) {
1355  return(false);
1356  }
1357  if (record->getOwnerProc() == fei::localProc(comm_)) {
1358  return(true);
1359  }
1360 
1361  return(false);
1362 }
1363 
1364 //----------------------------------------------------------------------------
1365 unsigned fei::VectorSpace::getFieldSize(int fieldID)
1366 {
1367  std::map<int,unsigned>::const_iterator
1368  f_iter = fieldDatabase_.find(fieldID);
1369 
1370  if (f_iter == fieldDatabase_.end()) {
1371  FEI_OSTRINGSTREAM osstr;
1372  osstr << "fei::VectorSpace";
1373  if (name_.length() > 0) {
1374  osstr << "(name: "<<name_<<")";
1375  }
1376  osstr << "::getFieldSize: fieldID " << fieldID << " not found.";
1377  throw std::runtime_error(osstr.str());
1378  }
1379 
1380  return((*f_iter).second);
1381 }
1382 
1383 //----------------------------------------------------------------------------
1384 void fei::VectorSpace::getFields(int idType, int ID, std::vector<int>& fieldIDs)
1385 {
1386  int idindex = fei::binarySearch(idType, idTypes_);
1387  if (idindex < 0) {
1388  fieldIDs.resize(0);
1389  return;
1390  }
1391 
1392  fei::Record<int>* record = recordCollections_[idindex]->getRecordWithID(ID);
1393  if (record == NULL) {
1394  voidERReturn;
1395  }
1396 
1397  fei::FieldMask* fieldMask = record->getFieldMask();
1398  std::vector<int>& maskFieldIDs = fieldMask->getFieldIDs();
1399  int numFields = maskFieldIDs.size();
1400  fieldIDs.resize(numFields);
1401  for(int i=0; i<numFields; ++i) {
1402  fieldIDs[i] = maskFieldIDs[i];
1403  }
1404 }
1405 
1406 //----------------------------------------------------------------------------
1407 void fei::VectorSpace::getGlobalIndexOffsets(std::vector<int>& globalOffsets) const
1408 {
1409  globalOffsets = globalOffsets_;
1410 }
1411 
1412 //----------------------------------------------------------------------------
1413 void fei::VectorSpace::getGlobalBlkIndexOffsets(std::vector<int>& globalBlkOffsets) const
1414 {
1415  globalBlkOffsets = globalIDOffsets_;
1416 }
1417 
1418 //----------------------------------------------------------------------------
1420 {
1421  if (globalIndex < 0) return(-1);
1422 
1423  unsigned len = globalOffsets_.size();
1424  for(int i=0; i<(int)(len-1); ++i) {
1425  if (globalIndex < globalOffsets_[i+1]) {
1426  return(i);
1427  }
1428  }
1429 
1430  return(-1);
1431 }
1432 
1433 //----------------------------------------------------------------------------
1435 {
1436  if (globalIndex < 0) return(-1);
1437 
1438  unsigned len = globalOffsets_.size();
1439  for(int i=0; i<(int)(len-1); ++i) {
1440  if (globalIndex < globalIDOffsets_[i+1]) {
1441  return(i);
1442  }
1443  }
1444 
1445  return(-1);
1446 }
1447 
1448 //----------------------------------------------------------------------------
1450 {
1451  int idx = fei::binarySearch(idType, idTypes_);
1452  if (idx < 0) return(0);
1453 
1454  return( recordCollections_[idx]->getNumRecords() );
1455 }
1456 
1457 //----------------------------------------------------------------------------
1459 {
1460  int idx = fei::binarySearch(idType, idTypes_);
1461  if (idx < 0) return(0);
1462 
1463  fei::RecordAttributeCounter attrCounter(fei::localProc(comm_));
1464  runRecords(attrCounter, idx);
1465 
1466  return( attrCounter.numLocallyOwnedIDs_ );
1467 }
1468 
1469 //----------------------------------------------------------------------------
1470 int fei::VectorSpace::getNumSharedIDs(int idType, int& numShared)
1471 {
1472  numShared = sharedIDTables_[idType].getSharedIDs().size();
1473 
1474  return(0);
1475 }
1476 
1477 //----------------------------------------------------------------------------
1479  int lenList,
1480  int* IDs,
1481  int& numLocalIDs)
1482 {
1483  int idx = fei::binarySearch(idType, idTypes_);
1484  if (idx < 0) return(-1);
1485 
1486  snl_fei::RecordCollection* records = recordCollections_[idx];
1487  std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
1488  std::map<int,int>::iterator
1489  it = global_to_local.begin(),
1490  end = global_to_local.end();
1491  numLocalIDs = records->getNumRecords();
1492  int len = numLocalIDs;
1493  if (lenList < len) len = lenList;
1494  int i=0;
1495  for(; it!=end; ++it) {
1496  int local_id = it->second;
1497  IDs[i++] = records->getRecordWithLocalID(local_id)->getID();
1498  if (i >= len) break;
1499  }
1500 
1501  return(0);
1502 }
1503 
1504 //----------------------------------------------------------------------------
1506  int lenList,
1507  int* IDs,
1508  int& numLocalIDs)
1509 {
1510  int idx = fei::binarySearch(idType, idTypes_);
1511  if (idx < 0) return(-1);
1512 
1513  snl_fei::RecordCollection* records = recordCollections_[idx];
1514 
1515  std::map<int,int>& global_to_local = records->getGlobalToLocalMap();
1516  std::map<int,int>::iterator
1517  it = global_to_local.begin(),
1518  end = global_to_local.end();
1519 
1520  numLocalIDs = 0;
1521 
1522  for(; it!=end; ++it) {
1523  fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
1524 
1525  if (thisrecord.getOwnerProc() == fei::localProc(comm_)) {
1526  if (numLocalIDs < lenList) {
1527  IDs[numLocalIDs] = thisrecord.getID();
1528  }
1529  ++numLocalIDs;
1530  }
1531  }
1532 
1533  return(0);
1534 }
1535 
1536 //----------------------------------------------------------------------------
1538 {
1539  return(eqnNumbers_.size());
1540 }
1541 
1542 //----------------------------------------------------------------------------
1543 int fei::VectorSpace::getIndices_SharedAndOwned(std::vector<int>& globalIndices) const
1544 {
1545  if (eqnNumbers_.size() == 0) {
1546  globalIndices.resize(0);
1547  return(0);
1548  }
1549 
1550  size_t numIndices = eqnNumbers_.size();
1551  const int* indicesPtr = &eqnNumbers_[0];
1552 
1553  globalIndices.resize(numIndices);
1554  for(size_t i=0; i<numIndices; ++i) {
1555  globalIndices[i] = indicesPtr[i];
1556  }
1557 
1558  return(0);
1559 }
1560 
1561 //----------------------------------------------------------------------------
1563 {
1564  numBlkIndices = 0;
1565  for(size_t i=0; i<recordCollections_.size(); ++i) {
1566  numBlkIndices += recordCollections_[i]->getNumRecords();
1567  }
1568 
1569  return(0);
1570 }
1571 
1572 //----------------------------------------------------------------------------
1574  int* globalBlkIndices,
1575  int* blkSizes,
1576  int& numBlkIndices)
1577 {
1578  if (!sharedRecordsSynchronized_) {
1579  numBlkIndices = 0;
1580  return(-1);
1581  }
1582 
1583  fei::BlkIndexAccessor blkIndAccessor(lenBlkIndices,
1584  globalBlkIndices, blkSizes);
1585  runRecords(blkIndAccessor);
1586 
1587  numBlkIndices = blkIndAccessor.numBlkIndices_;
1588 
1589  return(0);
1590 }
1591 
1592 //----------------------------------------------------------------------------
1594 {
1595  if (globalOffsets_.size() < 1) return(0);
1596  return(globalOffsets_[globalOffsets_.size()-1]);
1597 }
1598 
1599 //----------------------------------------------------------------------------
1601 {
1602  if (!sharedRecordsSynchronized_) {
1603  return(-1);
1604  }
1605 
1606  int localProc = fei::localProc(comm_);
1607  int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
1608 
1609  return(numIndices);
1610 }
1611 
1612 //----------------------------------------------------------------------------
1613 int fei::VectorSpace::getIndices_Owned(std::vector<int>& globalIndices) const
1614 {
1615  if (!sharedRecordsSynchronized_) {
1616  globalIndices.clear();
1617  return(-1);
1618  }
1619 
1620  int localProc = fei::localProc(comm_);
1621  int numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
1622 
1623  globalIndices.resize(numIndices);
1624 
1625  int firstOffset = globalOffsets_[localProc];
1626  for(int i=0; i<numIndices; ++i) {
1627  globalIndices[i] = firstOffset+i;
1628  }
1629 
1630  return(0);
1631 }
1632 
1633 //----------------------------------------------------------------------------
1634 int fei::VectorSpace::getIndices_Owned(int lenIndices, int* globalIndices, int& numIndices) const
1635 {
1636  if (!sharedRecordsSynchronized_) {
1637  numIndices = 0;
1638  return(-1);
1639  }
1640 
1641  int localProc = fei::localProc(comm_);
1642  numIndices = globalOffsets_[localProc+1]-globalOffsets_[localProc];
1643 
1644  int len = lenIndices >= numIndices ? numIndices : lenIndices;
1645 
1646  int firstOffset = globalOffsets_[localProc];
1647  for(int i=0; i<len; ++i) {
1648  globalIndices[i] = firstOffset+i;
1649  }
1650 
1651  return(0);
1652 }
1653 
1654 //----------------------------------------------------------------------------
1656 {
1657  if (!sharedRecordsSynchronized_) {
1658  return(-1);
1659  }
1660 
1661  int localProc = fei::localProc(comm_);
1662  int numBlkIndices = globalIDOffsets_[localProc+1]-globalIDOffsets_[localProc];
1663 
1664  return(numBlkIndices);
1665 }
1666 
1667 //----------------------------------------------------------------------------
1669 {
1670  int numBlkIndices = 0;
1671  unsigned len = globalIDOffsets_.size();
1672  if (len > 0) {
1673  numBlkIndices = globalIDOffsets_[len-1];
1674  }
1675 
1676  return(numBlkIndices);
1677 }
1678 
1679 //----------------------------------------------------------------------------
1681  int* globalBlkIndices,
1682  int* blkSizes,
1683  int& numBlkIndices)
1684 {
1685  if (!sharedRecordsSynchronized_) {
1686  numBlkIndices = 0;
1687  return(-1);
1688  }
1689 
1690  int localProc = fei::localProc(comm_);
1691  fei::BlkIndexAccessor blkIndAccessor(localProc, lenBlkIndices,
1692  globalBlkIndices, blkSizes);
1693  runRecords(blkIndAccessor);
1694 
1695  numBlkIndices = blkIndAccessor.numBlkIndices_;
1696 
1697  return(0);
1698 }
1699 
1700 //----------------------------------------------------------------------------
1702  snl_fei::RecordCollection*& records)
1703 {
1704  int idx = fei::binarySearch(idType, idTypes_);
1705  if (idx < 0) return(-1);
1706 
1707  records = recordCollections_[idx];
1708  return(0);
1709 }
1710 
1711 //----------------------------------------------------------------------------
1713  const snl_fei::RecordCollection*& records) const
1714 {
1715  int idx = fei::binarySearch(idType, idTypes_);
1716  if (idx < 0) return(-1);
1717 
1718  records = recordCollections_[idx];
1719  return(0);
1720 }
1721 
1722 //----------------------------------------------------------------------------
1723 void fei::VectorSpace::setOwners_shared()
1724 {
1725  //first, add localProc to each of the sharing-proc lists, in case it wasn't
1726  //included via initSharedIDs().
1727  int localProc = fei::localProc(comm_);
1728 
1729  std::map<int, fei::SharedIDs<int> >::iterator
1730  s_iter = sharedIDTables_.begin(), s_end = sharedIDTables_.end();
1731 
1732  for(; s_iter != s_end; ++s_iter) {
1733  fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
1734 
1736  t_iter = shid_table.begin(),
1737  t_end = shid_table.end();
1738 
1739  for(; t_iter != t_end; ++t_iter) {
1740  std::set<int>& shProcs = t_iter->second;
1741  shProcs.insert(localProc);
1742  }
1743  }
1744 
1745  //now set the owningProcs for the SharedIDs records, and the owning procs on
1746  //the appropriate records in the recordCollections. Set the owner to be the
1747  //lowest-numbered sharing proc in all cases.
1748  s_iter = sharedIDTables_.begin();
1749  for(; s_iter != s_end; ++s_iter) {
1750  fei::SharedIDs<int>::map_type& shid_table = s_iter->second.getSharedIDs();
1751 
1752  std::vector<int>& owningProcs = s_iter->second.getOwningProcs();
1753 
1754  int len = shid_table.size();
1755  owningProcs.resize(len);
1756  int j = 0;
1757 
1759  t_iter = shid_table.begin(),
1760  t_end = shid_table.end();
1761 
1762  for(; t_iter != t_end; ++t_iter, ++j) {
1763  std::set<int>& shProcs = (*t_iter).second;
1764  int lowest = *(shProcs.begin());
1765  owningProcs[j] = lowest;
1766  }
1767 
1768  int idx = fei::binarySearch(s_iter->first, idTypes_);
1769  if (idx < 0) {
1770  throw std::runtime_error("internal error in fei::VectorSapce::setOwners_lowestSharing");
1771  }
1772 
1773  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1774  recordCollections_[idx]->setDebugOutput(output_stream_);
1775  }
1776 
1777  recordCollections_[idx]->setOwners_lowestSharing(s_iter->second);
1778  }
1779 }
1780 
1781 //----------------------------------------------------------------------------
1782 int fei::VectorSpace::calculateGlobalIndices()
1783 {
1784  int localProc = fei::localProc(comm_);
1785  int numProcs = fei::numProcs(comm_);
1786  std::vector<int> localOffsets(numProcs +1, 0);
1787  std::vector<int> localIDOffsets(numProcs+1, 0);
1788  globalOffsets_.resize(numProcs + 1, 0);
1789 
1790  globalIDOffsets_.assign(globalIDOffsets_.size(), 0);
1791 
1792  //first we'll calculate the number of local degrees of freedom, and the
1793  //number of local identifiers.
1794  fei::RecordAttributeCounter counter(localProc);
1795  runRecords(counter);
1796 
1797  int numLocalDOF = counter.numLocalDOF_;
1798  int numLocalIDs = counter.numLocalIDs_;
1799  int numRemoteSharedDOF = counter.numRemoteSharedDOF_;
1800 
1801  eqnNumbers_.resize(numLocalDOF+numRemoteSharedDOF);
1802 
1803  localOffsets[localProc] = numLocalDOF;
1804  CHK_MPI( fei::GlobalMax(comm_, localOffsets, globalOffsets_) );
1805 
1806  localIDOffsets[localProc] = numLocalIDs;
1807  CHK_MPI( fei::GlobalMax(comm_, localIDOffsets, globalIDOffsets_) );
1808 
1809  //Now globalOffsets_ contains numLocalDOF for proc i, in the i-th position.
1810  //(and similarly for globalIDOffsets_)
1811  //So all we need to do is turn that data into global-offsets (i.e., the
1812  //starting global-offset for each processor).
1813  int localOffset = 0;
1814  int localIDOffset = 0;
1815  for(int p=0; p<numProcs; ++p) {
1816  numLocalDOF = globalOffsets_[p];
1817  globalOffsets_[p] = localOffset;
1818  localOffset += numLocalDOF;
1819  numLocalIDs = globalIDOffsets_[p];
1820  globalIDOffsets_[p] = localIDOffset;
1821  localIDOffset += numLocalIDs;
1822  }
1823  globalOffsets_[numProcs] = localOffset;
1824  globalIDOffsets_[numProcs] = localIDOffset;
1825 
1826  firstLocalOffset_ = globalOffsets_[localProc];
1827  lastLocalOffset_ = globalOffsets_[localProc+1] - 1;
1828 
1829  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
1830  FEI_OSTREAM& os = *output_stream_;
1831  os <<dbgprefix_<< " firstLocalOffset_: " << firstLocalOffset_ << ", lastLocalOffset_: "
1832  << lastLocalOffset_ << FEI_ENDL;
1833  }
1834 
1835  //Now we're ready to set the equation-numbers on all local records.
1836  CHK_ERR( setLocalEqnNumbers() );
1837 
1838  return(0);
1839 }
1840 
1841 //----------------------------------------------------------------------------
1842 int fei::VectorSpace::synchronizeSharedRecords()
1843 {
1844  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1845  FEI_OSTREAM& os = *output_stream_;
1846  os<<dbgprefix_<<"#synchronizeSharedRecords num-field-masks: "<<fieldMasks_.size()<<FEI_ENDL;
1847  for(unsigned fm=0; fm<fieldMasks_.size(); ++fm) {
1848  os << dbgprefix_<<"# maskID["<<fm<<"]: " << fieldMasks_[fm]->getMaskID() << FEI_ENDL;
1849  }
1850  }
1851 
1852  bool safetyCheck = checkSharedIDs_;
1853 
1854  int numProcs = fei::numProcs(comm_);
1855  int localProc = fei::localProc(comm_);
1856  int numShTables = sharedIDTables_.size();
1857 
1858  if (numProcs < 2) {
1859  sharedRecordsSynchronized_ = true;
1860  return(0);
1861  }
1862 
1863  if (safetyCheck) {
1864  if (localProc == 0 && output_level_ >= fei::BRIEF_LOGS) {
1865  FEI_COUT << "fei::VectorSpace: global consistency-check of shared ID"
1866  << " data (involves all-to-all communication). This is done "
1867  << "only if requested by parameter 'FEI_CHECK_SHARED_IDS true'."
1868  << FEI_ENDL;
1869  }
1870 
1871  int totalNumShTables = 0;
1872  CHK_ERR( fei::GlobalSum(comm_, numShTables, totalNumShTables) );
1873  if (totalNumShTables != numShTables * numProcs) {
1874  //not all processors have the same number of shared-id tables. This means
1875  //that one or more processors is 'empty', which may not be an error. But
1876  //it does mean that the safety check can't be performed because the safety
1877  //check involves all-to-all communication and if one or more processors
1878  //don't participate then we'll hang...
1879  safetyCheck = false;
1880  }
1881  }
1882 
1883  //Create a list of fei::comm_maps which will be the communication-patterns
1884  //for each of the sharedIDTables. A sharedIDTable maps lists of processors
1885  //to each shared ID. The communication-pattern will be a transpose of
1886  //that, mapping lists of IDs to owning or sharing processors.
1887  //
1888 
1889  int local_err = 0;
1890 
1891  for(size_t i=0; i<idTypes_.size(); ++i) {
1892 
1893  fei::SharedIDs<int>& shared = getSharedIDs(idTypes_[i]);
1894 
1895  //now create/initialize ownerPatterns which map owning processors to lists
1896  //of ids that are shared locally, and sharerPatterns which map sharing
1897  //processors to lists of ids that are owned locally.
1898  fei::comm_map* ownerPattern = new fei::comm_map(1, numProcs);
1899 
1900  std::map<int,fei::comm_map*>::iterator iter = ownerPatterns_.find(idTypes_[i]);
1901  if (iter == ownerPatterns_.end()) {
1902  ownerPatterns_.insert(std::make_pair(idTypes_[i], ownerPattern));
1903  }
1904  else {
1905  delete iter->second;
1906  iter->second = ownerPattern;
1907  }
1908 
1909  fei::comm_map* sharerPattern = new fei::comm_map(1, numProcs);
1910 
1911  iter = sharerPatterns_.find(idTypes_[i]);
1912  if (iter == sharerPatterns_.end()) {
1913  sharerPatterns_.insert(std::make_pair(idTypes_[i], sharerPattern));
1914  }
1915  else {
1916  delete iter->second;
1917  iter->second = sharerPattern;
1918  }
1919 
1920  std::vector<int>& owningProcs = shared.getOwningProcs();
1921 
1922  fei::SharedIDs<int>::map_type& shtable = shared.getSharedIDs();
1923 
1925  s_iter = shtable.begin(),
1926  s_end = shtable.end();
1927 
1928  int j = 0;
1929  for(; s_iter != s_end; ++s_iter, ++j) {
1930  int ID = (*s_iter).first;
1931  std::set<int>& shProcs = (*s_iter).second;
1932 
1933  int owner = owningProcs[j];
1934 
1935  if (owner == localProc) {
1936  std::set<int>::const_iterator
1937  p_iter = shProcs.begin(),
1938  p_end = shProcs.end();
1939  for(; p_iter != p_end; ++p_iter) {
1940  if (*p_iter != localProc) {
1941  sharerPattern->addIndices(*p_iter, 1, &ID);
1942  }
1943  }
1944  }
1945  else {
1946  ownerPattern->addIndices(owner, 1, &ID);
1947  }
1948  }
1949 
1950  if (safetyCheck) {
1951  fei::comm_map* checkPattern = NULL;
1952  CHK_ERR( fei::mirrorCommPattern(comm_, ownerPattern, checkPattern) );
1953  fei::Barrier(comm_);
1954  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1955  FEI_OSTREAM& os = *output_stream_;
1956  fei::comm_map::map_type& owner_map = ownerPattern->getMap();
1957  int numKeys = owner_map.size();
1958  fei::comm_map::map_type::const_iterator
1959  omap_iter = owner_map.begin(),
1960  omap_end = owner_map.end();
1961 
1962  os << dbgprefix_<<"# synchronizeSharedRecords" << FEI_ENDL
1963  << dbgprefix_<<"# ownerPattern, num-procs-to-send-to: " << numKeys << FEI_ENDL;
1964  for(int sk=0; omap_iter != omap_end; ++sk, ++omap_iter) {
1965  os << dbgprefix_<<"# sendProc["<<sk<<"]: " << omap_iter->first << " IDs: ";
1966  fei::comm_map::row_type::const_iterator
1967  val_iter = omap_iter->second->begin(),
1968  val_end = omap_iter->second->end();
1969  for(; val_iter != val_end; ++val_iter) {
1970  os << *val_iter << " ";
1971  }
1972  os << FEI_ENDL;
1973  }
1974 
1975  fei::comm_map::map_type& check_map = checkPattern->getMap();
1976  int numCheckKeys = check_map.size();
1977  fei::comm_map::map_type::const_iterator
1978  cmap_iter = check_map.begin(),
1979  cmap_end = check_map.end();
1980 
1981  os <<dbgprefix_<< "# synchronizeSharedRecords" << FEI_ENDL
1982  <<dbgprefix_<< "# checkPattern (send mirror), num-procs: "
1983  << numCheckKeys << FEI_ENDL;
1984  for(int sk=0; cmap_iter != cmap_end; ++sk, ++cmap_iter) {
1985  os <<dbgprefix_<< "# proc["<<sk<<"]: " << cmap_iter->first << " IDs: ";
1986  fei::comm_map::row_type::const_iterator
1987  val_iter = cmap_iter->second->begin(),
1988  val_end = cmap_iter->second->end();
1989  for(; val_iter != val_end; ++val_iter) {
1990  os << *val_iter << " ";
1991  }
1992  os << FEI_ENDL;
1993  }
1994  }
1995 
1996  int err = 0;
1997  bool quiet = false;
1998  if (!checkPattern->equal(*sharerPattern, quiet)) {
1999  //fei::console_out() << "shared-ID safety-check failed..."
2000  // << FEI_ENDL;
2001  err = -1;
2002  }
2003  int globalErr = 0;
2004  CHK_ERR( fei::GlobalSum(comm_, err, globalErr) );
2005 
2006  delete checkPattern;
2007 
2008  if (globalErr != 0) {
2009  return(globalErr);
2010  }
2011  }
2012 
2013  local_err += exchangeFieldInfo(ownerPattern, sharerPattern,
2014  recordCollections_[i], fieldMasks_);
2015  }
2016 
2017  int global_err = 0;
2018  CHK_ERR( fei::GlobalSum(comm_, local_err, global_err) );
2019 
2020  if (global_err != 0) {
2021  ERReturn(-1);
2022  }
2023 
2024  sharedRecordsSynchronized_ = true;
2025 
2026  return(0);
2027 }
2028 
2029 //----------------------------------------------------------------------------
2030 int fei::VectorSpace::exchangeGlobalIndices()
2031 {
2032  for(size_t i=0; i<idTypes_.size(); ++i) {
2033  snl_fei::RecordMsgHandler recmsghndlr(fei::localProc(comm_),
2034  recordCollections_[i], *ptBlkMap_,
2035  fieldMasks_, eqnNumbers_);
2036  recmsghndlr.setTask(snl_fei::RecordMsgHandler::_EqnNumbers_);
2037 
2038  recmsghndlr.setSendPattern(sharerPatterns_[idTypes_[i]]);
2039  recmsghndlr.setRecvPattern(ownerPatterns_[idTypes_[i]]);
2040  CHK_ERR( fei::exchange(comm_, &recmsghndlr) );
2041  }
2042 
2043  return(0);
2044 }
2045 
2046 //----------------------------------------------------------------------------
2047 void fei::VectorSpace::runRecords(fei::Record_Operator<int>& record_op)
2048 {
2049  for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
2050  snl_fei::RecordCollection* records = recordCollections_[rec];
2051  std::map<int,int>& g2l = records->getGlobalToLocalMap();
2052  std::map<int,int>::iterator
2053  it = g2l.begin(),
2054  end= g2l.end();
2055 
2056  for(; it!=end; ++it) {
2057  fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
2058 
2059  record_op(thisrecord);
2060  }
2061  }
2062 }
2063 
2064 //----------------------------------------------------------------------------
2065 void fei::VectorSpace::runRecords(fei::Record_Operator<int>& record_op, int recordIndex)
2066 {
2067  snl_fei::RecordCollection* records = recordCollections_[recordIndex];
2068  std::map<int,int>& g2l = records->getGlobalToLocalMap();
2069  std::map<int,int>::iterator
2070  it = g2l.begin(),
2071  end= g2l.end();
2072 
2073  for(; it!=end; ++it) {
2074  fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
2075 
2076  record_op(thisrecord);
2077  }
2078 }
2079 
2080 //----------------------------------------------------------------------------
2081 int fei::VectorSpace::setLocalEqnNumbers()
2082 {
2083  int proc = fei::localProc(comm_);
2084  int eqnNumber = firstLocalOffset_;
2085  int idNumber = globalIDOffsets_[proc];
2086 
2087  int numProcs = fei::numProcs(comm_);
2088  int localProc = fei::localProc(comm_);
2089 
2090  if (ptBlkMap_ != NULL) {
2091  delete ptBlkMap_;
2092  }
2093  ptBlkMap_ = new snl_fei::PointBlockMap;
2094 
2095  int maxNumIndices = 0;
2096  for(unsigned i=0; i<fieldMasks_.size(); ++i) {
2097  if (fieldMasks_[i]->getNumIndices() > maxNumIndices) {
2098  maxNumIndices = fieldMasks_[i]->getNumIndices();
2099  }
2100  }
2101 
2102  if (maxNumIndices == 1) {
2103  ptBlkMap_->setPtEqualBlk();
2104  }
2105 
2106  FEI_OSTREAM* id2eqnStream = NULL;
2107  if (output_level_ >= fei::BRIEF_LOGS) {
2108  std::string path = fei::LogManager::getLogManager().getOutputPath();
2109  if (path == "") path = ".";
2110  FEI_OSTRINGSTREAM osstr;
2111  osstr << path << "/fei_ID2Eqn";
2112  if (name_.size() > 0) osstr << "_" << name_;
2113  osstr << "." <<numProcs<< "." << localProc;
2114 
2115  id2eqnStream = new FEI_OFSTREAM(osstr.str().c_str(), IOS_OUT);
2116  FEI_OSTREAM& os = *id2eqnStream;
2117  os << "# Each line contains:\n# ID blk-eqn eqn" << FEI_ENDL;
2118  }
2119 
2120  int eqnNumberOffset = 0;
2121  int maxNumDOF = 0;
2122 
2123  for(size_t rec=0; rec<recordCollections_.size(); ++rec) {
2124  snl_fei::RecordCollection* records = recordCollections_[rec];
2125 
2126  const std::map<int,int>& rmap = records->getGlobalToLocalMap();
2127  std::map<int,int>::const_iterator
2128  it = rmap.begin(), it_end = rmap.end();
2129 
2130  int* eqnNumPtr = eqnNumbers_.empty() ? NULL : &eqnNumbers_[0];
2131 
2132  for(; it!=it_end; ++it) {
2133  fei::Record<int>& thisrecord = *records->getRecordWithLocalID(it->second);
2134 
2135  fei::FieldMask* mask = thisrecord.getFieldMask();
2136 
2137  thisrecord.setOffsetIntoEqnNumbers(eqnNumberOffset);
2138 
2139  int owner = thisrecord.getOwnerProc();
2140  if (owner == proc) {
2141  thisrecord.setNumber(idNumber++);
2142  }
2143 
2144  int* eqnNumbers = eqnNumPtr
2145  + thisrecord.getOffsetIntoEqnNumbers();
2146 
2147  int numDOF = mask->getNumIndices();
2148  eqnNumberOffset += numDOF;
2149 
2150  if (output_level_ >= fei::BRIEF_LOGS) {
2151  for(int ii=0; ii<numDOF; ++ii) {
2152  if (isLogEqn(eqnNumber+ii) && output_stream_ != NULL) {
2153  FEI_OSTREAM& os = *output_stream_;
2154  os <<dbgprefix_<< "setLocalEqnNumbers: ID "<<thisrecord.getID()
2155  << " <--> eqn " << eqnNumber+ii<<FEI_ENDL;
2156  }
2157  }
2158  }
2159 
2160  if (owner == proc) {
2161  int offset = 0;
2162  for(int n=0; n<numDOF; ++n) {
2163  eqnNumbers[offset++] = eqnNumber++;
2164  }
2165  }
2166 
2167  if (numDOF > maxNumDOF) maxNumDOF = numDOF;
2168 
2169  if (owner == proc) {
2170  int thiseqn = eqnNumber-numDOF;
2171  int thisrecordnumber = thisrecord.getNumber();
2172  if (maxNumIndices > 1) {
2173  CHK_ERR( ptBlkMap_->setEqn(thiseqn, thisrecordnumber, numDOF) );
2174  if (numDOF > 1) {
2175  for(int i=1; i<numDOF; ++i) {
2176  CHK_ERR( ptBlkMap_->setEqn(thiseqn+i, thisrecordnumber, numDOF) );
2177  }
2178  }
2179  }
2180  }
2181 
2182  if (id2eqnStream != NULL) {
2183  if (owner == proc) {
2184  FEI_OSTREAM& os = *id2eqnStream;
2185  for(int n=0; n<numDOF; ++n) {
2186  os << thisrecord.getID() << " " << thisrecord.getNumber() << " "
2187  << eqnNumber-numDOF+n<<FEI_ENDL;
2188  }
2189  }
2190  }
2191  }
2192 
2193  }
2194 
2195  ptBlkMap_->setMaxBlkEqnSize(maxNumDOF);
2196 
2197  int globalMaxNumDOF;
2198  CHK_ERR( fei::GlobalMax(comm_, maxNumDOF, globalMaxNumDOF) );
2199 
2200  if (globalMaxNumDOF == 1) {
2201  ptBlkMap_->setPtEqualBlk();
2202  }
2203 
2204  delete id2eqnStream;
2205 
2206  return(0);
2207 }
2208 
2209 //----------------------------------------------------------------------------
2210 int fei::VectorSpace::exchangeFieldInfo(fei::comm_map* ownerPattern,
2211  fei::comm_map* sharerPattern,
2212  snl_fei::RecordCollection* recordCollection,
2213  std::vector<fei::FieldMask*>& fieldMasks)
2214 {
2215  //ownerPattern maps owning processors to lists of IDs that we (the local
2216  //processor) share.
2217  //
2218  //sharerPattern maps sharing processors to lists of IDs that we own.
2219  //
2220  //In this function we need to perform these tasks:
2221  //1. processors exchange and combine their sets of field-masks so that all
2222  // processors will have the super-set of field-masks that they need.
2223  //2. sharing processors send maskIDs for their shared IDs to the owning
2224  // processors. The owning processors then combine masks as necessary to
2225  // make sure each shared ID has the union of field-masks that are held by
2226  // each of the sharing processors. all owning processors now send their
2227  // maskIDs out to the sharing processors. At this point all processors
2228  // should have the same view of the masks that belong on shared IDs.
2229  //3. exchange info describing inactive DOF offsets for shared records.
2230  //
2231 
2232  int numProcs = fei::numProcs(comm_);
2233  if (numProcs < 2) return(0);
2234 
2235  if (ptBlkMap_ == 0) ptBlkMap_ = new snl_fei::PointBlockMap;
2236 
2237  snl_fei::RecordMsgHandler recMsgHandler(fei::localProc(comm_),
2238  recordCollection,
2239  *ptBlkMap_,
2240  fieldMasks, eqnNumbers_);
2241 
2242  //Step 1a.
2243  recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
2244 
2245  recMsgHandler.setSendPattern(ownerPattern);
2246  recMsgHandler.setRecvPattern(sharerPattern);
2247  CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
2248 
2249  //Step 2a.
2250  recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
2251 
2252  recMsgHandler.setSendPattern(ownerPattern);
2253  recMsgHandler.setRecvPattern(sharerPattern);
2254  CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
2255 
2256  //Step 1b.
2257  recMsgHandler.setTask(snl_fei::RecordMsgHandler::_FieldMasks_);
2258 
2259  recMsgHandler.setSendPattern(sharerPattern);
2260  recMsgHandler.setRecvPattern(ownerPattern);
2261  CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
2262 
2263  //Step 2b.
2264  recMsgHandler.setTask(snl_fei::RecordMsgHandler::_MaskIDs_);
2265 
2266  recMsgHandler.setSendPattern(sharerPattern);
2267  recMsgHandler.setRecvPattern(ownerPattern);
2268  CHK_ERR( fei::exchange(comm_, &recMsgHandler) );
2269 
2270  return(0);
2271 }
2272 
2273 //----------------------------------------------------------------------------
2275 fei::VectorSpace::getFieldDofMap()
2276 {
2277  return fieldDofMap_;
2278 }
2279 
2280 //----------------------------------------------------------------------------
2281 void fei::VectorSpace::setName(const char* name)
2282 {
2283  if (name == NULL) return;
2284 
2285  if (name_ == name) return;
2286 
2287  name_ = name;
2288  dbgprefix_ = "VecSpc_"+name_+": ";
2289 }
2290 
int getGlobalNumBlkIndices() const
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
GlobalIDType getNumber() const
Definition: fei_Record.hpp:58
MPI_Comm getCommunicator() const
int sortedListInsert(const T &item, std::vector< T > &list)
int getGlobalIndicesLocalIDs(int numIDs, const int *localIDs, int idType, int fieldID, int *globalIndices)
int initSharedIDs(int numShared, int idType, const int *sharedIDs, const int *numSharingProcsPerID, const int *sharingProcs)
snl_fei::RecordCollection *const * getRecordCollections() const
Definition: fei_Pattern.hpp:70
int getBlkIndices_Owned(int lenBlkIndices, int *globalBlkIndices, int *blkSizes, int &numBlkIndices)
void setNumber(const GlobalIDType &num)
Definition: fei_Record.hpp:52
ParamType getType() const
Definition: fei_Param.hpp:98
VectorSpace(MPI_Comm comm, const char *name=NULL)
std::vector< int > & getOwningProcs()
int getGlobalNumIndices() const
int getIndices_SharedAndOwned(std::vector< int > &globalIndices) const
int getNumIndices() const
Definition: fei_Pattern.hpp:92
std::map< T, std::set< int > > map_type
const Param * get(const char *name) const
map_type & getSharedIDs()
const int * getNumFieldsPerID() const
Definition: fei_Pattern.hpp:73
const std::string & getOutputPath()
int mirrorCommPattern(MPI_Comm comm, comm_map *inPattern, comm_map *&outPattern)
fei::OutputLevel string_to_output_level(const std::string &str)
Definition: fei_utils.cpp:58
int getOwnedIDs(int idtype, int lenList, int *IDs, int &numLocalIDs)
bool getBoolValue() const
Definition: fei_Param.hpp:122
int getOwnerProc() const
Definition: fei_Record.hpp:94
const int * getFieldIDs() const
Definition: fei_Pattern.hpp:76
void defineFields(int numFields, const int *fieldIDs, const int *fieldSizes, const int *fieldTypes=NULL)
int getTotalNumFields() const
Definition: fei_Pattern.hpp:85
void set_shared_ids(MPI_Comm comm, const snl_fei::RecordCollection &records, fei::SharedIDs< int > &sharedIDs, int lowest_global_id, int highest_global_id)
int addDOFs(int fieldID, int idType, int numIDs, const int *IDs)
int getNumIndices() const
int getBlkIndices_SharedAndOwned(int lenBlkIndices, int *globalBlkIndices, int *blkSizes, int &numBlkIndices)
int getIntValue() const
Definition: fei_Param.hpp:116
void setParameters(const fei::ParameterSet &paramset)
virtual fei::SharedPtr< VectorSpace > createVectorSpace(MPI_Comm, const char *name)
int addVectorSpace(fei::VectorSpace *inputSpace)
int getGlobalIndex(int idType, int ID, int fieldID, int fieldOffset, int whichComponentOfField, int &globalIndex)
void getIDTypes(std::vector< int > &idTypes) const
int getGlobalIndices(int numIDs, const int *IDs, int idType, int fieldID, int *globalIndices)
int binarySearch(const T &item, const T *list, int len)
int getOwnedAndSharedIDs(int idtype, int lenList, int *IDs, int &numOwnedAndSharedIDs)
int getNumOwnedIDs(int idType)
int getOwnerProcBlkIndex(int globalIndex)
size_t getNumFields() const
void setOutputLevel(OutputLevel olevel)
int getNumBlkIndices_SharedAndOwned(int &numBlkIndices) const
void defineIDTypes(int numIDTypes, const int *idTypes)
bool isLocallyOwned(int idType, int ID)
int GlobalMin(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
PatternType getPatternType() const
Definition: fei_Pattern.hpp:61
int getOwnerProcPtIndex(int globalIndex)
void getGlobalBlkIndexOffsets(std::vector< int > &globalBlkOffsets) const
std::vector< int > & getFieldIDs()
int getNumSharedIDs(int idType, int &numShared)
void setOffsetIntoEqnNumbers(int offset)
Definition: fei_Record.hpp:119
int getGlobalBlkIndex(int idType, int ID, int &globalBlkIndex)
int getIndices_Owned(std::vector< int > &globalIndices) const
void getFields(std::vector< int > &fieldIDs)
void getGlobalIndexOffsets(std::vector< int > &globalOffsets) const
void setOwnerProc(int owner)
Definition: fei_Record.hpp:88
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int getNumIDs() const
Definition: fei_Pattern.hpp:64
const std::string & getStringValue() const
Definition: fei_Param.hpp:104
std::map< int, int > & getGlobalToLocalMap()
void addSharedID(const T &ID, size_t numSharingProcs, const int *sharingProcs)
int getNumBlkIndices_Owned() const
int getNumOwnedAndSharedIDs(int idType)
int localProc(MPI_Comm comm)
static LogManager & getLogManager()
int getNumIndices_SharedAndOwned() const
int getFieldEqnOffset(int fieldID, int &offset) const
GlobalIDType getID() const
Definition: fei_Record.hpp:46
unsigned getFieldSize(int fieldID)
void addIndices(int row, int numIndices, const int *indices)
bool equal(const RaggedTable< MAP_TYPE, SET_TYPE > &rhs, bool quiet=true) const
int getOffsetIntoEqnNumbers() const
Definition: fei_Record.hpp:126
int getNumDegreesOfFreedom(int idType, int ID)
fei::FieldMask * getFieldMask()
Definition: fei_Record.hpp:106
int getRecordCollection(int idType, snl_fei::RecordCollection *&records)
int getGlobalBlkIndices(int numIDs, const int *IDs, int idType, int *globalBlkIndices)
int numProcs(MPI_Comm comm)
int getNumIndices_Owned() const
bool isLocal(int idType, int ID)