FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SNL_FEI_Structure.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 <limits>
13 #include <cmath>
14 #include <assert.h>
15 
16 #include "fei_defs.h"
17 
18 #include "fei_TemplateUtils.hpp"
19 #include <fei_CommUtils.hpp>
20 #include "snl_fei_Constraint.hpp"
22 
23 #include "fei_NodeDescriptor.hpp"
24 #include "fei_NodeCommMgr.hpp"
25 #include "fei_NodeDatabase.hpp"
26 
27 #include "fei_SlaveVariable.hpp"
28 
29 #include "fei_BlockDescriptor.hpp"
30 
31 #include "snl_fei_PointBlockMap.hpp"
32 #include "fei_ProcEqns.hpp"
33 #include "fei_EqnBuffer.hpp"
34 #include <fei_FillableMat.hpp>
35 #include <fei_CSRMat.hpp>
36 #include <fei_CSVec.hpp>
37 #include "fei_EqnCommMgr.hpp"
38 
39 #include "fei_Lookup.hpp"
40 #include "fei_ConnectivityTable.hpp"
41 #include "snl_fei_Utils.hpp"
42 #include "SNL_FEI_Structure.hpp"
43 
44 #undef fei_file
45 #define fei_file "SNL_FEI_Structure.cpp"
46 #include "fei_ErrMacros.hpp"
47 
48 //-----Constructor-------------------------------------------------------------
50  : comm_(comm),
51  localProc_(0),
52  masterProc_(0),
53  numProcs_(1),
54  fieldIDs_(),
55  fieldSizes_(),
56  fieldDatabase_(new std::map<int,int>),
57  fieldDofMap_(),
58  workarray_(),
59  blockIDs_(),
60  blocks_(),
61  connTables_(),
62  nodeDatabase_(NULL),
63  activeNodesInitialized_(false),
64  globalNodeOffsets_(),
65  globalEqnOffsets_(),
66  globalBlkEqnOffsets_(),
67  slaveVars_(NULL),
68  slaveEqns_(NULL),
69  slvEqnNumbers_(NULL),
70  numSlvs_(0),
71  lowestSlv_(0),
72  highestSlv_(0),
73  slaveMatrix_(NULL),
74  globalNumNodesVanished_(),
75  localVanishedNodeNumbers_(),
76  nodeCommMgr_(NULL),
77  eqnCommMgr_(NULL),
78  slvCommMgr_(NULL),
79  numGlobalEqns_(0),
80  numLocalEqns_(0),
81  localStartRow_(0),
82  localEndRow_(0),
83  numLocalNodalEqns_(0),
84  numLocalElemDOF_(0),
85  numLocalMultCRs_(0),
86  reducedStartRow_(0),
87  reducedEndRow_(0),
88  numLocalReducedRows_(0),
89  Kid_(NULL),
90  Kdi_(NULL),
91  Kdd_(NULL),
92  tmpMat1_(),
93  tmpMat2_(),
94  reducedEqnCounter_(0),
95  reducedRHSCounter_(0),
96  rSlave_(),
97  cSlave_(),
98  work_nodePtrs_(),
99  structureFinalized_(false),
100  generateGraph_(true),
101  sysMatIndices_(NULL),
102  blockMatrix_(false),
103  numGlobalEqnBlks_(0),
104  numLocalEqnBlks_(0),
105  numLocalReducedEqnBlks_(0),
106  localBlkOffset_(0),
107  localReducedBlkOffset_(0),
108  globalMaxBlkSize_(0),
109  firstLocalNodeNumber_(-1),
110  numGlobalNodes_(0),
111  sysBlkMatIndices_(NULL),
112  matIndicesDestroyed_(false),
113  workSpace_(),
114  blkEqnMapper_(new snl_fei::PointBlockMap()),
115  multCRs_(),
116  penCRs_(),
117  checkSharedNodes_(false),
118  name_(),
119  outputLevel_(0),
120  debugOutput_(false),
121  dbgPath_(),
122  dbgOStreamPtr_(NULL),
123  setDbgOutCalled_(true)
124 {
125  numProcs_ = 1, localProc_ = 0, masterProc_ = 0;
126 
127  localProc_ = fei::localProc(comm_);
128  numProcs_ = fei::numProcs(comm_);
129  masterProc_ = 0;
130 
131  slaveVars_ = new std::vector<SlaveVariable*>;
132  slaveEqns_ = new EqnBuffer;
133 
134  nodeCommMgr_ = new NodeCommMgr(comm_, *this);
135 
136  eqnCommMgr_ = new EqnCommMgr(comm_);
137  eqnCommMgr_->setNumRHSs(1);
138 
139  nodeDatabase_ = new NodeDatabase(fieldDatabase_, nodeCommMgr_);
140 
141  Kid_ = new fei::FillableMat;
142  Kdi_ = new fei::FillableMat;
143  Kdd_ = new fei::FillableMat;
144 }
145 
146 //-----------------------------------------------------------------------------
147 int SNL_FEI_Structure::parameters(int numParams, const char*const* paramStrings)
148 {
149  const char* param = NULL;
150 
151  param = snl_fei::getParamValue("outputLevel",numParams,paramStrings);
152  if (param != NULL){
153  std::string str(param);
154  FEI_ISTRINGSTREAM isstr(str);
155  isstr >> outputLevel_;
156  }
157 
158  param = snl_fei::getParam("debugOutput",numParams,paramStrings);
159  if (param != NULL){
160  debugOutput_ = true;
161  }
162 
163  param = snl_fei::getParam("debugOutputOff",numParams,paramStrings);
164  if (param != NULL){
165  debugOutput_ = false;
166  }
167 
168  param = snl_fei::getParam("FEI_CHECK_SHARED_IDS", numParams,paramStrings);
169  if (param != NULL){
170  checkSharedNodes_ = true;
171  }
172 
173  param = snl_fei::getParamValue("sharedNodeOwnership",
174  numParams,paramStrings);
175  if (param != NULL){
176  if (!strcmp(param, "LowNumberedProc")) {
177  nodeCommMgr_->setSharedOwnershipRule(NodeCommMgr::STRICTLY_LOW_PROC);
178  }
179  if (!strcmp(param, "ProcWithLocalElem")) {
180  nodeCommMgr_->setSharedOwnershipRule(NodeCommMgr::PROC_WITH_LOCAL_ELEM);
181  }
182  if (!strcmp(param, "SierraSpecifies")) {
183  nodeCommMgr_->setSharedOwnershipRule(NodeCommMgr::CALLER_SPECIFIES);
184  }
185  }
186 
187  return(0);
188 }
189 
190 //-----Destructor--------------------------------------------------------------
192 {
193  for(size_t j=0; j<slaveVars_->size(); j++) {
194  delete (*slaveVars_)[j];
195  }
196  delete slaveVars_;
197 
198  delete slaveEqns_;
199  delete slaveMatrix_;
200 
201  destroyBlockRoster();
202  destroyConnectivityTables();
203 
204  delete nodeCommMgr_;
205  delete eqnCommMgr_;
206  delete blkEqnMapper_;
207 
208  destroyMatIndices();
209 
210  deleteMultCRs();
211 
212  int numPCRs = penCRs_.size();
213  if (numPCRs > 0) {
214  fei::destroyValues(penCRs_);
215  penCRs_.clear();
216  }
217 
218  delete nodeDatabase_;
219  delete fieldDatabase_;
220 
221  delete Kid_;
222  delete Kdi_;
223  delete Kdd_;
224 }
225 
226 //------------------------------------------------------------------------------
227 int SNL_FEI_Structure::setDbgOut(FEI_OSTREAM& ostr,
228  const char* path, const char* feiName)
229 {
230  dbgOStreamPtr_ = &ostr;
231  setDbgOutCalled_ = true;
232  if (path != NULL) {
233  dbgPath_ = path;
234  }
235  else {
236  dbgPath_ = ".";
237  }
238  if (feiName != NULL) {
239  name_ = feiName;
240  }
241 
242  debugOutput_ = true;
243  return(0);
244 }
245 
246 //------------------------------------------------------------------------------
247 void SNL_FEI_Structure::destroyBlockRoster()
248 {
249  for(size_t i=0; i<blockIDs_.size(); i++) delete blocks_[i];
250  blocks_.resize(0);
251 }
252 
253 //------------------------------------------------------------------------------
254 void SNL_FEI_Structure::destroyConnectivityTables()
255 {
256  for(size_t i=0; i<blockIDs_.size(); i++) {
257  delete connTables_[i];
258  }
259  connTables_.resize(0);
260 }
261 
262 //------------------------------------------------------------------------------
263 void SNL_FEI_Structure::destroyMatIndices()
264 {
265  if (matIndicesDestroyed_ == true) return;
266 
267  delete [] sysMatIndices_;
268  sysMatIndices_ = NULL;
269 
270  delete [] sysBlkMatIndices_;
271  sysBlkMatIndices_ = NULL;
272 
273  matIndicesDestroyed_ = true;
274 }
275 
276 //------------------------------------------------------------------------------
277 const int* SNL_FEI_Structure::getNumFieldsPerNode(GlobalID blockID)
278 {
279  BlockDescriptor* block = NULL;
280  int err = getBlockDescriptor(blockID, block);
281  if (err) return(NULL);
282 
283  return(block->fieldsPerNodePtr());
284 }
285 
286 //------------------------------------------------------------------------------
287 const int* const* SNL_FEI_Structure::getFieldIDsTable(GlobalID blockID)
288 {
289  BlockDescriptor* block = NULL;
290  int err = getBlockDescriptor(blockID, block);
291  if (err) return(NULL);
292 
293  return(block->fieldIDsTablePtr());
294 }
295 
296 //------------------------------------------------------------------------------
297 void SNL_FEI_Structure::getElemBlockInfo(GlobalID blockID,
298  int& interleaveStrategy, int& lumpingStrategy,
299  int& numElemDOF, int& numElements,
300  int& numNodesPerElem, int& numEqnsPerElem)
301 {
302  BlockDescriptor* block = NULL;
303  int err = getBlockDescriptor(blockID, block);
304  if (err) {
305  interleaveStrategy = -1;
306  lumpingStrategy = -1;
307  numElemDOF = -1;
308  numElements = -1;
309  numNodesPerElem = -1;
310  numEqnsPerElem = -1;
311  }
312 
313  interleaveStrategy = block->getInterleaveStrategy();
314  lumpingStrategy = block->getLumpingStrategy();
315  numElemDOF = block->getNumElemDOFPerElement();
316  numElements = block->getNumElements();
317  numNodesPerElem = block->getNumNodesPerElement();
318  numEqnsPerElem = block->getNumEqnsPerElement();
319 }
320 
321 //------------------------------------------------------------------------------
322 int SNL_FEI_Structure::getEqnNumber(int nodeNumber, int fieldID)
323 {
324  //this function is only used by clients of the Lookup interface, who are
325  //expecting equations in the 'reduced' equation space.
326 
327  int eqnNumber;
328 
329  const NodeDescriptor* node = NULL;
330  CHK_ERR( nodeDatabase_->getNodeWithNumber(nodeNumber, node) );
331 
332  bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
333  if (!hasField) {
334  fei::console_out() << "SNL_FEI_Structure::getEqnNumber: ERROR, node with nodeNumber "
335  << nodeNumber << " doesn't have fieldID " << fieldID << FEI_ENDL;
336  ERReturn(-1);
337  }
338 
339  int reducedEqn = -1;
340  bool isSlave = translateToReducedEqn(eqnNumber, reducedEqn);
341  if (isSlave) return(-1);
342 
343  return(reducedEqn);
344 }
345 
346 //------------------------------------------------------------------------------
348 {
349  if (eqn < 0) return(-1);
350 
351  int len = globalEqnOffsets_.size(); // len is numProcs+1...
352 
353  for(int i=0; i<len-1; i++) {
354  if (eqn >= globalEqnOffsets_[i] && eqn < globalEqnOffsets_[i+1]) return(i);
355  }
356 
357  return(-1);
358 }
359 
360 //------------------------------------------------------------------------------
361 int SNL_FEI_Structure::initFields(int numFields,
362  const int *fieldSizes,
363  const int *fieldIDs,
364  const int *fieldTypes)
365 {
366  // store the incoming solution fields
367  //
368  if (debugOutput_) {
369  FEI_OSTREAM& ostr = dbgOut();
370  ostr << "FEI: initFields" << FEI_ENDL
371  << "#num-fields" << FEI_ENDL << numFields << FEI_ENDL;
372  int nf;
373  ostr << "#field-sizes" << FEI_ENDL;
374  for(nf=0; nf<numFields; nf++) {
375  ostr <<fieldSizes[nf] << " ";
376  }
377  ostr << FEI_ENDL << "#field-ids" << FEI_ENDL;
378  for(nf=0; nf<numFields; nf++) {
379  ostr << fieldIDs[nf] << " ";
380  }
381  if (fieldTypes != NULL) {
382  ostr << FEI_ENDL << "#field-types" << FEI_ENDL;
383  for(nf=0; nf<numFields; nf++) {
384  ostr << fieldTypes[nf] << " ";
385  }
386  }
387  ostr<<FEI_ENDL;
388  }
389 
390  for (int i=0; i<numFields; i++) {
391  fieldDatabase_->insert(std::pair<int,int>(fieldIDs[i], fieldSizes[i]));
392  int offs = fei::sortedListInsert(fieldIDs[i], fieldIDs_);
393  if (offs >= 0) fieldSizes_.insert(fieldSizes_.begin()+offs, fieldSizes[i]);
394 
395  if (fieldIDs[i] >= 0) {
396  if (fieldTypes != NULL) {
397  fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i], fieldTypes[i]);
398  }
399  else {
400  fieldDofMap_.add_field(fieldIDs[i], fieldSizes[i]);
401  }
402  }
403  }
404 
405  return(FEI_SUCCESS);
406 }
407 
408 //------------------------------------------------------------------------------
409 int SNL_FEI_Structure::initElemBlock(GlobalID elemBlockID,
410  int numElements,
411  int numNodesPerElement,
412  const int* numFieldsPerNode,
413  const int* const* nodalFieldIDs,
414  int numElemDofFieldsPerElement,
415  const int* elemDofFieldIDs,
416  int interleaveStrategy)
417 {
418  int nn, nf;
419  if (debugOutput_) {
420  FEI_OSTREAM& os = dbgOut();
421  os << "FEI: initElemBlock" << FEI_ENDL << "#elemBlockID" << FEI_ENDL
422  << (int)elemBlockID << FEI_ENDL;
423  os << "#numElements"<< FEI_ENDL << numElements << FEI_ENDL;
424  os << "#numNodesPerElement"<< FEI_ENDL <<numNodesPerElement << FEI_ENDL;
425  os << "#numFieldsPerNode -- one entry per node " << FEI_ENDL;
426  for(nn=0; nn<numNodesPerElement; nn++) os << numFieldsPerNode[nn]<<" ";
427  os << FEI_ENDL << "#nodalFieldIDs -- one row per node" << FEI_ENDL;
428  for(nn=0; nn<numNodesPerElement; ++nn) {
429  for(nf=0; nf<numFieldsPerNode[nn]; ++nf) os << nodalFieldIDs[nn][nf] << " ";
430  os << FEI_ENDL;
431  }
432  os << "#numElemDofFieldsPerElement" << FEI_ENDL
433  << numElemDofFieldsPerElement<<FEI_ENDL;
434  os << "#elemDofFieldIDs -- 'numElemDofFieldsPerElement' entries" << FEI_ENDL;
435  for(nn=0; nn<numElemDofFieldsPerElement; ++nn) {
436  os << elemDofFieldIDs[nn] << " ";
437  }
438  if (numElemDofFieldsPerElement > 0) os << FEI_ENDL;
439  os << "#interleaveStrategy" << FEI_ENDL << interleaveStrategy << FEI_ENDL;
440  }
441  int j;
442 
443  CHK_ERR( addBlock(elemBlockID) );
444  BlockDescriptor* block = NULL;
445  CHK_ERR( getBlockDescriptor(elemBlockID, block) );
446 
447  CHK_ERR( block->setNumNodesPerElement(numNodesPerElement) );
448  block->setNumElements(numElements);
449  block->setElemDofFieldIDs(numElemDofFieldsPerElement, elemDofFieldIDs);
450  block->setInterleaveStrategy(interleaveStrategy);
451 
452  int *fieldsPerNodePtr = block->fieldsPerNodePtr();
453 
454 // construct the list of nodal solution cardinalities for this block
455 
456  int numNodalEqns = 0;
457  int countDOF;
458  std::vector<int> distinctFields(0);
459 
460  for(j=0; j<numNodesPerElement; j++) {
461 
462  countDOF = 0;
463  for(int k = 0; k < numFieldsPerNode[j]; k++) {
464  fei::sortedListInsert(nodalFieldIDs[j][k], distinctFields);
465 
466  int fieldSize = getFieldSize(nodalFieldIDs[j][k]);
467  if (fieldSize < 0) {
468  fei::console_out() << "SNL_FEI_Structure::initElemBlock ERROR: fieldID " <<
469  nodalFieldIDs[j][k] << " has negative size. " << FEI_ENDL;
470  ERReturn(-1);
471  }
472  countDOF += fieldSize;
473  }
474 
475  fieldsPerNodePtr[j] = numFieldsPerNode[j];
476  numNodalEqns += countDOF;
477  }
478 
479  block->setNumDistinctFields(distinctFields.size());
480 
481  int numElemDOFPerElement = 0;
482  for(j=0; j<numElemDofFieldsPerElement; j++) {
483  int fieldSize = getFieldSize(elemDofFieldIDs[j]);
484  if (fieldSize < 0) {
485  fei::console_out() << "SNL_FEI_Structure::initElemBlock ERROR: elemDoffieldID " <<
486  elemDofFieldIDs[j] << " has negative size. " << FEI_ENDL;
487  ERReturn(-1);
488  }
489  numElemDOFPerElement += fieldSize;
490  }
491 
492  block->setNumElemDOFPerElement(numElemDOFPerElement);
493  block->setNumEqnsPerElement(numNodalEqns + numElemDOFPerElement);
494  block->setNumBlkEqnsPerElement(numNodesPerElement + numElemDOFPerElement);
495 
496 // cache a copy of the element fields array for later use...
497 
498  CHK_ERR( block->allocateFieldIDsTable() );
499  int **fieldIDsTablePtr = block->fieldIDsTablePtr();
500 
501  for (j = 0; j < numNodesPerElement; j++) {
502  for(int k = 0; k < numFieldsPerNode[j]; k++) {
503  fieldIDsTablePtr[j][k] = nodalFieldIDs[j][k];
504  }
505  }
506 
507 // create data structures for storage of element ID and topology info
508 
509  if (numElements > 0) {
510  CHK_ERR( allocateBlockConnectivity(elemBlockID) );
511  }
512 
513  return(FEI_SUCCESS);
514 }
515 
516 //------------------------------------------------------------------------------
517 int SNL_FEI_Structure::initElem(GlobalID elemBlockID,
518  GlobalID elemID,
519  const GlobalID* elemConn)
520 {
521  if (debugOutput_ && outputLevel_ > 2) {
522  FEI_OSTREAM& os = dbgOut();
523  os << "FEI: initElem"<<FEI_ENDL;
524  os << "#blkID"<<FEI_ENDL<<(int)elemBlockID<<FEI_ENDL<<"#elmID"<<FEI_ENDL
525  <<(int)elemID<< FEI_ENDL;
526  }
527 
528  //first get the block-descriptor for this elemBlockID...
529 
530  BlockDescriptor* block = NULL;
531  CHK_ERR( getBlockDescriptor(elemBlockID, block) );
532 
533  ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
534 
535  std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
536  GlobalID* conn = &((*connTable.elem_conn_ids)[0]);
537 
538  int elemIndex = elemIDList.size();
539  std::map<GlobalID,int>::iterator
540  iter = elemIDList.find(elemID);
541 
542  bool redundantElem = false;
543 
544  if (iter != elemIDList.end()) {
545  elemIndex = iter->second;
546  redundantElem = true;
547  }
548  else {
549  elemIDList.insert(std::make_pair(elemID,elemIndex));
550  }
551 
552  int numNodes = block->getNumNodesPerElement();
553 
554  if (debugOutput_ && outputLevel_ > 2) {
555  FEI_OSTREAM& os = dbgOut();
556  os << "#n-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL<<"#nodeIDs"<<FEI_ENDL;
557  for(int nn=0; nn<numNodes; nn++) os << (int)elemConn[nn] << " ";
558  os << FEI_ENDL;
559  }
560 
561  if (redundantElem) {
562  //redundantElem == true means this element has been initialized before.
563  //So we'll simply make sure the connectivity is the same as it was, and
564  //if it isn't return -1.
565 
566  int offset = elemIndex*numNodes;
567  for(int j=0; j<numNodes; j++) {
568  if ( conn[offset+j] != elemConn[j]) {
569  fei::console_out() << "SNL_FEI_Structure::initElem ERROR, elemID " << (int)elemID
570  << " registered more than once, with differing connectivity."
571  << FEI_ENDL;
572  return(-1);
573  }
574  }
575  }
576  else {
577  int offset = elemIndex*numNodes;
578  for(int j = 0; j < numNodes; j++) {
579  conn[offset+j] = elemConn[j];
580  }
581 
582  CHK_ERR( nodeDatabase_->initNodeIDs(&(conn[offset]), numNodes) );
583  }
584 
585  return(FEI_SUCCESS);
586 }
587 
588 //------------------------------------------------------------------------------
589 int SNL_FEI_Structure::initSlaveVariable(GlobalID slaveNodeID,
590  int slaveFieldID,
591  int offsetIntoSlaveField,
592  int numMasterNodes,
593  const GlobalID* masterNodeIDs,
594  const int* masterFieldIDs,
595  const double* weights,
596  double rhsValue)
597 {
598  if (debugOutput_) {
599  FEI_OSTREAM& os = dbgOut();
600  os << "FEI: initSlaveVariable" << FEI_ENDL;
601  os << "#slvNodeID" << FEI_ENDL << (int)slaveNodeID << FEI_ENDL
602  << "#slvFieldID"<< FEI_ENDL << slaveFieldID << FEI_ENDL
603  << "#offsetIntoSlvField" << FEI_ENDL << offsetIntoSlaveField << FEI_ENDL
604  << "#numMasterNodes" << FEI_ENDL << numMasterNodes << FEI_ENDL
605  << "#masterNodeIDs" << FEI_ENDL;
606  int nn;
607  for(nn=0; nn<numMasterNodes; ++nn) {
608  os <<(int)masterNodeIDs[nn]<<" ";
609  }
610  os << FEI_ENDL << "#masterFieldIDs" << FEI_ENDL;
611  for(nn=0; nn<numMasterNodes; ++nn) {
612  os <<masterFieldIDs[nn] << " ";
613  }
614  os << FEI_ENDL << "#field-sizes" << FEI_ENDL;
615  for(nn=0; nn<numMasterNodes; ++nn) {
616  int size = getFieldSize(masterFieldIDs[nn]);
617  os << size << " ";
618  }
619  os << FEI_ENDL << "#weights" << FEI_ENDL;
620  int offset = 0;
621  for(nn=0; nn<numMasterNodes; ++nn) {
622  int size = getFieldSize(masterFieldIDs[nn]);
623  for(int j=0; j<size; ++j) {
624  os << weights[offset++] << " ";
625  }
626  }
627  os << FEI_ENDL << "#rhsValue" << FEI_ENDL << rhsValue << FEI_ENDL;
628  }
629 
630  //Create and initialize a slave-variable object with the incoming data,
631  //and add it to our list.
632 
633  SlaveVariable* svar = new SlaveVariable;
634  svar->setNodeID(slaveNodeID);
635  svar->setFieldID(slaveFieldID);
636  svar->setFieldOffset(offsetIntoSlaveField);
637 
638  int woffset = 0;
639 
640  for(int i=0; i<numMasterNodes; i++) {
641  svar->addMasterNodeID(masterNodeIDs[i]);
642  svar->addMasterField(masterFieldIDs[i]);
643  int fieldSize = getFieldSize(masterFieldIDs[i]);
644  if (fieldSize < 0) ERReturn(-1);
645 
646  for(int j=0; j<fieldSize; j++) {
647  svar->addWeight(weights[woffset++]);
648  }
649  }
650 
651  addSlaveVariable(svar);
652 
653  return(FEI_SUCCESS);
654 }
655 
656 //------------------------------------------------------------------------------
657 int SNL_FEI_Structure::deleteMultCRs()
658 {
659  int numMCRs = multCRs_.size();
660  if (numMCRs > 0) {
661  fei::destroyValues(multCRs_);
662  multCRs_.clear();
663  }
664  return(0);
665 }
666 
667 //------------------------------------------------------------------------------
668 int SNL_FEI_Structure::initSharedNodes(int numSharedNodes,
669  const GlobalID *sharedNodeIDs,
670  const int* numProcsPerNode,
671  const int *const *sharingProcIDs)
672 {
673  if (debugOutput_) {
674  FEI_OSTREAM& os = dbgOut();
675  os << "FEI: initSharedNodes"<<FEI_ENDL;
676  os << "#n-nodes"<<FEI_ENDL<<numSharedNodes<<FEI_ENDL;
677  os << "#num-procs-per-node"<<FEI_ENDL;
678  int nn;
679  for(nn=0; nn<numSharedNodes; ++nn) {
680  os << numProcsPerNode[nn] << " ";
681  }
682  os << FEI_ENDL << "#following lines: nodeID sharing-procs" << FEI_ENDL;
683  for(nn=0; nn<numSharedNodes; nn++) {
684  os << (int)sharedNodeIDs[nn]<<" ";
685  for(int np=0; np<numProcsPerNode[nn]; np++) os<<sharingProcIDs[nn][np]<<" ";
686  os << FEI_ENDL;
687  }
688  os << "#end shared nodes"<<FEI_ENDL;
689  }
690 
691  // In this function we simply accumulate the incoming data into the
692  // NodeCommMgr object.
693  //
694  CHK_ERR( nodeCommMgr_->addSharedNodes(sharedNodeIDs, numSharedNodes,
695  sharingProcIDs, numProcsPerNode) );
696 
697  if (activeNodesInitialized_) {
698  //if active nodes have already been initialized, then we won't be
699  //re-running the element-connectivities during the next call to
700  //initComplete(), which means we won't have an opportunity to call the
701  //NodeCommMgr::informLocal() method for nodes that reside locally. So we
702  //need to check now for any locally residing shared nodes and make that
703  //call. This is expensive, but it's a case that only arises when constraints
704  //change between solves and nothing else changes. In general, constraints
705  //are the only structural features allowed to change without requiring a
706  //complete destruction and re-calculation of the FEI structure.
707  //
708  for(int i=0; i<numSharedNodes; ++i) {
709  for(size_t nc=0; nc<connTables_.size(); ++nc) {
710  if (connTables_[nc]->elem_conn_ids == NULL) continue;
711  int len = connTables_[nc]->elem_conn_ids->size();
712  if (len < 1) continue;
713  GlobalID* conn_ids = &((*connTables_[nc]->elem_conn_ids)[0]);
714  NodeDescriptor** nodes = &((*connTables_[nc]->elem_conn_ptrs)[0]);
715 
716  for(int j=0; j<len; ++j) {
717  if (conn_ids[j] == sharedNodeIDs[i]) {
718  CHK_ERR( nodeCommMgr_->informLocal(*(nodes[j])));
719  break;
720  }
721  }
722  }
723  }
724  }
725 
726  return(FEI_SUCCESS);
727 }
728 
729 //------------------------------------------------------------------------------
730 int SNL_FEI_Structure::initCRMult(int numCRNodes,
731  const GlobalID* CRNodes,
732  const int *CRFields,
733  int& CRID)
734 {
735  if (debugOutput_) {
736  FEI_OSTREAM& os = dbgOut();
737  os << "FEI: initCRMult" << FEI_ENDL << "# numCRNodes: "<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
738  os << "#CRNodes:"<<FEI_ENDL;
739  int nn;
740  for(nn=0; nn<numCRNodes; nn++) {
741  os << (int)CRNodes[nn]<<" ";
742  }
743  os << FEI_ENDL<<"#fields:"<<FEI_ENDL;
744  for(nn=0; nn<numCRNodes; nn++) {
745  os << CRFields[nn]<<" ";
746  }
747  os<<FEI_ENDL;
748  }
749  // tasks: add local constraint equations, determine sparsity
750  // patterns for the new constraint relation
751  //
752 
753  CRID = localProc_*100000 + multCRs_.size();
754  ConstraintType* multCRPtr = NULL;
755  addCR(CRID, multCRPtr, multCRs_);
756 
757  ConstraintType& multCR = *multCRPtr;
758 
759  multCR.setConstraintID(CRID);
760 
761  std::vector<GlobalID>& CRNodeArray = multCR.getMasters();
762  std::vector<int>& CRFieldArray = multCR.getMasterFieldIDs();
763 
764  for(int j = 0; j < numCRNodes; j++) {
765  CRNodeArray.push_back(CRNodes[j]);
766  CRFieldArray.push_back(CRFields[j]);
767  }
768 
769  if (debugOutput_) dbgOut() << "#(output) CRID:"<<FEI_ENDL << CRID << FEI_ENDL;
770 
771  return(FEI_SUCCESS);
772 }
773 
774 //------------------------------------------------------------------------------
775 int SNL_FEI_Structure::initCRPen(int numCRNodes,
776  const GlobalID* CRNodes,
777  const int *CRFields,
778  int& CRID)
779 {
780  if (debugOutput_) {
781  FEI_OSTREAM& os = dbgOut();
782  os << "initCRPen, numCRNodes: "<<numCRNodes<<FEI_ENDL;
783  for(int nn=0; nn<numCRNodes; nn++) {
784  os << " crNodeID: "<<(int)CRNodes[nn]<<", field: "<<CRFields[nn]<<FEI_ENDL;
785  }
786  }
787 
788  CRID = localProc_*100000 + penCRs_.size();
789 
790  ConstraintType* penCRPtr = NULL;
791  addCR(CRID, penCRPtr, penCRs_);
792 
793  ConstraintType& penCR = *penCRPtr;
794 
795  penCR.setConstraintID(CRID);
796  penCR.setIsPenalty(true);
797 
798  std::vector<GlobalID>& CRNodesArray = penCR.getMasters();
799 
800  std::vector<int>& CRFieldArray = penCR.getMasterFieldIDs();
801 
802  for(int i = 0; i < numCRNodes; i++) {
803  CRNodesArray.push_back(CRNodes[i]);
804  CRFieldArray.push_back(CRFields[i]);
805  }
806 
807  return(FEI_SUCCESS);
808 }
809 
810 //------------------------------------------------------------------------------
811 bool SNL_FEI_Structure::isInLocalElement(int nodeNumber)
812 {
813  //I'm going to make an assumption: if we're running in serial, then there
814  //are no remote anythings, and the node in question must be in a local
815  //element.
816  if (numProcs_ < 2) return(true);
817 
818  const NodeDescriptor* node = NULL;
819  int err = nodeDatabase_->getNodeWithNumber(nodeNumber, node);
820  if (err != 0) return(false);
821 
822  GlobalID nodeID = node->getGlobalNodeID();
823 
824  //now let's find out if this node is a shared node.
825  int shIndx = nodeCommMgr_->getSharedNodeIndex(nodeID);
826  if (shIndx < 0) {
827  //if shIndx < 0, then this isn't a shared node. For now, we will assume
828  //that it must be a local node.
829  return(true);
830  }
831 
832  //If we reach this line, then the node is a shared node. Let's now ask if
833  //it appears locally...
834  std::vector<GlobalID>& localNodes = nodeCommMgr_->getLocalNodeIDs();
835  int index = fei::binarySearch(nodeID, &localNodes[0], localNodes.size());
836  if (index >= 0) return(true);
837 
838  //If we reach this line, then the node is shared, but doesn't appear in the
839  //"localNodeIDs" held by NodeCommMgr. This means it is not in a local element,
840  //so we'll return false.
841  return(false);
842 }
843 
844 //------------------------------------------------------------------------------
845 int SNL_FEI_Structure::initComplete(bool generateGraph)
846 {
847  //This is the most significant function in SNL_FEI_Structure. This is where
848  //the underlying matrix structure is calculated from all the data that has
849  //been provided since construction. Calculating the equation space and
850  //forming the matrix structure is a multi-step process, which proceeds like
851  //this:
852  //
853  //1. finalizeActiveNodes()
854  // - makes sure that all shared nodes have been identified to the
855  // NodeDatabase, then has the NodeDatabase allocate its internal list of
856  // NodeDescriptor objects.
857  // - runs the element-connectivity tables and sets the nodal field and block
858  // info on the NodeDescriptors in the NodeDatabase. IMPORTANT: At this
859  // point, the nodeIDs in the connectivitiy tables are replaced with
860  // indices into the nodeDatabase to speed future lookups.
861  // - As the connectivity tables are being run, the NodeCommMgr is given the
862  // nodeID of each node that is connected to a local element, and the node's
863  // owner is set to the initial value of localProc_.
864  //1a. finalizeNodeCommMgr()
865  // - NodeCommMgr::initComplete is called, at which time the ownership of all
866  // shared nodes is determined.
867  //
868  //2. Next, all local nodal degrees of freedom can be counted. This, together
869  // with the number of lagrange multipliers and element-centered DOFs are used
870  // by the function calcGlobalEqnInfo(), which determines global equation
871  // counts and offsets.
872  //
873  //3. setNodalEqnInfo() runs all locally active nodes and for each one that's
874  // locally owned, sets global equation-numbers on it, as well as its
875  // numNodalDOF and block-entry equation-number.
876  //
877  //4. setElemDOFEqnInfo() and setMultCREqnInfo() associate global equation-
878  // numbers with each of the lagrange multipliers and element-dofs.
879  //
880  //5. NodeDatabase::synchronize() is called, which associates nodeNumbers with
881  // each node, and also in turn calls NodeCommMgr::exchangeEqnInfo() to obtain
882  // equation-numbers and field/block info for remotely-owned shared nodes.
883  //
884  //6. setNumNodesAndEqnsPerBlock() then runs the active nodes again and counts
885  // active nodes and equations per element-block. This couldn't be done before
886  // because NodeCommMgr::exchangeEqnInfo() may have found out about new blocks
887  // that are associated with shared nodes on other processors.
888  //
889  //7. calculateSlaveEqns() associates equation-numbers with slave-nodes and
890  // master-nodes that have been identified locally, then gathers slave
891  // equation info from all processors onto all other processors, so that all
892  // processors know about all slave equations.
893  //
894  //8. initializeEqnCommMgr() runs the shared-nodes from NodeCommMgr and lets
895  // EqnCommMgr know which equations will be sent to other processors, and
896  // which will be received from other processors.
897  //
898  //9. translate localStartRow_ and localEndRow_ into the reduced equation
899  // space, where they are called reducedStartRow_ and reducedEndRow_, after
900  // which we know how many local equations there are not including any
901  // slave equations. At this point we can allocate the sysMatIndices_ array
902  // of arrays, which will hold the point-entry matrix structure.
903  //
904  //10. formMatrixStructure()
905  // - initElemBlockStructure() for each element in each element-block,
906  // obtains scatter-indices and makes the corresponding structurally
907  // symmetric contributions to the matrix structure. This process includes
908  // calculations associated with reducing out slave equations, and with
909  // putting off-processor contributions (element contributions with shared
910  // nodes) into the EqnCommMgr object.
911  // - initMultCRStructure() and initPenCRStructure() follow the same general
912  // routine as for the initElemBlockStructure() function.
913  // - EqnCommMgr::exchangeIndices() is called, which sends all shared
914  // contributions to the processors that own the corresponding equations,
915  // after which the receiving processors insert those contributions into
916  // the local matrix structure.
917  //
918  //11. initializeBlkEqnMapper() run all nodes, elem-dof, multiplier constraints
919  // and pass global point-equation-numbers with corresponding block-equation
920  // numbers to the snl_fei::PointBlockMap object which maintains a mapping
921  // between the point-entry and block-entry equation spaces.
922  //
923  //12. EqnCommMgr::exchangePtToBlkInfo() ensures that the point-to-block
924  // mapping data is globally consistent across all processors.
925  //
926  //13. The sysBlkMatIndices_ array is allocated, which is the array of
927  // arrays that will hold the block-entry matrix structure. This block-entry
928  // matrix structure is not filled, however, until later if/when the
929  // method getMatrixStructure is called specifically requesting block-entry
930  // structure data.
931  //
932  //14. Finally, if the debugOutputLevel_ is greater than 0, a file is written
933  // which contains a mapping from equations to nodes. This file is often very
934  // useful for post-mortem debugging operations, and is also necessary if the
935  // matrix produced from a serial run is to be compared with a matrix
936  // produced from a similar run on multiple processors. The equation-to-node
937  // mapping is the only way to reconcile the equation-orderings between the
938  // two matrices.
939  //
940 
941  if (debugOutput_) {
942  FEI_OSTREAM& os = dbgOut();
943  os << "FEI: initComplete" << FEI_ENDL;
944  }
945  // marshall all of the nodes we received in element connectivities into an
946  // active node list, and set the fields, blocks and owning processors
947  // associated with those nodes.
948 
949  //IMPORTANT!!! Note that at this point, the nodeIDs in the
950  //connectivitiy tables are being replaced with indices from the nodeDatabase.
951 
952  CHK_ERR( finalizeActiveNodes() );
953 
954  CHK_ERR( finalizeNodeCommMgr() );
955 
956  numLocalElemDOF_ = calcTotalNumElemDOF();
957  numLocalMultCRs_ = calcNumMultCREqns();
958 
959  // ok, now the active equation calculations -- we need to count how many
960  // local equations there are, then obtain the global equation info, then
961  // set equation numbers on all of our active nodes.
962 
963  int numLocalNodes = nodeDatabase_->countLocalNodeDescriptors(localProc_);
964  numLocalNodalEqns_ = nodeDatabase_->countLocalNodalEqns(localProc_);
965 
966  numLocalEqns_ = numLocalNodalEqns_ + numLocalElemDOF_ + numLocalMultCRs_;
967  numLocalEqnBlks_ = numLocalNodes + numLocalElemDOF_ + numLocalMultCRs_;
968 
969  if (debugOutput_) {
970  FEI_OSTREAM& os = dbgOut();
971  os << "# numMultCRs: " << numLocalMultCRs_ << ", numElemDOF: "
972  << numLocalElemDOF_ << ", numLocalNodalEqns: " <<numLocalNodalEqns_
973  << FEI_ENDL;
974  os << "# numLocalEqns_: " << numLocalEqns_ << FEI_ENDL;
975  }
976 
977  calcGlobalEqnInfo(numLocalNodes, numLocalEqns_, numLocalEqnBlks_);
978 
979  CHK_ERR( setNodalEqnInfo() );
980 
981  //now we need to set equation numbers for the element-DOF's, and for the
982  //lagrange-multipler constraints.
983  //(exactly where to put these element DOF is an interesting issue... here,
984  //just toss them at the end of the nodal active eqns, which may or may not
985  //be such a great choice.)
986 
987  setElemDOFEqnInfo();
988 
989  CHK_ERR( setMultCREqnInfo() );
990 
991  //now have the NodeDatabase run through the list of active nodes and
992  //do whatever final initializations it needs to do. This function also
993  //calls nodeCommMgr::exchangeEqnInfo.
994  CHK_ERR( nodeDatabase_->synchronize(firstLocalNodeNumber_, localStartRow_,
995  localProc_, comm_) );
996 
997  //now run the active nodes again and set the number of
998  //active nodes and equations per block. We couldn't do this
999  //before, because the nodeCommMgr::initComplete() and
1000  //nodeCommMgr::exchangeEqnInfo functions may both have found out about new
1001  //blocks that are associated with shared nodes.
1002 
1003  setNumNodesAndEqnsPerBlock();
1004 
1005  //at this point we'll run through any SlaveVariable records that were set,
1006  //and establish an EqnBuffer of slave equation numbers, and gather the slave
1007  //equations initialized on any other processors.
1008 
1009  slvCommMgr_ = new EqnCommMgr(comm_);
1010  slvCommMgr_->setNumRHSs(1);
1011 
1012  CHK_ERR( calculateSlaveEqns(comm_) );
1013 
1014  //From this point on, all computations involving equation numbers will assume
1015  //that those numbers are in the reduced equation space. So we'll start by
1016  //translating the contents of eqnCommMgr_ to the reduced space.
1017  //CHK_ERR( translateToReducedEqns(*eqnCommMgr_) );
1018 
1019  initializeEqnCommMgr();
1020 
1021  translateToReducedEqn(localEndRow_, reducedEndRow_);
1022  int startRow = localStartRow_;
1023  while(isSlaveEqn(startRow)) startRow++;
1024  translateToReducedEqn(startRow, reducedStartRow_);
1025  numLocalReducedRows_ = reducedEndRow_ - reducedStartRow_ + 1;
1026 
1027  generateGraph_ = generateGraph;
1028 
1029  if (generateGraph_) {
1030  sysMatIndices_ = new fei::ctg_set<int>[numLocalReducedRows_];
1031  }
1032 
1033  //Now we'll run through all the data structures that have been initialized
1034  //and form the matrix structure (lists of column-indices).
1035 
1036  CHK_ERR( formMatrixStructure() );
1037 
1038  CHK_ERR( initializeBlkEqnMapper() );
1039 
1040  if (globalMaxBlkSize_ > 1) {
1041  CHK_ERR( eqnCommMgr_->exchangePtToBlkInfo(*blkEqnMapper_) );
1042  if (numSlvs_ > 0) {
1043  try {
1044  CHK_ERR( slvCommMgr_->exchangeIndices() );
1045  }
1046  catch (std::runtime_error& exc) {
1047  fei::console_out() << exc.what() << FEI_ENDL;
1048  ERReturn(-1);
1049  }
1050 
1051  CHK_ERR( slvCommMgr_->exchangePtToBlkInfo(*blkEqnMapper_) );
1052  }
1053  }
1054 
1055  localReducedBlkOffset_ = ptEqnToBlkEqn(reducedStartRow_);
1056  int lastLocalReducedBlkEqn = ptEqnToBlkEqn(reducedEndRow_);
1057  numLocalReducedEqnBlks_ = lastLocalReducedBlkEqn - localReducedBlkOffset_ + 1;
1058 
1059  if (globalMaxBlkSize_ > 1 && generateGraph_) {
1060  sysBlkMatIndices_ = numLocalReducedEqnBlks_>0 ?
1061  new fei::ctg_set<int>[numLocalReducedEqnBlks_] : NULL;
1062  }
1063 
1064  delete slvCommMgr_;
1065 
1066  blockMatrix_ = true;
1067 
1068  structureFinalized_ = true;
1069  matIndicesDestroyed_ = false;
1070 
1071  //while we're here, let's write an equation -> nodeNumber map into a file
1072  //for possible debugging purposes... it will come in handy if we need to
1073  //compare matrices/vectors from different parallel runs.
1074  if (debugOutput_) writeEqn2NodeMap();
1075 
1076  if (debugOutput_) dbgOut() << "#leaving initComplete" << FEI_ENDL;
1077  return(FEI_SUCCESS);
1078 }
1079 
1080 //------------------------------------------------------------------------------
1081 int SNL_FEI_Structure::formMatrixStructure()
1082 {
1083  FEI_OSTREAM& os = dbgOut();
1084  if (debugOutput_) os << "# formMatrixStructure" << FEI_ENDL;
1085 
1086  //do our local equation profile calculations (i.e., determine
1087  //how long each row is). use the element scatter arrays to determine
1088  //the sparsity pattern
1089 
1090  CHK_ERR( initElemBlockStructure() );
1091 
1092  // next, handle the matrix structure imposed by the constraint relations...
1093  //
1094 
1095  CHK_ERR( initMultCRStructure() );
1096 
1097  // we also need to accomodate penalty constraints, so let's process them
1098  // now...
1099 
1100  CHK_ERR( initPenCRStructure() );
1101 
1102  //we've now processed all of the local data that produces matrix structure.
1103  //so now let's have the equation comm manager exchange indices with the
1104  //other processors so we can account for remote contributions to our
1105  //matrix structure.
1106  try {
1107  CHK_ERR( eqnCommMgr_->exchangeIndices(&os) );
1108  }
1109  catch(std::runtime_error& exc) {
1110  fei::console_out() << exc.what() << FEI_ENDL;
1111  ERReturn(-1);
1112  }
1113 
1114  //so now the remote contributions should be available, let's get them out
1115  //of the eqn comm mgr and put them into our local matrix structure.
1116 
1117  int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
1118  std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
1119  std::vector<fei::CSVec*>& recvEqns = eqnCommMgr_->localEqns();
1120  int i;
1121  if (debugOutput_) {
1122  os << "# after eqnCommMgr_->exchangeIndices, numRecvEqns: "
1123  << numRecvEqns << FEI_ENDL;
1124  }
1125 
1126  for(i=0; i<numRecvEqns; i++) {
1127  int eqn = recvEqnNumbers[i];
1128  if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
1129  fei::console_out() << "SNL_FEI_Structure::initComplete: ERROR, recvEqn " << eqn
1130  << " out of range. (reducedStartRow_: " << reducedStartRow_
1131  << ", reducedEndRow_: " << reducedEndRow_ << ", localProc_: "
1132  << localProc_ << ")" << FEI_ENDL;
1133  return(-1);
1134  }
1135 
1136  for(size_t j=0; j<recvEqns[i]->size(); j++) {
1137  CHK_ERR( createMatrixPosition(eqn, recvEqns[i]->indices()[j],
1138  "frmMatStr") );
1139  }
1140  }
1141 
1142  if (debugOutput_) {
1143  os << "# leaving formMatrixStructure" << FEI_ENDL;
1144  }
1145 
1146  return(0);
1147 }
1148 
1149 //------------------------------------------------------------------------------
1150 int SNL_FEI_Structure::initElemBlockStructure()
1151 {
1152  int numBlocks = getNumElemBlocks();
1153 
1154  FEI_OSTREAM& os = dbgOut();
1155  if (debugOutput_) {
1156  os << "# initElemBlockStructure" << FEI_ENDL;
1157  os << "# numElemBlocks: " << numBlocks << FEI_ENDL;
1158  }
1159 
1160  for (int bIndex = 0; bIndex < numBlocks; bIndex++) {
1161  BlockDescriptor* block = NULL;
1162  CHK_ERR( getBlockDescriptor_index(bIndex, block) );
1163 
1164  int numEqns = block->getNumEqnsPerElement();
1165  int interleave = block->getInterleaveStrategy();
1166  std::vector<int> scatterIndices(numEqns);
1167 
1168  GlobalID elemBlockID = block->getGlobalBlockID();
1169  ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
1170 
1171  std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
1172  block->setNumElements(elemIDList.size());
1173  int numBlockElems = block->getNumElements();
1174 
1175  //loop over all the elements, determining the elemental (both from nodes
1176  //and from element DOF) contributions to the sparse matrix structure
1177  if (debugOutput_) {
1178  os << "# block " << bIndex << ", numElems: " << numBlockElems<<FEI_ENDL;
1179  }
1180  for(int elemIndex = 0; elemIndex < numBlockElems; elemIndex++) {
1181 
1182  getScatterIndices_index(bIndex, elemIndex, interleave,
1183  &scatterIndices[0]);
1184 
1185  //now, store the structure that will arise from this contribution,
1186  //after first performing any manipulations associated with slave
1187  //reduction.
1188 
1189  CHK_ERR( createSymmEqnStructure(scatterIndices) );
1190  }
1191  }
1192 
1193  if (reducedEqnCounter_ > 0) CHK_ERR( assembleReducedStructure() );
1194 
1195  if (debugOutput_) os << "# leaving initElemBlockStructure" << FEI_ENDL;
1196  return(FEI_SUCCESS);
1197 }
1198 
1199 //------------------------------------------------------------------------------
1200 int SNL_FEI_Structure::getMatrixRowLengths(std::vector<int>& rowLengths)
1201 {
1202  if (!structureFinalized_) return(-1);
1203 
1204  rowLengths.resize(numLocalReducedRows_);
1205 
1206  for(int i=0; i<numLocalReducedRows_; i++) {
1207  rowLengths[i] = sysMatIndices_[i].size();
1208  }
1209  return(0);
1210 }
1211 
1212 //------------------------------------------------------------------------------
1213 int SNL_FEI_Structure::getMatrixStructure(int** indices,
1214  std::vector<int>& rowLengths)
1215 {
1216  if (!structureFinalized_) return(-1);
1217 
1218  rowLengths.resize(numLocalReducedRows_);
1219 
1220  for(int i=0; i<numLocalReducedRows_; i++) {
1221  rowLengths[i] = sysMatIndices_[i].size();
1222  fei::copySetToArray(sysMatIndices_[i], rowLengths[i], indices[i]);
1223  }
1224 
1225  return(0);
1226 }
1227 
1228 //------------------------------------------------------------------------------
1229 int SNL_FEI_Structure::getMatrixStructure(int** ptColIndices,
1230  std::vector<int>& ptRowLengths,
1231  int** blkColIndices,
1232  int* blkIndices_1D,
1233  std::vector<int>& blkRowLengths,
1234  std::vector<int>& numPtRowsPerBlkRow)
1235 {
1236  int err = getMatrixStructure(ptColIndices, ptRowLengths);
1237  if (err != 0) return(-1);
1238 
1239  if (globalMaxBlkSize_ == 1) {
1240  //No block-equations have more than one point-equation, so we'll just assign
1241  //the block-structure arrays to be the same as the point-structure arrays.
1242  int numRows = ptRowLengths.size();
1243  blkRowLengths.resize(numRows);
1244  numPtRowsPerBlkRow.assign(numRows, 1);
1245 
1246  blkRowLengths.resize(ptRowLengths.size());
1247  for(size_t ii=0; ii<ptRowLengths.size(); ++ii) {
1248  blkRowLengths[ii] = ptRowLengths[ii];
1249  }
1250 
1251  for(int i=0; i<numRows; i++) {
1252  blkColIndices[i] = ptColIndices[i];
1253  }
1254  }
1255  else {
1256  //There are block-equations with more than 1 point-equation, so let's
1257  //calculate the actual block-structure.
1258 
1259  std::map<int,int>* ptEqns = blkEqnMapper_->getPtEqns();
1260  int numPtEqns = ptEqns->size();
1261 
1262  std::map<int,int>::const_iterator
1263  pteq = ptEqns->begin();
1264 
1265  int lastBlkRow = -1;
1266 
1267  for(int jj=0; jj<numPtEqns; ++jj, ++pteq) {
1268  int ptEqn = (*pteq).first;
1269  int localPtEqn = ptEqn - reducedStartRow_;
1270  if (localPtEqn < 0 || localPtEqn >= numLocalReducedRows_) continue;
1271 
1272  int rowLength = sysMatIndices_[localPtEqn].size();
1273 
1274  int blkRow = blkEqnMapper_->eqnToBlkEqn(ptEqn);
1275  if (blkRow < 0) {
1276  ERReturn(-1);
1277  }
1278 
1279  if (blkRow == lastBlkRow) continue;
1280 
1281  int localBlkRow = blkRow - localReducedBlkOffset_;
1282  if (localBlkRow < 0 || localBlkRow >= numLocalReducedEqnBlks_) {
1283  ERReturn(-1);
1284  }
1285 
1286  fei::ctg_set<int>& sysBlkIndices = sysBlkMatIndices_[localBlkRow];
1287 
1288  int* theseColIndices = ptColIndices[localPtEqn];
1289 
1290  for(int colj=0; colj<rowLength; colj++) {
1291  int blkCol = blkEqnMapper_->eqnToBlkEqn(theseColIndices[colj]);
1292  if (blkCol < 0) {
1293  fei::console_out() << localProc_
1294  <<"SNL_FEI_Structure::getMatrixStructure ERROR pt row "
1295  << ptEqn << ", pt col "
1296  << ptColIndices[localPtEqn][colj]
1297  << " doesn't have a corresponding block" << FEI_ENDL;
1298  blkCol = blkEqnMapper_->eqnToBlkEqn(theseColIndices[colj]);
1299  std::abort();
1300  }
1301 
1302  sysBlkIndices.insert2(blkCol);
1303  }
1304 
1305  lastBlkRow = blkRow;
1306  }
1307 
1308  //now we're ready to set the block-sizes...
1309 
1310  numPtRowsPerBlkRow.resize(numLocalReducedEqnBlks_);
1311  blkRowLengths.resize(numLocalReducedEqnBlks_);
1312  int offset = 0;
1313  for(int i=0; i<numLocalReducedEqnBlks_; i++) {
1314  blkRowLengths[i] = sysBlkMatIndices_[i].size();
1315  fei::copySetToArray(sysBlkMatIndices_[i],blkRowLengths[i],
1316  &(blkIndices_1D[offset]));
1317  offset += blkRowLengths[i];
1318 
1319  int blkEqn = localReducedBlkOffset_ + i;
1320  numPtRowsPerBlkRow[i] = blkEqnMapper_->getBlkEqnSize(blkEqn);
1321  }
1322  }
1323 
1324  return(0);
1325 }
1326 
1327 //------------------------------------------------------------------------------
1328 bool SNL_FEI_Structure::nodalEqnsAllSlaves(const NodeDescriptor* node,
1329  std::vector<int>& slaveEqns)
1330 {
1331  int numFields = node->getNumFields();
1332  const int* fieldIDs = node->getFieldIDList();
1333  const int* fieldEqns= node->getFieldEqnNumbers();
1334 
1335  for(int i=0; i<numFields; ++i) {
1336  int fieldSize = getFieldSize(fieldIDs[i]);
1337  for(int eqn=0; eqn<fieldSize; ++eqn) {
1338  int thisEqn = fieldEqns[i] + eqn;
1339  if (fei::binarySearch(thisEqn, slaveEqns) < 0) {
1340  //found a nodal eqn that's not a slave eqn, so return false.
1341  return(false);
1342  }
1343  }
1344  }
1345 
1346  return(true);
1347 }
1348 
1349 //------------------------------------------------------------------------------
1350 NodeDescriptor& SNL_FEI_Structure::findNodeDescriptor(GlobalID nodeID)
1351 {
1352 //
1353 //This function returns a NodeDescriptor reference if nodeID is an active node.
1354 //
1355  NodeDescriptor* node = NULL;
1356  int err = nodeDatabase_->getNodeWithID(nodeID, node);
1357 
1358  if (err != 0) {
1359  fei::console_out() << "ERROR, findNodeDescriptor unable to find node " << (int)nodeID
1360  << FEI_ENDL;
1361  std::abort();
1362  }
1363 
1364  return( *node );
1365 }
1366 
1367 //------------------------------------------------------------------------------
1368 NodeDescriptor* SNL_FEI_Structure::findNode(GlobalID nodeID)
1369 {
1370  NodeDescriptor* node = NULL;
1371  nodeDatabase_->getNodeWithID(nodeID, node);
1372  return node;
1373 }
1374 
1375 //------------------------------------------------------------------------------
1376 int SNL_FEI_Structure::initMultCRStructure()
1377 {
1378  FEI_OSTREAM& os = dbgOut();
1379  if (debugOutput_) os << "# initMultCRStructure" << FEI_ENDL;
1380  //
1381  //Private SNL_FEI_Structure function, to be called from initComplete.
1382  //
1383  //Records the system matrix structure arising from Lagrange Multiplier
1384  //Constraint Relations.
1385  //
1386  // since at initialization all we are being passed is the
1387  // general form of the constraint equations, we can't check to see if any
1388  // of the weight terms are zeros at this stage of the game. Hence, we
1389  // have to reserve space for all the nodal weight vectors, even though
1390  // they might turn out to be zeros during the load step....
1391 
1392  std::map<GlobalID,ConstraintType*>::const_iterator
1393  cr_iter = multCRs_.begin(),
1394  cr_end = multCRs_.end();
1395 
1396  while(cr_iter != cr_end) {
1397  ConstraintType& multCR = *((*cr_iter).second);
1398 
1399  int lenList = multCR.getMasters().size();
1400 
1401  std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
1402  GlobalID *CRNodePtr = &CRNode_vec[0];
1403  std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
1404  int* CRFieldPtr = &CRField_vec[0];
1405 
1406  int crEqn = multCR.getEqnNumber();
1407  int reducedCrEqn = 0;
1408  translateToReducedEqn(crEqn, reducedCrEqn);
1409 
1410  createMatrixPosition(reducedCrEqn, reducedCrEqn, "initMCRStr");
1411 
1412  for(int j=0; j<lenList; j++) {
1413  GlobalID nodeID = CRNodePtr[j];
1414  int fieldID = CRFieldPtr[j];
1415 
1416  NodeDescriptor *nodePtr = findNode(nodeID);
1417  if(nodePtr==NULL) continue;
1418  NodeDescriptor& node = *nodePtr;
1419 
1420  //first, store the column indices associated with this node, in
1421  //the constraint's equation. i.e., position (crEqn, node)
1422 
1423  storeNodalColumnIndices(crEqn, node, fieldID);
1424 
1425  //now we need to store the transpose of the above contribution,
1426  //i.e., position(s) (node, crEqn)
1427 
1428  if (node.getOwnerProc() == localProc_) {
1429  //if we own this node, we will simply store its equation
1430  //numbers in local rows (equations) of the matrix.
1431 
1432  storeNodalRowIndices(node, fieldID, crEqn);
1433  }
1434  else {
1435  //if we don't own it, then we need to register with the
1436  //eqn comm mgr that we'll be sending contributions to
1437  //column crEqn of the remote equations associated with this
1438  //node.
1439 
1440  storeNodalSendIndex(node, fieldID, crEqn);
1441  }
1442  }
1443  ++cr_iter;
1444  }
1445  return(FEI_SUCCESS);
1446 }
1447 
1448 //------------------------------------------------------------------------------
1449 int SNL_FEI_Structure::initPenCRStructure()
1450 {
1451  FEI_OSTREAM& os = dbgOut();
1452  if (debugOutput_) os << "# initPenCRStructure" << FEI_ENDL;
1453 //
1454 //This function oversees the putting in of any matrix structure that results
1455 //from Penalty constraint relations.
1456 //
1457 // note that penalty constraints don't generate new equations
1458 // (unlike Lagrange constraints), but they do add terms to the system
1459 // stiffness matrix that couple all the nodes that contribute to the
1460 // penalty constraint. In addition, each submatrix is defined by the pair
1461 // of nodes that creates its contribution, hence a submatrix can be defined
1462 // in terms of two weight vectors (of length p and q) instead of the
1463 // generally larger product matrix (of size pq)
1464 
1465 // the additional terms take the form of little submatrices that look a lot
1466 // like an element stiffness and load matrix, where the nodes in the
1467 // constraint list take on the role of the nodes associated with an
1468 // element, and the individual matrix terms arise from the outer products
1469 // of the constraint nodal weight vectors
1470 
1471 // rather than get elegant and treat this task as such an elemental energy
1472 // term, we'll use some brute force to construct these penalty contributions
1473 // on the fly, primarily to simplify -reading- this thing, so that the
1474 // derivations in the annotated implementation document are more readily
1475 // followed...
1476  std::map<GlobalID,ConstraintType*>::const_iterator
1477  cr_iter = penCRs_.begin(),
1478  cr_end = penCRs_.end();
1479 
1480  while(cr_iter != cr_end) {
1481  ConstraintType& penCR = *((*cr_iter).second);
1482 
1483  int lenList = penCR.getMasters().size();
1484  std::vector<GlobalID>& CRNode_vec = penCR.getMasters();
1485  GlobalID* CRNodesPtr = &CRNode_vec[0];
1486 
1487  std::vector<int>& CRField_vec = penCR.getMasterFieldIDs();
1488  int* CRFieldPtr = &CRField_vec[0];
1489 
1490  // each constraint equation generates a set of nodal energy terms, so
1491  // we have to process a matrix of nodes for each constraint
1492 
1493  for(int i = 0; i < lenList; i++) {
1494  GlobalID iNodeID = CRNodesPtr[i];
1495  int iField = CRFieldPtr[i];
1496 
1497  NodeDescriptor* iNodePtr = findNode(iNodeID);
1498  if(!iNodePtr) continue;
1499  NodeDescriptor& iNode = *iNodePtr;
1500 
1501  for(int j = 0; j < lenList; j++) {
1502  GlobalID jNodeID = CRNodesPtr[j];
1503  int jField = CRFieldPtr[j];
1504 
1505  NodeDescriptor* jNodePtr = findNode(jNodeID);
1506  if(!jNodePtr) continue;
1507  NodeDescriptor &jNode = *jNodePtr;
1508 
1509  if (iNode.getOwnerProc() == localProc_) {
1510  //if iNode is local, we'll put equations into the local
1511  //matrix structure.
1512 
1513  storeLocalNodeIndices(iNode, iField, jNode, jField);
1514  }
1515  else {
1516  //if iNode is remotely owned, we'll be registering equations
1517  //to send to its owning processor.
1518 
1519  storeNodalSendIndices(iNode, iField, jNode, jField);
1520  }
1521  }
1522  } // end i loop
1523  ++cr_iter;
1524  } // end while loop
1525  return(FEI_SUCCESS);
1526 }
1527 
1528 //------------------------------------------------------------------------------
1529 void SNL_FEI_Structure::storeNodalSendIndex(NodeDescriptor& node, int fieldID,
1530  int col)
1531 {
1532  //
1533  //This is a private SNL_FEI_Structure function. We can safely assume that it
1534  //will only be called with a node that is not locally owned.
1535  //
1536  //This function tells the eqn comm mgr that we'll be sending contributions
1537  //to column 'col' for the equations associated with 'fieldID', on 'node', on
1538  //node's owning processor.
1539  //
1540  int proc = node.getOwnerProc();
1541 
1542  int eqnNumber = -1;
1543  if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
1544 
1545  int numEqns = getFieldSize(fieldID);
1546  if (numEqns < 1) {
1547  fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
1548  <<") with size "<<numEqns<<FEI_ENDL;
1549  voidERReturn;
1550  }
1551 
1552  for(int i=0; i<numEqns; i++) {
1553  eqnCommMgr_->addRemoteIndices(eqnNumber+i, proc, &col, 1);
1554  }
1555 }
1556 
1557 //------------------------------------------------------------------------------
1558 void SNL_FEI_Structure::storeNodalSendIndices(NodeDescriptor& iNode, int iField,
1559  NodeDescriptor& jNode, int jField){
1560 //
1561 //This function will register with the eqn comm mgr the equations associated
1562 //with iNode, field 'iField' having column indices that are the equations
1563 //associated with jNode, field 'jField', to be sent to the owner of iNode.
1564 //
1565  int proc = iNode.getOwnerProc();
1566 
1567  int iEqn = -1, jEqn = -1;
1568  if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
1569  if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
1570 
1571  int iNumParams = getFieldSize(iField);
1572  int jNumParams = getFieldSize(jField);
1573  if (iNumParams < 1 || jNumParams < 1) {
1574  fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
1575  << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
1576  << jNumParams<<FEI_ENDL;
1577  voidERReturn;
1578  }
1579 
1580  for(int i=0; i<iNumParams; i++) {
1581  int eqn = iEqn + i;
1582 
1583  for(int j=0; j<jNumParams; j++) {
1584  int col = jEqn + j;
1585  eqnCommMgr_->addRemoteIndices(eqn, proc, &col, 1);
1586  }
1587  }
1588 }
1589 
1590 //------------------------------------------------------------------------------
1591 void SNL_FEI_Structure::storeLocalNodeIndices(NodeDescriptor& iNode, int iField,
1592  NodeDescriptor& jNode, int jField)
1593 {
1594 //
1595 //This function will add to the local matrix structure the equations associated
1596 //with iNode at iField, having column indices that are the equations associated
1597 //with jNode at jField.
1598 //
1599  int iEqn = -1, jEqn = -1;
1600  if (!iNode.getFieldEqnNumber(iField, iEqn)) return; //voidERReturn;
1601  if (!jNode.getFieldEqnNumber(jField, jEqn)) return; //voidERReturn;
1602 
1603  int iNumParams = getFieldSize(iField);
1604  int jNumParams = getFieldSize(jField);
1605  if (iNumParams < 1 || jNumParams < 1) {
1606  fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
1607  << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
1608  << jNumParams<<FEI_ENDL;
1609  voidERReturn;
1610  }
1611 
1612  for(int i=0; i<iNumParams; i++) {
1613  int row = iEqn + i;
1614  int reducedRow = -1;
1615  bool isSlave = translateToReducedEqn(row, reducedRow);
1616  if (isSlave) continue;
1617 
1618  for(int j=0; j<jNumParams; j++) {
1619  int col = jEqn + j;
1620  int reducedCol = -1;
1621  isSlave = translateToReducedEqn(col, reducedCol);
1622  if (isSlave) continue;
1623 
1624  createMatrixPosition(reducedRow, reducedCol,
1625  "storeLocNdInd");
1626  }
1627  }
1628 }
1629 
1630 //------------------------------------------------------------------------------
1631 void SNL_FEI_Structure::storeNodalColumnIndices(int eqn, NodeDescriptor& node,
1632  int fieldID)
1633 {
1634  //
1635  //This function stores the equation numbers associated with 'node' at
1636  //'fieldID' as column indices in row 'eqn' of the system matrix structure.
1637  //
1638  if ((localStartRow_ > eqn) || (eqn > localEndRow_)) {
1639  return;
1640  }
1641 
1642  int colNumber = -1;
1643  if (!node.getFieldEqnNumber(fieldID, colNumber)) voidERReturn;
1644 
1645  int numParams = getFieldSize(fieldID);
1646  if (numParams < 1) {
1647  fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
1648  <<") with size "<<numParams<<FEI_ENDL;
1649  voidERReturn;
1650  }
1651 
1652  int reducedEqn = -1;
1653  bool isSlave = translateToReducedEqn(eqn, reducedEqn);
1654  if (isSlave) return;
1655 
1656  for(int j=0; j<numParams; j++) {
1657  int reducedCol = -1;
1658  isSlave = translateToReducedEqn(colNumber+j, reducedCol);
1659  if (isSlave) continue;
1660 
1661  createMatrixPosition(reducedEqn, reducedCol,
1662  "storeNdColInd");
1663  }
1664 }
1665 
1666 //------------------------------------------------------------------------------
1667 void SNL_FEI_Structure::storeNodalRowIndices(NodeDescriptor& node,
1668  int fieldID, int eqn)
1669 {
1670  //
1671  //This function stores column 'eqn' in equation numbers associated with
1672  //'node' at 'fieldID' in the system matrix structure.
1673  //
1674  int eqnNumber = -1;
1675  if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
1676 
1677  int numParams = getFieldSize(fieldID);
1678  if (numParams < 1) {
1679  fei::console_out() << "FEI error, attempt to store indices for field ("<<fieldID
1680  <<") with size "<<numParams<<FEI_ENDL;
1681  voidERReturn;
1682  }
1683 
1684  int reducedEqn = -1;
1685  bool isSlave = translateToReducedEqn(eqn, reducedEqn);
1686  if (isSlave) return;
1687 
1688  for(int j=0; j<numParams; j++) {
1689  int reducedRow = -1;
1690  isSlave = translateToReducedEqn(eqnNumber+j, reducedRow);
1691  if (isSlave) continue;
1692 
1693  createMatrixPosition(reducedRow, reducedEqn, "storeNdRowInd");
1694  }
1695 }
1696 
1697 //------------------------------------------------------------------------------
1698 int SNL_FEI_Structure::createSymmEqnStructure(std::vector<int>& scatterIndices)
1699 {
1700  //scatterIndices are both the rows and the columns for a structurally
1701  //symmetric 2-dimensional contribution to the matrix.
1702  //
1703 
1704  //if there aren't any slaves in the whole problem, store the scatter-indices
1705  //using a fast function (no translations-to-reduced-space) and return.
1706  if (numSlaveEquations() == 0) {
1707  CHK_ERR( storeElementScatterIndices_noSlaves(scatterIndices) );
1708  return(FEI_SUCCESS);
1709  }
1710 
1711  try {
1712 
1713  int len = scatterIndices.size();
1714  bool anySlaves = false;
1715  rSlave_.resize(len);
1716  for(int is=0; is<len; is++) {
1717  rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
1718  if (rSlave_[is] == 1) anySlaves = true;
1719  }
1720 
1721  //if there aren't any slaves in this contribution, then just store it
1722  //and return
1723  if (!anySlaves) {
1724  CHK_ERR( storeElementScatterIndices(scatterIndices) );
1725  return(FEI_SUCCESS);
1726  }
1727 
1728  int* scatterPtr = &scatterIndices[0];
1729 
1730  workSpace_.resize(len);
1731  for(int j=0; j<len; j++) {
1732  translateToReducedEqn(scatterPtr[j], workSpace_[j]);
1733  }
1734 
1735  for(int i=0; i<len; i++) {
1736  int row = scatterPtr[i];
1737  if (rSlave_[i] == 1) {
1738  reducedEqnCounter_++;
1739  //
1740  //'row' is a slave equation, so add this row to Kdi_. But as we do that,
1741  //watch for columns that are slave equations and add them to Kid_.
1742  //
1743  for(int jj=0; jj<len; jj++) {
1744  int col = scatterPtr[jj];
1745 
1746  if (rSlave_[jj] == 1) {
1747  //'col' is a slave equation, so add this column to Kid_.
1748  for(int ii=0; ii<len; ii++) {
1749  int rowi = scatterPtr[ii];
1750 
1751  //only add the non-slave rows for this column.
1752  if (rSlave_[ii] == 1) continue;
1753 
1754  Kid_->createPosition(rowi, col);
1755  }
1756  //now skip to the next column
1757  continue;
1758  }
1759 
1760  Kdi_->createPosition(row, col);
1761  }
1762 
1763  //finally, add all slave columns in this slave row to Kdd_.
1764  for(int kk=0; kk<len; kk++) {
1765  int colk = scatterPtr[kk];
1766 
1767  if (rSlave_[kk] != 1) continue;
1768 
1769  Kdd_->createPosition(row, colk);
1770  }
1771  }
1772  else {
1773  //this is not a slave row, so we will loop across it, creating matrix
1774  //positions for all non-slave columns in it.
1775  int reducedRow = -1;
1776  translateToReducedEqn(row, reducedRow);
1777 
1778  bool rowIsLocal = true;
1779  int owner = localProc_;
1780  if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
1781  rowIsLocal = false;
1782  owner = getOwnerProcForEqn(row);
1783  }
1784 
1785  for(int j=0; j<len; j++) {
1786  if (rSlave_[j] == 1) continue;
1787 
1788  int reducedCol = workSpace_[j];
1789 
1790  if (rowIsLocal) {
1791  CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
1792  "crtSymmEqnStr") );
1793  }
1794  else {
1795  eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
1796  }
1797  }
1798  }
1799  }
1800 
1801  if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
1802 
1803  }
1804  catch(std::runtime_error& exc) {
1805  fei::console_out() << exc.what() << FEI_ENDL;
1806  ERReturn(-1);
1807  }
1808 
1809  return(FEI_SUCCESS);
1810 }
1811 
1812 //------------------------------------------------------------------------------
1813 int SNL_FEI_Structure::createBlkSymmEqnStructure(std::vector<int>& scatterIndices)
1814 {
1815  //scatterIndices are both the rows and the columns for a structurally
1816  //symmetric 2-dimensional contribution to the matrix.
1817  //
1818 
1819  //if there aren't any slaves in the whole problem, store the scatter-indices
1820  //using a fast function (no translations-to-reduced-space) and return.
1821  if (numSlaveEquations() == 0) {
1822  return( storeElementScatterBlkIndices_noSlaves(scatterIndices) );
1823  }
1824 
1825  try {
1826 
1827  int len = scatterIndices.size();
1828  bool anySlaves = false;
1829  rSlave_.resize(len);
1830  for(int is=0; is<len; is++) {
1831  rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
1832  if (rSlave_[is] == 1) anySlaves = true;
1833  }
1834 
1835  //if there aren't any slaves in this contribution, then just store it
1836  //and return
1837  if (!anySlaves) {
1838  CHK_ERR( storeElementScatterIndices(scatterIndices) );
1839  return(FEI_SUCCESS);
1840  }
1841 
1842  int* scatterPtr = &scatterIndices[0];
1843 
1844  workSpace_.resize(len);
1845  for(int j=0; j<len; j++) {
1846  translateToReducedEqn(scatterPtr[j], workSpace_[j]);
1847  }
1848 
1849  for(int i=0; i<len; i++) {
1850  int row = scatterPtr[i];
1851  if (rSlave_[i] == 1) {
1852  reducedEqnCounter_++;
1853  //
1854  //'row' is a slave equation, so add this row to Kdi_. But as we do that,
1855  //watch for columns that are slave equations and add them to Kid_.
1856  //
1857  for(int jj=0; jj<len; jj++) {
1858  int col = scatterPtr[jj];
1859 
1860  if (rSlave_[jj] == 1) {
1861  //'col' is a slave equation, so add this column to Kid_.
1862  for(int ii=0; ii<len; ii++) {
1863  int rowi = scatterPtr[ii];
1864 
1865  //only add the non-slave rows for this column.
1866  if (rSlave_[ii] == 1) continue;
1867 
1868  Kid_->createPosition(rowi, col);
1869  }
1870  //now skip to the next column
1871  continue;
1872  }
1873 
1874  Kdi_->createPosition(row, col);
1875  }
1876 
1877  //finally, add all slave columns in this slave row to Kdd_.
1878  for(int kk=0; kk<len; kk++) {
1879  int colk = scatterPtr[kk];
1880 
1881  if (rSlave_[kk] != 1) continue;
1882 
1883  Kdd_->createPosition(row, colk);
1884  }
1885  }
1886  else {
1887  //this is not a slave row, so we will loop across it, creating matrix
1888  //positions for all non-slave columns in it.
1889  int reducedRow = -1;
1890  translateToReducedEqn(row, reducedRow);
1891 
1892  bool rowIsLocal = true;
1893  int owner = localProc_;
1894  if (reducedStartRow_ > reducedRow || reducedEndRow_ < reducedRow) {
1895  rowIsLocal = false;
1896  owner = getOwnerProcForEqn(row);
1897  }
1898 
1899  for(int j=0; j<len; j++) {
1900  if (rSlave_[j] == 1) continue;
1901 
1902  int reducedCol = workSpace_[j];
1903 
1904  if (rowIsLocal) {
1905  CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
1906  "crtSymmEqnStr") );
1907  }
1908  else {
1909  eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
1910  }
1911  }
1912  }
1913  }
1914 
1915  if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
1916 
1917  }
1918  catch(std::runtime_error& exc) {
1919  fei::console_out() << exc.what() << FEI_ENDL;
1920  ERReturn(-1);
1921  }
1922 
1923  return(FEI_SUCCESS);
1924 }
1925 
1926 //------------------------------------------------------------------------------
1927 int SNL_FEI_Structure::storeElementScatterIndices(std::vector<int>& indices)
1928 {
1929 //This function takes a list of equation numbers, as input, and for each
1930 //one, if it is a local equation number, goes to that row of the sysMatIndices
1931 //structure and stores the whole list of equation numbers as column indices
1932 //in that row of the matrix structure. If the equation number is not local,
1933 //then it is given to the EqnCommMgr.
1934 //
1935  int numIndices = indices.size();
1936  int* indPtr = &indices[0];
1937 
1938  workSpace_.resize(numIndices);
1939  int* workPtr = &workSpace_[0];
1940  int reducedEqn = -1;
1941  for(size_t j=0; j<indices.size(); j++) {
1942  bool isSlave = translateToReducedEqn(indPtr[j], reducedEqn);
1943  if (isSlave) ERReturn(-1);
1944  workPtr[j] = reducedEqn;
1945  }
1946 
1947  for(int i=0; i<numIndices; i++) {
1948  int row = indPtr[i];
1949  int rrow = workPtr[i];
1950 
1951  if ((localStartRow_ > row) || (row > localEndRow_)) {
1952 
1953  int owner = getOwnerProcForEqn(row);
1954  eqnCommMgr_->addRemoteIndices(rrow, owner, workPtr, numIndices);
1955  continue;
1956  }
1957 
1958  CHK_ERR( createMatrixPositions(rrow, numIndices, workPtr,
1959  "storeElemScttrInd") );
1960  }
1961 
1962  return(FEI_SUCCESS);
1963 }
1964 
1965 //------------------------------------------------------------------------------
1966 int SNL_FEI_Structure::storeElementScatterIndices_noSlaves(std::vector<int>& indices)
1967 {
1968  //This function takes a list of equation numbers as input and for each one,
1969  //if it is a local equation number, goes to that row of the sysMatIndices
1970  //structure and stores the whole list of equation numbers as column indices
1971  //in that row of the matrix structure. If the equation number is not local,
1972  //then it is given to the EqnCommMgr.
1973  //
1974  int i, numIndices = indices.size();
1975  int* indPtr = &indices[0];
1976 
1977  for(i=0; i<numIndices; i++) {
1978  int row = indPtr[i];
1979  if (row < localStartRow_ || row > localEndRow_) {
1980  int owner = getOwnerProcForEqn(row);
1981  eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
1982  continue;
1983  }
1984 
1985  if (generateGraph_) {
1986  fei::ctg_set<int>& thisRow =
1987  sysMatIndices_[row - reducedStartRow_];
1988 
1989  for(int j=0; j<numIndices; ++j) {
1990  if (debugOutput_ && outputLevel_ > 2) {
1991  dbgOut() << "# storeElemScttrInds_ns crtMatPoss("
1992  << row << "," << indPtr[j] << ")"<<FEI_ENDL;
1993  }
1994 
1995  thisRow.insert2(indPtr[j]);
1996  }
1997  }
1998  }
1999 
2000  return(FEI_SUCCESS);
2001 }
2002 
2003 //------------------------------------------------------------------------------
2004 int SNL_FEI_Structure::storeElementScatterBlkIndices_noSlaves(std::vector<int>& indices)
2005 {
2006  //This function takes a list of equation numbers as input and for each one,
2007  //if it is a local equation number, goes to that row of the sysMatIndices
2008  //structure and stores the whole list of equation numbers as column indices
2009  //in that row of the matrix structure. If the equation number is not local,
2010  //then it is given to the EqnCommMgr.
2011  //
2012  int i, numIndices = indices.size();
2013  int* indPtr = &indices[0];
2014 
2015  for(i=0; i<numIndices; i++) {
2016  int row = indPtr[i];
2017  if (row < localStartRow_ || row > localEndRow_) {
2018  int owner = getOwnerProcForEqn(row);
2019  eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
2020  continue;
2021  }
2022 
2023  if (generateGraph_) {
2024  fei::ctg_set<int>& thisRow =
2025  sysMatIndices_[row - reducedStartRow_];
2026 
2027  for(int j=0; j<numIndices; ++j) {
2028  if (debugOutput_ && outputLevel_ > 2) {
2029  dbgOut() << "# storeElemScttrInds_ns crtMatPoss("
2030  << row << "," << indPtr[j] << ")"<<FEI_ENDL;
2031  }
2032 
2033  thisRow.insert2(indPtr[j]);
2034  }
2035  }
2036  }
2037 
2038  return(FEI_SUCCESS);
2039 }
2040 
2041 //------------------------------------------------------------------------------
2042 int SNL_FEI_Structure::assembleReducedStructure()
2043 {
2044  //This function performs the appropriate manipulations (matrix-matrix-products,
2045  //matrix-vector-products) on the Kid,Kdi,Kdd matrices and then assembles the
2046  //results into the global matrix structure. This function also adds any
2047  //resulting contributions to off-processor portions of the matrix structure to
2048  //the EqnCommMgr.
2049  //
2050  fei::FillableMat* D = getSlaveDependencies();
2051 
2052  csrD = *D;
2053  csrKid = *Kid_;
2054  csrKdi = *Kdi_;
2055  csrKdd = *Kdd_;
2056 
2057  //form tmpMat1_ = Kid_*D
2058  fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
2059 
2060  //form tmpMat2_ = D^T*Kdi_
2061  fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
2062 
2063  CHK_ERR( translateMatToReducedEqns(tmpMat1_) );
2064  CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
2065  if (numProcs_ > 1) {
2066  CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat1_, true) );
2067  CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
2068  }
2069 
2070  CHK_ERR( createMatrixPositions(tmpMat1_) );
2071  CHK_ERR( createMatrixPositions(tmpMat2_) );
2072 
2073  //form tmpMat1_ = D^T*Kdd_
2074  fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
2075 
2076  //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
2077  fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
2078 
2079 
2080  CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
2081  if (numProcs_ > 1) {
2082  CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
2083  }
2084  CHK_ERR( createMatrixPositions(tmpMat2_) );
2085 
2086  Kdi_->clear();
2087  Kid_->clear();
2088  Kdd_->clear();
2089  reducedEqnCounter_ = 0;
2090 
2091  return(FEI_SUCCESS);
2092 }
2093 
2094 //------------------------------------------------------------------------------
2096 {
2097  if (numSlvs_ < 1) return(0);
2098 
2099  CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvProcEqns())) );
2100  CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getRecvEqns())) );
2101  CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendProcEqns())) );
2102  CHK_ERR( translateToReducedEqns(*(eqnCommMgr.getSendEqns())) );
2103 
2104  return(0);
2105 }
2106 
2107 //------------------------------------------------------------------------------
2109 {
2110  int numEqns = eqnBuf.getNumEqns();
2111  int* eqnNumbers = &(eqnBuf.eqnNumbers()[0]);
2112  std::vector<fei::CSVec*>& eqnArray = eqnBuf.eqns();
2113  for(int i=0; i<numEqns; ++i) {
2114  int reducedEqn;
2115  translateToReducedEqn(eqnNumbers[i], reducedEqn);
2116  eqnNumbers[i] = reducedEqn;
2117 
2118  int* indicesPtr = &(eqnArray[i]->indices()[0]);
2119  int numIndices = eqnArray[i]->size();
2120  for(int j=0; j<numIndices; ++j) {
2121  translateToReducedEqn(indicesPtr[j], reducedEqn);
2122  indicesPtr[j] = reducedEqn;
2123  }
2124  }
2125 
2126  return(0);
2127 }
2128 
2129 //------------------------------------------------------------------------------
2131 {
2132  std::vector<std::vector<int>*>& eqnTable = procEqns.procEqnNumbersPtr();
2133  int dim1 = eqnTable.size();
2134  for(int i=0; i<dim1; ++i) {
2135  int dim2 = eqnTable[i]->size();
2136  int* indicesPtr = &(*(eqnTable[i]))[0];
2137  for(int j=0; j<dim2; ++j) {
2138  int reducedEqn;
2139  translateToReducedEqn(indicesPtr[j], reducedEqn);
2140  indicesPtr[j] = reducedEqn;
2141  }
2142  }
2143 
2144  return(0);
2145 }
2146 
2147 //------------------------------------------------------------------------------
2149 {
2150  //Given a matrix in unreduced numbering, convert all of its contents to the
2151  //"reduced" equation space. If any of the row-numbers or column-indices in
2152  //this matrix object are slave equations, they will not be referenced. In
2153  //this case, a positive warning-code will be returned.
2154 
2155  int reducedEqn = -1;
2156  int foundSlave = 0;
2157 
2158  fei::SparseRowGraph& srg = mat.getGraph();
2159 
2160  std::vector<int>& rowNumbers = srg.rowNumbers;
2161  for(size_t i=0; i<rowNumbers.size(); ++i) {
2162  bool isSlave = translateToReducedEqn(rowNumbers[i], reducedEqn);
2163  if (isSlave) foundSlave = 1;
2164  else rowNumbers[i] = reducedEqn;
2165  }
2166 
2167  std::vector<int>& colIndices = srg.packedColumnIndices;
2168  for(size_t i=0; i<colIndices.size(); ++i) {
2169  bool isSlave = translateToReducedEqn(colIndices[i], reducedEqn);
2170  if (isSlave) foundSlave = 1;
2171  else colIndices[i] = reducedEqn;
2172  }
2173 
2174  return(foundSlave);
2175 }
2176 
2177 //------------------------------------------------------------------------------
2178 int SNL_FEI_Structure::createMatrixPositions(fei::CSRMat& mat)
2179 {
2180  if (!generateGraph_) {
2181  return(0);
2182  }
2183 
2184  //This function must be called with a CSRMat object that already has its
2185  //contents numbered in "reduced" equation-numbers.
2186  //
2187  //This function has one simple task -- for each row,col pair stored in 'mat',
2188  //call 'createMatrixPosition' to make an entry in the global matrix structure
2189  //if it doesn't already exist.
2190  //
2191  const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
2192  const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
2193  const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
2194 
2195  for(size_t i=0; i<rowNumbers.size(); ++i) {
2196  int row = rowNumbers[i];
2197  int offset = rowOffsets[i];
2198  int rowlen = rowOffsets[i+1]-offset;
2199  const int* indices = &pckColInds[offset];
2200 
2201  for(int j=0; j<rowlen; j++) {
2202  CHK_ERR( createMatrixPosition(row, indices[j], "crtMatPos(CSRMat)") );
2203  }
2204  }
2205 
2206  return(0);
2207 }
2208 
2209 //------------------------------------------------------------------------------
2210 int SNL_FEI_Structure::createMatrixPosition(int row, int col,
2211  const char* callingFunction)
2212 {
2213  if (!generateGraph_) {
2214  return(0);
2215  }
2216 
2217  //
2218  //This function inserts the column index 'col' in the system matrix structure
2219  //in 'row', if it isn't already there.
2220  //
2221  //This is a private SNL_FEI_Structure function. These row/col numbers must be
2222  //global, 0-based equation numbers in the "reduced" equation space..
2223  //
2224 
2225  //if the equation isn't local, just return. (The assumption being that
2226  //this call is also being made on the processor that owns the equation.)
2227  if (row < reducedStartRow_ || row > reducedEndRow_) {
2228  if (debugOutput_) {
2229  dbgOut() << " " << callingFunction << " passed " <<row<<","<<col
2230  << ", not local." << FEI_ENDL;
2231  }
2232 
2233  return(0);
2234  }
2235 
2236  if (debugOutput_ && outputLevel_ > 2) {
2237  dbgOut() << "# " << callingFunction << " crtMatPos("
2238  << row << "," << col << ")"<<FEI_ENDL;
2239  }
2240 
2241  fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
2242 
2243  thisRow.insert2(col);
2244 
2245  return(0);
2246 }
2247 
2248 //------------------------------------------------------------------------------
2249 int SNL_FEI_Structure::createMatrixPositions(int row, int numCols, int* cols,
2250  const char* callingFunction)
2251 {
2252  if (!generateGraph_) {
2253  return(0);
2254  }
2255 
2256  //
2257  //This function inserts the column indices 'cols' in the system matrix structure
2258  //in 'row', if they aren't already there.
2259  //
2260  //This is a private SNL_FEI_Structure function. These row/col numbers must be
2261  //global, 0-based equation numbers in the "reduced" equation space..
2262  //
2263 
2264  //if the equation isn't local, just return. (The assumption being that
2265  //this call is also being made on the processor that owns the equation.)
2266  if (row < reducedStartRow_ || row > reducedEndRow_) {
2267  if (debugOutput_) {
2268  dbgOut() << " " << callingFunction << " passed " <<row
2269  << ", not local." << FEI_ENDL;
2270  }
2271 
2272  return(0);
2273  }
2274 
2275  fei::ctg_set<int>& thisRow = sysMatIndices_[row - reducedStartRow_];
2276 
2277  for(int i=0; i<numCols; i++) {
2278  if (debugOutput_ && outputLevel_ > 2) {
2279  dbgOut() << "# " << callingFunction << " crtMatPoss("
2280  << row << "," << cols[i] << ")"<<FEI_ENDL;
2281  }
2282 
2283  thisRow.insert2(cols[i]);
2284  }
2285 
2286  return(0);
2287 }
2288 
2289 //------------------------------------------------------------------------------
2290 int SNL_FEI_Structure::initializeBlkEqnMapper()
2291 {
2292  //Run the nodes, elem-dof's and multiplier constraints, establishing the
2293  //point-equation to block-equation mappings. Note that this data is
2294  //initialized as non-reduced (i.e., not accounting for any slave equations)
2295  //equation numbers, then translated to the reduced space at the end of this
2296  //function.
2297 
2298  if (blkEqnMapper_ != NULL) delete blkEqnMapper_;
2299  blkEqnMapper_ = new snl_fei::PointBlockMap();
2300  snl_fei::PointBlockMap& blkEqnMapper = getBlkEqnMapper();
2301 
2302  //First, if there's only one (non-negative) fieldID and that fieldID's size
2303  //is 1, then we don't need to execute the rest of this function.
2304  int numFields = fieldDatabase_->size();
2305  const int* fieldIDs = getFieldIDsPtr();
2306  const int* fieldSizes = fieldIDs+numFields;
2307  int numNonNegativeFields = 0;
2308  int maxFieldSize = 0;
2309  for(int f=0; f<numFields; f++) {
2310  if (fieldIDs[f] >= 0) {
2311  numNonNegativeFields++;
2312  if (fieldSizes[f] > maxFieldSize) maxFieldSize = fieldSizes[f];
2313  }
2314  }
2315 
2316  if (numNonNegativeFields == 1 && maxFieldSize == 1) {
2317  globalMaxBlkSize_ = 1;
2318  blkEqnMapper.setPtEqualBlk();
2319  return(0);
2320  }
2321 
2322  int localVanishedNodeAdjustment = 0;
2323  for(int i=0; i<localProc_; ++i) {
2324  localVanishedNodeAdjustment += globalNumNodesVanished_[i];
2325  }
2326 
2327  int eqnNumber, blkEqnNumber;
2328 
2329  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
2330  std::map<GlobalID,int>::iterator
2331  iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
2332 
2333  for(; iter!=iter_end; ++iter) {
2334  NodeDescriptor* node = NULL;
2335  nodeDatabase_->getNodeAtIndex(iter->second, node);
2336 
2337  // If the node doesn't exist, skip.
2338  if (node==NULL) continue;
2339 
2340  // If the node doesn't have any fields, skip
2341  int numFields = node->getNumFields();
2342  if(numFields == 0) continue;
2343 
2344  int firstPtEqn = node->getFieldEqnNumbers()[0];
2345  int numNodalDOF = node->getNumNodalDOF();
2346  blkEqnNumber = node->getBlkEqnNumber() - localVanishedNodeAdjustment;
2347 
2348  int blkSize = numNodalDOF;
2349 
2350  for(int j=0; j<numNodalDOF; ++j) {
2351  bool isSlave = translateToReducedEqn(firstPtEqn+j, eqnNumber);
2352  if (isSlave) {
2353  --blkSize;
2354  }
2355  else {
2356  CHK_ERR( blkEqnMapper.setEqn(eqnNumber, blkEqnNumber) );
2357  }
2358  }
2359 
2360  if (blkSize > 0) {
2361  CHK_ERR( blkEqnMapper.setBlkEqnSize(blkEqnNumber, blkSize) );
2362  }
2363  else {
2364  ++localVanishedNodeAdjustment;
2365  }
2366  }
2367 
2368  //
2369  //Now the element-dofs...
2370  //
2371  int numBlocks = getNumElemBlocks();
2372  int numLocalNodes = globalNodeOffsets_[localProc_+1]-globalNodeOffsets_[localProc_];
2373  eqnNumber = localStartRow_ + numLocalNodalEqns_;
2374  blkEqnNumber = localBlkOffset_ + numLocalNodes;
2375 
2376  for(int i = 0; i < numBlocks; i++) {
2377  BlockDescriptor* block = NULL;
2378  CHK_ERR( getBlockDescriptor_index(i, block) );
2379 
2380  int numElemDOF = block->getNumElemDOFPerElement();
2381 
2382  if (numElemDOF > 0) {
2383  int numElems = block->getNumElements();
2384 
2385  for(int j=0; j<numElems; j++) {
2386 
2387  for(int ede=0; ede<numElemDOF; ede++) {
2388  blkEqnMapper.setEqn(eqnNumber+ede, blkEqnNumber+ede);
2389  blkEqnMapper.setBlkEqnSize(blkEqnNumber+ede, 1);
2390  }
2391 
2392  eqnNumber += numElemDOF;
2393  blkEqnNumber += numElemDOF;
2394  }
2395  }
2396  }
2397 
2398  //Now we'll set mappings for all local multiplier constraint-relations,
2399  //and also gather the equation numbers for other processors' multipliers
2400  //to record in our snl_fei::PointBlockMap object.
2401  //
2402  std::map<GlobalID,ConstraintType*>::const_iterator
2403  cr_iter = multCRs_.begin(),
2404  cr_end = multCRs_.end();
2405 
2406  std::vector<int> localMultEqns;
2407  localMultEqns.reserve(multCRs_.size()*2);
2408  while(cr_iter != cr_end) {
2409  ConstraintType* multCR = (*cr_iter).second;
2410  eqnNumber = multCR->getEqnNumber();
2411  blkEqnNumber = multCR->getBlkEqnNumber();
2412 
2413  CHK_ERR( blkEqnMapper_->setEqn(eqnNumber, blkEqnNumber) );
2414  CHK_ERR( blkEqnMapper_->setBlkEqnSize(blkEqnNumber, 1) );
2415 
2416  localMultEqns.push_back(eqnNumber);
2417  localMultEqns.push_back(blkEqnNumber);
2418  ++cr_iter;
2419  }
2420 
2421  //Now gather and store the equation-numbers arising from other
2422  //processors' multiplier constraints. (We obviously only need to do this if
2423  //there is more than one processor, and we can skip it if the problem
2424  //only contains 1 scalar equation for each block-equation.)
2425 
2426  int localMaxBlkEqnSize = blkEqnMapper_->getMaxBlkEqnSize();
2427  globalMaxBlkSize_ = 0;
2428  CHK_ERR( fei::GlobalMax(comm_, localMaxBlkEqnSize, globalMaxBlkSize_) );
2429 
2430  blkEqnMapper_->setMaxBlkEqnSize(globalMaxBlkSize_);
2431 
2432  if (globalMaxBlkSize_ == 1) {
2433  //If the globalMax block-size is 1 then tell the blkEqnMapper, and it will
2434  //reclaim the memory it allocated in arrays, etc.
2435  blkEqnMapper_->setPtEqualBlk();
2436  return(0);
2437  }
2438 
2439  if (numProcs_ > 1) {
2440  std::vector<int> recvLengths, globalMultEqns;
2441  CHK_ERR(fei::Allgatherv(comm_, localMultEqns, recvLengths, globalMultEqns));
2442 
2443  int offset = 0;
2444  for(int p=0; p<numProcs_; p++) {
2445  if (p == localProc_) { offset += recvLengths[p]; continue; }
2446 
2447  for(int j=0; j<recvLengths[p]/2; j++) {
2448  blkEqnMapper_->setEqn(globalMultEqns[offset], globalMultEqns[offset+1]);
2449  blkEqnMapper_->setBlkEqnSize(globalMultEqns[offset+1], 1);
2450  offset += 2;
2451  }
2452  }
2453  }
2454 
2455  return(0);
2456 }
2457 
2458 //------------------------------------------------------------------------------
2459 int SNL_FEI_Structure::setMultCREqnInfo()
2460 {
2461  //We'll set equation-numbers on all local multiplier constraint-relations,
2462  //and also gather the equation numbers for other processors' multipliers
2463  //to record in our snl_fei::PointBlockMap object.
2464  //
2465  std::map<GlobalID,ConstraintType*>::const_iterator
2466  cr_iter = multCRs_.begin(),
2467  cr_end = multCRs_.end();
2468 
2469  int eqnNumber = localStartRow_ + numLocalNodalEqns_ + numLocalElemDOF_;
2470  int blkEqnNum = localBlkOffset_ + numLocalEqnBlks_ - numLocalMultCRs_;
2471 
2472  while(cr_iter != cr_end) {
2473  ConstraintType* multCR = (*cr_iter).second;
2474 
2475  multCR->setEqnNumber(eqnNumber);
2476 
2477  multCR->setBlkEqnNumber(blkEqnNum);
2478 
2479  ++eqnNumber;
2480  ++blkEqnNum;
2481  ++cr_iter;
2482  }
2483 
2484  return(0);
2485 }
2486 
2487 //------------------------------------------------------------------------------
2488 int SNL_FEI_Structure::writeEqn2NodeMap()
2489 {
2490  FEI_OSTRINGSTREAM osstr;
2491  osstr << dbgPath_ << "/FEI" << name_ << "_nID_nNum_blkEq_ptEq."
2492  << numProcs_ << "." << localProc_;
2493 
2494  FEI_OFSTREAM e2nFile(osstr.str().c_str());
2495 
2496  if (!e2nFile) {
2497  fei::console_out() << "SNL_FEI_Structure::writeEqn2NodeMap: ERROR, couldn't open file "
2498  << osstr.str() << FEI_ENDL;
2499  ERReturn(-1);
2500  }
2501 
2502  std::map<GlobalID,ConstraintType*>::const_iterator
2503  cr_iter = multCRs_.begin(),
2504  cr_end = multCRs_.end();
2505 
2506  while(cr_iter != cr_end) {
2507  ConstraintType* multCR = (*cr_iter).second;
2508  int eqnNumber = multCR->getEqnNumber();
2509  int blkEqnNumber = multCR->getBlkEqnNumber();
2510 
2511  e2nFile << "-999 -999 " << blkEqnNumber << " " << eqnNumber << FEI_ENDL;
2512  ++cr_iter;
2513  }
2514 
2515  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
2516  std::map<GlobalID,int>::iterator
2517  iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
2518 
2519  for(; iter!=iter_end; ++iter) {
2520  NodeDescriptor* node = NULL;
2521  nodeDatabase_->getNodeAtIndex(iter->second, node);
2522 
2523  if (node==NULL || node->getOwnerProc() != localProc_) {
2524  continue;
2525  }
2526 
2527  const int* fieldIDList = node->getFieldIDList();
2528  const int* eqnNumbers = node->getFieldEqnNumbers();
2529  int nodeNum = node->getNodeNumber();
2530  int blkEqnNumber = node->getBlkEqnNumber();
2531  int numFields = node->getNumFields();
2532 
2533  for(int j=0; j<numFields; j++) {
2534  int numFieldParams = getFieldSize(fieldIDList[j]);
2535  assert( numFieldParams > 0 );
2536 
2537  for(int k=0; k<numFieldParams; k++) {
2538  e2nFile << (int)node->getGlobalNodeID() << " " << nodeNum
2539  << " " << blkEqnNumber << " " << eqnNumbers[j]+k<<FEI_ENDL;
2540  }
2541  }
2542  }
2543 
2544  return(FEI_SUCCESS);
2545 }
2546 
2547 //------------------------------------------------------------------------------
2548 int SNL_FEI_Structure::calcTotalNumElemDOF()
2549 {
2550  int numBlocks = getNumElemBlocks();
2551  int totalNumElemDOF = 0;
2552 
2553  for(int i = 0; i < numBlocks; i++) {
2554  BlockDescriptor* block = NULL;
2555  CHK_ERR( getBlockDescriptor_index(i, block) );
2556  int numElemDOF = block->getNumElemDOFPerElement();
2557  if (numElemDOF > 0) {
2558  int numElems = block->getNumElements();
2559 
2560  totalNumElemDOF += numElems*numElemDOF;
2561  }
2562  }
2563 
2564  return(totalNumElemDOF);
2565 }
2566 
2567 //------------------------------------------------------------------------------
2568 int SNL_FEI_Structure::calcNumMultCREqns()
2569 {
2570  int totalNumMultCRs = 0;
2571  int numMCRecords = getNumMultConstRecords();
2572 
2573  totalNumMultCRs = numMCRecords;
2574 
2575  return( totalNumMultCRs );
2576 }
2577 
2578 //------------------------------------------------------------------------------
2579 void SNL_FEI_Structure::calcGlobalEqnInfo(int numLocallyOwnedNodes,
2580  int numLocalEqns,
2581  int numLocalEqnBlks)
2582 {
2583  FEI_OSTREAM& os = dbgOut();
2584  if (debugOutput_) os << "# entering calcGlobalEqnInfo" << FEI_ENDL;
2585 
2586 //
2587 //This function does the inter-process communication necessary to obtain,
2588 //on each processor, the global number of equations, and the local starting
2589 //and ending equation numbers.
2590 //While we're here, we do the same thing for node-numbers.
2591 //
2592 
2593  //first, get each processor's local number of equations on the master proc.
2594 
2595  std::vector<int> numProcNodes(numProcs_);
2596  std::vector<int> numProcEqns(numProcs_);
2597  std::vector<int> numProcEqnBlks(numProcs_);
2598 
2599  if (numProcs_ > 1) {
2600 #ifndef FEI_SER
2601  int glist[3];
2602  std::vector<int> recvList(3*numProcs_);
2603 
2604  glist[0] = numLocallyOwnedNodes;
2605  glist[1] = numLocalEqns;
2606  glist[2] = numLocalEqnBlks;
2607  if (MPI_Gather(&glist[0], 3, MPI_INT, &recvList[0], 3,
2608  MPI_INT, masterProc_, comm_) != MPI_SUCCESS) {
2609  fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Gather" << FEI_ENDL;
2610  MPI_Abort(comm_, -1);
2611  }
2612 
2613  for(int p=0; p<numProcs_; p++) {
2614  numProcNodes[p] = recvList[p*3];
2615  numProcEqns[p] = recvList[p*3+1];
2616  numProcEqnBlks[p] = recvList[p*3+2];
2617  }
2618 #endif
2619  }
2620  else {
2621  numProcNodes[0] = numLocallyOwnedNodes;
2622  numProcEqns[0] = numLocalEqns;
2623  numProcEqnBlks[0] = numLocalEqnBlks;
2624  }
2625 
2626  //compute offsets for all processors (starting index for local equations)
2627 
2628  globalNodeOffsets_.resize(numProcs_+1);
2629  globalEqnOffsets_.resize(numProcs_+1);
2630  globalBlkEqnOffsets_.resize(numProcs_+1);
2631 
2632  globalNodeOffsets_[0] = 0; // 0-based node-numbers
2633  globalEqnOffsets_[0] = 0; // we're going to start rows & cols at 0 (global)
2634  globalBlkEqnOffsets_[0] = 0;
2635 
2636  if (localProc_ == masterProc_) {
2637  for(int i=1;i<numProcs_;i++) {
2638  globalNodeOffsets_[i] = globalNodeOffsets_[i-1] +numProcNodes[i-1];
2639  globalEqnOffsets_[i] = globalEqnOffsets_[i-1] + numProcEqns[i-1];
2640  globalBlkEqnOffsets_[i] = globalBlkEqnOffsets_[i-1] +
2641  numProcEqnBlks[i-1];
2642  }
2643 
2644  globalNodeOffsets_[numProcs_] = globalNodeOffsets_[numProcs_-1] +
2645  numProcNodes[numProcs_-1];
2646  globalEqnOffsets_[numProcs_] = globalEqnOffsets_[numProcs_-1] +
2647  numProcEqns[numProcs_-1];
2648  globalBlkEqnOffsets_[numProcs_] = globalBlkEqnOffsets_[numProcs_-1] +
2649  numProcEqnBlks[numProcs_-1];
2650  }
2651 
2652  //now, broadcast vector of eqnOffsets to all processors
2653 
2654  if (numProcs_ > 1) {
2655 #ifndef FEI_SER
2656  std::vector<int> blist(3*(numProcs_+1));
2657 
2658  if (localProc_ == masterProc_) {
2659  int offset = 0;
2660  for(int i=0; i<numProcs_+1; i++) {
2661  blist[offset++] = globalNodeOffsets_[i];
2662  blist[offset++] = globalEqnOffsets_[i];
2663  blist[offset++] = globalBlkEqnOffsets_[i];
2664  }
2665  }
2666 
2667  if (MPI_Bcast(&blist[0], 3*(numProcs_+1), MPI_INT,
2668  masterProc_, comm_) != MPI_SUCCESS) {
2669  fei::console_out() << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Bcast" << FEI_ENDL;
2670  MPI_Abort(comm_, -1);
2671  }
2672 
2673  if (localProc_ != masterProc_) {
2674  int offset = 0;
2675  for(int i=0; i<numProcs_+1; i++) {
2676  globalNodeOffsets_[i] = blist[offset++];
2677  globalEqnOffsets_[i] = blist[offset++];
2678  globalBlkEqnOffsets_[i] = blist[offset++];
2679  }
2680  }
2681 #endif
2682  }
2683 
2684  localStartRow_ = globalEqnOffsets_[localProc_];
2685  localEndRow_ = globalEqnOffsets_[localProc_+1]-1;
2686  numGlobalEqns_ = globalEqnOffsets_[numProcs_];
2687 
2688  firstLocalNodeNumber_ = globalNodeOffsets_[localProc_];
2689  numGlobalNodes_ = globalNodeOffsets_[numProcs_];
2690 
2691  localBlkOffset_ = globalBlkEqnOffsets_[localProc_];
2692  numGlobalEqnBlks_ = globalBlkEqnOffsets_[numProcs_];
2693 
2694  if (debugOutput_) {
2695  os << "# localStartRow_: " << localStartRow_ << FEI_ENDL;
2696  os << "# localEndRow_: " << localEndRow_ << FEI_ENDL;
2697  os << "# numGlobalEqns_: " << numGlobalEqns_ << FEI_ENDL;
2698  os << "# numGlobalEqnBlks_: " << numGlobalEqnBlks_ << FEI_ENDL;
2699  os << "# localBlkOffset_: " << localBlkOffset_ << FEI_ENDL;
2700  os << "# firstLocalNodeNumber_: " << firstLocalNodeNumber_ << FEI_ENDL;
2701  size_t n;
2702  os << "# numGlobalNodes_: " << numGlobalNodes_ << FEI_ENDL;
2703  os << "# " << globalNodeOffsets_.size() << " globalNodeOffsets_: ";
2704  for(size_t nn=0; nn<globalNodeOffsets_.size(); nn++)
2705  os << globalNodeOffsets_[nn] << " ";
2706  os << FEI_ENDL;
2707  os << "# " << globalEqnOffsets_.size() << " globalEqnOffsets_: ";
2708  for(n=0; n<globalEqnOffsets_.size(); n++)
2709  os << globalEqnOffsets_[n] << " ";
2710  os << FEI_ENDL;
2711  os << "# " << globalBlkEqnOffsets_.size() << " globalBlkEqnOffsets_: ";
2712  for(n=0; n<globalBlkEqnOffsets_.size(); n++)
2713  os << globalBlkEqnOffsets_[n] << " ";
2714  os << FEI_ENDL;
2715  os << "# leaving calcGlobalEqnInfo" << FEI_ENDL;
2716  }
2717 }
2718 
2719 //------------------------------------------------------------------------------
2720 int SNL_FEI_Structure::setNodalEqnInfo()
2721 {
2722 //
2723 //Private SNL_FEI_Structure function, to be called in initComplete, after
2724 //'calcGlobalEqnInfo()'.
2725 //
2726 //This function sets global equation numbers on the nodes.
2727 //
2728 //The return value is an error-code.
2729 //
2730  int eqn = localStartRow_;
2731 
2732  //localBlkOffset_ is 0-based, and so is blkEqnNumber.
2733  int blkEqnNumber = localBlkOffset_;
2734 
2735  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
2736  std::map<GlobalID,int>::iterator
2737  iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
2738 
2739  for(; iter!=iter_end; ++iter) {
2740  NodeDescriptor* node = NULL;
2741  nodeDatabase_->getNodeAtIndex(iter->second, node);
2742 
2743  if (node==NULL) continue;
2744 
2745  int numFields = node->getNumFields();
2746  const int* fieldIDs = node->getFieldIDList();
2747 
2748  if (node->getOwnerProc() == localProc_) {
2749  int numNodalDOF = 0;
2750 
2751  for(int j=0; j<numFields; j++) {
2752  node->setFieldEqnNumber(fieldIDs[j], eqn);
2753  int fieldSize = getFieldSize(fieldIDs[j]);
2754  eqn += fieldSize;
2755  numNodalDOF += fieldSize;
2756  }
2757 
2758  node->setNumNodalDOF(numNodalDOF);
2759  node->setBlkEqnNumber(blkEqnNumber++);
2760  }
2761  }
2762 
2763  return(0);
2764 }
2765 
2766 //------------------------------------------------------------------------------
2767 void SNL_FEI_Structure::setElemDOFEqnInfo()
2768 {
2769  //
2770  //Private SNL_FEI_Structure function, to be called from initComplete after
2771  //setBlkEqnInfo() has been called.
2772  //
2773  //Sets global equation numbers for all element-DOF.
2774  //
2775  int numBlocks = getNumElemBlocks();
2776 
2777  int eqnNumber = localStartRow_ + numLocalNodalEqns_;
2778 
2779  for(int i = 0; i < numBlocks; i++) {
2780  BlockDescriptor* block = NULL;
2781  int err = getBlockDescriptor_index(i, block);
2782  if (err) voidERReturn;
2783 
2784  int numElemDOF = block->getNumElemDOFPerElement();
2785 
2786  if (numElemDOF > 0) {
2787  std::vector<int>& elemDOFEqns = block->elemDOFEqnNumbers();
2788  std::map<GlobalID,int>& elemIDs = connTables_[i]->elemIDs;
2789 
2790  std::map<GlobalID,int>::const_iterator
2791  iter = elemIDs.begin(),
2792  iter_end = elemIDs.end();
2793 
2794  for(; iter != iter_end; ++iter) {
2795  elemDOFEqns[iter->second] = eqnNumber;
2796 
2797  eqnNumber += numElemDOF;
2798  }
2799  }
2800  }
2801 }
2802 
2803 //------------------------------------------------------------------------------
2804 int SNL_FEI_Structure::addBlock(GlobalID blockID) {
2805 //
2806 //Append a blockID to our (unsorted) list of blockIDs if it isn't
2807 //already present. Also, add a block-descriptor if blockID wasn't
2808 //already present, and a ConnectivityTable to go with it.
2809 //
2810  int insertPoint = -1;
2811  int found = fei::binarySearch(blockID, blockIDs_, insertPoint);
2812 
2813  if (found<0) {
2814  blockIDs_.insert(blockIDs_.begin()+insertPoint, blockID);
2815 
2816  BlockDescriptor* block = new BlockDescriptor;
2817  block->setGlobalBlockID(blockID);
2818 
2819  blocks_.insert(blocks_.begin()+insertPoint, block);
2820 
2821  ConnectivityTable* newConnTable = new ConnectivityTable;
2822  connTables_.insert(connTables_.begin()+insertPoint, newConnTable);
2823  }
2824 
2825  return(FEI_SUCCESS);
2826 }
2827 
2828 //------------------------------------------------------------------------------
2829 int SNL_FEI_Structure::getBlockDescriptor(GlobalID blockID,
2830  BlockDescriptor*& block)
2831 {
2832  int index = fei::binarySearch(blockID, blockIDs_);
2833 
2834  if (index < 0) {
2835  fei::console_out() << "SNL_FEI_Structure::getBlockDescriptor: ERROR, blockID "
2836  << (int)blockID << " not found." << FEI_ENDL;
2837  return(-1);
2838  }
2839 
2840  block = blocks_[index];
2841 
2842  return(FEI_SUCCESS);
2843 }
2844 
2845 //------------------------------------------------------------------------------
2846 int SNL_FEI_Structure::getIndexOfBlock(GlobalID blockID) const
2847 {
2848  int index = fei::binarySearch(blockID, blockIDs_);
2849  return(index);
2850 }
2851 
2852 //------------------------------------------------------------------------------
2854  BlockDescriptor*& block)
2855 {
2856  if (index < 0 || index >= (int)blockIDs_.size()) ERReturn(FEI_FATAL_ERROR);
2857 
2858  block = blocks_[index];
2859 
2860  return(FEI_SUCCESS);
2861 }
2862 
2863 //------------------------------------------------------------------------------
2864 int SNL_FEI_Structure::allocateBlockConnectivity(GlobalID blockID) {
2865 
2866  int index = fei::binarySearch(blockID, blockIDs_);
2867 
2868  if (index < 0) {
2869  fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, blockID "
2870  << (int)blockID << " not found. Aborting." << FEI_ENDL;
2871  MPI_Abort(comm_, -1);
2872  }
2873 
2874  connTables_[index]->numNodesPerElem =
2875  blocks_[index]->getNumNodesPerElement();
2876 
2877  int numRows = blocks_[index]->getNumElements();
2878  int numCols = connTables_[index]->numNodesPerElem;
2879 
2880  if ((numRows <= 0) || (numCols <= 0)) {
2881  fei::console_out() << "SNL_FEI_Structure::allocateConnectivityTable: ERROR, either "
2882  << "numElems or numNodesPerElem not yet set for blockID "
2883  << (int)blockID << ". Aborting." << FEI_ENDL;
2884  MPI_Abort(comm_, -1);
2885  }
2886 
2887  connTables_[index]->elemNumbers.resize(numRows);
2888 
2889  connTables_[index]->elem_conn_ids = new std::vector<GlobalID>(numRows*numCols);
2890 
2891  return(FEI_SUCCESS);
2892 }
2893 
2894 //------------------------------------------------------------------------------
2895 ConnectivityTable& SNL_FEI_Structure::getBlockConnectivity(GlobalID blockID) {
2896 
2897  int index = fei::binarySearch(blockID, blockIDs_);
2898 
2899  if (index < 0) {
2900  fei::console_out() << "SNL_FEI_Structure::getBlockConnectivity: ERROR, blockID "
2901  << (int)blockID << " not found. Aborting." << FEI_ENDL;
2902  MPI_Abort(comm_, -1);
2903  }
2904 
2905  return(*(connTables_[index]));
2906 }
2907 
2908 //------------------------------------------------------------------------------
2909 int SNL_FEI_Structure::finalizeActiveNodes()
2910 {
2911  FEI_OSTREAM& os = dbgOut();
2912  if (debugOutput_) {
2913  os << "# finalizeActiveNodes" << FEI_ENDL;
2914  }
2915  //First, make sure that the shared-node IDs are in the nodeDatabase, since
2916  //they might not have been added yet...
2917 
2918  std::vector<GlobalID>& shNodeIDs = nodeCommMgr_->getSharedNodeIDs();
2919  for(unsigned i=0; i<shNodeIDs.size(); i++) {
2920  CHK_ERR( nodeDatabase_->initNodeID(shNodeIDs[i]) );
2921  }
2922 
2923  if (activeNodesInitialized_) return(0);
2924  //
2925  //Now run through the connectivity tables and set the nodal field and block
2926  //info. Also, replace the nodeID in the element-connectivity lists with an
2927  //index from the NodeDatabase object. This will save us a lot of time when
2928  //looking up nodes later.
2929  //
2930  //One other thing we're going to do is give all elements an element-number.
2931  //This element-number will start at zero locally (meaning that it's not
2932  //globally unique).
2933  //
2934  int elemNumber = 0;
2935  int numBlocks = blockIDs_.size();
2936  for(int bIndex=0; bIndex<numBlocks; bIndex++) {
2937 
2938  GlobalID blockID = blocks_[bIndex]->getGlobalBlockID();
2939  ConnectivityTable& conn = *(connTables_[bIndex]);
2940 
2941  GlobalID* elemConn = NULL;
2942  NodeDescriptor** elemNodeDescPtrs = NULL;
2943 
2944  int numElems = conn.elemIDs.size();
2945  if (numElems > 0) {
2946  elemConn = &((*conn.elem_conn_ids)[0]);
2947  if (!activeNodesInitialized_) {
2948  int elemConnLen = conn.elem_conn_ids->size();
2949  conn.elem_conn_ptrs = new std::vector<NodeDescriptor*>(elemConnLen);
2950  }
2951  elemNodeDescPtrs = &((*conn.elem_conn_ptrs)[0]);
2952  }
2953  int nodesPerElem = conn.numNodesPerElem;
2954 
2955  int* fieldsPerNodePtr = blocks_[bIndex]->fieldsPerNodePtr();
2956  int** fieldIDsTablePtr = blocks_[bIndex]->fieldIDsTablePtr();
2957 
2958  //
2959  //Now we want to go through the connectivity table, and for each node,
2960  //add its fields and this block to the correct NodeDescriptor from the
2961  //nodeDatabase_.
2962  //(Note that the addField and addBlock functions only add the input if
2963  //it isn't already present in that NodeDescriptor.)
2964  //
2965  //We also need to set its owner proc to localProc_ as a default, and
2966  //inform the nodeCommMgr that the node appears in local connectivities.
2967  //
2968 
2969  conn.elemNumbers.resize(numElems);
2970 
2971  int numDistinctFields = blocks_[bIndex]->getNumDistinctFields();
2972  int fieldID = fieldIDsTablePtr[0][0];
2973 
2974  int elem;
2975  for(elem=0; elem<numElems; elem++) {
2976  conn.elemNumbers[elem] = elemNumber++;
2977  GlobalID* elemNodes = &(elemConn[elem*nodesPerElem]);
2978  NodeDescriptor** elemNodePtrs = &(elemNodeDescPtrs[elem*nodesPerElem]);
2979 
2980  for(int n=0; n<nodesPerElem; n++) {
2981  NodeDescriptor* node = NULL;
2982  int index = nodeDatabase_->getNodeWithID( elemNodes[n], node );
2983  if (index < 0) {
2984  fei::console_out() << "ERROR in SNL_FEI_Structure::initializeActiveNodes, "
2985  << FEI_ENDL << "failed to find node "
2986  << (int)(elemNodes[n]) << FEI_ENDL;
2987  }
2988 
2989  if (numDistinctFields == 1) {
2990  node->addField(fieldID);
2991  }
2992  else {
2993  for(int i=0; i<fieldsPerNodePtr[n]; i++) {
2994  node->addField(fieldIDsTablePtr[n][i]);
2995  }
2996  }
2997 
2998  int blk_idx = getIndexOfBlock(blockID);
2999  if (blk_idx >= 0) {
3000  node->addBlockIndex(blk_idx);
3001  }
3002 
3003  //now store this NodeDescriptor pointer, for fast future lookups
3004  elemNodePtrs[n] = node;
3005 
3006  node->setOwnerProc(localProc_);
3007  if (numProcs_ > 1) nodeCommMgr_->informLocal(*node);
3008  }
3009  }
3010  }
3011 
3012  //
3013  //we also need to run through the penalty constraints and inform the
3014  //nodeCommMgr that the constrained nodes are local.
3015  //
3016 
3017  if (numProcs_ > 1) {
3018  std::map<GlobalID,ConstraintType*>::const_iterator
3019  cr_iter = getPenConstRecords().begin(),
3020  cr_end = getPenConstRecords().end();
3021 
3022  while(cr_iter != cr_end) {
3023  ConstraintType& cr = *((*cr_iter).second);
3024  std::vector<GlobalID>& nodeID_vec = cr.getMasters();
3025  GlobalID* nodeIDs = &nodeID_vec[0];
3026  int numNodes = cr.getMasters().size();
3027 
3028  NodeDescriptor* node = NULL;
3029  for(int k=0; k<numNodes; ++k) {
3030  int err = nodeDatabase_->getNodeWithID(nodeIDs[k], node) ;
3031  if ( err != -1 ) {
3032  nodeCommMgr_->informLocal(*node);
3033  }
3034  //CHK_ERR( nodeDatabase_->getNodeWithID(nodeIDs[k], node) );
3035  //CHK_ERR( nodeCommMgr_->informLocal(*node) );
3036  }
3037  ++cr_iter;
3038  }
3039  }
3040 
3041  activeNodesInitialized_ = true;
3042 
3043  if (debugOutput_) os << "# leaving finalizeActiveNodes" << FEI_ENDL;
3044  return(FEI_SUCCESS);
3045 }
3046 
3047 //------------------------------------------------------------------------------
3048 int SNL_FEI_Structure::finalizeNodeCommMgr()
3049 {
3050  //call initComplete() on the nodeCommMgr, so that it can
3051  //finalize shared-node ownerships, etc.
3052 
3053  //request the safetyCheck if the user requested it via the
3054  //parameters function.
3055  bool safetyCheck = checkSharedNodes_;
3056 
3057  if (safetyCheck && localProc_==0 && numProcs_>1 && outputLevel_>0) {
3058  FEI_COUT << "FEI Info: A consistency-check of shared-node data will be "
3059  << "performed, which involves all-to-all communication. This check is "
3060  << "done only if explicitly requested by parameter "
3061  << "'FEI_CHECK_SHARED_IDS'."<<FEI_ENDL;
3062  }
3063 
3064  CHK_ERR( nodeCommMgr_->initComplete(*nodeDatabase_, safetyCheck) );
3065 
3066  if (debugOutput_) {
3067  FEI_OSTREAM& os = dbgOut();
3068  int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
3069  os << "# numSharedNodes: " << numSharedNodes << FEI_ENDL;
3070  for(int ns=0; ns<numSharedNodes; ns++) {
3071  NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(ns);
3072  GlobalID nodeID = node.getGlobalNodeID();
3073  int proc = node.getOwnerProc();
3074  os << "# shNodeID " << (int)nodeID << ", owned by proc "<<proc<<FEI_ENDL;
3075  }
3076  }
3077 
3078  return(0);
3079 }
3080 
3081 //------------------------------------------------------------------------------
3082 int SNL_FEI_Structure::setNumNodesAndEqnsPerBlock()
3083 {
3084  //This function will count how
3085  //many active nodes there are for each block.
3086 
3087  int numBlocks = blockIDs_.size();
3088  std::vector<int> nodesPerBlock(numBlocks);
3089  std::vector<int> eqnsPerBlock(numBlocks);
3090 
3091  int j;
3092  for(j=0; j<numBlocks; j++) {
3093  nodesPerBlock[j] = 0;
3094  eqnsPerBlock[j] = 0;
3095  }
3096 
3097  int numNodes = nodeDatabase_->getNumNodeDescriptors();
3098 
3099  for(j=0; j<numNodes; j++) {
3100  NodeDescriptor* node = NULL;
3101  nodeDatabase_->getNodeAtIndex(j, node);
3102  if (node==NULL) continue;
3103 
3104  const int* fieldIDList = node->getFieldIDList();
3105 
3106  int numFields = node->getNumFields();
3107 
3108  const std::vector<unsigned>& blkIndices = node->getBlockIndexList();
3109  for(std::vector<unsigned>::const_iterator b=blkIndices.begin(), bend=blkIndices.end();
3110  b!=bend; ++b) {
3111  nodesPerBlock[*b]++;
3112 
3113  for(int fld=0; fld<numFields; fld++) {
3114  if (blocks_[*b]->containsField(fieldIDList[fld])) {
3115  int fSize = getFieldSize(fieldIDList[fld]);
3116  assert(fSize >= 0);
3117  eqnsPerBlock[*b] += fSize;
3118  }
3119  }
3120  }
3121  }
3122 
3123  for(j=0; j<numBlocks; j++) {
3124  blocks_[j]->setNumActiveNodes(nodesPerBlock[j]);
3125  }
3126 
3127  //now we need to add the elem-dof to the eqn-count for each block,
3128  //and then set the total number of eqns on each block.
3129 
3130  for(j=0; j<numBlocks; j++) {
3131  eqnsPerBlock[j] += blocks_[j]->getNumElemDOFPerElement() *
3132  blocks_[j]->getNumElements();
3133 
3134  blocks_[j]->setTotalNumEqns(eqnsPerBlock[j]);
3135  }
3136  return(0);
3137 }
3138 
3139 //------------------------------------------------------------------------------
3140 void SNL_FEI_Structure::initializeEqnCommMgr()
3141 {
3142  FEI_OSTREAM& os = dbgOut();
3143  if (debugOutput_) os << "# initializeEqnCommMgr" << FEI_ENDL;
3144 
3145  //This function will take information from the (internal member) nodeCommMgr
3146  //and use it to tell the eqnCommMgr which equations we can expect other
3147  //processors to contribute to, and also which equations we'll be sending to
3148  //other processors.
3149  //
3150  //This function can not be called until after initComplete() has been called
3151  //on the nodeCommMgr.
3152  //
3153  int numSharedNodes = nodeCommMgr_->getNumSharedNodes();
3154 
3155  for(int i=0; i<numSharedNodes; i++) {
3156  NodeDescriptor& node = nodeCommMgr_->getSharedNodeAtIndex(i);
3157 
3158  int numFields = node.getNumFields();
3159  const int* nodeFieldIDList = node.getFieldIDList();
3160  const int* nodeEqnNums = node.getFieldEqnNumbers();
3161 
3162  int owner = node.getOwnerProc();
3163  if (owner == localProc_) {
3164  std::vector<int>& proc = nodeCommMgr_->getSharedNodeProcs(i);
3165  int numProcs = proc.size();
3166 
3167  for(int p=0; p<numProcs; p++) {
3168 
3169  if (proc[p] == localProc_) continue;
3170 
3171  for(int j=0; j<numFields; j++) {
3172  int numEqns = getFieldSize(nodeFieldIDList[j]);
3173  assert(numEqns >= 0);
3174 
3175  int eqn;
3176  for(eqn=0; eqn<numEqns; eqn++) {
3177  int reducedEqn = -1;
3178  bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn,
3179  reducedEqn);
3180  if (isSlave) continue;
3181 
3182  eqnCommMgr_->addLocalEqn(reducedEqn, proc[p]);
3183  }
3184  }
3185  }
3186  }
3187  else {
3188  for(int j=0; j<numFields; j++) {
3189  int numEqns = getFieldSize(nodeFieldIDList[j]);
3190  assert(numEqns >= 0);
3191  for(int eqn=0; eqn<numEqns; eqn++) {
3192  int reducedEqn = -1;
3193  bool isSlave = translateToReducedEqn(nodeEqnNums[j]+eqn, reducedEqn);
3194  if (isSlave) continue;
3195 
3196  eqnCommMgr_->addRemoteIndices(reducedEqn, owner, NULL, 0);
3197  }
3198  }
3199  }
3200  }
3201 
3202  if (debugOutput_) os << "# leaving initializeEqnCommMgr" << FEI_ENDL;
3203 }
3204 
3205 //------------------------------------------------------------------------------
3206 void SNL_FEI_Structure::getEqnInfo(int& numGlobalEqns, int& numLocalEqns,
3207  int& localStartRow, int& localEndRow){
3208 
3209  numGlobalEqns = numGlobalEqns_;
3210  numLocalEqns = numLocalEqns_;
3211  localStartRow = localStartRow_;
3212  localEndRow = localEndRow_;
3213 }
3214 
3215 //------------------------------------------------------------------------------
3216 int SNL_FEI_Structure::getEqnNumbers(GlobalID ID,
3217  int idType, int fieldID,
3218  int& numEqns, int* eqnNumbers)
3219 {
3220  //for now, only allow node ID. allowance for element ID is coming soon!!!
3221  if (idType != FEI_NODE) ERReturn(-1);
3222 
3223  NodeDescriptor* node = NULL;
3224  CHK_ERR( nodeDatabase_->getNodeWithID(ID, node) );
3225 
3226  numEqns = getFieldSize(fieldID);
3227  if (numEqns < 0) ERReturn(-1);
3228 
3229  int nodeFieldEqnNumber = -1;
3230  if (!node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber)) {
3231  ERReturn(-1);
3232  }
3233 
3234  for(int i=0; i<numEqns; i++) eqnNumbers[i] = nodeFieldEqnNumber + i;
3235 
3236  return(FEI_SUCCESS);
3237 }
3238 
3239 //------------------------------------------------------------------------------
3240 int SNL_FEI_Structure::getEqnNumbers(int numIDs,
3241  const GlobalID* IDs,
3242  int idType, int fieldID,
3243  int& numEqns, int* eqnNumbers)
3244 {
3245  //for now, only allow node ID. allowance for element ID is coming soon!!!
3246  if (idType != FEI_NODE) ERReturn(-1);
3247 
3248  int fieldSize = getFieldSize(fieldID);
3249  if (fieldSize < 0) {
3250  ERReturn(-1);
3251  }
3252 
3253  int offset = 0;
3254  for(int i=0; i<numIDs; ++i) {
3255  NodeDescriptor* node = NULL;
3256 
3257  if ( nodeDatabase_->getNodeWithID(IDs[i], node) != 0 ) {
3258  // fei::console_out() << "SNL_FEI_Structure::getEqnNumbers: ERROR getting node " << IDs[i] << FEI_ENDL;
3259  for(int ii=0; ii<fieldSize; ii++) {
3260  eqnNumbers[offset++] = -1;
3261  }
3262  continue;
3263  // ERReturn(-1);
3264  }
3265 
3266  int nodeFieldEqnNumber = -1;
3267 
3268  if ( !node->getFieldEqnNumber(fieldID, nodeFieldEqnNumber) ) {
3269  //if we didn't find a field eqn-number, set eqnNumbers to -1. These
3270  //positions will then be skipped by code that passes the data on down to
3271  //underlying solver objects.
3272  for(int ii=0; ii<fieldSize; ii++) {
3273  eqnNumbers[offset++] = -1;
3274  }
3275  }
3276  else {
3277  //else we did find valid eqn-numbers...
3278  for(int ii=0; ii<fieldSize; ii++) {
3279  eqnNumbers[offset++] = nodeFieldEqnNumber + ii;
3280  }
3281  }
3282  }
3283 
3284  numEqns = fieldSize*numIDs;
3285 
3286  return(FEI_SUCCESS);
3287 }
3288 
3289 //------------------------------------------------------------------------------
3290 int SNL_FEI_Structure::translateToReducedNodeNumber(int nodeNumber, int proc)
3291 {
3292  if (proc != localProc_) {
3293  return(nodeNumber - globalNumNodesVanished_[proc]);
3294  }
3295 
3296  int insertPoint = -1;
3297  int index =
3298  fei::binarySearch(nodeNumber, localVanishedNodeNumbers_, insertPoint);
3299 
3300  int localAdjustment = index < 0 ? insertPoint : index + 1;
3301 
3302  return(nodeNumber - localAdjustment - globalNumNodesVanished_[proc]);
3303 }
3304 
3305 //------------------------------------------------------------------------------
3306 void SNL_FEI_Structure::getEqnBlkInfo(int& numGlobalEqnBlks,
3307  int& numLocalEqnBlks,
3308  int& localBlkOffset) {
3309 
3310  numGlobalEqnBlks = numGlobalEqnBlks_;
3311  numLocalEqnBlks = numLocalEqnBlks_;
3312  localBlkOffset = localBlkOffset_;
3313 }
3314 
3315 //------------------------------------------------------------------------------
3316 int SNL_FEI_Structure::calculateSlaveEqns(MPI_Comm comm)
3317 {
3318  FEI_OSTREAM& os = dbgOut();
3319  if (debugOutput_) os << "# calculateSlaveEqns" << FEI_ENDL;
3320 
3321  if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
3322  eqnCommMgr_ = new EqnCommMgr(comm_);
3323  eqnCommMgr_->setNumRHSs(1);
3324 
3325  std::vector<int> eqns;
3326  std::vector<int> mEqns;
3327  std::vector<double> mCoefs;
3328 
3329  for(size_t i=0; i<slaveVars_->size(); i++) {
3330  int numEqns;
3331  SlaveVariable* svar = (*slaveVars_)[i];
3332 
3333  eqns.resize( getFieldSize(svar->getFieldID()));
3334  CHK_ERR( getEqnNumbers(svar->getNodeID(), FEI_NODE, svar->getFieldID(),
3335  numEqns, &eqns[0]));
3336 
3337  int slaveEqn = eqns[svar->getFieldOffset()];
3338 
3339  const std::vector<GlobalID>* mNodes = svar->getMasterNodeIDs();
3340  const std::vector<int>* mFields = svar->getMasterFields();
3341  const std::vector<double>* mWeights = svar->getWeights();
3342  const std::vector<double>& mWeightsRef = *mWeights;
3343  int mwOffset = 0;
3344 
3345  for(size_t j=0; j<mNodes->size(); j++) {
3346  int mfSize = getFieldSize((*mFields)[j]);
3347 
3348  eqns.resize(mfSize);
3349  GlobalID mNodeID = (*mNodes)[j];
3350  int mFieldID = (*mFields)[j];
3351  CHK_ERR( getEqnNumbers( mNodeID, FEI_NODE, mFieldID,
3352  mfSize, &eqns[0]));
3353 
3354  double fei_eps = 1.e-49;
3355 
3356  for(int k=0; k<mfSize; k++) {
3357  if (std::abs(mWeightsRef[mwOffset++]) > fei_eps) {
3358  mEqns.push_back(eqns[k]);
3359  mCoefs.push_back(mWeightsRef[mwOffset-1]);
3360  }
3361  }
3362  }
3363 
3364  CHK_ERR( slaveEqns_->addEqn(slaveEqn, &mCoefs[0], &mEqns[0],
3365  mEqns.size(), false) );
3366  mEqns.resize(0);
3367  mCoefs.resize(0);
3368  }
3369 
3370 #ifndef FEI_SER
3371  int numLocalSlaves = slaveVars_->size();
3372  int globalMaxSlaves = 0;
3373  CHK_ERR( fei::GlobalMax(comm_, numLocalSlaves, globalMaxSlaves) );
3374 
3375  if (globalMaxSlaves > 0) {
3376  CHK_ERR( gatherSlaveEqns(comm, eqnCommMgr_, slaveEqns_) );
3377  }
3378 #endif
3379 
3380  globalNumNodesVanished_.resize(numProcs_+1, 0);
3381 
3382  slvEqnNumbers_ = &(slaveEqns_->eqnNumbers());
3383  numSlvs_ = slvEqnNumbers_->size();
3384  if (numSlvs_ > 0) {
3385  //first, let's remove any 'couplings' among the slave equations. A coupling
3386  //is where a slave depends on a master which is also a slave that depends on
3387  //other masters.
3388 
3389  if (debugOutput_) {
3390  os << "# slave-equations:" << FEI_ENDL;
3391  os << *slaveEqns_;
3392  os << "# leaving calculateSlaveEqns" << FEI_ENDL;
3393  }
3394 
3395  int levelsOfCoupling;
3396  CHK_ERR( removeCouplings(*slaveEqns_, levelsOfCoupling) );
3397 
3398  if (debugOutput_) {
3399  os << "# SNL_FEI_Structure: " << levelsOfCoupling
3400  << " level(s) of coupling removed: " << FEI_ENDL;
3401  }
3402 
3403  lowestSlv_ = (*slvEqnNumbers_)[0];
3404  highestSlv_ = (*slvEqnNumbers_)[numSlvs_-1];
3405 
3406  //For each slave equation, we need to find out if we (the local proc) either
3407  //own or share the node from which that equation arises. If we own or share
3408  //a slave node, then we will need access to the solution for each of the
3409  //associated master equations, and all other processors that share the slave
3410  //will also need access to all of the master solutions.
3411  //So,
3412  // 1. for each slave node that we share but don't own, we need to add the
3413  // master equations to the eqnCommMgr_ object using addRemoteIndices, if
3414  // they're not local.
3415  // 2. for each slave node that we own, we need to add the master equations
3416  // to the eqnCommMgr_ object using addLocalEqn for each processor that
3417  // shares the slave node.
3418 
3419  std::vector<int>& slvEqns = *slvEqnNumbers_;
3420  std::vector<fei::CSVec*>& mstrEqns = slaveEqns_->eqns();
3421 
3422  //keep track of the number of locally owned nodes that vanish due to the
3423  //fact that all equations at that node are slave equations.
3424  int numLocalNodesVanished = 0;
3425 
3426  GlobalID lastNodeID = -1;
3427 
3428  for(int i=0; i<numSlvs_; i++) {
3429  const NodeDescriptor* node = NULL;
3430  int reducedSlaveEqn;
3431  translateToReducedEqn(slvEqns[i], reducedSlaveEqn);
3432  int numMasters = mstrEqns[i]->size();
3433 
3434  int err = nodeDatabase_->getNodeWithEqn(slvEqns[i], node);
3435  if (err != 0) {
3436  if (debugOutput_) {
3437  os << "# no local node for slave eqn " << slvEqns[i] << FEI_ENDL;
3438  }
3439 
3440  continue;
3441  }
3442 
3443  if (node->getGlobalNodeID() != lastNodeID &&
3444  node->getOwnerProc() == localProc_) {
3445  if (nodalEqnsAllSlaves(node, slvEqns)) {
3446  numLocalNodesVanished++;
3447  if (fei::sortedListInsert(node->getNodeNumber(), localVanishedNodeNumbers_)
3448  == -2) {
3449  ERReturn(-1);
3450  }
3451  }
3452  lastNodeID = node->getGlobalNodeID();
3453  }
3454 
3455  GlobalID slvNodeID = node->getGlobalNodeID();
3456  int shIndex = nodeCommMgr_->getSharedNodeIndex(slvNodeID);
3457  if (shIndex < 0) {
3458  continue;
3459  }
3460 
3461  std::vector<int>& sharingProcs = nodeCommMgr_->getSharedNodeProcs(shIndex);
3462 
3463  for(int j=0; j<numMasters; j++) {
3464  int masterEquation = mstrEqns[i]->indices()[j];
3465  int owner = getOwnerProcForEqn( masterEquation );
3466 
3467  int reducedMasterEqn;
3468  translateToReducedEqn(masterEquation, reducedMasterEqn);
3469 
3470  if (owner == localProc_) {
3471  int numSharing = sharingProcs.size();
3472  for(int sp=0; sp<numSharing; sp++) {
3473  if (sharingProcs[sp] == localProc_) continue;
3474 
3475  if (debugOutput_) {
3476  os << "# slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
3477  << ", master eqn " << masterEquation << " eqnCommMgr_->addLocal "
3478  << reducedMasterEqn << ", proc " << sharingProcs[sp] << FEI_ENDL;
3479  }
3480  eqnCommMgr_->addLocalEqn(reducedMasterEqn, sharingProcs[sp]);
3481  slvCommMgr_->addRemoteIndices(reducedSlaveEqn, sharingProcs[sp],
3482  &reducedMasterEqn, 1);
3483  }
3484  }
3485  else {
3486  if (debugOutput_) {
3487  os << "# slave node " << (int)slvNodeID << ",eqn "<<slvEqns[i]
3488  << ", master eqn " << masterEquation << " eqnCommMgr_->addRemote "
3489  << reducedMasterEqn << ", proc " << owner << FEI_ENDL;
3490  }
3491  eqnCommMgr_->addRemoteIndices(reducedMasterEqn, owner,
3492  &reducedMasterEqn, 1);
3493  slvCommMgr_->addLocalEqn(reducedSlaveEqn, owner);
3494  }
3495  }
3496  }
3497 
3498  std::vector<int> tmp(1), tmp2(numProcs_);
3499  tmp[0] = numLocalNodesVanished;
3500  CHK_ERR( fei::Allgatherv(comm_, tmp, tmp2, globalNumNodesVanished_) );
3501 
3502  if ((int)globalNumNodesVanished_.size() != numProcs_) {
3503  ERReturn(-1);
3504  }
3505 
3506  globalNumNodesVanished_.resize(numProcs_+1);
3507  int tmpNumVanished = 0;
3508  for(int proc=0; proc<numProcs_; ++proc) {
3509  int temporary = tmpNumVanished;
3510  tmpNumVanished += globalNumNodesVanished_[proc];
3511  globalNumNodesVanished_[proc] = temporary;
3512  }
3513  globalNumNodesVanished_[numProcs_] = tmpNumVanished;
3514  }
3515 
3516  if (slaveMatrix_ != NULL) delete slaveMatrix_;
3517  slaveMatrix_ = new fei::FillableMat(*slaveEqns_);
3518 
3519  if (debugOutput_) {
3520  os << "# slave-equations:" << FEI_ENDL;
3521  os << *slaveEqns_;
3522  os << "# leaving calculateSlaveEqns" << FEI_ENDL;
3523  }
3524 
3525  return(0);
3526 }
3527 
3528 //------------------------------------------------------------------------------
3529 int SNL_FEI_Structure::removeCouplings(EqnBuffer& eqnbuf, int& levelsOfCoupling)
3530 {
3531  std::vector<int>& eqnNumbers = eqnbuf.eqnNumbers();
3532  std::vector<fei::CSVec*>& eqns = eqnbuf.eqns();
3533 
3534  std::vector<double> tempCoefs;
3535  std::vector<int> tempIndices;
3536 
3537  levelsOfCoupling = 0;
3538  bool finished = false;
3539  while(!finished) {
3540  bool foundCoupling = false;
3541  for(size_t i=0; i<eqnNumbers.size(); i++) {
3542  int rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
3543 
3544  if (rowIndex == (int)i) {
3545  fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
3546  << " illegal master-slave constraint coupling. Eqn "
3547  << eqnNumbers[i] << " is both master and slave. " << FEI_ENDL;
3548  ERReturn(-1);
3549  }
3550 
3551  while(rowIndex >= 0) {
3552  foundCoupling = true;
3553  double coef = 0.0;
3554  CHK_ERR( eqnbuf.getCoefAndRemoveIndex( eqnNumbers[rowIndex],
3555  eqnNumbers[i], coef) );
3556 
3557  std::vector<int>& indicesRef = eqns[i]->indices();
3558  std::vector<double>& coefsRef = eqns[i]->coefs();
3559 
3560  int len = indicesRef.size();
3561  tempCoefs.resize(len);
3562  tempIndices.resize(len);
3563 
3564  double* tempCoefsPtr = &tempCoefs[0];
3565  int* tempIndicesPtr = &tempIndices[0];
3566  double* coefsPtr = &coefsRef[0];
3567  int* indicesPtr = &indicesRef[0];
3568 
3569  for(int j=0; j<len; ++j) {
3570  tempIndicesPtr[j] = indicesPtr[j];
3571  tempCoefsPtr[j] = coef*coefsPtr[j];
3572  }
3573 
3574  CHK_ERR( eqnbuf.addEqn(eqnNumbers[rowIndex], tempCoefsPtr,
3575  tempIndicesPtr, len, true) );
3576 
3577  rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
3578  }
3579  }
3580  if (foundCoupling) ++levelsOfCoupling;
3581  else finished = true;
3582 
3583  if (levelsOfCoupling>1 && !finished) {
3584  fei::console_out() <<" SNL_FEI_Structure::removeCouplings ERROR,"
3585  << " too many (>1) levels of master-slave constraint coupling. "
3586  << "Hint: coupling is considered infinite if two slaves depend on "
3587  << "each other. This may or may not be the case here." << FEI_ENDL;
3588  ERReturn(-1);
3589  }
3590  }
3591 
3592  return(0);
3593 }
3594 
3595 #ifndef FEI_SER
3596 //------------------------------------------------------------------------------
3597 int SNL_FEI_Structure::gatherSlaveEqns(MPI_Comm comm,
3598  EqnCommMgr* eqnCommMgr,
3599  EqnBuffer* slaveEqns)
3600 {
3601  int numProcs = 0;
3602  if (MPI_Comm_size(comm, &numProcs) != MPI_SUCCESS) return(-1);
3603  if (numProcs == 1) return(0);
3604  int localProc;
3605  if (MPI_Comm_rank(comm, &localProc) != MPI_SUCCESS) return(-1);
3606 
3607  //We're going to send all of our slave equations to all other processors, and
3608  //receive the slave equations from all other processors.
3609  //So we'll first fill a ProcEqns object with all of our eqn/proc pairs,
3610  //then use the regular exchange functions from EqnCommMgr.
3611  ProcEqns localProcEqns, remoteProcEqns;
3612  std::vector<int>& slvEqnNums = slaveEqns->eqnNumbers();
3613  fei::CSVec** slvEqnsPtr = NULL;
3614  if (slaveEqns->eqns().size() > 0) slvEqnsPtr = &(slaveEqns->eqns()[0]);
3615 
3616  for(size_t i=0; i<slvEqnNums.size(); i++) {
3617  for(int p=0; p<numProcs; p++) {
3618  if (p == localProc) continue;
3619 
3620  localProcEqns.addEqn(slvEqnNums[i], slvEqnsPtr[i]->size(), p);
3621  }
3622  }
3623 
3624  CHK_ERR( eqnCommMgr->mirrorProcEqns(localProcEqns, remoteProcEqns) );
3625 
3626  CHK_ERR( eqnCommMgr->mirrorProcEqnLengths(localProcEqns,
3627  remoteProcEqns));
3628 
3629  EqnBuffer remoteEqns;
3630  CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm, &localProcEqns, slaveEqns,
3631  &remoteProcEqns, &remoteEqns,
3632  false) );
3633 
3634  //Now the remoteEqns buffer should hold the slave equations from all other
3635  //processors, so let's add them to the ones we already have.
3636  CHK_ERR( slaveEqns->addEqns(remoteEqns, false) );
3637 
3638  return(0);
3639 }
3640 #endif
3641 
3642 //------------------------------------------------------------------------------
3643 bool SNL_FEI_Structure::isSlaveEqn(int eqn)
3644 {
3645  if (numSlvs_ == 0) return(false);
3646 
3647  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
3648  int insertPoint = -1;
3649  int foundOffset = fei::binarySearch(eqn, slvEqns, insertPoint);
3650 
3651  if (foundOffset >= 0) return(true);
3652  else return(false);
3653 }
3654 
3655 //------------------------------------------------------------------------------
3656 bool SNL_FEI_Structure::translateToReducedEqn(int eqn, int& reducedEqn)
3657 {
3658  if (numSlvs_ == 0) { reducedEqn = eqn; return(false); }
3659 
3660  if (eqn < lowestSlv_) {reducedEqn = eqn; return(false); }
3661  if (eqn > highestSlv_) {reducedEqn = eqn - numSlvs_; return(false); }
3662 
3663  int index = 0;
3664  int foundOffset = fei::binarySearch(eqn, *slvEqnNumbers_, index);
3665 
3666  bool isSlave = false;
3667 
3668  if (foundOffset < 0) {
3669  reducedEqn = eqn - index;
3670  }
3671  else {
3672  isSlave = true; reducedEqn = eqn - (foundOffset+1);
3673  }
3674 
3675  return(isSlave);
3676 }
3677 
3678 //------------------------------------------------------------------------------
3680 {
3681  int numSlvs = slaveEqns_->getNumEqns();
3682 
3683  if (numSlvs == 0) return(reducedEqn);
3684 
3685  const int* slvEqns = &(slaveEqns_->eqnNumbers()[0]);
3686 
3687  if (reducedEqn < slvEqns[0]) return(reducedEqn);
3688 
3689  int eqn = reducedEqn;
3690 
3691  for(int i=0; i<numSlvs; i++) {
3692  if (eqn < slvEqns[i]) return(eqn);
3693  eqn++;
3694  }
3695 
3696  return(eqn);
3697 }
3698 
3699 //------------------------------------------------------------------------------
3701  std::vector<int>*& masterEqns)
3702 {
3703  if (slaveEqns_->getNumEqns() == 0) {
3704  masterEqns = NULL;
3705  return(0);
3706  }
3707 
3708  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
3709  int index = 0;
3710  int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
3711 
3712  if (foundOffset >= 0) {
3713  masterEqns = &(slaveEqns_->eqns()[foundOffset]->indices());
3714  }
3715  else {
3716  masterEqns = NULL;
3717  }
3718 
3719  return(0);
3720 }
3721 
3722 //------------------------------------------------------------------------------
3723 int SNL_FEI_Structure::getMasterEqnCoefs(int slaveEqn,
3724  std::vector<double>*& masterCoefs)
3725 {
3726  if (slaveEqns_->getNumEqns() == 0) {
3727  masterCoefs = NULL;
3728  return(0);
3729  }
3730 
3731  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
3732  int index = 0;
3733  int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
3734 
3735  if (foundOffset >= 0) {
3736  masterCoefs = &(slaveEqns_->eqns()[foundOffset]->coefs());
3737  }
3738  else {
3739  masterCoefs = NULL;
3740  }
3741 
3742  return(0);
3743 }
3744 
3745 //------------------------------------------------------------------------------
3746 int SNL_FEI_Structure::getMasterEqnRHS(int slaveEqn,
3747  double& rhsValue)
3748 {
3749  if (slaveEqns_->getNumEqns() == 0) {
3750  return(0);
3751  }
3752 
3753  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
3754  int index = 0;
3755  int foundOffset = fei::binarySearch(slaveEqn, slvEqns, index);
3756 
3757  if (foundOffset >= 0) {
3758  std::vector<double>* rhsCoefsPtr = (*(slaveEqns_->rhsCoefsPtr()))[foundOffset];
3759  rhsValue = (*rhsCoefsPtr)[0];
3760  }
3761 
3762  return(0);
3763 }
3764 
3765 //------------------------------------------------------------------------------
3766 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID,
3767  int interleaveStrategy,
3768  int* scatterIndices)
3769 {
3770  int index = fei::binarySearch(blockID, blockIDs_);
3771 
3772  if (index < 0) {
3773  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
3774  << (int)blockID << " not found. Aborting." << FEI_ENDL;
3775  std::abort();
3776  }
3777 
3778  std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
3779 
3780  std::map<GlobalID,int>::const_iterator
3781  iter = elemIDs.find(elemID);
3782 
3783  if (iter == elemIDs.end()) {
3784  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: "
3785  << (int)blockID << ", elemID "
3786  << (int)elemID << " not found. Aborting." << FEI_ENDL;
3787  std::abort();
3788  }
3789 
3790  int elemIndex = iter->second;
3791 
3792  getScatterIndices_index(index, elemIndex,
3793  interleaveStrategy, scatterIndices);
3794 }
3795 
3796 //------------------------------------------------------------------------------
3797 void SNL_FEI_Structure::getScatterIndices_ID(GlobalID blockID, GlobalID elemID,
3798  int interleaveStrategy,
3799  int* scatterIndices,
3800  int* blkScatterIndices,
3801  int* blkSizes)
3802 {
3803  int index = fei::binarySearch(blockID, blockIDs_);
3804 
3805  if (index < 0) {
3806  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID "
3807  << (int)blockID << " not found. Aborting." << FEI_ENDL;
3808  std::abort();
3809  }
3810 
3811  std::map<GlobalID,int>& elemIDs = connTables_[index]->elemIDs;
3812 
3813  std::map<GlobalID,int>::const_iterator
3814  iter = elemIDs.find(elemID);
3815 
3816  if (iter == elemIDs.end()) {
3817  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_ID: ERROR, blockID: "
3818  << (int)blockID << ", elemID "
3819  << (int)elemID << " not found. Aborting." << FEI_ENDL;
3820  std::abort();
3821  }
3822 
3823  int elemIndex = iter->second;
3824 
3825  getScatterIndices_index(index, elemIndex,
3826  interleaveStrategy, scatterIndices,
3827  blkScatterIndices, blkSizes);
3828 }
3829 
3830 //------------------------------------------------------------------------------
3831 int SNL_FEI_Structure::getBlkScatterIndices_index(int blockIndex,
3832  int elemIndex,
3833  int* scatterIndices)
3834 {
3835  BlockDescriptor& block = *(blocks_[blockIndex]);
3836  int numNodes = block.getNumNodesPerElement();
3837  work_nodePtrs_.resize(numNodes);
3838  NodeDescriptor** nodes = &work_nodePtrs_[0];
3839  int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
3840  if (err) {
3841  fei::console_out() << "SNL_FEI_Structure::getBlkScatterIndices_index: ERROR getting"
3842  << " node descriptors." << FEI_ENDL;
3843  ERReturn(-1);
3844  }
3845 
3846  int offset = 0;
3847  return( getNodeBlkIndices(nodes, numNodes, scatterIndices, offset) );
3848 }
3849 
3850 //------------------------------------------------------------------------------
3851 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
3852  int interleaveStrategy,
3853  int* scatterIndices)
3854 {
3855 //On input, scatterIndices, is assumed to be allocated by the calling code,
3856 // and be of length the number of equations per element.
3857 //
3858  BlockDescriptor& block = *(blocks_[blockIndex]);
3859  int numNodes = block.getNumNodesPerElement();
3860  int* fieldsPerNode = block.fieldsPerNodePtr();
3861  int** fieldIDs = block.fieldIDsTablePtr();
3862 
3863  work_nodePtrs_.resize(numNodes);
3864  NodeDescriptor** nodes = &work_nodePtrs_[0];
3865 
3866  int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
3867  if (err) {
3868  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
3869  << " node descriptors." << FEI_ENDL;
3870  std::abort();
3871  }
3872 
3873  int offset = 0;
3874  if (fieldDatabase_->size() == 1) {
3875  err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
3876  scatterIndices, offset);
3877  if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
3878  }
3879  else {
3880  switch (interleaveStrategy) {
3881  case 0:
3882  err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
3883  scatterIndices, offset);
3884  if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
3885  break;
3886 
3887  case 1:
3888  err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
3889  scatterIndices, offset);
3890  if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
3891  break;
3892 
3893  default:
3894  fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
3895  break;
3896  }
3897  }
3898 
3899  //now the element-DOF.
3900  int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
3901  std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
3902 
3903  for(int i=0; i<numElemDOF; i++) {
3904  scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
3905  }
3906 }
3907 
3908 //------------------------------------------------------------------------------
3909 void SNL_FEI_Structure::getScatterIndices_index(int blockIndex, int elemIndex,
3910  int interleaveStrategy,
3911  int* scatterIndices,
3912  int* blkScatterIndices,
3913  int* blkSizes)
3914 {
3915 //On input, scatterIndices, is assumed to be allocated by the calling code,
3916 // and be of length the number of equations per element.
3917 //
3918  BlockDescriptor& block = *(blocks_[blockIndex]);
3919  int numNodes = block.getNumNodesPerElement();
3920  int* fieldsPerNode = block.fieldsPerNodePtr();
3921  int** fieldIDs = block.fieldIDsTablePtr();
3922 
3923  work_nodePtrs_.resize(numNodes);
3924  NodeDescriptor** nodes = &work_nodePtrs_[0];
3925 
3926  int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
3927  if (err) {
3928  fei::console_out() << "SNL_FEI_Structure::getScatterIndices_index: ERROR getting"
3929  << " node descriptors." << FEI_ENDL;
3930  std::abort();
3931  }
3932 
3933  int offset = 0, blkOffset = 0;
3934  if (fieldDatabase_->size() == 1) {
3935  err = getNodeIndices_simple(nodes, numNodes, fieldIDs[0][0],
3936  scatterIndices, offset,
3937  blkScatterIndices, blkSizes, blkOffset);
3938  if (err) fei::console_out() << "ERROR in getNodeIndices_simple." << FEI_ENDL;
3939  }
3940  else {
3941  switch (interleaveStrategy) {
3942  case 0:
3943  err = getNodeMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
3944  scatterIndices, offset,
3945  blkScatterIndices, blkSizes, blkOffset);
3946  if (err) fei::console_out() << "ERROR in getNodeMajorIndices." << FEI_ENDL;
3947  break;
3948 
3949  case 1:
3950  err = getFieldMajorIndices(nodes, numNodes, fieldIDs, fieldsPerNode,
3951  scatterIndices, offset);
3952  if (err) fei::console_out() << "ERROR in getFieldMajorIndices." << FEI_ENDL;
3953  break;
3954 
3955  default:
3956  fei::console_out() << "ERROR, unrecognized interleaveStrategy." << FEI_ENDL;
3957  break;
3958  }
3959  }
3960 
3961  //now the element-DOF.
3962  int numElemDOF = blocks_[blockIndex]->getNumElemDOFPerElement();
3963  std::vector<int>& elemDOFEqns = blocks_[blockIndex]->elemDOFEqnNumbers();
3964 
3965  for(int i=0; i<numElemDOF; i++) {
3966  scatterIndices[offset++] = elemDOFEqns[elemIndex] + i;
3967  if (interleaveStrategy == 0) {
3968  blkSizes[blkOffset] = 1;
3969  blkScatterIndices[blkOffset++] = elemDOFEqns[elemIndex] + i;
3970  }
3971  }
3972 }
3973 
3974 //------------------------------------------------------------------------------
3975 int SNL_FEI_Structure::getElemNodeDescriptors(int blockIndex, int elemIndex,
3976  NodeDescriptor** nodes)
3977 {
3978  //Private function, called by 'getScatterIndices_index'. We can safely
3979  //assume that blockIndex and elemIndex are valid in-range indices.
3980  //
3981  //This function's task is to obtain the NodeDescriptor objects, from the
3982  //nodeDatabase, that are connected to the specified element.
3983  //
3984  int numNodes = connTables_[blockIndex]->numNodesPerElem;
3985 
3986  if (activeNodesInitialized_) {
3987  NodeDescriptor** elemNodePtrs =
3988  &((*connTables_[blockIndex]->elem_conn_ptrs)[elemIndex*numNodes]);
3989  for(int i=0; i<numNodes; ++i) nodes[i] = elemNodePtrs[i];
3990  }
3991  else {
3992  const GlobalID* elemConn =
3993  &((*connTables_[blockIndex]->elem_conn_ids)[elemIndex*numNodes]);
3994  for(int i=0; i<numNodes; ++i) {
3995  CHK_ERR( nodeDatabase_->getNodeWithID(elemConn[i], nodes[i]));
3996  }
3997  }
3998 
3999  return(FEI_SUCCESS);
4000 }
4001 
4002 //------------------------------------------------------------------------------
4003 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
4004  int numNodes,
4005  int fieldID,
4006  int* scatterIndices,
4007  int& offset)
4008 {
4009  int fieldSize = getFieldSize(fieldID);
4010 
4011  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4012  NodeDescriptor& node = *(nodes[nodeIndex]);
4013  const int* eqnNumbers = node.getFieldEqnNumbers();
4014  int eqn = eqnNumbers[0];
4015  scatterIndices[offset++] = eqn;
4016  if (fieldSize > 1) {
4017  for(int i=1; i<fieldSize; i++) {
4018  scatterIndices[offset++] = eqn+i;
4019  }
4020  }
4021  }
4022  return(0);
4023 }
4024 
4025 //------------------------------------------------------------------------------
4026 int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
4027  int numNodes,
4028  int fieldID,
4029  int* scatterIndices,
4030  int& offset,
4031  int* blkScatterIndices,
4032  int* blkSizes,
4033  int& blkOffset)
4034 {
4035  int fieldSize = getFieldSize(fieldID);
4036 
4037  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4038  NodeDescriptor& node = *(nodes[nodeIndex]);
4039  const int* eqnNumbers = node.getFieldEqnNumbers();
4040  int eqn = eqnNumbers[0];
4041  scatterIndices[offset++] = eqn;
4042  if (fieldSize > 1) {
4043  for(int i=1; i<fieldSize; i++) {
4044  scatterIndices[offset++] = eqn+i;
4045  }
4046  }
4047  blkSizes[blkOffset] = node.getNumNodalDOF();
4048  blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
4049  }
4050  return(0);
4051 }
4052 
4053 //------------------------------------------------------------------------------
4054 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
4055  int** fieldIDs, int* fieldsPerNode,
4056  int* scatterIndices, int& offset)
4057 {
4058  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4059 
4060  NodeDescriptor& node = *(nodes[nodeIndex]);
4061 
4062  const int* nodeFieldIDList = node.getFieldIDList();
4063  const int* nodeEqnNums = node.getFieldEqnNumbers();
4064  int numFields = node.getNumFields();
4065 
4066  int* fieldID_ind = fieldIDs[nodeIndex];
4067 
4068  for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
4069  int numEqns = getFieldSize(fieldID_ind[j]);
4070  assert(numEqns >= 0);
4071 
4072  int insert = -1;
4073  int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
4074  numFields, insert);
4075 
4076  if (nind >= 0) {
4077  int eqn = nodeEqnNums[nind];
4078 
4079  if (eqn < 0) {
4080  FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
4081  << (int)node.getGlobalNodeID()
4082  << ", field " << nodeFieldIDList[nind]
4083  << " has equation number " << eqn << FEI_ENDL;
4084  std::abort();
4085  }
4086 
4087  for(int jj=0; jj<numEqns; jj++) {
4088  scatterIndices[offset++] = eqn+jj;
4089  }
4090  }
4091  else {
4092  if (outputLevel_ > 2) {
4093  fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
4094  << " not found for node "
4095  << (int)(node.getGlobalNodeID()) << FEI_ENDL;
4096  }
4097  }
4098  }
4099  }
4100 
4101  return(FEI_SUCCESS);
4102 }
4103 
4104 //------------------------------------------------------------------------------
4105 int SNL_FEI_Structure::getNodeBlkIndices(NodeDescriptor** nodes,
4106  int numNodes,
4107  int* scatterIndices,
4108  int& offset)
4109 {
4110  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4111  NodeDescriptor* node = nodes[nodeIndex];
4112  scatterIndices[offset++] = node->getBlkEqnNumber();
4113  }
4114  return(0);
4115 }
4116 
4117 //------------------------------------------------------------------------------
4118 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
4119  int** fieldIDs, int* fieldsPerNode,
4120  int* scatterIndices, int& offset,
4121  int* blkScatterIndices,
4122  int* blkSizes,
4123  int& blkOffset)
4124 {
4125  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4126 
4127  NodeDescriptor& node = *(nodes[nodeIndex]);
4128 
4129  const int* nodeFieldIDList = node.getFieldIDList();
4130  const int* nodeEqnNums = node.getFieldEqnNumbers();
4131  int numFields = node.getNumFields();
4132 
4133  blkSizes[blkOffset] = node.getNumNodalDOF();
4134  blkScatterIndices[blkOffset++] = node.getBlkEqnNumber();
4135 
4136  int* fieldID_ind = fieldIDs[nodeIndex];
4137 
4138  for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
4139  int numEqns = getFieldSize(fieldID_ind[j]);
4140  assert(numEqns >= 0);
4141 
4142  int insert = -1;
4143  int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
4144  numFields, insert);
4145 
4146  if (nind >= 0) {
4147  int eqn = nodeEqnNums[nind];
4148  if (eqn < 0) {
4149  FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
4150  << (int)node.getGlobalNodeID()
4151  << ", field " << nodeFieldIDList[nind]
4152  << " has equation number " << eqn << FEI_ENDL;
4153  std::abort();
4154  }
4155 
4156  for(int jj=0; jj<numEqns; jj++) {
4157  scatterIndices[offset++] = eqn+jj;
4158  }
4159  }
4160  else {
4161  if (outputLevel_ > 2) {
4162  fei::console_out() << "WARNING, field ID " << fieldIDs[nodeIndex][j]
4163  << " not found for node "
4164  << (int)(node.getGlobalNodeID()) << FEI_ENDL;
4165  }
4166  }
4167  }
4168  }
4169 
4170  return(FEI_SUCCESS);
4171 }
4172 
4173 //------------------------------------------------------------------------------
4174 int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
4175  std::vector<int>* fieldIDs,
4176  std::vector<int>& fieldsPerNode,
4177  std::vector<int>& scatterIndices)
4178 {
4179  int offset = 0;
4180  scatterIndices.resize(0);
4181 
4182  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4183 
4184  NodeDescriptor& node = *(nodes[nodeIndex]);
4185 
4186  const int* nodeFieldIDList = node.getFieldIDList();
4187  const int* nodeEqnNums = node.getFieldEqnNumbers();
4188  int numFields = node.getNumFields();
4189 
4190  int* fieldID_ind = &(fieldIDs[nodeIndex][0]);
4191 
4192  for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
4193  int numEqns = getFieldSize(fieldID_ind[j]);
4194  assert(numEqns >= 0);
4195 
4196  int insert = -1;
4197  int nind = fei::binarySearch(fieldID_ind[j], nodeFieldIDList,
4198  numFields, insert);
4199 
4200  if (nind >= 0) {
4201  int eqn = nodeEqnNums[nind];
4202 
4203  if (eqn < 0) {
4204  FEI_COUT << "SNL_FEI_Structure::getNodeMajorIndices: ERROR, node "
4205  << (int)node.getGlobalNodeID()
4206  << ", field " << nodeFieldIDList[nind]
4207  << " has equation number " << eqn << FEI_ENDL;
4208  MPI_Abort(comm_, -1);
4209  }
4210 
4211  scatterIndices.resize(offset+numEqns);
4212  int* inds = &scatterIndices[0];
4213 
4214  for(int jj=0; jj<numEqns; jj++) {
4215  inds[offset+jj] = eqn+jj;
4216  }
4217  offset += numEqns;
4218  }
4219  else {
4220  if (outputLevel_ > 2) {
4221  fei::console_out() << "WARNING, field ID " << fieldID_ind[j]
4222  << " not found for node "
4223  << (int)node.getGlobalNodeID() << FEI_ENDL;
4224  }
4225  }
4226  }
4227  }
4228 
4229  return(FEI_SUCCESS);
4230 }
4231 
4232 //------------------------------------------------------------------------------
4233 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
4234  int** fieldIDs, int* fieldsPerNode,
4235  int* scatterIndices, int& offset)
4236 {
4237  //In this function we want to group equations by field, but
4238  //in what order should the fields be?
4239  //Let's just run through the fieldIDs table, and add the fields to a
4240  //flat list, in the order we encounter them, but making sure no fieldID
4241  //gets added more than once.
4242 
4243  std::vector<int> fields;
4244  for(int ii=0; ii<numNodes; ii++) {
4245  for(int j=0; j<fieldsPerNode[ii]; j++) {
4246  if (std::find(fields.begin(), fields.end(), fieldIDs[ii][j]) == fields.end()) {
4247  fields.push_back(fieldIDs[ii][j]);
4248  }
4249  }
4250  }
4251 
4252  int* fieldsPtr = fields.size()>0 ? &fields[0] : NULL;
4253 
4254  //ok, we've got our flat list of fields, so let's proceed to get the
4255  //scatter indices.
4256 
4257  for(size_t i=0; i<fields.size(); i++) {
4258  int field = fieldsPtr[i];
4259 
4260  for(int nodeIndex = 0; nodeIndex < numNodes; ++nodeIndex) {
4261  int fidx = fei::searchList(field, fieldIDs[nodeIndex],
4262  fieldsPerNode[nodeIndex]);
4263  if (fidx < 0) {
4264  continue;
4265  }
4266 
4267  NodeDescriptor* node = nodes[nodeIndex];
4268 
4269  const int* nodeFieldIDList = node->getFieldIDList();
4270  const int* nodeEqnNums = node->getFieldEqnNumbers();
4271  int numFields = node->getNumFields();
4272 
4273  int numEqns = getFieldSize(field);
4274  assert(numEqns >= 0);
4275 
4276  int insert = -1;
4277  int nind = fei::binarySearch(field, nodeFieldIDList,
4278  numFields, insert);
4279 
4280  if (nind > -1) {
4281  for(int jj=0; jj<numEqns; ++jj) {
4282  scatterIndices[offset++] = nodeEqnNums[nind]+jj;
4283  }
4284  }
4285  else {
4286  ERReturn(-1);
4287  }
4288  }
4289  }
4290 
4291  return(FEI_SUCCESS);
4292 }
4293 
4294 //------------------------------------------------------------------------------
4295 int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
4296  std::vector<int>* fieldIDs,
4297  std::vector<int>& fieldsPerNode,
4298  std::vector<int>& scatterIndices)
4299 {
4300  //In this function we want to group equations by field, but
4301  //in what order should the fields be?
4302  //Let's just run through the fieldIDs table, and add the fields to a
4303  //flat list, in the order we encounter them, but making sure no fieldID
4304  //gets added more than once.
4305 
4306  std::vector<int> fields;
4307  for(int ii=0; ii<numNodes; ii++) {
4308  for(int j=0; j<fieldsPerNode[ii]; j++) {
4309  std::vector<int>& fldIDs = fieldIDs[ii];
4310  if (std::find(fields.begin(), fields.end(), fldIDs[j]) == fields.end()) {
4311  fields.push_back(fldIDs[j]);
4312  }
4313  }
4314  }
4315 
4316  //ok, we've got our flat list of fields, so let's proceed to get the
4317  //scatter indices.
4318 
4319  int offset = 0;
4320  scatterIndices.resize(0);
4321 
4322  for(size_t i=0; i<fields.size(); i++) {
4323  for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
4324 
4325  const int* nodeFieldIDList = nodes[nodeIndex]->getFieldIDList();
4326  const int* nodeEqnNums = nodes[nodeIndex]->getFieldEqnNumbers();
4327  int numFields = nodes[nodeIndex]->getNumFields();
4328 
4329  int numEqns = getFieldSize(fields[i]);
4330  assert(numEqns >= 0);
4331 
4332  int insert = -1;
4333  int nind = fei::binarySearch(fields[i], nodeFieldIDList,
4334  numFields, insert);
4335 
4336  if (nind >= 0) {
4337  for(int jj=0; jj<numEqns; jj++) {
4338  scatterIndices.push_back(nodeEqnNums[nind]+jj);
4339  offset++;
4340  }
4341  }
4342  else {
4343  if (outputLevel_ > 2) {
4344  fei::console_out() << "WARNING, field ID " << fields[i]
4345  << " not found for node "
4346  << (int)nodes[nodeIndex]->getGlobalNodeID() << FEI_ENDL;
4347  }
4348  }
4349  }
4350  }
4351 
4352  return(FEI_SUCCESS);
4353 }
4354 
4355 //------------------------------------------------------------------------------
4356 void SNL_FEI_Structure::addCR(int CRID,
4358  std::map<GlobalID,snl_fei::Constraint<GlobalID>* >& crDB)
4359 {
4360  std::map<GlobalID,snl_fei::Constraint<GlobalID>* >::iterator
4361  cr_iter = crDB.find(CRID);
4362 
4363  if (cr_iter == crDB.end()) {
4364  cr = new ConstraintType;
4365  crDB.insert(std::pair<GlobalID,snl_fei::Constraint<GlobalID>*>(CRID, cr));
4366  }
4367 }
int isInIndices(int eqn)
int sortedListInsert(const T &item, std::vector< T > &list)
int translateMatToReducedEqns(fei::CSRMat &mat)
int getFieldSize(int fieldID)
const int * getFieldIDsPtr()
int Allgatherv(MPI_Comm comm, std::vector< T > &sendbuf, std::vector< int > &recvLengths, std::vector< T > &recvbuf)
int countLocalNodeDescriptors(int localRank)
std::vector< int > & getMasterFieldIDs()
int getBlockDescriptor_index(int index, BlockDescriptor *&block)
std::map< GlobalID, int > & getNodeIDs()
int getNodeWithNumber(int nodeNumber, const NodeDescriptor *&node) const
int addEqn(int eqnNumber, const double *coefs, const int *indices, int len, bool accumulate, bool create_indices_union=false)
SNL_FEI_Structure(MPI_Comm comm)
void multiply_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
Definition: fei_CSRMat.cpp:189
int getMasterEqnRHS(int slaveEqn, double &rhsValue)
int initNodeID(GlobalID nodeID)
int size() const
int ptEqnToBlkEqn(int ptEqn)
std::vector< fei::CSVec * > & eqns()
int getMasterEqnCoefs(int slaveEqn, std::vector< double > *&masterCoefs)
int mirrorProcEqns(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
void addEqn(int eqnNumber, int proc)
std::vector< int > rowNumbers
std::vector< int > & eqnNumbers()
int getEqnNumber() const
void copySetToArray(const SET_TYPE &set_obj, int lenList, int *list)
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
void setEqnNumber(int eqn)
void setConstraintID(int id)
int getMasterEqnNumbers(int slaveEqn, std::vector< int > *&masterEqns)
std::vector< int > packedColumnIndices
std::vector< std::vector< double > * > * rhsCoefsPtr()
std::vector< int > rowOffsets
void setBlkEqnNumber(int blkEqn)
void addLocalEqn(int eqnNumber, int srcProc)
int ** fieldIDsTablePtr()
int binarySearch(const T &item, const T *list, int len)
int getNodeWithEqn(int eqnNumber, const NodeDescriptor *&node) const
int synchronize(int firstLocalNodeNumber, int firstLocalEqn, int localRank, MPI_Comm comm)
std::vector< std::vector< int > * > & procEqnNumbersPtr()
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int setBlkEqnSize(int blkEqn, int size)
void insert2(const T &item)
int translateToReducedEqns(EqnCommMgr &eqnCommMgr)
int getCoefAndRemoveIndex(int eqnNumber, int colIndex, double &coef)
bool isInLocalElement(int nodeNumber)
int setEqn(int ptEqn, int blkEqn)
const int * getNumFieldsPerNode(GlobalID blockID)
int eqnToBlkEqn(int eqn) const
int countLocalNodalEqns(int localRank)
int addEqns(EqnBuffer &inputEqns, bool accumulate)
bool translateToReducedEqn(int eqn, int &reducedEqn)
std::ostream & console_out()
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int getOwnerProcForEqn(int eqn)
void setMaxBlkEqnSize(int sz)
void setIsPenalty(bool isPenaltyConstr)
void multiply_trans_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
Definition: fei_CSRMat.cpp:268
int localProc(MPI_Comm comm)
int mirrorProcEqnLengths(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
std::map< int, int > * getPtEqns()
int getBlkEqnNumber() const
std::vector< int > & getMasters()
void destroyValues(MAP_TYPE &map_obj)
int getNumNodeDescriptors() const
int translateFromReducedEqn(int reducedEqn)
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
int parameters(int numParams, const char *const *paramStrings)
void getElemBlockInfo(GlobalID blockID, int &interleaveStrategy, int &lumpingStrategy, int &numElemDOF, int &numElements, int &numNodesPerElem, int &numEqnsPerElem)
int searchList(const T &item, const T *list, int len)
int setNumNodesPerElement(int numNodes)
int numProcs(MPI_Comm comm)
int getIndexOfBlock(GlobalID blockID) const
const int *const * getFieldIDsTable(GlobalID blockID)
int getEqnNumber(int nodeNumber, int fieldID)
int initNodeIDs(GlobalID *nodeIDs, int numNodes)
int getNumEqns()