← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/fei/base/SNL_FEI_Structure.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/fei/base/SNL_FEI_Structure.cpp
15
15
 
16
16
#include "fei_defs.h"
17
17
 
18
 
#include "feiArray.hpp"
19
18
#include "fei_TemplateUtils.hpp"
20
 
#include "snl_fei_CommUtils.hpp"
 
19
#include <fei_CommUtils.hpp>
21
20
#include "snl_fei_Constraint.hpp"
22
21
typedef snl_fei::Constraint<GlobalID> ConstraintType;
23
22
 
27
26
 
28
27
#include "fei_SlaveVariable.hpp"
29
28
 
30
 
#include "fei_PatternDescriptor.hpp"
31
29
#include "fei_BlockDescriptor.hpp"
32
30
 
33
31
#include "snl_fei_PointBlockMap.hpp"
34
32
#include "fei_ProcEqns.hpp"
35
33
#include "fei_EqnBuffer.hpp"
36
 
#include "fei_SSVec.hpp"
37
 
#include "fei_SSMat.hpp"
38
 
#include "fei_SSGraph.hpp"
 
34
#include <fei_FillableMat.hpp>
 
35
#include <fei_FillableVec.hpp>
 
36
#include <fei_CSRMat.hpp>
 
37
#include <fei_CSVec.hpp>
39
38
#include "fei_EqnCommMgr.hpp"
40
39
 
41
40
#include "fei_Lookup.hpp"
55
54
   numProcs_(1),
56
55
   fieldDatabase_(new std::map<int,int>),
57
56
   workarray_(),
58
 
   blockIDs_(0, 8),
59
 
   blocks_(0, 8),
60
 
   connTables_(0, 8),
61
 
   patternIDs_(0, 4),
62
 
   patterns_(NULL),
63
 
   patternConn_(NULL),
 
57
   blockIDs_(),
 
58
   blocks_(),
 
59
   connTables_(),
64
60
   nodeDatabase_(NULL),
65
61
   activeNodesInitialized_(false),
66
62
   globalNodeOffsets_(),
74
70
   highestSlv_(0),
75
71
   slaveMatrix_(NULL),
76
72
   globalNumNodesVanished_(),
77
 
   localVanishedNodeNumbers_(0, 32),
78
 
   commUtilsInt_(NULL),
79
 
   commUtilsDouble_(NULL),
 
73
   localVanishedNodeNumbers_(),
80
74
   nodeCommMgr_(NULL),
81
75
   eqnCommMgr_(NULL),
82
76
   slvCommMgr_(NULL),
93
87
   Kid_(NULL),
94
88
   Kdi_(NULL),
95
89
   Kdd_(NULL),
96
 
   tmpMat1_(NULL),
97
 
   tmpMat2_(NULL),
 
90
   tmpMat1_(),
 
91
   tmpMat2_(),
98
92
   reducedEqnCounter_(0),
99
93
   reducedRHSCounter_(0),
100
94
   rSlave_(),
101
95
   cSlave_(),
102
 
   work_nodePtrs_(0),
 
96
   work_nodePtrs_(),
103
97
   structureFinalized_(false),
104
98
   generateGraph_(true),
105
99
   sysMatIndices_(NULL),
114
108
   numGlobalNodes_(0),
115
109
   sysBlkMatIndices_(NULL),
116
110
   matIndicesDestroyed_(false),
117
 
   workSpace_(0, 8),
118
 
   workSpace2_(0, 8),
 
111
   workSpace_(),
119
112
   blkEqnMapper_(new snl_fei::PointBlockMap()),
120
113
   multCRs_(),
121
114
   penCRs_(),
129
122
{
130
123
  numProcs_ = 1, localProc_ = 0, masterProc_ = 0;
131
124
 
132
 
#ifndef FEI_SER
133
 
  MPI_Comm_rank(comm_, &localProc_);
134
 
  MPI_Comm_size(comm_, &numProcs_);
 
125
  localProc_ = fei::localProc(comm_);
 
126
  numProcs_ = fei::numProcs(comm_);
135
127
  masterProc_ = 0;
136
 
#endif
137
128
 
138
 
  slaveVars_ = new feiArray<SlaveVariable*>;
 
129
  slaveVars_ = new std::vector<SlaveVariable*>;
139
130
  slaveEqns_ = new EqnBuffer;
140
131
 
141
 
  commUtilsInt_ = new snl_fei::CommUtils<int>(comm_);
142
 
  commUtilsDouble_ = new snl_fei::CommUtils<double>(comm_);
143
 
 
144
132
  nodeCommMgr_ = new NodeCommMgr(comm_);
145
133
 
146
 
  eqnCommMgr_ = new EqnCommMgr(*commUtilsInt_);
 
134
  eqnCommMgr_ = new EqnCommMgr(comm_);
147
135
  eqnCommMgr_->setNumRHSs(1);
148
136
 
149
137
  nodeDatabase_ = new NodeDatabase(fieldDatabase_, nodeCommMgr_);
150
138
 
151
 
  Kid_ = new SSMat;
152
 
  Kdi_ = new SSMat;
153
 
  Kdd_ = new SSMat;
154
 
  tmpMat1_ = new SSMat;
155
 
  tmpMat2_ = new SSMat;
 
139
  Kid_ = new fei::FillableMat;
 
140
  Kdi_ = new fei::FillableMat;
 
141
  Kdd_ = new fei::FillableMat;
156
142
}
157
143
 
158
144
//-----------------------------------------------------------------------------
199
185
//-----Destructor--------------------------------------------------------------
200
186
SNL_FEI_Structure::~SNL_FEI_Structure()
201
187
{
202
 
  int j;
203
 
  for(j=0; j<slaveVars_->length(); j++) {
 
188
  for(size_t j=0; j<slaveVars_->size(); j++) {
204
189
    delete (*slaveVars_)[j];
205
190
  }
206
191
  delete slaveVars_;
213
198
 
214
199
  delete nodeCommMgr_;
215
200
  delete eqnCommMgr_;
216
 
  delete commUtilsInt_;
217
 
  delete commUtilsDouble_;
218
201
  delete blkEqnMapper_;
219
202
 
220
203
  destroyMatIndices();
227
210
    penCRs_.clear();
228
211
  }
229
212
 
230
 
  for(int i=0; i<patternIDs_.length(); i++) {
231
 
    delete patterns_[i];
232
 
    delete patternConn_[i];
233
 
  }
234
 
  delete [] patterns_;
235
 
  delete [] patternConn_;
236
 
 
237
213
  delete nodeDatabase_;
238
214
  delete fieldDatabase_;
239
215
 
240
216
  delete Kid_;
241
217
  delete Kdi_;
242
218
  delete Kdd_;
243
 
  delete tmpMat1_;
244
 
  delete tmpMat2_;
245
219
}
246
220
 
247
221
//------------------------------------------------------------------------------
267
241
//------------------------------------------------------------------------------
268
242
void SNL_FEI_Structure::destroyBlockRoster()
269
243
{
270
 
  for(int i=0; i<blockIDs_.length(); i++) delete blocks_[i];
 
244
  for(size_t i=0; i<blockIDs_.size(); i++) delete blocks_[i];
271
245
  blocks_.resize(0);
272
246
}
273
247
 
