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