274
248
//------------------------------------------------------------------------------
275
249
void SNL_FEI_Structure::destroyConnectivityTables()
276
250
{
277
 
  for(int i=0; i<blockIDs_.length(); i++) {
 
251
  for(size_t i=0; i<blockIDs_.size(); i++) {
278
252
    delete connTables_[i];
279
253
  }
280
254
  connTables_.resize(0);
335
309
  lumpingStrategy = block->getLumpingStrategy();
336
310
  numElemDOF = block->getNumElemDOFPerElement();
337
311
  numElements = block->getNumElements();
338
 
  numNodesPerElem = block->numNodesPerElement;
 
312
  numNodesPerElem = block->getNumNodesPerElement();
339
313
  numEqnsPerElem = block->getNumEqnsPerElement();
340
314
}
341
315
 
347
321
 
348
322
  int eqnNumber;
349
323
 
350
 
  NodeDescriptor* node = NULL;
 
324
  const NodeDescriptor* node = NULL;
351
325
  CHK_ERR( nodeDatabase_->getNodeWithNumber(nodeNumber, node) );
352
326
 
353
327
  bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
369
343
{
370
344
  if (eqn < 0) return(-1);
371
345
 
372
 
  int len = globalEqnOffsets_.length(); // len is numProcs+1...
 
346
  int len = globalEqnOffsets_.size(); // len is numProcs+1...
373
347
 
374
348
  for(int i=0; i<len-1; i++) {
375
349
    if (eqn >= globalEqnOffsets_[i] && eqn < globalEqnOffsets_[i+1]) return(i);
458
432
 
459
433
   int numNodalEqns = 0;
460
434
   int countDOF;
461
 
   feiArray<int> distinctFields(0);
 
435
   std::vector<int> distinctFields(0);
462
436
 
463
437
   for(j=0; j<numNodesPerElement; j++) {
464
438
 
479
453
      numNodalEqns += countDOF;
480
454
   }
481
455
 
482
 
   block->setNumDistinctFields(distinctFields.length());
 
456
   block->setNumDistinctFields(distinctFields.size());
483
457
 
484
458
   int numElemDOFPerElement = 0;
485
459
   for(j=0; j<numElemDofFieldsPerElement; j++) {
536
510
  ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
537
511
 
538
512
  std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
539
 
  GlobalID* conn = connTable.elem_conn_ids->dataPtr();
 
513
  GlobalID* conn = &((*connTable.elem_conn_ids)[0]);
540
514
 
541
515
  int elemIndex = elemIDList.size();
542
516
  std::map<GlobalID,int>::iterator
552
526
    elemIDList.insert(std::make_pair(elemID,elemIndex));
553
527
  }
554
528
 
555
 
  int numNodes = block->numNodesPerElement;
 
529
  int numNodes = block->getNumNodesPerElement();
556
530
 
557
531
  if (debugOutput_ && outputLevel_ > 2) {
558
532
    FEI_OSTREAM& os = dbgOut();
589
563
}
590
564
 
591
565
//------------------------------------------------------------------------------
592
 
int SNL_FEI_Structure::initCoefAccessPattern(int patternID,
593
 
                                            int numRowIDs,
594
 
                                            const int* numFieldsPerRow,
595
 
                                            const int* const* rowFieldIDs,
596
 
                                            int numColIDsPerRow,
597
 
                                            const int* numFieldsPerCol,
598
 
                                            const int* const* colFieldIDs,
599
 
                                            int interleaveStrategy)
600
 
{
601
 
//
602
 
//This is a simple function -- just stores this incoming data in a new
603
 
//PatternDescriptor, if this patternID hasn't already been used. If it
604
 
//has already been used, return an error (nonzero).
605
 
//
606
 
   CHK_ERR( addPattern(patternID) )
607
 
   PatternDescriptor* pattern = NULL;
608
 
   CHK_ERR( getPatternDescriptor(patternID, pattern) )
609
 
 
610
 
   if (numRowIDs <= 0 || numColIDsPerRow <= 0) {
611
 
      FEI_CERR << "SNL_FEI_Structure::initCoefAccessPattern: numRowIDs: " << numRowIDs
612
 
          << ", numColIDsPerRow: " << numColIDsPerRow
613
 
          << ". These values must both be" << " > 0." << FEI_ENDL;
614
 
      ERReturn(-1);
615
 
   }
616
 
 
617
 
   int i;
618
 
 
619
 
   CHK_ERR( pattern->setNumRowIDs(numRowIDs) )
620
 
   feiArray<int>& fieldsPerRow = pattern->getNumFieldsPerRow();
621
 
   feiArray<int>* rFieldIDs = pattern->getRowFieldIDs();
622
 
 
623
 
   for(i=0; i<numRowIDs; i++) {
624
 
      fieldsPerRow[i] = numFieldsPerRow[i];
625
 
 
626
 
      for(int j=0; j<fieldsPerRow[i]; j++) {
627
 
        rFieldIDs[i].append(rowFieldIDs[i][j]);
628
 
      }
629
 
   }
630
 
 
631
 
 
632
 
   CHK_ERR( pattern->setNumColIDsPerRow(numColIDsPerRow) )
633
 
   feiArray<int>& fieldsPerCol = pattern->getNumFieldsPerCol();
634
 
   feiArray<int>* cFieldIDs = pattern->getColFieldIDs();
635
 
 
636
 
   for(i=0; i<numColIDsPerRow; i++) {
637
 
      fieldsPerCol[i] = numFieldsPerCol[i];
638
 
 
639
 
      for(int j=0; j<numFieldsPerCol[i]; j++) {
640
 
        cFieldIDs[i].append(colFieldIDs[i][j]);
641
 
      }
642
 
   }
643
 
 
644
 
   pattern->setInterleaveStrategy(interleaveStrategy);
645
 
 
646
 
   return(FEI_SUCCESS);
647
 
}
648
 
 
649
 
//------------------------------------------------------------------------------
650
 
int SNL_FEI_Structure::initCoefAccess(int patternID,
651
 
                                     const int* rowIDTypes,
652
 
                                     const GlobalID* rowIDs,
653
 
                                     const int* colIDTypes,
654
 
                                     const GlobalID* colIDs)
655
 
{
656
 
  PatternDescriptor* pattern = NULL;
657
 
  CHK_ERR( getPatternDescriptor(patternID, pattern) );
658
 
 
659
 
  int numRowIDs = pattern->getNumRowIDs();
660
 
  int numColIDsPerRow = pattern->getNumColIDsPerRow();
661
 
 
662
 
  CHK_ERR( appendPatternConnectivity(patternID,
663
 
                                     numRowIDs, rowIDTypes, rowIDs,
664
 
                                     numColIDsPerRow,
665
 
                                     colIDTypes, colIDs) );
666
 
 
667
 
  for(int i=0; i<numRowIDs; i++){
668
 
    if (rowIDTypes[i] == FEI_NODE) {
669
 
      CHK_ERR( nodeDatabase_->initNodeID(rowIDs[i]) );
670
 
    }
671
 
  }
672
 
 
673
 
   for(int j=0; j<numRowIDs*numColIDsPerRow; j++) {
674
 
      if (colIDTypes[j] == FEI_NODE) {
675
 
        CHK_ERR( nodeDatabase_->initNodeID(colIDs[j]) );
676
 
      }
677
 
   }
678
 
 
679
 
   return(FEI_SUCCESS);
680
 
}
681
 
 
682
 
//------------------------------------------------------------------------------
683
566
int SNL_FEI_Structure::initSlaveVariable(GlobalID slaveNodeID, 
684
567
                                         int slaveFieldID,
685
568
                                         int offsetIntoSlaveField,
732
615
  int woffset = 0;
733
616
 
734
617
  for(int i=0; i<numMasterNodes; i++) {
735
 
    CHK_ERR( svar->addMasterNodeID(masterNodeIDs[i]) );
736
 
    CHK_ERR( svar->addMasterField(masterFieldIDs[i]) );
 
618
    svar->addMasterNodeID(masterNodeIDs[i]);
 
619
    svar->addMasterField(masterFieldIDs[i]);
737
620
    int fieldSize = getFieldSize(masterFieldIDs[i]);
738
621
    if (fieldSize < 0) ERReturn(-1);
739
622
 
740
623
    for(int j=0; j<fieldSize; j++) {
741
 
      CHK_ERR( svar->addWeight(weights[woffset++]) );
 
624
      svar->addWeight(weights[woffset++]);
742
625
    }
743
626
  }
744
627
 
745
 
  CHK_ERR( addSlaveVariable(svar) );
 
628
  addSlaveVariable(svar);
746
629
 
747
630
  return(FEI_SUCCESS);
748
631
}
800
683
    //complete destruction and re-calculation of the FEI structure.
801
684
    //
802
685
    for(int i=0; i<numSharedNodes; ++i) {
803
 
      for(int nc=0; nc<connTables_.length(); ++nc) {
 
686
      for(size_t nc=0; nc<connTables_.size(); ++nc) {
804
687
        if (connTables_[nc]->elem_conn_ids == NULL) continue;
805
 
        int len = connTables_[nc]->elem_conn_ids->length();
 
688
        int len = connTables_[nc]->elem_conn_ids->size();
806
689
        if (len < 1) continue;
807
 
        GlobalID* conn_ids = connTables_[nc]->elem_conn_ids->dataPtr();
808
 
        NodeDescriptor** nodes = connTables_[nc]->elem_conn_ptrs->dataPtr();
 
690
        GlobalID* conn_ids = &((*connTables_[nc]->elem_conn_ids)[0]);
 
691
        NodeDescriptor** nodes = &((*connTables_[nc]->elem_conn_ptrs)[0]);
809
692
 
810
693
        for(int j=0; j<len; ++j) {
811
694
          if (conn_ids[j] == sharedNodeIDs[i]) {
851
734
  ConstraintType& multCR = *multCRPtr;
852
735
 
853
736
  multCR.setConstraintID(CRID);
854
 
  CHK_ERR( multCR.allocate() );
855
737
 
856
 
  feiArray<GlobalID>* CRNodeArray = multCR.getMasters();
857
 
  feiArray<int>* CRFieldArray = multCR.getMasterFieldIDs();
 
738
  std::vector<GlobalID>& CRNodeArray = multCR.getMasters();
 
739
  std::vector<int>& CRFieldArray = multCR.getMasterFieldIDs();
858
740
 
859
741
  for(int j = 0; j < numCRNodes; j++) {
860
 
    CRNodeArray->append(CRNodes[j]);
861
 
    CRFieldArray->append(CRFields[j]);
 
742
    CRNodeArray.push_back(CRNodes[j]);
 
743
    CRFieldArray.push_back(CRFields[j]);
862
744
  }
863
745
 
864
746
  if (debugOutput_) dbgOut() << "#(output) CRID:"<<FEI_ENDL << CRID << FEI_ENDL;
890
772
  penCR.setConstraintID(CRID);
891
773
  penCR.setIsPenalty(true);
892
774
 
893
 
  CHK_ERR( penCR.allocate() );
894
 
 
895
 
  feiArray<GlobalID>* CRNodesArray = penCR.getMasters();
896
 
 
897
 
  feiArray<int>* CRFieldArray = penCR.getMasterFieldIDs();
 
775
  std::vector<GlobalID>& CRNodesArray = penCR.getMasters();
 
776
 
 
777
  std::vector<int>& CRFieldArray = penCR.getMasterFieldIDs();
898
778
 
899
779
  for(int i = 0; i < numCRNodes; i++) {
900
 
    CRNodesArray->append(CRNodes[i]);
901
 
    CRFieldArray->append(CRFields[i]);
 
780
    CRNodesArray.push_back(CRNodes[i]);
 
781
    CRFieldArray.push_back(CRFields[i]);
902
782
  }
903
783
 
904
784
  return(FEI_SUCCESS);
912
792
  //element.
913
793
  if (numProcs_ < 2) return(true);
914
794
 
915
 
  NodeDescriptor* node = NULL;
 
795
  const NodeDescriptor* node = NULL;
916
796
  int err = nodeDatabase_->getNodeWithNumber(nodeNumber, node);
917
797
  if (err != 0) return(false);
918
798
 
958
838
  //  - As the connectivity tables are being run, the NodeCommMgr is given the
959
839
  //    nodeID of each node that is connected to a local element, and the node's
960
840
  //    owner is set to the initial value of localProc_.
961
 
  //  - Any Patterns are also run, and nodal field and ownership information is
962
 
  //    initialized similarly to what's done for the element connectivity
963
 
  //    tables.
964
841
  //1a. finalizeNodeCommMgr()
965
842
  //  - NodeCommMgr::initComplete is called, at which time the ownership of all
966
843
  //    shared nodes is determined.
1008
885
  //    calculations associated with reducing out slave equations, and with
1009
886
  //    putting off-processor contributions (element contributions with shared
1010
887
  //    nodes) into the EqnCommMgr object.
1011
 
  //  - initAccessPatternStructure(), initMultCRStructure() and
1012
 
  //    initPenCRStructure() follow the same general routine as for the
1013
 
  //    initElemBlockStructure() function.
 
888
  //  - initMultCRStructure() and initPenCRStructure() follow the same general
 
889
  //    routine as for the initElemBlockStructure() function.
1014
890
  //  - EqnCommMgr::exchangeIndices() is called, which sends all shared
1015
891
  //    contributions to the processors that own the corresponding equations,
1016
892
  //    after which the receiving processors insert those contributions into
1107
983
  //and establish an EqnBuffer of slave equation numbers, and gather the slave
1108
984
  //equations initialized on any other processors.
1109
985
 
1110
 
  slvCommMgr_ = new EqnCommMgr(*commUtilsInt_);
 
986
  slvCommMgr_ = new EqnCommMgr(comm_);
1111
987
  slvCommMgr_->setNumRHSs(1);
1112
988
 
1113
989
  CHK_ERR( calculateSlaveEqns(comm_) );
1144
1020
      try {
1145
1021
      CHK_ERR( slvCommMgr_->exchangeIndices() );
1146
1022
      }
1147
 
      catch (fei::Exception& exc) {
 
1023
      catch (std::runtime_error& exc) {
1148
1024
        FEI_CERR << exc.what() << FEI_ENDL;
1149
1025
        ERReturn(-1);
1150
1026
      }
1190
1066
 
1191
1067
  CHK_ERR( initElemBlockStructure() );
1192
1068
 
1193
 
  // now we'll add the connectivities for any 'access patterns' that were
1194
 
  // initialized...
1195
 
  CHK_ERR( initAccessPatternStructure() );
1196
 
 
1197
1069
  // next, handle the matrix structure imposed by the constraint relations...
1198
1070
  //
1199
1071
 
1211
1083
  try {
1212
1084
  CHK_ERR( eqnCommMgr_->exchangeIndices(&os) );
1213
1085
  }
1214
 
  catch(fei::Exception& exc) {
 
1086
  catch(std::runtime_error& exc) {
1215
1087
    FEI_CERR << exc.what() << FEI_ENDL;
1216
1088
    ERReturn(-1);
1217
1089
  }
1220
1092
  //of the eqn comm mgr and put them into our local matrix structure.
1221
1093
 
1222
1094
  int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
1223
 
  feiArray<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbersPtr();
1224
 
  feiArray<SSVec*>& recvEqns = eqnCommMgr_->localEqns();
 
1095
  std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
 
1096
  std::vector<fei::CSVec*>& recvEqns = eqnCommMgr_->localEqns();
1225
1097
  int i;
1226
1098
  if (debugOutput_) {
1227
1099
    os << "#     after eqnCommMgr_->exchangeIndices, numRecvEqns: "
1238
1110
      return(-1);
1239
1111
    }
1240
1112
 
1241
 
    for(int j=0; j<recvEqns[i]->length(); j++) {
 
1113
    for(size_t j=0; j<recvEqns[i]->size(); j++) {
1242
1114
      CHK_ERR( createMatrixPosition(eqn, recvEqns[i]->indices()[j],
1243
1115
                                    "frmMatStr") );
1244
1116
    }
1268
1140
 
1269
1141
    int numEqns = block->getNumEqnsPerElement();
1270
1142
    int interleave = block->getInterleaveStrategy();
1271
 
    feiArray<int> scatterIndices(numEqns);
 
1143
    std::vector<int> scatterIndices(numEqns);
1272
1144
 
1273
1145
    GlobalID elemBlockID = block->getGlobalBlockID();
1274
1146
    ConnectivityTable& connTable = getBlockConnectivity(elemBlockID);
1285
1157
    for(int elemIndex = 0; elemIndex < numBlockElems; elemIndex++) {
1286
1158
 
1287
1159
      getScatterIndices_index(bIndex, elemIndex, interleave,
1288
 
                              scatterIndices.dataPtr());
 
1160
                              &scatterIndices[0]);
1289
1161
 
1290
1162
      //now, store the structure that will arise from this contribution,
1291
1163
      //after first performing any manipulations associated with slave
1302
1174
}
1303
1175
 
1304
1176
//------------------------------------------------------------------------------
1305
 
int SNL_FEI_Structure::getMatrixRowLengths(feiArray<int>& rowLengths)
 
1177
int SNL_FEI_Structure::getMatrixRowLengths(std::vector<int>& rowLengths)
1306
1178
{
1307
1179
  if (!structureFinalized_) return(-1);
1308
1180
 
1309
1181
  rowLengths.resize(numLocalReducedRows_);
1310
1182
 
1311
 
  int* rowLenPtr = rowLengths.dataPtr();
1312
 
 
1313
1183
  for(int i=0; i<numLocalReducedRows_; i++) {
1314
 
    rowLenPtr[i] = sysMatIndices_[i].size();
 
1184
    rowLengths[i] = sysMatIndices_[i].size();
1315
1185
  }
1316
1186
  return(0);
1317
1187
}
1318
1188
 
1319
1189
//------------------------------------------------------------------------------
1320
1190
int SNL_FEI_Structure::getMatrixStructure(int** indices,
1321
 
                                         feiArray<int>& rowLengths)
 
1191
                                         std::vector<int>& rowLengths)
1322
1192
{
1323
1193
  if (!structureFinalized_) return(-1);
1324
1194
 
1325
1195
  rowLengths.resize(numLocalReducedRows_);
1326
1196
 
1327
 
  int* rowLenPtr = rowLengths.dataPtr();
1328
 
 
1329
1197
  for(int i=0; i<numLocalReducedRows_; i++) {
1330
 
    rowLenPtr[i] = sysMatIndices_[i].size();
1331
 
    fei::copySetToArray(sysMatIndices_[i], rowLenPtr[i], indices[i]);
 
1198
    rowLengths[i] = sysMatIndices_[i].size();
 
1199
    fei::copySetToArray(sysMatIndices_[i], rowLengths[i], indices[i]);
1332
1200
  }
1333
1201
 
1334
1202
  return(0);
1336
1204
 
1337
1205
//------------------------------------------------------------------------------
1338
1206
int SNL_FEI_Structure::getMatrixStructure(int** ptColIndices,
1339
 
                                         feiArray<int>& ptRowLengths,
 
1207
                                         std::vector<int>& ptRowLengths,
1340
1208
                                         int** blkColIndices,
1341
1209
                                          int* blkIndices_1D,
1342
 
                                         feiArray<int>& blkRowLengths,
1343
 
                                         feiArray<int>& numPtRowsPerBlkRow)
 
1210
                                         std::vector<int>& blkRowLengths,
 
1211
                                         std::vector<int>& numPtRowsPerBlkRow)
1344
1212
{
1345
1213
  int err = getMatrixStructure(ptColIndices, ptRowLengths);
1346
1214
  if (err != 0) return(-1);
1348
1216
  if (globalMaxBlkSize_ == 1) {
1349
1217
    //No block-equations have more than one point-equation, so we'll just assign
1350
1218
    //the block-structure arrays to be the same as the point-structure arrays.
1351
 
    int numRows = ptRowLengths.length();
 
1219
    int numRows = ptRowLengths.size();
1352
1220
    blkRowLengths.resize(numRows);
1353
 
    numPtRowsPerBlkRow.resize(numRows);
 
1221
    numPtRowsPerBlkRow.assign(numRows, 1);
1354
1222
 
1355
 
    //feiArray::operator=
1356
 
    blkRowLengths = ptRowLengths;
1357
 
    numPtRowsPerBlkRow = 1;
 
1223
    blkRowLengths.resize(ptRowLengths.size());
 
1224
    for(size_t ii=0; ii<ptRowLengths.size(); ++ii) {
 
1225
      blkRowLengths[ii] = ptRowLengths[ii];
 
1226
    }
1358
1227
 
1359
1228
    for(int i=0; i<numRows; i++) {
1360
1229
      blkColIndices[i] = ptColIndices[i];
1434
1303
}
1435
1304
 
1436
1305
//------------------------------------------------------------------------------
1437
 
bool SNL_FEI_Structure::nodalEqnsAllSlaves(NodeDescriptor* node,
1438
 
                                           feiArray<int>& slaveEqns)
 
1306
bool SNL_FEI_Structure::nodalEqnsAllSlaves(const NodeDescriptor* node,
 
1307
                                           std::vector<int>& slaveEqns)
1439
1308
{
1440
1309
  int numFields = node->getNumFields();
1441
1310
  const int* fieldIDs = node->getFieldIDList();
1456
1325
}
1457
1326
 
1458
1327
//------------------------------------------------------------------------------
1459
 
int SNL_FEI_Structure::initAccessPatternStructure()
1460
 
{
1461
 
  FEI_OSTREAM& os = dbgOut();
1462
 
  if (debugOutput_) os << "#   initAccessPatternStructure" << FEI_ENDL;
1463
 
 
1464
 
   int numPatterns = getNumPatterns();
1465
 
 
1466
 
   for(int i=0; i<numPatterns; i++) {
1467
 
     PatternDescriptor* pattern = NULL;
1468
 
     CHK_ERR( getPatternDescriptor_index(i, pattern));
1469
 
 
1470
 
      int patternID = pattern->getPatternID();
1471
 
      ConnectivityTable& conn = getPatternConnectivity(patternID);
1472
 
 
1473
 
      //The ConnectivityTable struct is meant to hold element-connectivities,
1474
 
      //but I'm also using it to hold pattern connectivities.
1475
 
 
1476
 
      feiArray<int> rowIndices(0, 8), colIndices(0,8);
1477
 
 
1478
 
      //The rows of the 'connectivities' table alternate, holding rowNodes and
1479
 
      //then colNodes.  ... confusing, yes. But...
1480
 
 
1481
 
      int offset = 0, loopCount = conn.numRows/2;
1482
 
      for(int row=0; row<loopCount; row++) {
1483
 
 
1484
 
         feiArray<GlobalID>& rownodes = *(conn.connectivities[offset++]);
1485
 
         feiArray<GlobalID>& colnodes = *(conn.connectivities[offset++]);
1486
 
 
1487
 
         feiArray<int> rowColOffsets(0, rownodes.length());
1488
 
         CHK_ERR( rowColOffsets.reAllocate(rownodes.length()) )
1489
 
 
1490
 
         int numColsPerRow;
1491
 
 
1492
 
         CHK_ERR( getPatternScatterIndices(patternID, rownodes.dataPtr(),
1493
 
                                           colnodes.dataPtr(), rowIndices,
1494
 
                                           rowColOffsets,numColsPerRow,
1495
 
                                           colIndices));
1496
 
 
1497
 
         if (rownodes.length() == 0 || colnodes.length() == 0) ERReturn(-1)
1498
 
 
1499
 
         int colIndicesPerRow = colIndices.length()/rownodes.length();
1500
 
 
1501
 
         SSGraph ssgraph(rowIndices.length(),
1502
 
                         rowIndices.dataPtr(),
1503
 
                         colIndicesPerRow,
1504
 
                         rowColOffsets.dataPtr(),
1505
 
                         colIndices.dataPtr());
1506
 
 
1507
 
         CHK_ERR( createEqnStructure(ssgraph) );
1508
 
      }
1509
 
   }
1510
 
 
1511
 
   return(FEI_SUCCESS);
1512
 
}
1513
 
 
1514
 
//------------------------------------------------------------------------------
1515
1328
NodeDescriptor& SNL_FEI_Structure::findNodeDescriptor(GlobalID nodeID)
1516
1329
{
1517
1330
//
1553
1366
  while(cr_iter != cr_end) {
1554
1367
    ConstraintType& multCR = *((*cr_iter).second);
1555
1368
 
1556
 
    int lenList = multCR.getMasters()->length();
 
1369
    int lenList = multCR.getMasters().size();
1557
1370
 
1558
 
    GlobalID *CRNodePtr = multCR.getMasters()->dataPtr();
1559
 
    int* CRFieldPtr = multCR.getMasterFieldIDs()->dataPtr();
 
1371
    std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
 
1372
    GlobalID *CRNodePtr = &CRNode_vec[0];
 
1373
    std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
 
1374
    int* CRFieldPtr = &CRField_vec[0];
1560
1375
 
1561
1376
    int crEqn = multCR.getEqnNumber();
1562
1377
    int reducedCrEqn = 0;
1633
1448
  while(cr_iter != cr_end) {
1634
1449
    ConstraintType& penCR = *((*cr_iter).second);
1635
1450
 
1636
 
    int lenList = penCR.getMasters()->length();
1637
 
    GlobalID* CRNodesPtr = penCR.getMasters()->dataPtr();
 
1451
    int lenList = penCR.getMasters().size();
 
1452
    std::vector<GlobalID>& CRNode_vec = penCR.getMasters();
 
1453
    GlobalID* CRNodesPtr = &CRNode_vec[0];
1638
1454
 
1639
 
    int* CRFieldPtr = penCR.getMasterFieldIDs()->dataPtr();
 
1455
    std::vector<int>& CRField_vec = penCR.getMasterFieldIDs();
 
1456
    int* CRFieldPtr = &CRField_vec[0];
1640
1457
 
1641
1458
    // each constraint equation generates a set of nodal energy terms, so
1642
1459
    // we have to process a matrix of nodes for each constraint
1842
1659
}
1843
1660
 
1844
1661
//------------------------------------------------------------------------------
1845
 
int SNL_FEI_Structure::storePatternScatterIndices(SSGraph& mat)
1846
 
{
1847
 
  //
1848
 
  //This function takes a "small" matrix-graph structure, and puts the indices
1849
 
  //either into the global sparse matrix structure, or into the EqnCommMgr if
1850
 
  //they represent an off-processor matrix location.
1851
 
  //
1852
 
 
1853
 
  int numRows = mat.getRows().length();
1854
 
  int* rows = mat.getRows().dataPtr();
1855
 
 
1856
 
  if (numRows == 0) return(FEI_SUCCESS);
1857
 
 
1858
 
  feiArray<feiArray<int>*>& indices = mat.getIndices();
1859
 
 
1860
 
  for(int i=0; i<numRows; i++) {
1861
 
    int row = rows[i];
1862
 
    int reducedRow = -1;
1863
 
    bool isSlave = translateToReducedEqn(row, reducedRow);
1864
 
    if (isSlave) ERReturn(-1);
1865
 
 
1866
 
    feiArray<int>& colIndices = *(indices[i]);
1867
 
 
1868
 
    workSpace_.resize(colIndices.length());
1869
 
    for(int jj=0; jj<colIndices.length(); jj++) {
1870
 
      isSlave = translateToReducedEqn(colIndices[jj], workSpace_[jj]);
1871
 
      if (isSlave) ERReturn(-1);
1872
 
    }
1873
 
 
1874
 
    if ((localStartRow_ > row) || (row > localEndRow_)) {
1875
 
      int owner = getOwnerProcForEqn(row);
1876
 
      eqnCommMgr_->addRemoteIndices(reducedRow, owner,
1877
 
                                    workSpace_.dataPtr(),
1878
 
                                    workSpace_.length());
1879
 
      continue;
1880
 
    }
1881
 
 
1882
 
    for(int j=0; j<workSpace_.length(); j++) {
1883
 
      CHK_ERR( createMatrixPosition(reducedRow, workSpace_[j],
1884
 
                                    "storePatScttrInd") );
1885
 
    }
1886
 
  }
1887
 
  return(FEI_SUCCESS);
1888
 
}
1889
 
 
1890
 
//------------------------------------------------------------------------------
1891
 
int SNL_FEI_Structure::storePatternScatterIndices_noSlaves(SSGraph& mat)
1892
 
{
1893
 
  //
1894
 
  //This function takes a "small" matrix-graph structure, and puts the indices
1895
 
  //either into the global sparse matrix structure, or into the EqnCommMgr if
1896
 
  //they represent an off-processor matrix location.
1897
 
  //
1898
 
 
1899
 
  int numRows = mat.getRows().length();
1900
 
  int* rows = mat.getRows().dataPtr();
1901
 
 
1902
 
  if (numRows == 0) return(FEI_SUCCESS);
1903
 
 
1904
 
  feiArray<feiArray<int>*>& indices = mat.getIndices();
1905
 
 
1906
 
  for(int i=0; i<numRows; i++) {
1907
 
    int row = rows[i];
1908
 
 
1909
 
    feiArray<int>& colIndices = *(indices[i]);
1910
 
 
1911
 
    if ((localStartRow_ > row) || (row > localEndRow_)) {
1912
 
      int owner = getOwnerProcForEqn(row);
1913
 
      eqnCommMgr_->addRemoteIndices(row, owner,
1914
 
                                    colIndices.dataPtr(),
1915
 
                                    colIndices.length());
1916
 
      continue;
1917
 
    }
1918
 
 
1919
 
    int* colIndPtr = colIndices.dataPtr();
1920
 
    for(int j=0; j<colIndices.length(); j++) {
1921
 
      CHK_ERR( createMatrixPosition(row, colIndPtr[j],
1922
 
                                    "storePatScttrInd_noSlvs") );
1923
 
    }
1924
 
  }
1925
 
 
1926
 
  return(FEI_SUCCESS);
1927
 
}
1928
 
 
1929
 
//------------------------------------------------------------------------------
1930
 
int SNL_FEI_Structure::createSymmEqnStructure(feiArray<int>& scatterIndices)
 
1662
int SNL_FEI_Structure::createSymmEqnStructure(std::vector<int>& scatterIndices)
1931
1663
{
1932
1664
  //scatterIndices are both the rows and the columns for a structurally
1933
1665
  //symmetric 2-dimensional contribution to the matrix.
1942
1674
 
1943
1675
  try {
1944
1676
 
1945
 
  int len = scatterIndices.length();
 
1677
  int len = scatterIndices.size();
1946
1678
  bool anySlaves = false;
1947
1679
  rSlave_.resize(len);
1948
1680
  for(int is=0; is<len; is++) { 
1949
 
    rSlave_[is] = isSlaveEqn(scatterIndices[is]);
1950
 
    if (rSlave_[is]) anySlaves = true;
 
1681
    rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
 
1682
    if (rSlave_[is] == 1) anySlaves = true;
1951
1683
  }
1952
1684
 
1953
1685
  //if there aren't any slaves in this contribution, then just store it
1957
1689
    return(FEI_SUCCESS);
1958
1690
  }
1959
1691
 
1960
 
  int* scatterPtr = scatterIndices.dataPtr();
 
1692
  int* scatterPtr = &scatterIndices[0];
1961
1693
 
1962
1694
  workSpace_.resize(len);
1963
1695
  for(int j=0; j<len; j++) {
1966
1698
  
1967
1699
  for(int i=0; i<len; i++) {
1968
1700
    int row = scatterPtr[i];
1969
 
    if (rSlave_[i]) {
 
1701
    if (rSlave_[i] == 1) {
1970
1702
      reducedEqnCounter_++;
1971
1703
      //
1972
1704
      //'row' is a slave equation, so add this row to Kdi_. But as we do that,
1975
1707
      for(int jj=0; jj<len; jj++) {
1976
1708
        int col = scatterPtr[jj];
1977
1709
 
1978
 
        if (rSlave_[jj]) {
 
1710
        if (rSlave_[jj] == 1) {
1979
1711
          //'col' is a slave equation, so add this column to Kid_.
1980
1712
          for(int ii=0; ii<len; ii++) {
1981
1713
            int rowi = scatterPtr[ii];
1982
1714
 
1983
1715
            //only add the non-slave rows for this column.
1984
 
            if (rSlave_[ii]) continue;
 
1716
            if (rSlave_[ii] == 1) continue;
1985
1717
 
1986
1718
            Kid_->createPosition(rowi, col);
1987
1719
          }
1996
1728
      for(int kk=0; kk<len; kk++) {
1997
1729
        int colk = scatterPtr[kk];
1998
1730
 
1999
 
        if (!rSlave_[kk]) continue;
 
1731
        if (rSlave_[kk] != 1) continue;
2000
1732
 
2001
1733
        Kdd_->createPosition(row, colk);
2002
1734
      }
2015
1747
      }
2016
1748
 
2017
1749
      for(int j=0; j<len; j++) {
2018
 
        if (rSlave_[j]) continue;
 
1750
        if (rSlave_[j] == 1) continue;
2019
1751
 
2020
1752
        int reducedCol = workSpace_[j];
2021
1753
 
2033
1765
  if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
2034
1766
 
2035
1767
  }
2036
 
  catch(fei::Exception& exc) {
 
1768
  catch(std::runtime_error& exc) {
2037
1769
    FEI_CERR << exc.what() << FEI_ENDL;
2038
1770
    ERReturn(-1);
2039
1771
  }
2042
1774
}
2043
1775
 
2044
1776
//------------------------------------------------------------------------------
2045
 
int SNL_FEI_Structure::createBlkSymmEqnStructure(feiArray<int>& scatterIndices)
 
1777
int SNL_FEI_Structure::createBlkSymmEqnStructure(std::vector<int>& scatterIndices)
2046
1778
{
2047
1779
  //scatterIndices are both the rows and the columns for a structurally
2048
1780
  //symmetric 2-dimensional contribution to the matrix.
2056
1788
 
2057
1789
  try {
2058
1790
 
2059
 
  int len = scatterIndices.length();
 
1791
  int len = scatterIndices.size();
2060
1792
  bool anySlaves = false;
2061
1793
  rSlave_.resize(len);
2062
1794
  for(int is=0; is<len; is++) { 
2063
 
    rSlave_[is] = isSlaveEqn(scatterIndices[is]);
2064
 
    if (rSlave_[is]) anySlaves = true;
 
1795
    rSlave_[is] = isSlaveEqn(scatterIndices[is]) ? 1 : 0;
 
1796
    if (rSlave_[is] == 1) anySlaves = true;
2065
1797
  }
2066
1798
 
2067
1799
  //if there aren't any slaves in this contribution, then just store it
2071
1803
    return(FEI_SUCCESS);
2072
1804
  }
2073
1805
 
2074
 
  int* scatterPtr = scatterIndices.dataPtr();
 
1806
  int* scatterPtr = &scatterIndices[0];
2075
1807
 
2076
1808
  workSpace_.resize(len);
2077
1809
  for(int j=0; j<len; j++) {
2080
1812
  
2081
1813
  for(int i=0; i<len; i++) {
2082
1814
    int row = scatterPtr[i];
2083
 
    if (rSlave_[i]) {
 
1815
    if (rSlave_[i] == 1) {
2084
1816
      reducedEqnCounter_++;
2085
1817
      //
2086
1818
      //'row' is a slave equation, so add this row to Kdi_. But as we do that,
2089
1821
      for(int jj=0; jj<len; jj++) {
2090
1822
        int col = scatterPtr[jj];
2091
1823
 
2092
 
        if (rSlave_[jj]) {
 
1824
        if (rSlave_[jj] == 1) {
2093
1825
          //'col' is a slave equation, so add this column to Kid_.
2094
1826
          for(int ii=0; ii<len; ii++) {
2095
1827
            int rowi = scatterPtr[ii];
2096
1828
 
2097
1829
            //only add the non-slave rows for this column.
2098
 
            if (rSlave_[ii]) continue;
 
1830
            if (rSlave_[ii] == 1) continue;
2099
1831
 
2100
1832
            Kid_->createPosition(rowi, col);
2101
1833
          }
2110
1842
      for(int kk=0; kk<len; kk++) {
2111
1843
        int colk = scatterPtr[kk];
2112
1844
 
2113
 
        if (!rSlave_[kk]) continue;
 
1845
        if (rSlave_[kk] != 1) continue;
2114
1846
 
2115
1847
        Kdd_->createPosition(row, colk);
2116
1848
      }
2129
1861
      }
2130
1862
 
2131
1863
      for(int j=0; j<len; j++) {
2132
 
        if (rSlave_[j]) continue;
 
1864
        if (rSlave_[j] == 1) continue;
2133
1865
 
2134
1866
        int reducedCol = workSpace_[j];
2135
1867
 
2147
1879
  if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
2148
1880
 
2149
1881
  }
2150
 
  catch(fei::Exception& exc) {
 
1882
  catch(std::runtime_error& exc) {
2151
1883
    FEI_CERR << exc.what() << FEI_ENDL;
2152
1884
    ERReturn(-1);
2153
1885
  }
2156
1888
}
2157
1889
 
2158
1890
//------------------------------------------------------------------------------
2159
 
int SNL_FEI_Structure::storeElementScatterIndices(feiArray<int>& indices)
 
1891
int SNL_FEI_Structure::storeElementScatterIndices(std::vector<int>& indices)
2160
1892
{
2161
1893
//This function takes a list of equation numbers, as input, and for each
2162
1894
//one, if it is a local equation number, goes to that row of the sysMatIndices
2164
1896
//in that row of the matrix structure. If the equation number is not local,
2165
1897
//then it is given to the EqnCommMgr.
2166
1898
//
2167
 
   int numIndices = indices.length();
2168
 
   int* indPtr = indices.dataPtr();
 
1899
   int numIndices = indices.size();
 
1900
   int* indPtr = &indices[0];
2169
1901
 
2170
1902
   workSpace_.resize(numIndices);
2171
 
   int* workPtr = workSpace_.dataPtr();
 
1903
   int* workPtr = &workSpace_[0];
2172
1904
   int reducedEqn = -1;
2173
 
   for(int j=0; j<indices.length(); j++) {
 
1905
   for(size_t j=0; j<indices.size(); j++) {
2174
1906
     bool isSlave = translateToReducedEqn(indPtr[j], reducedEqn);
2175
1907
     if (isSlave) ERReturn(-1);
2176
1908
     workPtr[j] = reducedEqn;
2195
1927
}
2196
1928
 
2197
1929
//------------------------------------------------------------------------------
2198
 
int SNL_FEI_Structure::storeElementScatterIndices_noSlaves(feiArray<int>& indices)
2199
 
{
2200
 
  //This function takes a list of equation numbers as input and for each one,
2201
 
  //if it is a local equation number, goes to that row of the sysMatIndices
2202
 
  //structure and stores the whole list of equation numbers as column indices
2203
 
  //in that row of the matrix structure. If the equation number is not local,
2204
 
  //then it is given to the EqnCommMgr.
2205
 
  //
2206
 
  int i, numIndices = indices.length();
2207
 
  int* indPtr = indices.dataPtr();
2208
 
 
2209
 
  for(i=0; i<numIndices; i++) {
2210
 
    int row = indPtr[i];
2211
 
    if (row < localStartRow_ || row > localEndRow_) {
2212
 
      int owner = getOwnerProcForEqn(row);
2213
 
      eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
2214
 
      continue;
2215
 
    }
2216
 
 
2217
 
    if (generateGraph_) {
2218
 
      fei::ctg_set<int>& thisRow =
2219
 
        sysMatIndices_[row - reducedStartRow_];
2220
 
 
2221
 
      for(int j=0; j<numIndices; ++j) {
2222
 
        if (debugOutput_ && outputLevel_ > 2) {
2223
 
          dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
2224
 
                   << row << "," << indPtr[j] << ")"<<FEI_ENDL;
2225
 
        }
2226
 
 
2227
 
        thisRow.insert2(indPtr[j]);
2228
 
      }
2229
 
    }
2230
 
  }
2231
 
 
2232
 
  return(FEI_SUCCESS);
2233
 
}
2234
 
 
2235
 
//------------------------------------------------------------------------------
2236
 
int SNL_FEI_Structure::storeElementScatterBlkIndices_noSlaves(feiArray<int>& indices)
2237
 
{
2238
 
  //This function takes a list of equation numbers as input and for each one,
2239
 
  //if it is a local equation number, goes to that row of the sysMatIndices
2240
 
  //structure and stores the whole list of equation numbers as column indices
2241
 
  //in that row of the matrix structure. If the equation number is not local,
2242
 
  //then it is given to the EqnCommMgr.
2243
 
  //
2244
 
  int i, numIndices = indices.length();
2245
 
  int* indPtr = indices.dataPtr();
2246
 
 
2247
 
  for(i=0; i<numIndices; i++) {
2248
 
    int row = indPtr[i];
2249
 
    if (row < localStartRow_ || row > localEndRow_) {
2250
 
      int owner = getOwnerProcForEqn(row);
2251
 
      eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
2252
 
      continue;
2253
 
    }
2254
 
 
2255
 
    if (generateGraph_) {
2256
 
      fei::ctg_set<int>& thisRow =
2257
 
        sysMatIndices_[row - reducedStartRow_];
2258
 
 
2259
 
      for(int j=0; j<numIndices; ++j) {
2260
 
        if (debugOutput_ && outputLevel_ > 2) {
2261
 
          dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
2262
 
                   << row << "," << indPtr[j] << ")"<<FEI_ENDL;
2263
 
        }
2264
 
 
2265
 
        thisRow.insert2(indPtr[j]);
2266
 
      }
2267
 
    }
2268
 
  }
2269
 
 
2270
 
  return(FEI_SUCCESS);
2271
 
}
2272
 
 
2273
 
//------------------------------------------------------------------------------
2274
 
int SNL_FEI_Structure::createEqnStructure(SSGraph& mat)
2275
 
{
2276
 
  //This function creates the global matrix structure positions associated with
2277
 
  //the "small" matrix-structure in 'mat'.
2278
 
  //If there are no slave equations in either the whole problem or in 'mat',
2279
 
  //then mat is simply handed to 'storePatternScatterIndices' which either
2280
 
  //creates local matrix positions or makes insertions to EqnCommMgr, based on
2281
 
  //whether the locations represented in 'mat' are local or not.
2282
 
  //
2283
 
  //If there are slave equations present in 'mat', then its contents are
2284
 
  //distributed into the locally-held submatrices Kdi, Kid, Kdd which will
2285
 
  //later be assembled into the global matrix structure by the function
2286
 
  //'assembleReducedStructure'.
2287
 
  //
2288
 
  int numRows = mat.getRows().length();
2289
 
  int* rows = mat.getRows().dataPtr();
2290
 
 
2291
 
  if (numRows == 0) return(FEI_SUCCESS);
2292
 
 
2293
 
  if (numSlaveEquations() == 0) {
2294
 
    CHK_ERR( storePatternScatterIndices_noSlaves(mat) );
2295
 
    return(FEI_SUCCESS);
2296
 
  }
2297
 
 
2298
 
  feiArray<feiArray<int>*>& indices = mat.getIndices();
2299
 
 
2300
 
  rSlave_.resize(numRows);
2301
 
  cSlave_.resize(0);
2302
 
  bool anySlaves = false;
2303
 
  for(int r=0; r<numRows; r++) {
2304
 
    rSlave_[r] = isSlaveEqn(rows[r]);
2305
 
    const int* indPtr = indices[r]->dataPtr();
2306
 
    if (rSlave_[r]) anySlaves = true;
2307
 
 
2308
 
    for(int j=0; j<indices[r]->length(); j++) {
2309
 
      cSlave_.append(isSlaveEqn(indPtr[j]));
2310
 
      if (cSlave_[cSlave_.length()-1]) anySlaves = true;
2311
 
    }
2312
 
  }
2313
 
 
2314
 
  if (!anySlaves) {
2315
 
    CHK_ERR( storePatternScatterIndices(mat) );
2316
 
    return(FEI_SUCCESS);
2317
 
  }
2318
 
 
2319
 
  int offset = 0;
2320
 
  for(int i=0; i<numRows; i++) {
2321
 
    int row = rows[i];
2322
 
 
2323
 
    int numCols = indices[i]->length();
2324
 
    const int* indicesPtr = indices[i]->dataPtr();
2325
 
    bool* colSlave = cSlave_.dataPtr() + offset;
2326
 
    offset += numCols;
2327
 
 
2328
 
    if (rSlave_[i]) {
2329
 
      //Since this is a slave equation, the non-slave columns of this row go
2330
 
      //into 'Kdi_', and the slave columns go into 'Kdd_'.
2331
 
      for(int jj=0; jj<numCols; jj++) {
2332
 
        int col = indicesPtr[jj];
2333
 
        if (colSlave[jj]) {
2334
 
          Kdd_->createPosition(row, col);
2335
 
        }
2336
 
        else {
2337
 
          Kdi_->createPosition(row, col);
2338
 
        }
2339
 
      }
2340
 
 
2341
 
      //We also need to put the non-slave rows of column 'row' into 'K_id'.
2342
 
      for(int ii=0; ii<numRows; ii++) {
2343
 
        int rowi = rows[ii];
2344
 
        if (rSlave_[ii]) continue;
2345
 
 
2346
 
        int index = indices[ii]->find(row);
2347
 
        if (index < 0) continue;
2348
 
 
2349
 
        Kid_->createPosition(rowi, row);
2350
 
      }
2351
 
 
2352
 
      reducedEqnCounter_++;
2353
 
 
2354
 
      continue;
2355
 
    }
2356
 
    else {//row is not a slave eqn...
2357
 
      int reducedRow = -1;
2358
 
      bool isSlave = translateToReducedEqn(row, reducedRow);
2359
 
      if (isSlave) ERReturn(-1);
2360
 
 
2361
 
      bool rowIsLocal = true;
2362
 
      if (localStartRow_ > row || row > localEndRow_) {
2363
 
        rowIsLocal = false;
2364
 
      }
2365
 
 
2366
 
      //put all non-slave columns from this row into the assembled matrix.
2367
 
      for(int j=0; j<numCols; j++) {
2368
 
        int col = indicesPtr[j];
2369
 
 
2370
 
        if (colSlave[j]) continue;
2371
 
 
2372
 
        int reducedCol = -1;
2373
 
        isSlave = translateToReducedEqn(col, reducedCol);
2374
 
        if (isSlave) ERReturn(-1);
2375
 
 
2376
 
        if (rowIsLocal) {
2377
 
          CHK_ERR( createMatrixPosition(reducedRow, reducedCol,
2378
 
                                        "crtEqnStr") );
2379
 
        }
2380
 
        else {
2381
 
          int owner = getOwnerProcForEqn(row);
2382
 
          if (owner < 0) {
2383
 
            FEI_CERR << "SNL_FEI_Structure ERROR, owner-proc not found for row "
2384
 
                 << row << FEI_ENDL;
2385
 
            ERReturn(-1);
2386
 
          }
2387
 
 
2388
 
          eqnCommMgr_->addRemoteIndices(reducedRow, owner, &reducedCol, 1);
2389
 
        }
2390
 
      }
2391
 
    }
2392
 
  }
2393
 
 
2394
 
  if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedStructure() );
 
1930
int SNL_FEI_Structure::storeElementScatterIndices_noSlaves(std::vector<int>& indices)
 
1931
{
 
1932
  //This function takes a list of equation numbers as input and for each one,
 
1933
  //if it is a local equation number, goes to that row of the sysMatIndices
 
1934
  //structure and stores the whole list of equation numbers as column indices
 
1935
  //in that row of the matrix structure. If the equation number is not local,
 
1936
  //then it is given to the EqnCommMgr.
 
1937
  //
 
1938
  int i, numIndices = indices.size();
 
1939
  int* indPtr = &indices[0];
 
1940
 
 
1941
  for(i=0; i<numIndices; i++) {
 
1942
    int row = indPtr[i];
 
1943
    if (row < localStartRow_ || row > localEndRow_) {
 
1944
      int owner = getOwnerProcForEqn(row);
 
1945
      eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
 
1946
      continue;
 
1947
    }
 
1948
 
 
1949
    if (generateGraph_) {
 
1950
      fei::ctg_set<int>& thisRow =
 
1951
        sysMatIndices_[row - reducedStartRow_];
 
1952
 
 
1953
      for(int j=0; j<numIndices; ++j) {
 
1954
        if (debugOutput_ && outputLevel_ > 2) {
 
1955
          dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
 
1956
                   << row << "," << indPtr[j] << ")"<<FEI_ENDL;
 
1957
        }
 
1958
 
 
1959
        thisRow.insert2(indPtr[j]);
 
1960
      }
 
1961
    }
 
1962
  }
 
1963
 
 
1964
  return(FEI_SUCCESS);
 
1965
}
 
1966
 
 
1967
//------------------------------------------------------------------------------
 
1968
int SNL_FEI_Structure::storeElementScatterBlkIndices_noSlaves(std::vector<int>& indices)
 
1969
{
 
1970
  //This function takes a list of equation numbers as input and for each one,
 
1971
  //if it is a local equation number, goes to that row of the sysMatIndices
 
1972
  //structure and stores the whole list of equation numbers as column indices
 
1973
  //in that row of the matrix structure. If the equation number is not local,
 
1974
  //then it is given to the EqnCommMgr.
 
1975
  //
 
1976
  int i, numIndices = indices.size();
 
1977
  int* indPtr = &indices[0];
 
1978
 
 
1979
  for(i=0; i<numIndices; i++) {
 
1980
    int row = indPtr[i];
 
1981
    if (row < localStartRow_ || row > localEndRow_) {
 
1982
      int owner = getOwnerProcForEqn(row);
 
1983
      eqnCommMgr_->addRemoteIndices(row, owner, indPtr, numIndices);
 
1984
      continue;
 
1985
    }
 
1986
 
 
1987
    if (generateGraph_) {
 
1988
      fei::ctg_set<int>& thisRow =
 
1989
        sysMatIndices_[row - reducedStartRow_];
 
1990
 
 
1991
      for(int j=0; j<numIndices; ++j) {
 
1992
        if (debugOutput_ && outputLevel_ > 2) {
 
1993
          dbgOut() << "#       storeElemScttrInds_ns crtMatPoss("
 
1994
                   << row << "," << indPtr[j] << ")"<<FEI_ENDL;
 
1995
        }
 
1996
 
 
1997
        thisRow.insert2(indPtr[j]);
 
1998
      }
 
1999
    }
 
2000
  }
2395
2001
 
2396
2002
  return(FEI_SUCCESS);
2397
2003
}
2405
2011
  //resulting contributions to off-processor portions of the matrix structure to
2406
2012
  //the EqnCommMgr.
2407
2013
  //
2408
 
  SSMat* D = getSlaveDependencies();
2409
 
 
2410
 
  //form tmpMat1_ = Kid*D
2411
 
  CHK_ERR( Kid_->matMat(*D, *tmpMat1_) );
2412
 
 
2413
 
  //form tmpMat2_ = D^T*Kdi
2414
 
  CHK_ERR( D->matTransMat(*Kdi_, *tmpMat2_) );
2415
 
 
2416
 
  CHK_ERR( translateMatToReducedEqns(*tmpMat1_) );
2417
 
  CHK_ERR( translateMatToReducedEqns(*tmpMat2_) );
2418
 
  if (numProcs_ > 1) {
2419
 
    CHK_ERR( eqnCommMgr_->addRemoteEqns(*tmpMat1_, true) );
2420
 
    CHK_ERR( eqnCommMgr_->addRemoteEqns(*tmpMat2_, true) );
2421
 
  }
2422
 
 
2423
 
  CHK_ERR( createMatrixPositions(*tmpMat1_) );
2424
 
  CHK_ERR( createMatrixPositions(*tmpMat2_) );
2425
 
 
2426
 
  //form tmpMat1_ = D^T*Kdd
2427
 
  CHK_ERR( D->matTransMat(*Kdd_, *tmpMat1_) );
2428
 
 
2429
 
  //form tmpMat2_ = tpmMat1_*D = D^T*Kdd*D
2430
 
  CHK_ERR( tmpMat1_->matMat(*D, *tmpMat2_) );
2431
 
 
2432
 
  CHK_ERR( translateMatToReducedEqns(*tmpMat2_) );
2433
 
  if (numProcs_ > 1) {
2434
 
    CHK_ERR( eqnCommMgr_->addRemoteEqns(*tmpMat2_, true) );
2435
 
  }
2436
 
  CHK_ERR( createMatrixPositions(*tmpMat2_) );
2437
 
 
2438
 
  Kdi_->logicalClear();
2439
 
  Kid_->logicalClear();
2440
 
  Kdd_->logicalClear();
 
2014
  fei::FillableMat* D = getSlaveDependencies();
 
2015
 
 
2016
  csrD = *D;
 
2017
  csrKid = *Kid_;
 
2018
  csrKdi = *Kdi_;
 
2019
  csrKdd = *Kdd_;
 
2020
 
 
2021
  //form tmpMat1_ = Kid_*D
 
2022
  fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
 
2023
 
 
2024
  //form tmpMat2_ = D^T*Kdi_
 
2025
  fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
 
2026
 
 
2027
  CHK_ERR( translateMatToReducedEqns(tmpMat1_) );
 
2028
  CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
 
2029
  if (numProcs_ > 1) {
 
2030
    CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat1_, true) );
 
2031
    CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
 
2032
  }
 
2033
 
 
2034
  CHK_ERR( createMatrixPositions(tmpMat1_) );
 
2035
  CHK_ERR( createMatrixPositions(tmpMat2_) );
 
2036
 
 
2037
  //form tmpMat1_ = D^T*Kdd_
 
2038
  fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
 
2039
 
 
2040
  //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
 
2041
  fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
 
2042
 
 
2043
 
 
2044
  CHK_ERR( translateMatToReducedEqns(tmpMat2_) );
 
2045
  if (numProcs_ > 1) {
 
2046
    CHK_ERR( eqnCommMgr_->addRemoteEqns(tmpMat2_, true) );
 
2047
  }
 
2048
  CHK_ERR( createMatrixPositions(tmpMat2_) );
 
2049
 
 
2050
  Kdi_->clear();
 
2051
  Kid_->clear();
 
2052
  Kdd_->clear();
2441
2053
  reducedEqnCounter_ = 0;
2442
2054
 
2443
2055
  return(FEI_SUCCESS);
2460
2072
int SNL_FEI_Structure::translateToReducedEqns(EqnBuffer& eqnBuf)
2461
2073
{
2462
2074
  int numEqns = eqnBuf.getNumEqns();
2463
 
  int* eqnNumbers = eqnBuf.eqnNumbersPtr().dataPtr();
2464
 
  feiArray<SSVec*>& eqnArray = eqnBuf.eqns();
 
2075
  int* eqnNumbers = &(eqnBuf.eqnNumbers()[0]);
 
2076
  std::vector<fei::CSVec*>& eqnArray = eqnBuf.eqns();
2465
2077
  for(int i=0; i<numEqns; ++i) {
2466
2078
    int reducedEqn;
2467
2079
    translateToReducedEqn(eqnNumbers[i], reducedEqn);
2468
2080
    eqnNumbers[i] = reducedEqn;
2469
2081
 
2470
 
    int* indicesPtr = eqnArray[i]->indices().dataPtr();
2471
 
    int numIndices = eqnArray[i]->length();
 
2082
    int* indicesPtr = &(eqnArray[i]->indices()[0]);
 
2083
    int numIndices = eqnArray[i]->size();
2472
2084
    for(int j=0; j<numIndices; ++j) {
2473
2085
      translateToReducedEqn(indicesPtr[j], reducedEqn);
2474
2086
      indicesPtr[j] = reducedEqn;
2497
2109
}
2498
2110
 
2499
2111
//------------------------------------------------------------------------------
2500
 
int SNL_FEI_Structure::translateMatToReducedEqns(SSMat& mat)
 
2112
int SNL_FEI_Structure::translateMatToReducedEqns(fei::CSRMat& mat)
2501
2113
{
2502
 
  //Given a matrix in global numbering, convert all of its contents to the
 
2114
  //Given a matrix in unreduced numbering, convert all of its contents to the
2503
2115
  //"reduced" equation space. If any of the row-numbers or column-indices in
2504
2116
  //this matrix object are slave equations, they will not be referenced. In
2505
2117
  //this case, a positive warning-code will be returned.
2506
2118
 
2507
 
  feiArray<int>& rowNumbers = mat.getRowNumbers();
2508
 
  feiArray<SSVec*>& rows = mat.getRows();
2509
 
 
2510
2119
  int reducedEqn = -1;
2511
2120
  int foundSlave = 0;
2512
2121
 
2513
 
  for(int i=0; i<rowNumbers.length(); i++) {
 
2122
  fei::SparseRowGraph& srg = mat.getGraph();
 
2123
 
 
2124
  std::vector<int>& rowNumbers = srg.rowNumbers;
 
2125
  for(size_t i=0; i<rowNumbers.size(); ++i) {
2514
2126
    bool isSlave = translateToReducedEqn(rowNumbers[i], reducedEqn);
2515
2127
    if (isSlave) foundSlave = 1;
2516
2128
    else rowNumbers[i] = reducedEqn;
 
2129
  }
2517
2130
 
2518
 
    feiArray<int>& colInds = rows[i]->indices();
2519
 
    for(int j=0; j<colInds.length(); j++) {
2520
 
      bool isSlave = translateToReducedEqn(colInds[j], reducedEqn);
2521
 
      if (isSlave) foundSlave = 1;
2522
 
      else colInds[j] = reducedEqn;
2523
 
    }
 
2131
  std::vector<int>& colIndices = srg.packedColumnIndices;
 
2132
  for(size_t i=0; i<colIndices.size(); ++i) {
 
2133
    bool isSlave = translateToReducedEqn(colIndices[i], reducedEqn);
 
2134
    if (isSlave) foundSlave = 1;
 
2135
    else colIndices[i] = reducedEqn;
2524
2136
  }
2525
2137
 
2526
2138
  return(foundSlave);
2527
2139
}
2528
2140
 
2529
2141
//------------------------------------------------------------------------------
2530
 
int SNL_FEI_Structure::createMatrixPositions(SSMat& mat)
 
2142
int SNL_FEI_Structure::createMatrixPositions(fei::CSRMat& mat)
2531
2143
{
2532
2144
  if (!generateGraph_) {
2533
2145
    return(0);
2534
2146
  }
2535
2147
 
2536
 
  //This function must be called with a SSMat object that already has its
 
2148
  //This function must be called with a CSRMat object that already has its
2537
2149
  //contents numbered in "reduced" equation-numbers.
2538
2150
  //
2539
2151
  //This function has one simple task -- for each row,col pair stored in 'mat',
2540
2152
  //call 'createMatrixPosition' to make an entry in the global matrix structure
2541
2153
  //if it doesn't already exist.
2542
2154
  //
2543
 
  int numRows = mat.getRowNumbers().length();
2544
 
  int* rowNumbers = mat.getRowNumbers().dataPtr();
2545
 
  feiArray<SSVec*>& rows = mat.getRows();
2546
 
 
2547
 
  for(int i=0; i<numRows; i++) {
2548
 
    int* indicesRow = rows[i]->indices().dataPtr();
2549
 
 
2550
 
    for(int j=0; j<rows[i]->length(); j++) {
2551
 
      CHK_ERR( createMatrixPosition(rowNumbers[i], indicesRow[j],
2552
 
                                    "crtMatPos(SSMat)") );
 
2155
  const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
 
2156
  const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
 
2157
  const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
 
2158
 
 
2159
  for(size_t i=0; i<rowNumbers.size(); ++i) {
 
2160
    int row = rowNumbers[i];
 
2161
    int offset = rowOffsets[i];
 
2162
    int rowlen = rowOffsets[i+1]-offset;
 
2163
    const int* indices = &pckColInds[offset];
 
2164
 
 
2165
    for(int j=0; j<rowlen; j++) {
 
2166
      CHK_ERR( createMatrixPosition(row, indices[j], "crtMatPos(CSRMat)") );
2553
2167
    }
2554
2168
  }
2555
2169
 
2669
2283
    return(0);
2670
2284
  }
2671
2285
 
2672
 
  int i, numNodes = getNumActiveNodes();
2673
 
 
2674
2286
  int localVanishedNodeAdjustment = 0;
2675
 
  for(i=0; i<localProc_; ++i) {
 
2287
  for(int i=0; i<localProc_; ++i) {
2676
2288
    localVanishedNodeAdjustment += globalNumNodesVanished_[i];
2677
2289
  }
2678
2290
 
2679
2291
  int eqnNumber, blkEqnNumber;
2680
2292
 
2681
 
  for(i=0; i<numNodes; i++) {
 
2293
  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
 
2294
  std::map<GlobalID,int>::iterator
 
2295
    iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
 
2296
 
 
2297
  for(; iter!=iter_end; ++iter) {
2682
2298
    NodeDescriptor* node = NULL;
2683
 
    CHK_ERR( nodeDatabase_->getNodeAtIndex(i, node) );
 
2299
    CHK_ERR( nodeDatabase_->getNodeAtIndex(iter->second, node) );
2684
2300
 
2685
2301
    int firstPtEqn = node->getFieldEqnNumbers()[0];
2686
2302
    int numNodalDOF = node->getNumNodalDOF();
2710
2326
  //Now the element-dofs...
2711
2327
  //
2712
2328
  int numBlocks = getNumElemBlocks();
2713
 
  int localProc = commUtilsInt_->localProc();
2714
 
  int numLocalNodes = globalNodeOffsets_[localProc+1]-globalNodeOffsets_[localProc];
 
2329
  int numLocalNodes = globalNodeOffsets_[localProc_+1]-globalNodeOffsets_[localProc_];
2715
2330
  eqnNumber = localStartRow_ + numLocalNodalEqns_;
2716
2331
  blkEqnNumber = localBlkOffset_ + numLocalNodes;
2717
2332
 
2718
 
  for(i = 0; i < numBlocks; i++) {
 
2333
  for(int i = 0; i < numBlocks; i++) {
2719
2334
    BlockDescriptor* block = NULL;
2720
2335
    CHK_ERR( getBlockDescriptor_index(i, block) );
2721
2336
 
2767
2382
 
2768
2383
  int localMaxBlkEqnSize = blkEqnMapper_->getMaxBlkEqnSize();
2769
2384
  globalMaxBlkSize_ = 0;
2770
 
  CHK_ERR( commUtilsInt_->GlobalMax(localMaxBlkEqnSize, globalMaxBlkSize_) );
 
2385
  CHK_ERR( fei::GlobalMax(comm_, localMaxBlkEqnSize, globalMaxBlkSize_) );
2771
2386
 
2772
2387
  blkEqnMapper_->setMaxBlkEqnSize(globalMaxBlkSize_);
2773
2388
 
2780
2395
 
2781
2396
  if (numProcs_ > 1) {
2782
2397
    std::vector<int> recvLengths, globalMultEqns;
2783
 
    CHK_ERR(commUtilsInt_->Allgatherv(localMultEqns, recvLengths,
2784
 
                                      globalMultEqns));
 
2398
    CHK_ERR(fei::Allgatherv(comm_, localMultEqns, recvLengths, globalMultEqns));
2785
2399
 
2786
2400
    int offset = 0;
2787
2401
    for(int p=0; p<numProcs_; p++) {
2855
2469
    ++cr_iter;
2856
2470
  }
2857
2471
 
2858
 
  int numNodes = getNumActiveNodes();
 
2472
  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
 
2473
  std::map<GlobalID,int>::iterator
 
2474
    iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
2859
2475
 
2860
 
  for(int i=0; i<numNodes; i++) {
 
2476
  for(; iter!=iter_end; ++iter) {
2861
2477
    NodeDescriptor* node = NULL;
2862
 
    CHK_ERR( nodeDatabase_->getNodeAtIndex(i, node) );
 
2478
    CHK_ERR( nodeDatabase_->getNodeAtIndex(iter->second, node) );
2863
2479
 
2864
2480
    if (node->getOwnerProc() != localProc_) {
2865
2481
      continue;
2933
2549
 
2934
2550
  //first, get each processor's local number of equations on the master proc.
2935
2551
 
2936
 
  feiArray<int> numProcNodes(numProcs_);
2937
 
  feiArray<int> numProcEqns(numProcs_);
2938
 
  feiArray<int> numProcEqnBlks(numProcs_);
 
2552
  std::vector<int> numProcNodes(numProcs_);
 
2553
  std::vector<int> numProcEqns(numProcs_);
 
2554
  std::vector<int> numProcEqnBlks(numProcs_);
2939
2555
 
2940
2556
  if (numProcs_ > 1) {
2941
2557
#ifndef FEI_SER
2942
 
    feiArray<int> glist(3);
2943
 
    feiArray<int> recvList(3*numProcs_);
 
2558
    int glist[3];
 
2559
    std::vector<int> recvList(3*numProcs_);
2944
2560
 
2945
2561
    glist[0] = numLocallyOwnedNodes;
2946
2562
    glist[1] = numLocalEqns;
2947
2563
    glist[2] = numLocalEqnBlks;
2948
 
    if (MPI_Gather(glist.dataPtr(), 3, MPI_INT, recvList.dataPtr(), 3,
 
2564
    if (MPI_Gather(&glist[0], 3, MPI_INT, &recvList[0], 3,
2949
2565
                   MPI_INT, masterProc_, comm_) != MPI_SUCCESS) {
2950
2566
      FEI_CERR << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Gather" << FEI_ENDL;
2951
2567
      MPI_Abort(comm_, -1);
2994
2610
 
2995
2611
  if (numProcs_ > 1) {
2996
2612
#ifndef FEI_SER
2997
 
    feiArray<int> blist(3*(numProcs_+1));
 
2613
    std::vector<int> blist(3*(numProcs_+1));
2998
2614
 
2999
2615
    if (localProc_ == masterProc_) {
3000
2616
      int offset = 0;
3005
2621
      }
3006
2622
    }
3007
2623
 
3008
 
    if (MPI_Bcast(blist.dataPtr(), 3*(numProcs_+1), MPI_INT,
 
2624
    if (MPI_Bcast(&blist[0], 3*(numProcs_+1), MPI_INT,
3009
2625
                  masterProc_, comm_) != MPI_SUCCESS) {
3010
2626
      FEI_CERR << "SNL_FEI_Structure::calcGlobalEqnInfo: ERROR in MPI_Bcast" << FEI_ENDL;
3011
2627
      MPI_Abort(comm_, -1);
3039
2655
    os << "#     numGlobalEqnBlks_: " << numGlobalEqnBlks_ << FEI_ENDL;
3040
2656
    os << "#     localBlkOffset_: " << localBlkOffset_ << FEI_ENDL;
3041
2657
    os << "#     firstLocalNodeNumber_: " << firstLocalNodeNumber_ << FEI_ENDL;
3042
 
    int n;
 
2658
    size_t n;
3043
2659
    os << "#     numGlobalNodes_: " << numGlobalNodes_ << FEI_ENDL;
3044
 
    os << "#     " << globalNodeOffsets_.length() << " globalNodeOffsets_: ";
3045
 
    for(n=0; n<globalNodeOffsets_.length(); n++) 
3046
 
      os << globalNodeOffsets_[n] << " ";
 
2660
    os << "#     " << globalNodeOffsets_.size() << " globalNodeOffsets_: ";
 
2661
    for(size_t nn=0; nn<globalNodeOffsets_.size(); nn++) 
 
2662
      os << globalNodeOffsets_[nn] << " ";
3047
2663
    os << FEI_ENDL;
3048
 
    os << "#     " << globalEqnOffsets_.length() << " globalEqnOffsets_: ";
3049
 
    for(n=0; n<globalEqnOffsets_.length(); n++) 
 
2664
    os << "#     " << globalEqnOffsets_.size() << " globalEqnOffsets_: ";
 
2665
    for(n=0; n<globalEqnOffsets_.size(); n++) 
3050
2666
      os << globalEqnOffsets_[n] << " ";
3051
2667
    os << FEI_ENDL;
3052
 
    os << "#     " << globalBlkEqnOffsets_.length() << " globalBlkEqnOffsets_: ";
3053
 
    for(n=0; n<globalBlkEqnOffsets_.length(); n++) 
 
2668
    os << "#     " << globalBlkEqnOffsets_.size() << " globalBlkEqnOffsets_: ";
 
2669
    for(n=0; n<globalBlkEqnOffsets_.size(); n++) 
3054
2670
      os << globalBlkEqnOffsets_[n] << " ";
3055
2671
    os << FEI_ENDL;
3056
2672
    os << "#  leaving calcGlobalEqnInfo" << FEI_ENDL;
3068
2684
//
3069
2685
//The return value is an error-code.
3070
2686
//
3071
 
  int numNodes = getNumActiveNodes();
3072
2687
  int eqn = localStartRow_;
3073
2688
 
3074
2689
  //localBlkOffset_ is 0-based, and so is blkEqnNumber.
3075
2690
  int blkEqnNumber = localBlkOffset_;
3076
2691
 
3077
 
  for(int i=0; i<numNodes; i++) {
 
2692
  std::map<GlobalID,int>& nodeIDmap = nodeDatabase_->getNodeIDs();
 
2693
  std::map<GlobalID,int>::iterator
 
2694
    iter = nodeIDmap.begin(), iter_end = nodeIDmap.end();
 
2695
 
 
2696
  for(; iter!=iter_end; ++iter) {
3078
2697
    NodeDescriptor* node = NULL;
3079
 
    CHK_ERR( nodeDatabase_->getNodeAtIndex(i, node) );
 
2698
    CHK_ERR( nodeDatabase_->getNodeAtIndex(iter->second, node) );
3080
2699
 
3081
2700
    int numFields = node->getNumFields();
3082
2701
    const int* fieldIDs = node->getFieldIDList();
3147
2766
  int found = snl_fei::binarySearch(blockID, blockIDs_, insertPoint);
3148
2767
 
3149
2768
   if (found<0) {
3150
 
      blockIDs_.insert(blockID, insertPoint);
 
2769
      blockIDs_.insert(blockIDs_.begin()+insertPoint, blockID);
3151
2770
 
3152
2771
      BlockDescriptor* block = new BlockDescriptor;
3153
2772
      block->setGlobalBlockID(blockID);
3154
2773
 
3155
 
      blocks_.insert(block, insertPoint);
 
2774
      blocks_.insert(blocks_.begin()+insertPoint, block);
3156
2775
 
3157
2776
      ConnectivityTable* newConnTable = new ConnectivityTable;
3158
 
      connTables_.insert(newConnTable, insertPoint);
3159
 
   }
3160
 
 
3161
 
   return(FEI_SUCCESS);
3162
 
}
3163
 
 
3164
 
//------------------------------------------------------------------------------
3165
 
int SNL_FEI_Structure::addPattern(int patternID) {
3166
 
//
3167
 
//Append a patternID to our (unsorted) list of patternIDs if it isn't
3168
 
//already present. Also, add a pattern-descriptor if patternID wasn't
3169
 
//already present.
3170
 
//
3171
 
   int found = patternIDs_.find(patternID);
3172
 
 
3173
 
   if (found<0) {
3174
 
      patternIDs_.append(patternID);
3175
 
      int len = patternIDs_.length();
3176
 
 
3177
 
      PatternDescriptor** newPttrns = new PatternDescriptor*[len];
3178
 
      ConnectivityTable** newConn = new ConnectivityTable*[len];
3179
 
 
3180
 
      for(int i=0; i<len-1; i++) {
3181
 
         newPttrns[i] = patterns_[i];
3182
 
         newConn[i] = patternConn_[i];
3183
 
      }
3184
 
      newPttrns[len-1] = new PatternDescriptor;
3185
 
      newConn[len-1] = new ConnectivityTable;
3186
 
 
3187
 
      newPttrns[len-1]->setPatternID(patternID);
3188
 
 
3189
 
      delete [] patterns_;
3190
 
      patterns_ = newPttrns;
3191
 
 
3192
 
      delete [] patternConn_;
3193
 
      patternConn_ = newConn;
 
2777
      connTables_.insert(connTables_.begin()+insertPoint, newConnTable);
3194
2778
   }
3195
2779
 
3196
2780
   return(FEI_SUCCESS);
3224
2808
int SNL_FEI_Structure::getBlockDescriptor_index(int index,
3225
2809
                                               BlockDescriptor*& block)
3226
2810
{
3227
 
  if (index < 0 || index >= blockIDs_.length()) ERReturn(FEI_FATAL_ERROR);
 
2811
  if (index < 0 || index >= (int)blockIDs_.size()) ERReturn(FEI_FATAL_ERROR);
3228
2812
 
3229
2813
  block = blocks_[index];
3230
2814
 
3232
2816
}
3233
2817
 
3234
2818
//------------------------------------------------------------------------------
3235
 
int SNL_FEI_Structure::getPatternDescriptor(int patternID,
3236
 
                                         PatternDescriptor*& pattern)
3237
 
{
3238
 
   int index = patternIDs_.find(patternID);
3239
 
 
3240
 
   if (index < 0) {
3241
 
      FEI_CERR << "SNL_FEI_Structure::getPatternDescriptor: ERROR, patternID "
3242
 
           << (int)patternID << " not found." << FEI_ENDL;
3243
 
      return(-1);
3244
 
   }
3245
 
 
3246
 
   pattern = &(*(patterns_[index]));
3247
 
 
3248
 
   return(FEI_SUCCESS);
3249
 
}
3250
 
 
3251
 
//------------------------------------------------------------------------------
3252
 
int SNL_FEI_Structure::getPatternDescriptor_index(int index,
3253
 
                                         PatternDescriptor*& pattern)
3254
 
{
3255
 
  if (index < 0 || index >= patternIDs_.length()) return(FEI_FATAL_ERROR);
3256
 
 
3257
 
  pattern = patterns_[index];
3258
 
 
3259
 
  return(FEI_SUCCESS);
3260
 
}
3261
 
 
3262
 
//------------------------------------------------------------------------------
3263
2819
int SNL_FEI_Structure::allocateBlockConnectivity(GlobalID blockID) {
3264
2820
 
3265
2821
   int index = snl_fei::binarySearch(blockID, blockIDs_);
3271
2827
   }
3272
2828
 
3273
2829
   connTables_[index]->numNodesPerElem = 
3274
 
                         blocks_[index]->numNodesPerElement;
 
2830
                         blocks_[index]->getNumNodesPerElement();
3275
2831
 
3276
2832
   int numRows = blocks_[index]->getNumElements();
3277
2833
   int numCols = connTables_[index]->numNodesPerElem;
3283
2839
      MPI_Abort(comm_, -1);
3284
2840
   }
3285
2841
 
3286
 
   int err = connTables_[index]->elemNumbers.reAllocate(numRows);
3287
 
   if (err) return(-1);
3288
 
 
3289
 
   connTables_[index]->elem_conn_ids = new feiArray<GlobalID>(numRows*numCols);
3290
 
 
3291
 
   return(FEI_SUCCESS);
3292
 
}
3293
 
 
3294
 
//------------------------------------------------------------------------------
3295
 
int SNL_FEI_Structure::appendPatternConnectivity(int patternID,
3296
 
                                                int numRowIDs,
3297
 
                                                const int* rowIDTypes,
3298
 
                                                const GlobalID* rowIDs,
3299
 
                                                int numColIDsPerRow,
3300
 
                                                const int* colIDTypes,
3301
 
                                                const GlobalID* colIDs)
3302
 
{
3303
 
   int index = patternIDs_.find(patternID);
3304
 
 
3305
 
   if (index < 0) {
3306
 
      FEI_CERR << "SNL_FEI_Structure::appendPatternConnectivity: ERROR, patternID "
3307
 
           << (int)patternID << " not found. Aborting." << FEI_ENDL;
3308
 
      MPI_Abort(comm_, -1);
3309
 
   }
3310
 
 
3311
 
   if (numColIDsPerRow <= 0) return(FEI_SUCCESS);
3312
 
 
3313
 
   int& numRows = patternConn_[index]->numRows;
3314
 
   feiArray<GlobalID>**& conn = patternConn_[index]->connectivities;
3315
 
 
3316
 
   numRows += 2;
3317
 
 
3318
 
   int len = numRows;
3319
 
 
3320
 
   feiArray<GlobalID>** newConn = new feiArray<GlobalID>*[len];
3321
 
 
3322
 
   for(int i=0; i<len-2; i++) newConn[i] = conn[i];
3323
 
 
3324
 
   newConn[len-2] = new feiArray<GlobalID>(numRowIDs, numRowIDs);
3325
 
   newConn[len-1] = new feiArray<GlobalID>(numRowIDs*numColIDsPerRow,numRowIDs);
3326
 
 
3327
 
   feiArray<GlobalID>& rowConn = *(newConn[len-2]);
3328
 
   for(int j=0; j<numRowIDs; j++) rowConn[j] = rowIDs[j];
3329
 
 
3330
 
   feiArray<GlobalID>& colConn = *(newConn[len-1]);
3331
 
   for(int k=0; k<numRowIDs*numColIDsPerRow; k++) colConn[k] = colIDs[k];
3332
 
 
3333
 
   delete [] conn;
3334
 
   conn = newConn;
 
2842
   connTables_[index]->elemNumbers.resize(numRows);
 
2843
 
 
2844
   connTables_[index]->elem_conn_ids = new std::vector<GlobalID>(numRows*numCols);
3335
2845
 
3336
2846
   return(FEI_SUCCESS);
3337
2847
}
3351
2861
}
3352
2862
 
3353
2863
//------------------------------------------------------------------------------
3354
 
ConnectivityTable& SNL_FEI_Structure::getPatternConnectivity(int patternID) {
3355
 
 
3356
 
   int index = patternIDs_.find(patternID);
3357
 
 
3358
 
   if (index < 0) {
3359
 
      FEI_CERR << "SNL_FEI_Structure::getPatternConnectivity: ERROR, patternID "
3360
 
           << (int)patternID << " not found. Aborting." << FEI_ENDL;
3361
 
      MPI_Abort(comm_, -1);
3362
 
   }
3363
 
 
3364
 
   return(*(patternConn_[index]));
3365
 
}
3366
 
 
3367
 
//------------------------------------------------------------------------------
3368
2864
int SNL_FEI_Structure::finalizeActiveNodes()
3369
2865
{
3370
2866
  FEI_OSTREAM& os = dbgOut();
3379
2875
    CHK_ERR( nodeDatabase_->initNodeID(shNodeIDs[i]) );
3380
2876
  }
3381
2877
 
3382
 
  //That should be all of the nodeIDs, so allocate the list of NodeDescriptors.
3383
 
  //
3384
 
  CHK_ERR( nodeDatabase_->allocateNodeDescriptors() );
3385
 
 
3386
2878
  if (activeNodesInitialized_) return(0);
3387
2879
  //
3388
2880
  //Now run through the connectivity tables and set the nodal field and block
3395
2887
  //globally unique).
3396
2888
  //
3397
2889
  int elemNumber = 0;
3398
 
  int numBlocks = blockIDs_.length();
 
2890
  int numBlocks = blockIDs_.size();
3399
2891
  for(int bIndex=0; bIndex<numBlocks; bIndex++) {
3400
2892
 
3401
2893
    GlobalID blockID = blocks_[bIndex]->getGlobalBlockID();
3406
2898
 
3407
2899
    int numElems = conn.elemIDs.size();
3408
2900
    if (numElems > 0) {
3409
 
      elemConn = conn.elem_conn_ids->dataPtr();
 
2901
      elemConn = &((*conn.elem_conn_ids)[0]);
3410
2902
      if (!activeNodesInitialized_) {
3411
 
        int elemConnLen = conn.elem_conn_ids->length();
3412
 
        conn.elem_conn_ptrs = new feiArray<NodeDescriptor*>(elemConnLen);
 
2903
        int elemConnLen = conn.elem_conn_ids->size();
 
2904
        conn.elem_conn_ptrs = new std::vector<NodeDescriptor*>(elemConnLen);
3413
2905
      }
3414
 
      elemNodeDescPtrs = conn.elem_conn_ptrs->dataPtr();
 
2906
      elemNodeDescPtrs = &((*conn.elem_conn_ptrs)[0]);
3415
2907
    }
3416
2908
    int nodesPerElem = conn.numNodesPerElem;
3417
2909
 
3481
2973
 
3482
2974
    while(cr_iter != cr_end) {
3483
2975
      ConstraintType& cr = *((*cr_iter).second);
3484
 
      GlobalID* nodeIDs = cr.getMasters()->dataPtr();
3485
 
      int numNodes = cr.getMasters()->length();
 
2976
      std::vector<GlobalID>& nodeID_vec = cr.getMasters();
 
2977
      GlobalID* nodeIDs = &nodeID_vec[0];
 
2978
      int numNodes = cr.getMasters().size();
3486
2979
 
3487
2980
      NodeDescriptor* node = NULL;
3488
2981
      for(int k=0; k<numNodes; ++k) {
3493
2986
    }
3494
2987
  }
3495
2988
 
3496
 
  //
3497
 
  //we also need to run through the patterns and their connectivities, and
3498
 
  //add the fields for each pattern to the nodes that are accessed using that
3499
 
  //pattern.
3500
 
  //
3501
 
 
3502
 
  int numPatterns = patternIDs_.length();
3503
 
  for(int i=0; i<numPatterns; i++) {
3504
 
    PatternDescriptor& pattern = *(patterns_[i]);
3505
 
 
3506
 
    feiArray<int>& fieldsPerRow = pattern.getNumFieldsPerRow();
3507
 
    feiArray<int>* rowFields = pattern.getRowFieldIDs();
3508
 
    feiArray<int>& fieldsPerCol = pattern.getNumFieldsPerCol();
3509
 
    feiArray<int>* colFields = pattern.getColFieldIDs();
3510
 
    int numColIDsPerRow = pattern.getNumColIDsPerRow();
3511
 
 
3512
 
    int n,numConnRows = patternConn_[i]->numRows;
3513
 
    int offset = 0, loopCount = numConnRows/2;
3514
 
    for(int r=0; r<loopCount; r++) {
3515
 
      feiArray<GlobalID>& rowIDs = *(patternConn_[i]->connectivities[offset++]);
3516
 
      feiArray<GlobalID>& colIDs = *(patternConn_[i]->connectivities[offset++]);
3517
 
 
3518
 
      for(n=0; n<rowIDs.length(); n++) {
3519
 
        NodeDescriptor* node = NULL;
3520
 
        CHK_ERR( nodeDatabase_->getNodeWithID(rowIDs[n], node) );
3521
 
 
3522
 
        for(int f=0; f<fieldsPerRow[n]; f++) {
3523
 
          node->addField((rowFields[n])[f]);
3524
 
        }
3525
 
 
3526
 
        node->setOwnerProc(localProc_);
3527
 
        nodeCommMgr_->informLocal(*node);
3528
 
 
3529
 
        for(int c=0; c<fieldsPerCol.length(); c++) {
3530
 
          CHK_ERR(nodeDatabase_->getNodeWithID(colIDs[n*numColIDsPerRow + c],
3531
 
                                               node) );
3532
 
 
3533
 
          for(int f=0; f<fieldsPerCol[c]; f++) {
3534
 
            node->addField((colFields[c])[f]);
3535
 
          }
3536
 
 
3537
 
          node->setOwnerProc(localProc_);
3538
 
          //      nodeCommMgr_->informLocal(*node);       
3539
 
        }
3540
 
      }
3541
 
    }
3542
 
  }
3543
 
 
3544
2989
  activeNodesInitialized_ = true;
3545
2990
 
3546
2991
  if (debugOutput_) os << "#  leaving finalizeActiveNodes" << FEI_ENDL;
3587
3032
  //This function will count how 
3588
3033
  //many active nodes there are for each block.
3589
3034
 
3590
 
   int numBlocks = blockIDs_.length();
3591
 
   feiArray<int> nodesPerBlock(numBlocks);
3592
 
   feiArray<int> eqnsPerBlock(numBlocks);
3593
 
   int* nodesPerBlockPtr = nodesPerBlock.dataPtr();
3594
 
   int* eqnsPerBlockPtr = eqnsPerBlock.dataPtr();
3595
 
   GlobalID* blockIDsPtr = blockIDs_.dataPtr();
 
3035
   int numBlocks = blockIDs_.size();
 
3036
   std::vector<int> nodesPerBlock(numBlocks);
 
3037
   std::vector<int> eqnsPerBlock(numBlocks);
 
3038
   GlobalID* blockIDsPtr = &blockIDs_[0];
3596
3039
 
3597
3040
   int j;
3598
3041
   for(j=0; j<numBlocks; j++) {
3599
 
      nodesPerBlockPtr[j] = 0;
3600
 
      eqnsPerBlockPtr[j] = 0;
 
3042
      nodesPerBlock[j] = 0;
 
3043
      eqnsPerBlock[j] = 0;
3601
3044
   }
3602
3045
 
3603
3046
   int numNodes = nodeDatabase_->getNumNodeDescriptors();
3614
3057
       GlobalID blockID = blockIDsPtr[blk];
3615
3058
 
3616
3059
       if (node->containedInBlock(blockID)) {
3617
 
         nodesPerBlockPtr[blk]++;
 
3060
         nodesPerBlock[blk]++;
3618
3061
       }
3619
3062
       else {
3620
3063
         continue;
3624
3067
         if (blocks_[blk]->containsField(fieldIDList[fld])) {
3625
3068
           int fSize = getFieldSize(fieldIDList[fld]);
3626
3069
           assert(fSize >= 0);
3627
 
           eqnsPerBlockPtr[blk] += fSize;
 
3070
           eqnsPerBlock[blk] += fSize;
3628
3071
         }
3629
3072
       }
3630
3073
     }
3631
3074
   }
3632
3075
 
3633
3076
   for(j=0; j<numBlocks; j++) {
3634
 
     blocks_[j]->setNumActiveNodes(nodesPerBlockPtr[j]);
 
3077
     blocks_[j]->setNumActiveNodes(nodesPerBlock[j]);
3635
3078
   }
3636
3079
 
3637
3080
   //now we need to add the elem-dof to the eqn-count for each block,
3638
3081
   //and then set the total number of eqns on each block.
3639
3082
 
3640
3083
   for(j=0; j<numBlocks; j++) {
3641
 
     eqnsPerBlockPtr[j] += blocks_[j]->getNumElemDOFPerElement() *
 
3084
     eqnsPerBlock[j] += blocks_[j]->getNumElemDOFPerElement() *
3642
3085
                        blocks_[j]->getNumElements();
3643
3086
 
3644
 
     blocks_[j]->setTotalNumEqns(eqnsPerBlockPtr[j]);
 
3087
     blocks_[j]->setTotalNumEqns(eqnsPerBlock[j]);
3645
3088
   }
3646
3089
   return(0);
3647
3090
}
3766
3209
 
3767
3210
        if ( nodeDatabase_->getNodeWithID(IDs[i], node) != 0 ) {
3768
3211
            // FEI_CERR << "SNL_FEI_Structure::getEqnNumbers: ERROR getting node " << IDs[i] << FEI_ENDL;
3769
 
            for(int i=0; i<fieldSize; i++) {
 
3212
            for(int ii=0; ii<fieldSize; ii++) {
3770
3213
                eqnNumbers[offset++] = -1;
3771
3214
            }
3772
3215
            continue;
3779
3222
            ERReturn(-1);
3780
3223
        }
3781
3224
 
3782
 
        for(int i=0; i<fieldSize; i++) {
3783
 
            eqnNumbers[offset++] = nodeFieldEqnNumber + i;
 
3225
        for(int ii=0; ii<fieldSize; ii++) {
 
3226
            eqnNumbers[offset++] = nodeFieldEqnNumber + ii;
3784
3227
        }
3785
3228
    }
3786
3229
 
3822
3265
  if (debugOutput_) os << "#  calculateSlaveEqns" << FEI_ENDL;
3823
3266
 
3824
3267
  if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
3825
 
  eqnCommMgr_ = new EqnCommMgr(*commUtilsInt_);
 
3268
  eqnCommMgr_ = new EqnCommMgr(comm_);
3826
3269
  eqnCommMgr_->setNumRHSs(1);
3827
3270
 
3828
 
  int i;
3829
 
 
3830
 
  feiArray<int> eqns;
3831
 
  feiArray<int> mEqns;
3832
 
  feiArray<double> mCoefs;
3833
 
 
3834
 
  for(i=0; i<slaveVars_->length(); i++) {
 
3271
  std::vector<int> eqns;
 
3272
  std::vector<int> mEqns;
 
3273
  std::vector<double> mCoefs;
 
3274
 
 
3275
  for(size_t i=0; i<slaveVars_->size(); i++) {
3835
3276
    int numEqns;
3836
3277
    SlaveVariable* svar = (*slaveVars_)[i];
3837
3278
 
3838
3279
    eqns.resize( getFieldSize(svar->getFieldID()));
3839
3280
    CHK_ERR( getEqnNumbers(svar->getNodeID(), FEI_NODE, svar->getFieldID(),
3840
 
                           numEqns, eqns.dataPtr()));
 
3281
                           numEqns, &eqns[0]));
3841
3282
 
3842
3283
    int slaveEqn = eqns[svar->getFieldOffset()];
3843
3284
 
3844
 
    const feiArray<GlobalID>* mNodes = svar->getMasterNodeIDs();
3845
 
    const feiArray<int>* mFields = svar->getMasterFields();
3846
 
    const feiArray<double>* mWeights = svar->getWeights();
3847
 
    const feiArray<double>& mWeightsRef = *mWeights;
 
3285
    const std::vector<GlobalID>* mNodes = svar->getMasterNodeIDs();
 
3286
    const std::vector<int>* mFields = svar->getMasterFields();
 
3287
    const std::vector<double>* mWeights = svar->getWeights();
 
3288
    const std::vector<double>& mWeightsRef = *mWeights;
3848
3289
    int mwOffset = 0;
3849
3290
 
3850
 
    for(int j=0; j<mNodes->length(); j++) {
 
3291
    for(size_t j=0; j<mNodes->size(); j++) {
3851
3292
      int mfSize = getFieldSize((*mFields)[j]);
3852
3293
 
3853
3294
      eqns.resize(mfSize);
3854
3295
      GlobalID mNodeID = (*mNodes)[j];
3855
3296
      int mFieldID = (*mFields)[j];
3856
3297
      CHK_ERR( getEqnNumbers( mNodeID, FEI_NODE, mFieldID,
3857
 
                              mfSize, eqns.dataPtr()));
 
3298
                              mfSize, &eqns[0]));
3858
3299
 
3859
3300
      double fei_eps = 1.e-49;
3860
3301
 
3861
3302
      for(int k=0; k<mfSize; k++) {
3862
3303
        if (std::abs(mWeightsRef[mwOffset++]) > fei_eps) {
3863
 
          mEqns.append(eqns[k]);
3864
 
          mCoefs.append(mWeightsRef[mwOffset-1]);
 
3304
          mEqns.push_back(eqns[k]);
 
3305
          mCoefs.push_back(mWeightsRef[mwOffset-1]);
3865
3306
        }
3866
3307
      }
3867
3308
    }
3868
3309
 
3869
 
    CHK_ERR( slaveEqns_->addEqn(slaveEqn, mCoefs.dataPtr(), mEqns.dataPtr(),
3870
 
                       mEqns.length(), false) );
 
3310
    CHK_ERR( slaveEqns_->addEqn(slaveEqn, &mCoefs[0], &mEqns[0],
 
3311
                       mEqns.size(), false) );
3871
3312
    mEqns.resize(0);
3872
3313
    mCoefs.resize(0);
3873
3314
  }
3874
3315
 
3875
3316
#ifndef FEI_SER
3876
 
  int numLocalSlaves = slaveVars_->length();
 
3317
  int numLocalSlaves = slaveVars_->size();
3877
3318
  int globalMaxSlaves = 0;
3878
 
  CHK_ERR( commUtilsInt_->GlobalMax(numLocalSlaves, globalMaxSlaves) );
 
3319
  CHK_ERR( fei::GlobalMax(comm_, numLocalSlaves, globalMaxSlaves) );
3879
3320
 
3880
3321
  if (globalMaxSlaves > 0) {
3881
3322
    CHK_ERR( gatherSlaveEqns(comm, eqnCommMgr_, slaveEqns_) );
3884
3325
 
3885
3326
  globalNumNodesVanished_.resize(numProcs_+1, 0);
3886
3327
 
3887
 
  slvEqnNumbers_ = &(slaveEqns_->eqnNumbersPtr());
3888
 
  numSlvs_ = slvEqnNumbers_->length();
 
3328
  slvEqnNumbers_ = &(slaveEqns_->eqnNumbers());
 
3329
  numSlvs_ = slvEqnNumbers_->size();
3889
3330
  if (numSlvs_ > 0) {
3890
3331
    //first, let's remove any 'couplings' among the slave equations. A coupling
3891
3332
    //is where a slave depends on a master which is also a slave that depends on
3901
3342
    CHK_ERR( removeCouplings(*slaveEqns_, levelsOfCoupling) );
3902
3343
 
3903
3344
    if (debugOutput_) {
3904
 
      FEI_OSTREAM& os = dbgOut();
3905
3345
      os << "#     SNL_FEI_Structure: " << levelsOfCoupling
3906
3346
         << " level(s) of coupling removed: " << FEI_ENDL;
3907
3347
    }
3922
3362
    // to the eqnCommMgr_ object using addLocalEqn for each processor that
3923
3363
    // shares the slave node.
3924
3364
 
3925
 
    feiArray<int>& slvEqns = *slvEqnNumbers_;
3926
 
    feiArray<SSVec*>& mstrEqns = slaveEqns_->eqns();
 
3365
    std::vector<int>& slvEqns = *slvEqnNumbers_;
 
3366
    std::vector<fei::CSVec*>& mstrEqns = slaveEqns_->eqns();
3927
3367
 
3928
3368
    //keep track of the number of locally owned nodes that vanish due to the
3929
3369
    //fact that all equations at that node are slave equations.
3931
3371
 
3932
3372
    GlobalID lastNodeID = -1;
3933
3373
 
3934
 
    for(i=0; i<numSlvs_; i++) {
3935
 
      NodeDescriptor* node = NULL;
 
3374
    for(int i=0; i<numSlvs_; i++) {
 
3375
      const NodeDescriptor* node = NULL;
3936
3376
      int reducedSlaveEqn;
3937
3377
      translateToReducedEqn(slvEqns[i], reducedSlaveEqn);
3938
 
      int numMasters = mstrEqns[i]->length();
 
3378
      int numMasters = mstrEqns[i]->size();
3939
3379
 
3940
3380
      int err = nodeDatabase_->getNodeWithEqn(slvEqns[i], node);
3941
3381
      if (err != 0) {
4003
3443
 
4004
3444
    std::vector<int> tmp(1), tmp2(numProcs_);
4005
3445
    tmp[0] = numLocalNodesVanished;
4006
 
    CHK_ERR( commUtilsInt_->Allgatherv(tmp, tmp2, globalNumNodesVanished_) );
 
3446
    CHK_ERR( fei::Allgatherv(comm_, tmp, tmp2, globalNumNodesVanished_) );
4007
3447
 
4008
3448
    if ((int)globalNumNodesVanished_.size() != numProcs_) {
4009
3449
      ERReturn(-1);
4020
3460
  }
4021
3461
 
4022
3462
  if (slaveMatrix_ != NULL) delete slaveMatrix_;
4023
 
  slaveMatrix_ = new SSMat(*slaveEqns_);
 
3463
  slaveMatrix_ = new fei::FillableMat(*slaveEqns_);
4024
3464
 
4025
3465
  if (debugOutput_) {
4026
3466
    os << "#  slave-equations:" << FEI_ENDL;
4034
3474
//------------------------------------------------------------------------------
4035
3475
int SNL_FEI_Structure::removeCouplings(EqnBuffer& eqnbuf, int& levelsOfCoupling)
4036
3476
{
4037
 
  feiArray<int>& eqnNumbers = eqnbuf.eqnNumbersPtr();
4038
 
  feiArray<SSVec*>& eqns = eqnbuf.eqns();
 
3477
  std::vector<int>& eqnNumbers = eqnbuf.eqnNumbers();
 
3478
  std::vector<fei::CSVec*>& eqns = eqnbuf.eqns();
4039
3479
 
4040
 
  feiArray<double> tempCoefs;
4041
 
  feiArray<int> tempIndices;
 
3480
  std::vector<double> tempCoefs;
 
3481
  std::vector<int> tempIndices;
4042
3482
 
4043
3483
  levelsOfCoupling = 0;
4044
3484
  bool finished = false;
4045
3485
  while(!finished) {
4046
3486
    bool foundCoupling = false;
4047
 
    for(int i=0; i<eqnNumbers.length(); i++) {
 
3487
    for(size_t i=0; i<eqnNumbers.size(); i++) {
4048
3488
      int rowIndex = eqnbuf.isInIndices(eqnNumbers[i]);
4049
3489
 
4050
 
      if (rowIndex == i) {
 
3490
      if (rowIndex == (int)i) {
4051
3491
        FEI_CERR <<" SNL_FEI_Structure::removeCouplings ERROR,"
4052
3492
             << " illegal master-slave constraint coupling. Eqn "
4053
3493
             << eqnNumbers[i] << " is both master and slave. " << FEI_ENDL;
4060
3500
        CHK_ERR( eqnbuf.getCoefAndRemoveIndex( eqnNumbers[rowIndex],
4061
3501
                                               eqnNumbers[i], coef) );
4062
3502
 
4063
 
        feiArray<int>& indicesRef = eqns[i]->indices();
4064
 
        feiArray<double>& coefsRef = eqns[i]->coefs();
 
3503
        std::vector<int>& indicesRef = eqns[i]->indices();
 
3504
        std::vector<double>& coefsRef = eqns[i]->coefs();
4065
3505
 
4066
 
        int len = indicesRef.length();
 
3506
        int len = indicesRef.size();
4067
3507
        tempCoefs.resize(len);
4068
3508
        tempIndices.resize(len);
4069
3509
 
4070
 
        double* tempCoefsPtr = tempCoefs.dataPtr();
4071
 
        int* tempIndicesPtr = tempIndices.dataPtr();
4072
 
        double* coefsPtr = coefsRef.dataPtr();
4073
 
        int* indicesPtr = indicesRef.dataPtr();
 
3510
        double* tempCoefsPtr = &tempCoefs[0];
 
3511
        int* tempIndicesPtr = &tempIndices[0];
 
3512
        double* coefsPtr = &coefsRef[0];
 
3513
        int* indicesPtr = &indicesRef[0];
4074
3514
 
4075
3515
        for(int j=0; j<len; ++j) {
4076
3516
          tempIndicesPtr[j] = indicesPtr[j];
4117
3557
  //the most efficient way to do it, but it involves the least amount of new
4118
3558
  //code.)
4119
3559
  ProcEqns localProcEqns, remoteProcEqns;
4120
 
  feiArray<int>& slvEqnNums = slaveEqns->eqnNumbersPtr();
4121
 
  SSVec** slvEqnsPtr = slaveEqns->eqns().dataPtr();
 
3560
  std::vector<int>& slvEqnNums = slaveEqns->eqnNumbers();
 
3561
  fei::CSVec** slvEqnsPtr = &(slaveEqns->eqns()[0]);
4122
3562
 
4123
 
  for(int i=0; i<slvEqnNums.length(); i++) {
 
3563
  for(size_t i=0; i<slvEqnNums.size(); i++) {
4124
3564
    for(int p=0; p<numProcs; p++) {
4125
3565
      if (p == localProc) continue;
4126
3566
 
4127
 
      localProcEqns.addEqn(slvEqnNums[i], slvEqnsPtr[i]->length(), p);
 
3567
      localProcEqns.addEqn(slvEqnNums[i], slvEqnsPtr[i]->size(), p);
4128
3568
    }
4129
3569
  }
4130
3570
 
4151
3591
{
4152
3592
  if (numSlvs_ == 0) return(false);
4153
3593
 
4154
 
  feiArray<int>& slvEqns = slaveEqns_->eqnNumbersPtr();
 
3594
  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
4155
3595
  int insertPoint = -1;
4156
 
  int foundOffset = snl_fei::binarySearch(eqn, slvEqns.dataPtr(),
4157
 
                                                 slvEqns.length(), insertPoint);
 
3596
  int foundOffset = snl_fei::binarySearch(eqn, slvEqns, insertPoint);
4158
3597
 
4159
3598
  if (foundOffset >= 0) return(true);
4160
3599
  else return(false);
4169
3608
  if (eqn > highestSlv_) {reducedEqn = eqn - numSlvs_; return(false); }
4170
3609
 
4171
3610
  int index = 0;
4172
 
  int foundOffset = snl_fei::binarySearch(eqn, slvEqnNumbers_->dataPtr(),
4173
 
                                        slvEqnNumbers_->length(), index);
 
3611
  int foundOffset = snl_fei::binarySearch(eqn, *slvEqnNumbers_, index);
4174
3612
 
4175
3613
  bool isSlave = false;
4176
3614
 
4191
3629
 
4192
3630
  if (numSlvs == 0) return(reducedEqn);
4193
3631
 
4194
 
  const int* slvEqns = slaveEqns_->eqnNumbersPtr().dataPtr();
 
3632
  const int* slvEqns = &(slaveEqns_->eqnNumbers()[0]);
4195
3633
 
4196
3634
  if (reducedEqn < slvEqns[0]) return(reducedEqn);
4197
3635
 
4207
3645
 
4208
3646
//------------------------------------------------------------------------------
4209
3647
int SNL_FEI_Structure::getMasterEqnNumbers(int slaveEqn,
4210
 
                                          feiArray<int>*& masterEqns)
 
3648
                                          std::vector<int>*& masterEqns)
4211
3649
{
4212
3650
  if (slaveEqns_->getNumEqns() == 0) {
4213
3651
    masterEqns = NULL;
4214
3652
    return(0);
4215
3653
  }
4216
3654
 
4217
 
  feiArray<int>& slvEqns = slaveEqns_->eqnNumbersPtr();
 
3655
  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
4218
3656
  int index = 0;
4219
 
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns.dataPtr(),
4220
 
                                        slvEqns.length(), index);
 
3657
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns, index);
4221
3658
 
4222
3659
  if (foundOffset >= 0) {
4223
3660
    masterEqns = &(slaveEqns_->eqns()[foundOffset]->indices());
4231
3668
 
4232
3669
//------------------------------------------------------------------------------
4233
3670
int SNL_FEI_Structure::getMasterEqnCoefs(int slaveEqn,
4234
 
                                        feiArray<double>*& masterCoefs)
 
3671
                                        std::vector<double>*& masterCoefs)
4235
3672
{
4236
3673
  if (slaveEqns_->getNumEqns() == 0) {
4237
3674
    masterCoefs = NULL;
4238
3675
    return(0);
4239
3676
  }
4240
3677
 
4241
 
  feiArray<int>& slvEqns = slaveEqns_->eqnNumbersPtr();
 
3678
  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
4242
3679
  int index = 0;
4243
 
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns.dataPtr(),
4244
 
                                        slvEqns.length(), index);
 
3680
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns, index);
4245
3681
 
4246
3682
  if (foundOffset >= 0) {
4247
3683
    masterCoefs = &(slaveEqns_->eqns()[foundOffset]->coefs());
4261
3697
    return(0);
4262
3698
  }
4263
3699
 
4264
 
  feiArray<int>& slvEqns = slaveEqns_->eqnNumbersPtr();
 
3700
  std::vector<int>& slvEqns = slaveEqns_->eqnNumbers();
4265
3701
  int index = 0;
4266
 
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns.dataPtr(),
4267
 
                                        slvEqns.length(), index);
 
3702
  int foundOffset = snl_fei::binarySearch(slaveEqn, slvEqns, index);
4268
3703
 
4269
3704
  if (foundOffset >= 0) {
4270
 
    feiArray<double>* rhsCoefsPtr = (*(slaveEqns_->rhsCoefsPtr()))[foundOffset];
 
3705
    std::vector<double>* rhsCoefsPtr = (*(slaveEqns_->rhsCoefsPtr()))[foundOffset];
4271
3706
    rhsValue = (*rhsCoefsPtr)[0];
4272
3707
  }
4273
3708
 
4345
3780
                                                  int* scatterIndices)
4346
3781
{
4347
3782
  BlockDescriptor& block = *(blocks_[blockIndex]);
4348
 
  int numNodes = block.numNodesPerElement;
 
3783
  int numNodes = block.getNumNodesPerElement();
4349
3784
  work_nodePtrs_.resize(numNodes);
4350
 
  NodeDescriptor** nodes = work_nodePtrs_.dataPtr();
 
3785
  NodeDescriptor** nodes = &work_nodePtrs_[0];
4351
3786
  int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
4352
3787
  if (err) {
4353
3788
    FEI_CERR << "SNL_FEI_Structure::getBlkScatterIndices_index: ERROR getting"
4368
3803
// and be of length the number of equations per element.
4369
3804
//
4370
3805
   BlockDescriptor& block = *(blocks_[blockIndex]);
4371
 
   int numNodes = block.numNodesPerElement;
 
3806
   int numNodes = block.getNumNodesPerElement();
4372
3807
   int* fieldsPerNode = block.fieldsPerNodePtr();
4373
3808
   int** fieldIDs = block.fieldIDsTablePtr();
4374
3809
 
4375
3810
   work_nodePtrs_.resize(numNodes);
4376
 
   NodeDescriptor** nodes = work_nodePtrs_.dataPtr();
 
3811
   NodeDescriptor** nodes = &work_nodePtrs_[0];
4377
3812
 
4378
3813
   int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
4379
3814
   if (err) {
4428
3863
// and be of length the number of equations per element.
4429
3864
//
4430
3865
   BlockDescriptor& block = *(blocks_[blockIndex]);
4431
 
   int numNodes = block.numNodesPerElement;
 
3866
   int numNodes = block.getNumNodesPerElement();
4432
3867
   int* fieldsPerNode = block.fieldsPerNodePtr();
4433
3868
   int** fieldIDs = block.fieldIDsTablePtr();
4434
3869
 
4435
3870
   work_nodePtrs_.resize(numNodes);
4436
 
   NodeDescriptor** nodes = work_nodePtrs_.dataPtr();
 
3871
   NodeDescriptor** nodes = &work_nodePtrs_[0];
4437
3872
 
4438
3873
   int err = getElemNodeDescriptors(blockIndex, elemIndex, nodes);
4439
3874
   if (err) {
4484
3919
}
4485
3920
 
4486
3921
//------------------------------------------------------------------------------
4487
 
int SNL_FEI_Structure::getPatternScatterIndices(int patternID,
4488
 
                                               const GlobalID* rowNodes,
4489
 
                                               const GlobalID* colNodes,
4490
 
                                               feiArray<int>& rowIndices,
4491
 
                                               feiArray<int>& rowColOffsets,
4492
 
                                               int& colIndicesPerRow,
4493
 
                                               feiArray<int>& colIndices)
4494
 
{
4495
 
   int index = patternIDs_.find(patternID);
4496
 
 
4497
 
   if (index < 0) {
4498
 
      FEI_CERR << "SNL_FEI_Structure::getPatternScatterIndices: ERROR, patternID "
4499
 
           << (int)patternID << " not found." << FEI_ENDL;
4500
 
      return(-1);
4501
 
   }
4502
 
 
4503
 
   PatternDescriptor& pattern = *(patterns_[index]);
4504
 
 
4505
 
   NodeDescriptor** rNodeDesc = NULL;
4506
 
   NodeDescriptor** cNodeDesc = NULL;
4507
 
 
4508
 
   CHK_ERR( getPatternNodeDescriptors(pattern, rowNodes, colNodes,
4509
 
                                      rNodeDesc, cNodeDesc) );
4510
 
 
4511
 
   int numRows = pattern.getNumRowIDs();
4512
 
   int numColsPerRow = pattern.getNumColIDsPerRow();
4513
 
   int interleave = pattern.getInterleaveStrategy();
4514
 
 
4515
 
   feiArray<int>& fieldsPerRow = pattern.getNumFieldsPerRow();
4516
 
   feiArray<int>* rowFieldIDs = pattern.getRowFieldIDs();
4517
 
 
4518
 
   feiArray<int>& fieldsPerCol = pattern.getNumFieldsPerCol();
4519
 
   feiArray<int>* colFieldIDs = pattern.getColFieldIDs();
4520
 
   int** cFieldIDs = new int*[fieldsPerCol.length()];
4521
 
   for(int c=0; c<fieldsPerCol.length(); c++)
4522
 
     cFieldIDs[c] = colFieldIDs[c].dataPtr();
4523
 
 
4524
 
   int i;
4525
 
   colIndicesPerRow = 0;
4526
 
 
4527
 
   for(i=0; i<fieldsPerCol.length(); i++) {
4528
 
     for(int j=0; j<fieldsPerCol[i]; j++) 
4529
 
       colIndicesPerRow += getFieldSize(colFieldIDs[i][j]);
4530
 
   }
4531
 
 
4532
 
   colIndices.resize(numRows*colIndicesPerRow);
4533
 
 
4534
 
   CHK_ERR( setPatternRowColOffsets( pattern, colIndicesPerRow, rowColOffsets) );
4535
 
 
4536
 
   for(i=0; i<numRows; i++) {
4537
 
     int dummy2 = 0;
4538
 
 
4539
 
     if (interleave == FEI_NODE_MAJOR) {
4540
 
       CHK_ERR( getNodeMajorIndices(&(cNodeDesc[i*numColsPerRow]),
4541
 
                                    numColsPerRow,
4542
 
                                    cFieldIDs, fieldsPerCol.dataPtr(),
4543
 
                                    &(colIndices[i*colIndicesPerRow]), dummy2) )
4544
 
     }
4545
 
     else {
4546
 
       CHK_ERR(getFieldMajorIndices(&(cNodeDesc[i*numColsPerRow]),
4547
 
                                    numColsPerRow,
4548
 
                                    cFieldIDs, fieldsPerCol.dataPtr(),
4549
 
                                    &(colIndices[i*colIndicesPerRow]), dummy2) )
4550
 
     }
4551
 
   }
4552
 
 
4553
 
   delete [] cFieldIDs;
4554
 
 
4555
 
   if (interleave == FEI_NODE_MAJOR) {
4556
 
     CHK_ERR( getNodeMajorIndices(rNodeDesc, numRows, 
4557
 
                                  rowFieldIDs, fieldsPerRow, rowIndices ) );
4558
 
   }
4559
 
   else {
4560
 
     CHK_ERR( getFieldMajorIndices(rNodeDesc, numRows, 
4561
 
                                  rowFieldIDs, fieldsPerRow, rowIndices) );
4562
 
   }
4563
 
 
4564
 
   return(FEI_SUCCESS);
4565
 
}
4566
 
 
4567
 
//------------------------------------------------------------------------------
4568
 
int SNL_FEI_Structure::setPatternRowColOffsets(PatternDescriptor& pattern,
4569
 
                                              int numColIndices,
4570
 
                                              feiArray<int>& rowColOffsets)
4571
 
{
4572
 
  //This function sets the rowColOffsets array. The length of this array will
4573
 
  //be the number of row-indices associated with the pattern. The contents of
4574
 
  //this array will be offsets into an array of column-indices, which is of
4575
 
  //length num-row-IDs * numColIndices. The tricky thing here is that the
4576
 
  //order of those offsets varies depending on whether the pattern's interleave
4577
 
  //strategy is FEI_NODE_MAJOR or FEI_FIELD_MAJOR.
4578
 
 
4579
 
  rowColOffsets.resize(0);
4580
 
 
4581
 
  //First, get the row-field info and num-row-ids info.
4582
 
  int numRowIDs = pattern.getNumRowIDs();
4583
 
  feiArray<int>& numFieldsPerRow = pattern.getNumFieldsPerRow();
4584
 
  feiArray<int>* rowFieldIDs = pattern.getRowFieldIDs();
4585
 
 
4586
 
  //numFieldsPerRow is a list of length numRowIDs.
4587
 
  //rowFieldIDs is a list of feiArrays. There are 'numRowIDs' feiArrays, and
4588
 
  //the i-th feiArray is of length numFieldsPerRow[i].
4589
 
 
4590
 
  int interleave = pattern.getInterleaveStrategy();
4591
 
 
4592
 
 
4593
 
  //if interleave == FEI_NODE_MAJOR then we simply loop over the rowFieldIDs
4594
 
  //and increment the col-offset by numColIndices once each row-ID.
4595
 
  //
4596
 
  //if interleave == FEI_FIELD_MAJOR then we need to first create a flat list
4597
 
  //of the fieldIDs (each appearing only once) and then loop over the
4598
 
  //rowFieldIDs table once for each unique fieldID, and whenever that fieldID is
4599
 
  //encountered, append an offset to the list. The offset to be appended will be
4600
 
  //'current-row-id'*numColIndices.
4601
 
 
4602
 
  if (interleave == FEI_NODE_MAJOR) {
4603
 
    int colOffset = 0;
4604
 
    for(int rowID=0; rowID<numRowIDs; rowID++) {
4605
 
      for(int field=0; field<numFieldsPerRow[rowID]; field++) {
4606
 
        int numEqns = getFieldSize(rowFieldIDs[rowID][field]);
4607
 
        if (numEqns < 0) ERReturn(-1);
4608
 
 
4609
 
        for(int eqn=0; eqn<numEqns; eqn++) {
4610
 
          rowColOffsets.append(colOffset);
4611
 
        }
4612
 
      }
4613
 
      colOffset += numColIndices;
4614
 
    }
4615
 
  }
4616
 
  else if (interleave == FEI_FIELD_MAJOR) {
4617
 
    int i;
4618
 
    feiArray<int> fields;
4619
 
    for(i=0; i<numRowIDs; i++) {
4620
 
      for(int j=0; j<numFieldsPerRow[i]; j++) {
4621
 
        if (fields.find(rowFieldIDs[i][j]) < 0)
4622
 
          fields.append(rowFieldIDs[i][j]);
4623
 
      }
4624
 
    }
4625
 
 
4626
 
    for(i=0; i<fields.length(); i++) {
4627
 
      for(int rowID=0; rowID<numRowIDs; rowID++) {
4628
 
        for(int field=0; field<numFieldsPerRow[rowID]; field++) {
4629
 
          if (fields[i] == rowFieldIDs[rowID][field]) {
4630
 
            int numEqns = getFieldSize(fields[i]);
4631
 
            if (numEqns < 0) ERReturn(-1);
4632
 
 
4633
 
            for(int eqn=0; eqn<numEqns; eqn++) {
4634
 
              rowColOffsets.append(rowID*numColIndices);
4635
 
            }
4636
 
          }
4637
 
        }
4638
 
      }
4639
 
    }
4640
 
  }
4641
 
 
4642
 
  return(FEI_SUCCESS);
4643
 
}
4644
 
 
4645
 
//------------------------------------------------------------------------------
4646
 
int SNL_FEI_Structure::getPatternScatterIndices(int patternID,
4647
 
                                               const GlobalID* rowNodes,
4648
 
                                               feiArray<int>& rowIndices)
4649
 
{
4650
 
   int index = patternIDs_.find(patternID);
4651
 
 
4652
 
   if (index < 0) {
4653
 
      FEI_CERR << "SNL_FEI_Structure::getPatternScatterIndices: ERROR, patternID "
4654
 
           << (int)patternID << " not found." << FEI_ENDL;
4655
 
      return(-1);
4656
 
   }
4657
 
 
4658
 
   PatternDescriptor& pattern = *(patterns_[index]);
4659
 
 
4660
 
   int interleave = pattern.getInterleaveStrategy();
4661
 
 
4662
 
   int numRows = pattern.getNumRowIDs();
4663
 
   feiArray<int>& fieldsPerRow = pattern.getNumFieldsPerRow();
4664
 
   feiArray<int>* rowFieldIDs = pattern.getRowFieldIDs();
4665
 
 
4666
 
   NodeDescriptor** rNodeDesc = NULL;
4667
 
   NodeDescriptor** cNodeDesc = NULL;
4668
 
 
4669
 
   int err = getPatternNodeDescriptors(pattern, rowNodes, NULL,
4670
 
                                      rNodeDesc, cNodeDesc);
4671
 
   if (err) return(err);
4672
 
 
4673
 
   if (interleave == FEI_NODE_MAJOR) {
4674
 
     err = getNodeMajorIndices(rNodeDesc, numRows, 
4675
 
                                  rowFieldIDs, fieldsPerRow, rowIndices);
4676
 
     if (err) return(err);
4677
 
   }
4678
 
   else {
4679
 
     err = getFieldMajorIndices(rNodeDesc, numRows, 
4680
 
                                  rowFieldIDs, fieldsPerRow, rowIndices);
4681
 
     if (err) return(err);
4682
 
   }
4683
 
 
4684
 
   return(FEI_SUCCESS);
4685
 
}
4686
 
 
4687
 
//------------------------------------------------------------------------------
4688
3922
int SNL_FEI_Structure::getElemNodeDescriptors(int blockIndex, int elemIndex,
4689
3923
                                             NodeDescriptor** nodes)
4690
3924
{
4698
3932
 
4699
3933
  if (activeNodesInitialized_) {
4700
3934
    NodeDescriptor** elemNodePtrs =
4701
 
      connTables_[blockIndex]->elem_conn_ptrs->dataPtr() + elemIndex*numNodes;
 
3935
      &((*connTables_[blockIndex]->elem_conn_ptrs)[elemIndex*numNodes]);
4702
3936
    for(int i=0; i<numNodes; ++i) nodes[i] = elemNodePtrs[i];
4703
3937
  }
4704
3938
  else {
4705
3939
    const GlobalID* elemConn = 
4706
 
      connTables_[blockIndex]->elem_conn_ids->dataPtr() + elemIndex*numNodes;
 
3940
      &((*connTables_[blockIndex]->elem_conn_ids)[elemIndex*numNodes]);
4707
3941
    for(int i=0; i<numNodes; ++i) {
4708
3942
      CHK_ERR( nodeDatabase_->getNodeWithID(elemConn[i], nodes[i]));
4709
3943
    }
4713
3947
}
4714
3948
 
4715
3949
//------------------------------------------------------------------------------
4716
 
int SNL_FEI_Structure::getPatternNodeDescriptors(PatternDescriptor& pattern,
4717
 
                                                 const GlobalID* rowNodes,
4718
 
                                                 const GlobalID* colNodes,
4719
 
                                                 NodeDescriptor**& rNodeDesc,
4720
 
                                                 NodeDescriptor**& cNodeDesc)
4721
 
{
4722
 
  //This function's task is to obtain, from the nodeDatabase, the 
4723
 
  //NodeDescriptor objects for the specified "row-nodes" and "column-nodes".
4724
 
  //
4725
 
  //The PatternDescriptor is queried for the lengths of the 
4726
 
  //incoming array arguments.
4727
 
  //
4728
 
   int numRowNodes = pattern.getNumRowIDs();
4729
 
   int numColNodes = numRowNodes * pattern.getNumColIDsPerRow();
4730
 
   if (colNodes == NULL) numColNodes = 0;
4731
 
   int i, err;
4732
 
 
4733
 
   work_nodePtrs_.resize(numRowNodes+numColNodes);
4734
 
   work_nodePtrs_ = NULL;//feiArray::operator=
4735
 
   rNodeDesc = work_nodePtrs_.dataPtr();
4736
 
 
4737
 
   if (colNodes != NULL) {
4738
 
      cNodeDesc = rNodeDesc + numRowNodes;
4739
 
   }
4740
 
 
4741
 
   for(i=0; i<numRowNodes; i++) {
4742
 
     err = nodeDatabase_->getNodeWithID(rowNodes[i], rNodeDesc[i]);
4743
 
     if (err) return(err);
4744
 
   }
4745
 
 
4746
 
   for(i=0; i<numColNodes; i++) {
4747
 
     err = nodeDatabase_->getNodeWithID(colNodes[i], cNodeDesc[i]);
4748
 
     if (err) return(err);
4749
 
   }
4750
 
 
4751
 
   return(FEI_SUCCESS);
4752
 
}
4753
 
 
4754
 
//------------------------------------------------------------------------------
4755
3950
int SNL_FEI_Structure::getNodeIndices_simple(NodeDescriptor** nodes,
4756
3951
                                             int numNodes,
4757
3952
                                             int fieldID,
4924
4119
 
4925
4120
//------------------------------------------------------------------------------
4926
4121
int SNL_FEI_Structure::getNodeMajorIndices(NodeDescriptor** nodes, int numNodes,
4927
 
                                          feiArray<int>* fieldIDs,
4928
 
                                          feiArray<int>& fieldsPerNode,
4929
 
                                          feiArray<int>& scatterIndices)
 
4122
                                          std::vector<int>* fieldIDs,
 
4123
                                          std::vector<int>& fieldsPerNode,
 
4124
                                          std::vector<int>& scatterIndices)
4930
4125
{
4931
4126
   int offset = 0;
4932
4127
   scatterIndices.resize(0);
4939
4134
      const int* nodeEqnNums = node.getFieldEqnNumbers();
4940
4135
      int numFields = node.getNumFields();
4941
4136
 
4942
 
      int* fieldID_ind = fieldIDs[nodeIndex].dataPtr();
 
4137
      int* fieldID_ind = &(fieldIDs[nodeIndex][0]);
4943
4138
 
4944
4139
      for(int j=0; j<fieldsPerNode[nodeIndex]; j++) {
4945
4140
         int numEqns = getFieldSize(fieldID_ind[j]);
4961
4156
           }
4962
4157
 
4963
4158
           scatterIndices.resize(offset+numEqns);
4964
 
           int* inds = scatterIndices.dataPtr();
 
4159
           int* inds = &scatterIndices[0];
4965
4160
 
4966
4161
           for(int jj=0; jj<numEqns; jj++) {
4967
4162
             inds[offset+jj] = eqn+jj;
4992
4187
  //flat list, in the order we encounter them, but making sure no fieldID
4993
4188
  //gets added more than once.
4994
4189
 
4995
 
  int i;
4996
 
  feiArray<int> fields;
4997
 
  for(i=0; i<numNodes; i++) {
4998
 
    for(int j=0; j<fieldsPerNode[i]; j++) {
4999
 
      if (fields.find(fieldIDs[i][j]) < 0) fields.append(fieldIDs[i][j]);
 
4190
  std::vector<int> fields;
 
4191
  for(int ii=0; ii<numNodes; ii++) {
 
4192
    for(int j=0; j<fieldsPerNode[ii]; j++) {
 
4193
      if (std::find(fields.begin(), fields.end(), fieldIDs[ii][j]) == fields.end()) {
 
4194
        fields.push_back(fieldIDs[ii][j]);
 
4195
      }
5000
4196
    }
5001
4197
  }
5002
4198
 
5003
 
  int* fieldsPtr = fields.dataPtr();
 
4199
  int* fieldsPtr = fields.size()>0 ? &fields[0] : NULL;
5004
4200
 
5005
4201
  //ok, we've got our flat list of fields, so let's proceed to get the
5006
4202
  //scatter indices.
5007
4203
 
5008
 
  for(i=0; i<fields.length(); i++) {
 
4204
  for(size_t i=0; i<fields.size(); i++) {
5009
4205
    int field = fieldsPtr[i];
5010
4206
 
5011
4207
    for(int nodeIndex = 0; nodeIndex < numNodes; ++nodeIndex) {
5044
4240
 
5045
4241
//------------------------------------------------------------------------------
5046
4242
int SNL_FEI_Structure::getFieldMajorIndices(NodeDescriptor** nodes, int numNodes,
5047
 
                                          feiArray<int>* fieldIDs,
5048
 
                                          feiArray<int>& fieldsPerNode,
5049
 
                                          feiArray<int>& scatterIndices)
 
4243
                                          std::vector<int>* fieldIDs,
 
4244
                                          std::vector<int>& fieldsPerNode,
 
4245
                                          std::vector<int>& scatterIndices)
5050
4246
{
5051
4247
   //In this function we want to group equations by field, but
5052
4248
   //in what order should the fields be?
5054
4250
   //flat list, in the order we encounter them, but making sure no fieldID
5055
4251
   //gets added more than once.
5056
4252
 
5057
 
   int i;
5058
 
   feiArray<int> fields;
5059
 
   for(i=0; i<numNodes; i++) {
5060
 
      for(int j=0; j<fieldsPerNode[i]; j++) {
5061
 
         feiArray<int>& fldIDs = fieldIDs[i];
5062
 
         if (fields.find(fldIDs[j]) < 0) fields.append(fldIDs[j]);
 
4253
   std::vector<int> fields;
 
4254
   for(int ii=0; ii<numNodes; ii++) {
 
4255
      for(int j=0; j<fieldsPerNode[ii]; j++) {
 
4256
         std::vector<int>& fldIDs = fieldIDs[ii];
 
4257
         if (std::find(fields.begin(), fields.end(), fldIDs[j]) == fields.end()) {
 
4258
           fields.push_back(fldIDs[j]);
 
4259
         }
5063
4260
      }
5064
4261
   }
5065
4262
 
5069
4266
   int offset = 0;
5070
4267
   scatterIndices.resize(0);
5071
4268
 
5072
 
   for(i=0; i<fields.length(); i++) {
 
4269
   for(size_t i=0; i<fields.size(); i++) {
5073
4270
      for(int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
5074
4271
 
5075
4272
         const int* nodeFieldIDList = nodes[nodeIndex]->getFieldIDList();
5085
4282
 
5086
4283
         if (nind >= 0) {
5087
4284
            for(int jj=0; jj<numEqns; jj++) {
5088
 
               scatterIndices.append(nodeEqnNums[nind]+jj);
 
4285
               scatterIndices.push_back(nodeEqnNums[nind]+jj);
5089
4286
               offset++;
5090
4287
            }
5091
4288
         }