FEI Version of the Day
Loading...
Searching...
No Matches
fei_LinSysCoreFilter.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
11#include <limits>
12#include <cmath>
13#include <assert.h>
14
15#include "fei_utils.hpp"
16#include <fei_impl_utils.hpp>
17
18#include "fei_defs.h"
19
20#include <fei_ostream_ops.hpp>
21#include "fei_CommUtils.hpp"
22#include "fei_TemplateUtils.hpp"
23#include "snl_fei_Constraint.hpp"
25
26#include "fei_DirichletBCManager.hpp"
27#include "fei_FillableMat.hpp"
28#include "fei_CSVec.hpp"
29#include "FEI_Implementation.hpp"
30#include "fei_EqnCommMgr.hpp"
31#include "fei_ConnectivityTable.hpp"
32#include "fei_NodeDescriptor.hpp"
33#include "fei_NodeDatabase.hpp"
34#include "fei_BlockDescriptor.hpp"
35#include "SNL_FEI_Structure.hpp"
36#include "snl_fei_Utils.hpp"
37#include "fei_Data.hpp"
38#include "fei_LinearSystemCore.hpp"
39#include "fei_LinSysCore_flexible.hpp"
40
41#include "fei_LinSysCoreFilter.hpp"
42
43#undef fei_file
44#define fei_file "fei_LinSysCoreFilter.cpp"
45#include "fei_ErrMacros.hpp"
46
47#define ASSEMBLE_PUT 0
48#define ASSEMBLE_SUM 1
49
50
51//------------------------------------------------------------------------------
52LinSysCoreFilter::LinSysCoreFilter(FEI_Implementation* owner,
53 MPI_Comm comm,
54 SNL_FEI_Structure* probStruct,
56 int masterRank)
57 : Filter(probStruct),
58 timesInitializeCalled_(0),
59 lsc_(lsc),
60 useLookup_(true),
61 newMatrixData_(false),
62 newVectorData_(false),
63 newConstraintData_(false),
64 newBCData_(false),
65 connectivitiesInitialized_(false),
66 firstRemEqnExchange_(true),
67 needToCallMatrixLoadComplete_(false),
68 resolveConflictRequested_(false),
69 localStartRow_(0), //
70 localEndRow_(0), //Initialize all private variables here,
71 numLocalEqns_(0), //in the order that they are declared.
72 numGlobalEqns_(0), //(Don't bother initializing those that will
73 numLocallyOwnedNodes_(0), //be initialized in the body of the
74 numGlobalNodes_(0), //constructor below.)
75 firstLocalNodeNumber_(-1), //
76 blockMatrix_(false),
77 tooLateToChooseBlock_(false),
78 numLocalEqnBlks_(0),
79 localReducedBlkOffset_(0),
80 numLocalReducedEqnBlks_(0),
81 iterations_(0),
82 currentRHS_(0),
83 rhsIDs_(),
84 outputLevel_(0),
85 comm_(comm),
86 masterRank_(masterRank),
87 problemStructure_(probStruct),
88 matrixAllocated_(false),
89 rowIndices_(),
90 rowColOffsets_(0),
91 colIndices_(0),
92 nodeIDType_(0),
93 eqnCommMgr_(NULL),
94 eqnCommMgr_put_(NULL),
95 maxElemRows_(0),
96 scatterIndices_(),
97 blkScatterIndices_(),
98 iworkSpace_(),
99 iworkSpace2_(),
100 dworkSpace_(),
101 dworkSpace2_(),
102 eStiff_(NULL),
103 eStiff1D_(NULL),
104 eLoad_(NULL)
105{
106 localRank_ = fei::localProc(comm_);
107 numProcs_ = fei::numProcs(comm_);
108
109 internalFei_ = 0;
110
111 numRHSs_ = 1;
112 rhsIDs_.resize(numRHSs_);
113 rhsIDs_[0] = 0;
114
115 bcManager_ = new fei::DirichletBCManager(probStruct);
116 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
117 int err = createEqnCommMgr_put();
118
119 //We need to get the parameters from the owning FEI_Implementation...
120 int numParams = -1;
121 char** paramStrings = NULL;
122 err = owner->getParameters(numParams, paramStrings);
123
124 //Now let's pass them into our own parameter-handling mechanism.
125 err = parameters(numParams, paramStrings);
126 if (err != 0) {
127 fei::console_out() << "LinSysCoreFilter::LinSysCoreFilter ERROR, parameters failed." << FEI_ENDL;
128 MPI_Abort(comm_, -1);
129 }
130
131 Kid_ = new fei::FillableMat;
132 Kdi_ = new fei::FillableMat;
133 Kdd_ = new fei::FillableMat;
134 reducedEqnCounter_ = 0;
135 reducedRHSCounter_ = 0;
136
137 return;
138}
139
140//------------------------------------------------------------------------------
141LinSysCoreFilter::~LinSysCoreFilter() {
142//
143// Destructor function. Free allocated memory, etc.
144//
145 numRHSs_ = 0;
146
147 delete bcManager_;
148 delete eqnCommMgr_;
149 delete eqnCommMgr_put_;
150
151 delete [] eStiff_;
152 delete [] eStiff1D_;
153 delete [] eLoad_;
154
155 delete Kid_;
156 delete Kdi_;
157 delete Kdd_;
158}
159
160//------------------------------------------------------------------------------
161int LinSysCoreFilter::initialize()
162{
163// Determine final sparsity pattern for setting the structure of the
164// underlying sparse matrix.
165//
166 debugOutput("#LinSysCoreFilter::initialize");
167
168 numLocalEqns_ = problemStructure_->getNumLocalEqns();
169 numLocalEqnBlks_ = problemStructure_->getNumLocalEqnBlks();
170 numLocalReducedEqnBlks_ = problemStructure_->getNumLocalReducedEqnBlks();
171
172 // now, obtain the global equation info, such as how many equations there
173 // are globally, and what the local starting and ending row-numbers are.
174
175 // let's also get the number of global nodes, and a first-local-node-number.
176 // node-number is a globally 0-based number we are assigning to nodes.
177 // node-numbers are contiguous on a processor -- i.e., a processor owns a
178 // contiguous block of node-numbers. This provides an easier-to-work-with
179 // node numbering than the application-supplied nodeIDs which may not be
180 // assumed to be contiguous or 0-based, or anything else.
181
182 std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
183 localStartRow_ = eqnOffsets[localRank_];
184 localEndRow_ = eqnOffsets[localRank_+1]-1;
185 numGlobalEqns_ = eqnOffsets[numProcs_];
186
187 std::vector<int>& nodeOffsets = problemStructure_->getGlobalNodeOffsets();
188 firstLocalNodeNumber_ = nodeOffsets[localRank_];
189 numGlobalNodes_ = nodeOffsets[numProcs_];
190
191 //--------------------------------------------------------------------------
192 // ----- end active equation calculations -----
193
194 if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
195 eqnCommMgr_ = NULL;
196 if (eqnCommMgr_put_ != NULL) delete eqnCommMgr_put_;
197 eqnCommMgr_put_ = NULL;
198
199 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
200 if (eqnCommMgr_ == NULL) ERReturn(-1);
201
202 int err = createEqnCommMgr_put();
203 if (err != 0) ERReturn(err);
204
205 //(we need to set the number of RHSs in the eqn comm manager)
206 eqnCommMgr_->setNumRHSs(numRHSs_);
207
208 //let's let the underlying linear system know what the global offsets are.
209 //While we're dealing with global offsets, we'll also calculate the starting
210 //and ending 'reduced' rows, etc.
211 CHK_ERR( initLinSysCore() );
212
213 setLinSysCoreCREqns();
214
215 if (timesInitializeCalled_ == 0) {
216 //
217 // let's prepare some arrays for handing the matrix structure to
218 // the linear system.
219
220 std::vector<int> rowLengths;
221 CHK_ERR( problemStructure_->getMatrixRowLengths(rowLengths) );
222
223 int numReducedEqns = problemStructure_->getNumLocalReducedEqns();
224 int maxBlkSize = problemStructure_->getGlobalMaxBlkSize();
225 std::vector<int> blkSizes(numLocalReducedEqnBlks_, 1);
226
227 int numNonzeros = 0;
228 for(size_t ii=0; ii<rowLengths.size(); ++ii) {
229 numNonzeros += rowLengths[ii];
230 }
231
232 std::vector<int> indices_1D(numNonzeros);
233 std::vector<int*> indices(numReducedEqns);
234
235 int offset = 0;
236 for(size_t ii=0; ii<rowLengths.size(); ++ii) {
237 indices[ii] = &(indices_1D[offset]);
238 offset += rowLengths[ii];
239 }
240
241 if (maxBlkSize == 0) ERReturn(-1);
242
243 if (maxBlkSize == 1) {
244 CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths) );
245
246 debugOutput("#LinSysCoreFilter calling point lsc_->setMatrixStructure");
247 CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
248 &indices[0], &rowLengths[0], &blkSizes[0]) );
249 }
250 else {
251 std::vector<int> blkRowLengths;
252 int* blkIndices_1D = numNonzeros > 0 ? new int[numNonzeros] : NULL;
253 int** blkIndices = numLocalReducedEqnBlks_ ?
254 new int*[numLocalReducedEqnBlks_] : NULL;
255 if (blkIndices == NULL && numLocalReducedEqnBlks_ != 0) ERReturn(-1);
256
257 CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths,
258 blkIndices, blkIndices_1D,
259 blkRowLengths, blkSizes) );
260
261 offset = 0;
262 for(int ii=0; ii<numLocalReducedEqnBlks_; ++ii) {
263 blkIndices[ii] = &(blkIndices_1D[offset]);
264 offset += blkRowLengths[ii];
265 }
266
267 debugOutput("#LinSysCoreFilter calling block lsc_->setMatrixStructure");
268 CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
269 blkIndices, &blkRowLengths[0], &blkSizes[0]) );
270
271 if (numLocalReducedEqnBlks_ != 0) delete [] blkIndices;
272 if (numNonzeros > 0) delete [] blkIndices_1D;
273 }
274 }
275
276 matrixAllocated_ = true;
277
278 debugOutput("#leaving LinSysCoreFilter::initialize");
279 ++timesInitializeCalled_;
280 return(FEI_SUCCESS);
281}
282
283//------------------------------------------------------------------------------
284int LinSysCoreFilter::createEqnCommMgr_put()
285{
286 if (eqnCommMgr_put_ != NULL) return(0);
287
288 eqnCommMgr_put_ = eqnCommMgr_->deepCopy();
289 if (eqnCommMgr_put_ == NULL) ERReturn(-1);
290
291 eqnCommMgr_put_->resetCoefs();
292 eqnCommMgr_put_->accumulate_ = false;
293 return(0);
294}
295
296//==============================================================================
297int LinSysCoreFilter::initLinSysCore()
298{
299 debugOutput("#LinSysCoreFilter calling lsc_->setLookup");
300 int err = lsc_->setLookup(*problemStructure_);
301
302 if (err != 0) {
303 useLookup_ = false;
304 }
305
306 std::vector<int>& globalNodeOffsets = problemStructure_->getGlobalNodeOffsets();
307 std::vector<int>& globalEqnOffsets = problemStructure_->getGlobalEqnOffsets();
308 std::vector<int>& globalBlkEqnOffsets =
309 problemStructure_->getGlobalBlkEqnOffsets();
310
311 int startRow = localStartRow_;
312 while(problemStructure_->isSlaveEqn(startRow)) startRow++;
313 int endRow = localEndRow_;
314 while(problemStructure_->isSlaveEqn(endRow)) endRow--;
315
316 problemStructure_->translateToReducedEqn(startRow, reducedStartRow_);
317 problemStructure_->translateToReducedEqn(endRow, reducedEndRow_);
318 numReducedRows_ = reducedEndRow_ - reducedStartRow_ + 1;
319
320 std::vector<int> reducedEqnOffsets(globalEqnOffsets.size());
321 std::vector<int> reducedBlkEqnOffsets(globalBlkEqnOffsets.size());
322 std::vector<int> reducedNodeOffsets(globalNodeOffsets.size());
323
324 int numSlaveEqns = problemStructure_->numSlaveEquations();
325
326 int reducedNodeNum =
327 problemStructure_->getAssociatedNodeNumber(reducedStartRow_);
328
329 std::vector<int> tmpSend(2), tmpRecv, tmpRecvLengths;
330 tmpSend[0] = reducedStartRow_;
331 tmpSend[1] = reducedNodeNum;
332 CHK_ERR( fei::Allgatherv(comm_, tmpSend, tmpRecvLengths, tmpRecv) );
333
334 for(size_t ii=0; ii<globalEqnOffsets.size()-1; ++ii) {
335 reducedEqnOffsets[ii] = tmpRecv[ii*2];
336 reducedNodeOffsets[ii] = tmpRecv[ii*2+1];
337
338 //Major assumption: we're assuming that if there are slave equations, then
339 //there are no lagrange-multiplier constraints. Because if there are no
340 //lagrange-multiplier constraints, then blkEqn == nodeNum.
341 if (numSlaveEqns > 0) {
342 reducedBlkEqnOffsets[ii] = reducedNodeOffsets[ii];
343 }
344 else {
345 reducedBlkEqnOffsets[ii] = globalBlkEqnOffsets[ii];
346 }
347 }
348
349 if (localRank_ == numProcs_-1) {
350 reducedNodeNum = problemStructure_->
351 translateToReducedNodeNumber(globalNodeOffsets[numProcs_]-1, localRank_);
352 int blkEqn = globalBlkEqnOffsets[numProcs_]-1;
353 if (numSlaveEqns > 0) {
354 blkEqn = reducedNodeNum;
355 }
356
357 tmpSend.resize(3);
358 tmpSend[0] = reducedEndRow_;
359 tmpSend[1] = reducedNodeNum;
360 tmpSend[2] = blkEqn;
361 }
362 else {
363 tmpSend.resize(3);
364 }
365
366 CHK_ERR( fei::Bcast(comm_, tmpSend, numProcs_-1) );
367 reducedEqnOffsets[numProcs_] = tmpSend[0]+1;
368 reducedNodeOffsets[numProcs_] = tmpSend[1]+1;
369 reducedBlkEqnOffsets[numProcs_] = tmpSend[2]+1;
370
371 debugOutput("#LinSysCoreFilter calling lsc_->setGlobalOffsets");
372 CHK_ERR( lsc_->setGlobalOffsets(numProcs_+1,
373 &reducedNodeOffsets[0],
374 &reducedEqnOffsets[0],
375 &reducedBlkEqnOffsets[0]) );
376
377 if (connectivitiesInitialized_) return(0);
378
379 int numBlocks = problemStructure_->getNumElemBlocks();
380 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
381 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
382
383 int numNodes = nodeDB.getNumNodeDescriptors();
384 int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
385 nodeCommMgr.getLocalNodeIDs().size();
386 numNodes -= numRemoteNodes;
387
388 std::vector<int> numElemsPerBlock(numBlocks);
389 std::vector<int> numNodesPerElem(numBlocks);
390 std::vector<int> numDofPerNode(0,1);
391
392 for(int blk=0; blk<numBlocks; blk++) {
393 BlockDescriptor* block = NULL;
394 CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
395
396 numElemsPerBlock[blk] = block->getNumElements();
397 numNodesPerElem[blk] = block->getNumNodesPerElement();
398
399 int* fieldsPerNode = block->fieldsPerNodePtr();
400 int** fieldIDsTable = block->fieldIDsTablePtr();
401
402 for(int nn=0; nn<numNodesPerElem[blk]; nn++) {
403 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
404 numDofPerNode.push_back(0);
405 int indx = numDofPerNode.size()-1;
406
407 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
408 numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
409 }
410 }
411 }
412
413 for(int i=0; i<numBlocks; i++) {
414 BlockDescriptor* block = NULL;
415 CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
416
417 if (block->getNumElements() == 0) continue;
418
419 ConnectivityTable& ctbl =
420 problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
421
422 std::vector<int> cNodeList(block->getNumNodesPerElement());
423
424 int* fieldsPerNode = block->fieldsPerNodePtr();
425 int** fieldIDsTable = block->fieldIDsTablePtr();
426
427 numDofPerNode.resize(0);
428 for(int nn=0; nn<numNodesPerElem[i]; nn++) {
429 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
430 numDofPerNode.push_back(0);
431 int indx = numDofPerNode.size()-1;
432
433 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
434 numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
435 }
436 }
437
438 int nodesPerElement = block->getNumNodesPerElement();
439 NodeDescriptor** elemNodePtrs = &((*ctbl.elem_conn_ptrs)[0]);
440 int offset = 0;
441 for(int j=0; j<block->getNumElements(); j++) {
442
443 for(int k=0; k<nodesPerElement; k++) {
444 NodeDescriptor* node = elemNodePtrs[offset++];
445 cNodeList[k] = node->getNodeNumber();
446 }
447
448 int elemID = ctbl.elemNumbers[j];
449 int* nodeNumbers = &cNodeList[0];
450 debugOutput("#LinSysCoreFilter calling lsc->setConnectivities");
451 CHK_ERR( lsc_->setConnectivities(i, 1, nodesPerElement,
452 &elemID, &nodeNumbers) );
453 }
454 }
455
456 connectivitiesInitialized_ = true;
457
458 return(FEI_SUCCESS);
459}
460
461//==============================================================================
462void LinSysCoreFilter::setLinSysCoreCREqns()
463{
464 int err, i=0;
465
466 std::vector<int> iwork;
467
468 std::map<GlobalID,ConstraintType*>::const_iterator
469 cr_iter = problemStructure_->getMultConstRecords().begin(),
470 cr_end = problemStructure_->getMultConstRecords().end();
471
472 while(cr_iter != cr_end) {
473 ConstraintType& constraint = *((*cr_iter).second);
474 int numNodesPerCR = constraint.getMasters().size();
475 int meqn = constraint.getEqnNumber();
476
477 std::vector<GlobalID>& nodeID_vec = constraint.getMasters();
478 GlobalID* nodeIDPtr = &nodeID_vec[0];
479
480 if ((int)iwork.size() < 2*numNodesPerCR) {
481 iwork.resize(2*numNodesPerCR);
482 }
483
484 int* nodeList = &(iwork[0]);
485
486 int* eqnList = nodeList+numNodesPerCR;
487
488 std::vector<int>& fieldIDs_vec = constraint.getMasterFieldIDs();
489 int* fieldIDs = &fieldIDs_vec[0];
490
491 for(int k=0; k<numNodesPerCR; k++) {
492 const NodeDescriptor *node = Filter::findNode(nodeIDPtr[k]);
493 if(node == NULL)
494 {
495 nodeList[k] = -1; // Indicates that the node wasn't found
496 }
497 else
498 {
499 nodeList[k] = node->getNodeNumber();
500 }
501
502 int eqn = -1; // Indicates that the equation wasn't found.
503 if ( node ) {
504 node->getFieldEqnNumber(fieldIDs[k], eqn);
505 }
506 eqnList[k] = eqn;
507 }
508
509 int crMultID = constraint.getConstraintID() + i++;
510 if (Filter::logStream() != NULL) {
511 FEI_OSTREAM& os = *logStream();
512 os << "#LinSysCoreFilter calling lsc_->setMultCREqns"<<FEI_ENDL;
513 os << "# multiplier eqn: " << meqn << ", columns: ";
514 for(int j=0; j<numNodesPerCR; ++j) os << eqnList[j] << " ";
515 os << FEI_ENDL;
516 }
517
518 err = lsc_->setMultCREqns(crMultID, 1, numNodesPerCR,
519 &nodeList, &eqnList,
520 fieldIDs, &meqn);
521 if (err) voidERReturn;
522 ++cr_iter;
523 }
524
525 LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
526 if (lscf != NULL) {
527 debugOutput("LinSysCoreFilter calling lscf->setMultCRComplete");
528 err = lscf->setMultCRComplete();
529 if (err != 0) {
530 fei::console_out() << "LinSysCoreFilter::setLinSysCoreCREqns ERROR returned from "
531 << "lscf->setMultCRComplete()" << FEI_ENDL;
532 }
533 }
534
535 //
536 //Now the penalty CRs...
537 //
538 cr_iter = problemStructure_->getPenConstRecords().begin();
539 cr_end = problemStructure_->getPenConstRecords().end();
540
541 while(cr_iter != cr_end) {
542 ConstraintType& crset = *((*cr_iter).second);
543 int numNodesPerCR = crset.getMasters().size();
544
545 std::vector<GlobalID>& nodeIDs_vec = crset.getMasters();
546 GlobalID* nodeIDsPtr = &nodeIDs_vec[0];
547
548 if ((int)iwork.size() < 2*numNodesPerCR) {
549 iwork.resize(2*numNodesPerCR);
550 }
551
552 int* nodeList = &(iwork[0]);
553
554 int* eqnList = nodeList+numNodesPerCR;
555
556 std::vector<int>& fieldIDs_vec = crset.getMasterFieldIDs();
557 int* fieldIDs = &fieldIDs_vec[0];
558
559 for(int k=0; k<numNodesPerCR; k++) {
560 const NodeDescriptor& node = Filter::findNodeDescriptor(nodeIDsPtr[k]);
561 nodeList[k] = node.getNodeNumber();
562
563 int eqn = -1;
564 if (!node.getFieldEqnNumber(fieldIDs[k], eqn)) voidERReturn;
565
566 eqnList[k] = eqn;
567 }
568
569 int crPenID = crset.getConstraintID() + i;
570 err = lsc_->setPenCREqns(crPenID, 1, numNodesPerCR,
571 &nodeList, &eqnList,
572 fieldIDs);
573 if (err) voidERReturn;
574 ++cr_iter;
575 }
576}
577
578//==============================================================================
579int LinSysCoreFilter::storeNodalColumnCoefs(int eqn, const NodeDescriptor& node,
580 int fieldID, int fieldSize,
581 double* coefs)
582{
583 //
584 //This function stores the coeficients for 'node' at 'fieldID' at the correct
585 //column indices in row 'eqn' of the system matrix.
586 //
587 if ((localStartRow_ > eqn) || (eqn > localEndRow_)) return(0);
588
589 int eqnNumber = -1;
590 if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
591
592 int numParams = fieldSize;
593 if (numParams < 1) {
594 return(0);
595 }
596
597 if ((int)iworkSpace2_.size() < numParams) {
598 iworkSpace2_.resize(numParams);
599 }
600
601 int* cols = &iworkSpace2_[0];
602
603 for(int j=0; j<numParams; j++) {
604 cols[j] = eqnNumber + j;
605 }
606
607 CHK_ERR( giveToMatrix(1, &eqn, numParams, cols, &coefs, ASSEMBLE_SUM) );
608
609 return(FEI_SUCCESS);
610}
611
612//==============================================================================
613int LinSysCoreFilter::storeNodalRowCoefs(const NodeDescriptor& node,
614 int fieldID, int fieldSize,
615 double* coefs, int eqn)
616{
617 //
618 //This function stores coeficients in the equations for 'node', 'fieldID' at
619 //column index 'eqn' of the system matrix.
620 //
621 int eqnNumber = -1;
622 if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
623
624 if ((localStartRow_ > eqnNumber) || (eqnNumber > localEndRow_)) return(0);
625
626 int numParams = fieldSize;
627
628 if (numParams < 1) {
629 return(0);
630 }
631
632 if ((int)iworkSpace2_.size() < numParams) {
633 iworkSpace2_.resize(numParams);
634 }
635
636 int* ptRows = &iworkSpace2_[0];
637
638 if ((int)dworkSpace2_.size() < numParams) {
639 dworkSpace2_.resize(numParams);
640 }
641
642 const double* * values = &dworkSpace2_[0];
643
644 for(int j=0; j<numParams; j++) {
645 ptRows[j] = eqnNumber + j;
646 values[j] = &(coefs[j]);
647 }
648
649 CHK_ERR( giveToMatrix(numParams, ptRows, 1, &eqn, values, ASSEMBLE_SUM) );
650
651 return(FEI_SUCCESS);
652}
653
654//==============================================================================
655void LinSysCoreFilter::storeNodalSendEqn(const NodeDescriptor& node,
656 int fieldID, int col,
657 double* coefs)
658{
659 //
660 //This is a private LinSysCoreFilter function. We can safely assume that
661 //it will only be called with a node that is not locally owned.
662 //
663 //This function tells the eqn comm mgr that we'll be sending contributions
664 //to column 'col' for the equations associated with 'fieldID', on 'node', on
665 //node's owning processor.
666 //
667 int proc = node.getOwnerProc();
668
669 int eqnNumber = -1;
670 if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
671
672 int numEqns = problemStructure_->getFieldSize(fieldID);
673
674 for(int i=0; i<numEqns; i++) {
675 eqnCommMgr_->addRemoteEqn(eqnNumber+i, proc, &coefs[i], &col, 1);
676 }
677}
678
679//==============================================================================
680void LinSysCoreFilter::storePenNodeSendData(const NodeDescriptor& iNode,
681 int iField, int iFieldSize,
682 double* iCoefs,
683 const NodeDescriptor& jNode,
684 int jField, int jFieldSize,
685 double* jCoefs,
686 double penValue, double CRValue)
687{
688//
689//This function will register with the eqn comm mgr the equations associated
690//with iNode, field 'iField' having column indices that are the equations
691//associated with jNode, field 'jField', to be sent to the owner of iNode.
692//
693 int proc = iNode.getOwnerProc();
694
695 int iEqn = -1, jEqn = -1;
696 if (!iNode.getFieldEqnNumber(iField, iEqn)) voidERReturn;
697 if (!jNode.getFieldEqnNumber(jField, jEqn)) voidERReturn;
698
699 int iNumParams = iFieldSize;
700 int jNumParams = jFieldSize;
701 if (iNumParams < 1 || jNumParams < 1) {
702 fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
703 << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
704 << jNumParams<<FEI_ENDL;
705 voidERReturn;
706 }
707
708 if ((int)dworkSpace_.size() < jNumParams) {
709 dworkSpace_.resize(jNumParams);
710 }
711
712 double* coefs = &dworkSpace_[0];
713
714 if ((int)iworkSpace2_.size() < jNumParams) {
715 iworkSpace2_.resize(jNumParams);
716 }
717
718 int* cols = &iworkSpace2_[0];
719
720 for(int i=0; i<iNumParams; i++) {
721 for(int j=0; j<jNumParams; j++) {
722 cols[j] = jEqn + j;
723 coefs[j] = penValue*iCoefs[i]*jCoefs[j];
724 }
725
726 int row = iEqn + i;
727 eqnCommMgr_->addRemoteEqn(row, proc, coefs, cols, jNumParams);
728
729 double rhsValue = penValue*iCoefs[i]*CRValue;
730 eqnCommMgr_->addRemoteRHS(row, proc, currentRHS_, rhsValue);
731 }
732}
733
734//==============================================================================
735int LinSysCoreFilter::storePenNodeData(const NodeDescriptor& iNode,
736 int iField, int iFieldSize,
737 double* iCoefs,
738 const NodeDescriptor& jNode,
739 int jField, int jFieldSize,
740 double* jCoefs,
741 double penValue, double CRValue){
742//
743//This function will add to the local matrix the penalty constraint equations
744//associated with iNode at iField, having column indices that are the
745//equations associated with jNode at jField.
746//Also, add the penalty contribution to the RHS vector.
747//
748 int iEqn = -1, jEqn = -1;
749 if (!iNode.getFieldEqnNumber(iField, iEqn)) ERReturn(FEI_ID_NOT_FOUND);
750 if (!jNode.getFieldEqnNumber(jField, jEqn)) ERReturn(FEI_ID_NOT_FOUND);
751
752 int iNumParams = iFieldSize;
753 int jNumParams = jFieldSize;
754 if (iNumParams < 1 || jNumParams < 1) {
755 fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
756 << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
757 << jNumParams<<FEI_ENDL;
758 ERReturn(-1);
759 }
760
761 if ((int)dworkSpace2_.size() < iNumParams) {
762 dworkSpace2_.resize(iNumParams);
763 }
764
765 const double* * coefs = &dworkSpace2_[0];
766
767 if ((int)dworkSpace_.size() < iNumParams * jNumParams) {
768 dworkSpace_.resize(iNumParams * jNumParams);
769 }
770
771 double* coefList = &dworkSpace_[0];
772
773 if ((int)iworkSpace2_.size() < jNumParams+iNumParams) {
774 iworkSpace2_.resize(jNumParams+iNumParams);
775 }
776
777 int* cols = &iworkSpace2_[0];
778 int* rows = &iworkSpace2_[jNumParams];
779
780
781 for(int i=0; i<iNumParams; i++) {
782 double* coefPtr = coefList + i*jNumParams;
783 coefs[i] = coefPtr;
784
785 rows[i] = iEqn + i;
786
787 for(int j=0; j<jNumParams; j++) {
788 if (i==0) cols[j] = jEqn + j;
789 coefPtr[j] = penValue*iCoefs[i]*jCoefs[j];
790 }
791
792 double rhsValue = penValue*iCoefs[i]*CRValue;
793 CHK_ERR( giveToRHS(1, &rhsValue, &rows[i], ASSEMBLE_SUM) );
794 }
795
796 CHK_ERR( giveToMatrix(iNumParams, rows,
797 jNumParams, cols, coefs, ASSEMBLE_SUM) );
798
799 return(FEI_SUCCESS);
800}
801
802//------------------------------------------------------------------------------
803int LinSysCoreFilter::resetSystem(double s)
804{
805 //
806 // This puts the value s throughout both the matrix and the vector.
807 //
808 if (Filter::logStream() != NULL) {
809 (*logStream()) << "FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
810 }
811
812 CHK_ERR( resetTheMatrix(s) );
813 CHK_ERR( resetTheRHSVector(s) );
814
815 debugOutput("#LinSysCoreFilter leaving resetSystem");
816
817 return(FEI_SUCCESS);
818}
819
820//------------------------------------------------------------------------------
821int LinSysCoreFilter::deleteMultCRs()
822{
823 debugOutput("#LinSysCoreFilter::deleteMultCRs");
824
825 LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
826 if (lscf == NULL) {
827// fei::console_out() << "FEI::LinSysCoreFilter: ERROR deleteMultCRs requested, but "
828// << "underlying solver doesn't support this operation." << FEI_ENDL;
829 return(-1);
830 }
831
832 int err = lscf->resetConstraints(0.0);
833
834 debugOutput("#LinSysCoreFilter leaving deleteMultCRs");
835
836 return(err);
837}
838
839//------------------------------------------------------------------------------
840int LinSysCoreFilter::resetTheMatrix(double s)
841{
842 CHK_ERR( lsc_->resetMatrix(s) );
843
844 return(FEI_SUCCESS);
845}
846
847//------------------------------------------------------------------------------
848int LinSysCoreFilter::resetTheRHSVector(double s)
849{
850 CHK_ERR( lsc_->resetRHSVector(s) );
851
852 return(FEI_SUCCESS);
853}
854
855//------------------------------------------------------------------------------
856int LinSysCoreFilter::resetMatrix(double s) {
857//
858// This puts the value s throughout both the matrix and the vector.
859//
860
861 debugOutput("FEI: resetMatrix");
862
863 CHK_ERR( resetTheMatrix(s) )
864
865 //clear away any boundary condition data.
866 bcManager_->clearAllBCs();
867
868 eqnCommMgr_->resetCoefs();
869
870 debugOutput("#LinSysCoreFilter leaving resetMatrix");
871
872 return(FEI_SUCCESS);
873}
874
875//------------------------------------------------------------------------------
876int LinSysCoreFilter::resetRHSVector(double s) {
877//
878// This puts the value s throughout the rhs vector.
879//
880
881 debugOutput("FEI: resetRHSVector");
882
883 CHK_ERR( resetTheRHSVector(s) )
884
885 //clear away any boundary condition data.
886 bcManager_->clearAllBCs();
887
888 eqnCommMgr_->resetCoefs();
889
890 debugOutput("# LinSysCoreFilter leaving resetRHSVector");
891
892 return(FEI_SUCCESS);
893}
894
895//------------------------------------------------------------------------------
896int LinSysCoreFilter::resetInitialGuess(double s) {
897//
898// This puts the value s throughout the initial guess (solution) vector.
899//
900 if (Filter::logStream() != NULL) {
901 FEI_OSTREAM& os = *logStream();
902 os << "FEI: resetInitialGuess" << FEI_ENDL;
903 os << "#value to which initial guess is to be set" << FEI_ENDL;
904 os << s << FEI_ENDL;
905 }
906
907 int* eqns = new int[numReducedRows_];
908 double* values = new double[numReducedRows_];
909 if (eqns == NULL || values == NULL) return(-1);
910
911 for(int i=0; i<numReducedRows_; ++i) {
912 eqns[i] = reducedStartRow_ + i;
913 values[i] = s;
914 }
915
916 CHK_ERR( lsc_->putInitialGuess(eqns, values, numReducedRows_) );
917
918 delete [] eqns;
919 delete [] values;
920
921 debugOutput("# LinSysCoreFilter leaving resetInitialGuess");
922
923 return(FEI_SUCCESS);
924}
925
926//------------------------------------------------------------------------------
927int LinSysCoreFilter::loadNodeBCs(int numNodes,
928 const GlobalID *nodeIDs,
929 int fieldID,
930 const int* offsetsIntoField,
931 const double* prescribedValues)
932{
933 //
934 // load boundary condition information for a given set of nodes
935 //
936 int size = problemStructure_->getFieldSize(fieldID);
937 if (size < 1) {
938 fei::console_out() << "FEI Warning: loadNodeBCs called for fieldID "<<fieldID
939 <<", which was defined with size "<<size<<" (should be positive)."<<FEI_ENDL;
940 return(0);
941 }
942
943 if (Filter::logStream() != NULL) {
944 (*logStream())<<"FEI: loadNodeBCs"<<FEI_ENDL
945 <<"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
946 <<"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
947 <<"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
948 (*logStream())<<"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
949
950 for(int j=0; j<numNodes; j++) {
951 GlobalID nodeID = nodeIDs[j];
952 (*logStream())<<static_cast<int>(nodeID)<<" ";
953 (*logStream())<<offsetsIntoField[j]<<" "<<prescribedValues[j];
954 (*logStream())<<FEI_ENDL;
955 }
956 }
957
958 if (numNodes > 0) newBCData_ = true;
959
960 bcManager_->addBCRecords(numNodes, nodeIDType_, fieldID, nodeIDs,
961 offsetsIntoField, prescribedValues);
962
963 return(FEI_SUCCESS);
964}
965
966//------------------------------------------------------------------------------
967int LinSysCoreFilter::loadElemBCs(int numElems,
968 const GlobalID *elemIDs,
969 int fieldID,
970 const double *const *alpha,
971 const double *const *beta,
972 const double *const *gamma)
973{
974 return(-1);
975}
976
977//------------------------------------------------------------------------------
978void LinSysCoreFilter::allocElemStuff()
979{
980 int nb = problemStructure_->getNumElemBlocks();
981
982 for(int i=0; i<nb; i++) {
983 BlockDescriptor* block = NULL;
984 int err = problemStructure_->getBlockDescriptor_index(i, block);
985 if (err) voidERReturn;
986
987 int numEqns = block->getNumEqnsPerElement();
988 if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
989 }
990
991 scatterIndices_.resize(maxElemRows_);
992
993 if (eStiff_ != NULL) delete [] eStiff_;
994 if (eStiff1D_ != NULL) delete [] eStiff1D_;
995 if (eLoad_ != NULL) delete [] eLoad_;
996
997 eStiff_ = new double*[maxElemRows_];
998 eStiff1D_ = new double[maxElemRows_*maxElemRows_];
999
1000 if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
1001
1002 for(int r=0; r<maxElemRows_; r++) {
1003 eStiff_[r] = eStiff1D_ + r*maxElemRows_;
1004 }
1005
1006 eLoad_ = new double[maxElemRows_];
1007
1008 if (eLoad_ == NULL) voidERReturn
1009}
1010
1011//------------------------------------------------------------------------------
1012int LinSysCoreFilter::sumInElem(GlobalID elemBlockID,
1013 GlobalID elemID,
1014 const GlobalID* elemConn,
1015 const double* const* elemStiffness,
1016 const double* elemLoad,
1017 int elemFormat)
1018{
1019 if (Filter::logStream() != NULL && outputLevel_ > 2) {
1020 (*logStream()) << "FEI: sumInElem" << FEI_ENDL <<"#blkID" << FEI_ENDL
1021 << static_cast<int>(elemBlockID) << FEI_ENDL
1022 << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
1023 BlockDescriptor* block = NULL;
1024 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1025 int numNodes = block->getNumNodesPerElement();
1026 (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
1027 (*logStream()) << "#nodes" << FEI_ENDL;
1028 for(int i=0; i<numNodes; ++i) {
1029 GlobalID nodeID = elemConn[i];
1030 (*logStream())<<static_cast<int>(nodeID)<<" ";
1031 }
1032 (*logStream())<<FEI_ENDL;
1033 }
1034
1035 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
1036 elemLoad, elemFormat));
1037}
1038
1039//------------------------------------------------------------------------------
1040int LinSysCoreFilter::sumInElemMatrix(GlobalID elemBlockID,
1041 GlobalID elemID,
1042 const GlobalID* elemConn,
1043 const double* const* elemStiffness,
1044 int elemFormat)
1045{
1046 if (Filter::logStream() != NULL && outputLevel_ > 2) {
1047 (*logStream()) << "FEI: sumInElemMatrix"<<FEI_ENDL
1048 << "#blkID" << FEI_ENDL << static_cast<int>(elemBlockID) << FEI_ENDL
1049 << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
1050 BlockDescriptor* block = NULL;
1051 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1052 int numNodes = block->getNumNodesPerElement();
1053 (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
1054 (*logStream()) << "#nodes" << FEI_ENDL;
1055 for(int i=0; i<numNodes; ++i) {
1056 GlobalID nodeID = elemConn[i];
1057 (*logStream())<<static_cast<int>(nodeID)<<" ";
1058 }
1059 (*logStream())<<FEI_ENDL;
1060 }
1061
1062 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
1063 NULL, elemFormat));
1064}
1065
1066//------------------------------------------------------------------------------
1067int LinSysCoreFilter::sumInElemRHS(GlobalID elemBlockID,
1068 GlobalID elemID,
1069 const GlobalID* elemConn,
1070 const double* elemLoad)
1071{
1072 if (Filter::logStream() != NULL && outputLevel_ > 2) {
1073 (*logStream()) << "FEI: sumInElemRHS"<<FEI_ENDL<<"#blID" << FEI_ENDL
1074 <<(int)elemBlockID << FEI_ENDL
1075 << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
1076 BlockDescriptor* block = NULL;
1077 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1078 int numNodes = block->getNumNodesPerElement();
1079 (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
1080 (*logStream()) << "#nodes" << FEI_ENDL;
1081 for(int i=0; i<numNodes; ++i) {
1082 GlobalID nodeID = elemConn[i];
1083 (*logStream())<<static_cast<int>(nodeID)<<" ";
1084 }
1085 (*logStream())<<FEI_ENDL;
1086 }
1087
1088 return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
1089 elemLoad, -1));
1090}
1091
1092//------------------------------------------------------------------------------
1093int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
1094 GlobalID elemID,
1095 const GlobalID* elemConn,
1096 const double* const* elemStiffness,
1097 const double* elemLoad,
1098 int elemFormat)
1099{
1100 (void)elemConn;
1101 return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
1102 elemFormat) );
1103}
1104
1105//------------------------------------------------------------------------------
1106int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
1107 GlobalID elemID,
1108 const double* const* elemStiffness,
1109 const double* elemLoad,
1110 int elemFormat)
1111{
1112 //first get the block-descriptor for this elemBlockID...
1113
1114 BlockDescriptor* block = NULL;
1115 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1116
1117 //now allocate our local stiffness/load copy if we haven't already.
1118
1119 if (maxElemRows_ <= 0) allocElemStuff();
1120
1121 int numElemRows = block->getNumEqnsPerElement();
1122 int numBlkElemRows = block->getNumBlkEqnsPerElement();
1123 int interleave = block->getInterleaveStrategy();
1124
1125 //an std::vector.resize operation is free if the size is either shrinking or
1126 //staying the same.
1127
1128 scatterIndices_.resize(numElemRows);
1129 blkScatterIndices_.resize(numBlkElemRows*2);
1130
1131 const double* const* stiff = NULL;
1132 if (elemStiffness != NULL) stiff = elemStiffness;
1133
1134 const double* load = NULL;
1135 if (elemLoad != NULL) load = elemLoad;
1136
1137 //we'll make a local dense copy of the element stiffness array
1138 //if the stiffness array was passed in using elemFormat != FEI_DENSE_ROW
1139 //AND if the stiffness array is non-null.
1140 //Note that this is not optimal if elemFormat is one of the
1141 //diagonal or block-diagonal formats.
1142 if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
1143 if (elemFormat == FEI_BLOCK_DIAGONAL_ROW ||
1144 elemFormat == FEI_BLOCK_DIAGONAL_COL) {
1145 fei::console_out() << "LinSysCoreFilter::generalElemInput ERROR, elemFormat="
1146 << elemFormat << " not supported."<<FEI_ENDL;
1147 ERReturn(-1);
1148 }
1149
1150 Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
1151 stiff = eStiff_;
1152 }
1153
1154 if (stiff != NULL) newMatrixData_ = true;
1155 if (load != NULL) newVectorData_ = true;
1156
1157 if (Filter::logStream() != NULL && outputLevel_ > 2) {
1158 FEI_OSTREAM& os = *logStream();
1159
1160 if (stiff != NULL || load != NULL) {
1161 os << "#numRows"<< FEI_ENDL << numElemRows << FEI_ENDL;
1162 }
1163
1164 if (stiff != NULL) {
1165 os << "#elem-stiff (after copy into dense-row format)" << FEI_ENDL;
1166 for(int i=0; i<numElemRows; i++) {
1167 const double* stiff_i = stiff[i];
1168 for(int j=0; j<numElemRows; j++) {
1169 os << stiff_i[j] << " ";
1170 }
1171 os << FEI_ENDL;
1172 }
1173 }
1174
1175 if (load != NULL) {
1176 os << "#elem-load" << FEI_ENDL;
1177 for(int i=0; i<numElemRows; i++) {
1178 os << load[i] << " ";
1179 }
1180 os<<FEI_ENDL;
1181 }
1182
1183 if (stiff != NULL) {
1184 os << "#elemformat" << FEI_ENDL << elemFormat << FEI_ENDL;
1185 }
1186 }
1187
1188 //now let's obtain the scatter indices for assembling the equations
1189 //into their appropriate places in the global matrix and rhs vectors
1190
1191 int* indPtr = &scatterIndices_[0];
1192 int* blkIndPtr = &blkScatterIndices_[0];
1193 int* blkSizesPtr = blkIndPtr + numBlkElemRows;
1194
1195 bool useBlkEqns = false;
1196 if (interleave == 0) {
1197 //interleave==0 is node-major, so we'll get the block-indices too.
1198 problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
1199 interleave, indPtr,
1200 blkIndPtr, blkSizesPtr);
1201 int sumBlkSizes = 0;
1202 for(int ii=0; ii<numBlkElemRows; ++ii) {
1203 sumBlkSizes += blkSizesPtr[ii];
1204 }
1205 if (sumBlkSizes == numElemRows) {
1206 useBlkEqns = true;
1207 }
1208 }
1209 else {
1210 //interleave!=0 is field-major, and we'll only bother with point-indices.
1211 problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
1212 interleave, indPtr);
1213 }
1214
1215 if (stiff != NULL) {
1216 if (problemStructure_->numSlaveEquations() == 0) {
1217 //I'm not checking the return-value (error-code) on this call, because
1218 //I wasn't even calling it until recently, and I'm not sure if all
1219 //LinearSystemCore implementations even have it implemented.
1220 lsc_->setStiffnessMatrices(elemBlockID, 1, &elemID,
1221 &stiff, scatterIndices_.size(),
1222 &indPtr);
1223 }
1224
1225 int len = scatterIndices_.size();
1226 if (interleave == 0) {
1227 CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
1228 numBlkElemRows, blkIndPtr,
1229 blkSizesPtr, useBlkEqns, ASSEMBLE_SUM ) );
1230 }
1231 else {
1232 CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
1233 numBlkElemRows, blkIndPtr,
1234 blkSizesPtr, false, ASSEMBLE_SUM ) );
1235 }
1236 }
1237
1238 if (load != NULL) {
1239 if (problemStructure_->numSlaveEquations() == 0) {
1240 //I'm not checking the return-value (error-code) on this call, because
1241 //I wasn't even calling it until recently, and I'm not sure if all
1242 //LinearSystemCore implementations even have it implemented.
1243 lsc_->setLoadVectors(elemBlockID, 1, &elemID,
1244 &load, scatterIndices_.size(),
1245 &indPtr);
1246 }
1247
1248 CHK_ERR( assembleRHS(scatterIndices_.size(), indPtr, load, ASSEMBLE_SUM) );
1249 }
1250
1251 return(FEI_SUCCESS);
1252}
1253
1254//------------------------------------------------------------------------------
1255int LinSysCoreFilter::putIntoRHS(int IDType,
1256 int fieldID,
1257 int numIDs,
1258 const GlobalID* IDs,
1259 const double* rhsEntries)
1260{
1261 int fieldSize = problemStructure_->getFieldSize(fieldID);
1262
1263 rowIndices_.resize(fieldSize*numIDs);
1264 int checkNumEqns;
1265
1266 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1267 checkNumEqns,
1268 &rowIndices_[0]));
1269 if (checkNumEqns != numIDs*fieldSize) {
1270 ERReturn(-1);
1271 }
1272
1273 CHK_ERR( exchangeRemoteEquations() );
1274
1275 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
1276
1277 return(0);
1278}
1279
1280//------------------------------------------------------------------------------
1281int LinSysCoreFilter::sumIntoRHS(int IDType,
1282 int fieldID,
1283 int numIDs,
1284 const GlobalID* IDs,
1285 const double* rhsEntries)
1286{
1287 int fieldSize = problemStructure_->getFieldSize(fieldID);
1288
1289 rowIndices_.resize(fieldSize*numIDs);
1290 int checkNumEqns;
1291
1292 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1293 checkNumEqns,
1294 &rowIndices_[0]));
1295 if (checkNumEqns != numIDs*fieldSize) {
1296 ERReturn(-1);
1297 }
1298
1299 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
1300
1301 return(0);
1302}
1303
1304//------------------------------------------------------------------------------
1305int LinSysCoreFilter::implementAllBCs() {
1306//
1307// This function will handle the modifications to the stiffness and load
1308// necessary to enforce nodal boundary conditions.
1309//
1310 debugOutput("# implementAllBCs");
1311
1312 std::vector<int> essEqns;
1313 std::vector<double> essGamma;
1314
1315 EqnBuffer bcEqns;
1316
1317 CHK_ERR( bcManager_->finalizeBCEqns(bcEqns) );
1318
1319 if (resolveConflictRequested_) {
1320 CHK_ERR( resolveConflictingCRs(bcEqns) );
1321 }
1322
1323 CHK_ERR( eqnCommMgr_->gatherSharedBCs(bcEqns) );
1324
1325 //now separate the boundary-condition equations into arrays
1326 fei::FillableMat bcEqns_mat(bcEqns);
1327 fei::impl_utils::separate_BC_eqns(bcEqns_mat, essEqns, essGamma);
1328
1329 std::vector<double> essAlpha(essEqns.size(), 1);
1330
1331 exchangeRemoteBCs(essEqns, essAlpha, essGamma);
1332
1333 if (essEqns.size() > 0) {
1334 CHK_ERR( enforceEssentialBCs(&essEqns[0],
1335 &essAlpha[0],
1336 &essGamma[0], essEqns.size()) );
1337 }
1338
1339 debugOutput("#LinSysCoreFilter leaving implementAllBCs");
1340 return(FEI_SUCCESS);
1341}
1342
1343//------------------------------------------------------------------------------
1344int LinSysCoreFilter::enforceEssentialBCs(const int* eqns,
1345 const double* alpha,
1346 const double* gamma,
1347 int numEqns)
1348{
1349 int* cc_eqns = const_cast<int*>(eqns);
1350 double* cc_alpha = const_cast<double*>(alpha);
1351 double* cc_gamma = const_cast<double*>(gamma);
1352
1353 if (problemStructure_->numSlaveEquations() == 0) {
1354 CHK_ERR( lsc_->enforceEssentialBC(cc_eqns,
1355 cc_alpha, cc_gamma,
1356 numEqns) );
1357 }
1358 else {
1359 std::vector<int> reducedEqns(numEqns);
1360 for(int i=0; i<numEqns; i++) {
1361 problemStructure_->translateToReducedEqn(eqns[i], reducedEqns[i]);
1362 }
1363
1364 CHK_ERR( lsc_->enforceEssentialBC(&reducedEqns[0], cc_alpha, cc_gamma, numEqns) );
1365 }
1366
1367 return(FEI_SUCCESS);
1368}
1369
1370//------------------------------------------------------------------------------
1371int LinSysCoreFilter::enforceRemoteEssBCs(int numEqns, const int* eqns,
1372 const int* const* colIndices, const int* colIndLens,
1373 const double* const* BCcoefs)
1374{
1375 //These eqn-numbers were reduced to the slave-eqn space by
1376 //LinSysCore::exchangeRemoteBCs, which is the only function that calls THIS
1377 //function, so we can simply pass them straight on in to LinearSystemCore.
1378 //
1379
1380 int* cc_eqns = const_cast<int*>(eqns);
1381 int** cc_colIndices = const_cast<int**>(colIndices);
1382 int* cc_colIndLens = const_cast<int*>(colIndLens);
1383 double** cc_BCcoefs = const_cast<double**>(BCcoefs);
1384
1385 CHK_ERR( lsc_->enforceRemoteEssBCs(numEqns, cc_eqns, cc_colIndices,
1386 cc_colIndLens, cc_BCcoefs) );
1387
1388 return(FEI_SUCCESS);
1389}
1390
1391//------------------------------------------------------------------------------
1392int LinSysCoreFilter::resolveConflictingCRs(EqnBuffer& bcEqns)
1393{
1394 int numMultCRs = problemStructure_->getNumMultConstRecords();
1395 if (numMultCRs < 1) {
1396 return(0);
1397 }
1398
1399 std::map<GlobalID,ConstraintType*>::const_iterator
1400 cr_iter = problemStructure_->getMultConstRecords().begin(),
1401 cr_end = problemStructure_->getMultConstRecords().end();
1402
1403 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1404
1405 std::vector<int>& bcEqnNumbers = bcEqns.eqnNumbers();
1406
1407 double coefs[3];
1408 int indices[3];
1409 indices[0] = 0;
1410 indices[1] = 1;
1411 indices[2] = 2;
1412
1413 double fei_eps = 1.e-49;
1414
1415 while(cr_iter != cr_end) {
1416 ConstraintType& multCR = *((*cr_iter).second);
1417
1418 int lenList = multCR.getMasters().size();
1419
1420 std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
1421 GlobalID *CRNodePtr = &CRNode_vec[0];
1422 std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
1423 int* CRFieldPtr = &CRField_vec[0];
1424 std::vector<double>& weights_vec = multCR.getMasterWeights();
1425 double* weights = &weights_vec[0];
1426
1427 int offset = 0;
1428 for(int j=0; j<lenList; ++j) {
1429 int fieldSize = problemStructure_->getFieldSize(CRFieldPtr[j]);
1430 for(int k=0; k<fieldSize; ++k) {
1431 if (std::abs(weights[offset++] + 1.0) < fei_eps) {
1432 const NodeDescriptor* node = NULL;
1433 CHK_ERR( nodeDB.getNodeWithID(CRNodePtr[j], node) );
1434 int eqn = 0;
1435 node->getFieldEqnNumber(CRFieldPtr[j], eqn);
1436 eqn += k;
1437
1438 if (fei::binarySearch(eqn, bcEqnNumbers) >= 0) {
1439 coefs[0] = 1.0;
1440 coefs[1] = 0.0;
1441 coefs[2] = 1.0;
1442 int crEqn = multCR.getEqnNumber();
1443 CHK_ERR( bcEqns.addEqn(crEqn, coefs, indices, 3, false) );
1444 }
1445 }
1446 }
1447 }
1448 ++cr_iter;
1449 }
1450
1451 return(0);
1452}
1453
1454//------------------------------------------------------------------------------
1455int LinSysCoreFilter::exchangeRemoteEquations()
1456{
1457 //
1458 // This function is where processors send local contributions to remote
1459 // equations to the owners of those equations, and receive remote
1460 // contributions to local equations.
1461 //
1462 debugOutput("#LinSysCoreFilter::exchangeRemoteEquations");
1463
1464 if (reducedEqnCounter_ > 0) CHK_ERR( assembleReducedEqns() );
1465
1466 if (reducedRHSCounter_ > 0) CHK_ERR( assembleReducedRHS() );
1467
1468 int len = 4;
1469 std::vector<int> flags(len), globalFlags(len);
1470 flags[0] = newMatrixData_ ? 1 : 0;
1471 flags[1] = newVectorData_ ? 1 : 0;
1472 flags[2] = newConstraintData_ ? 1 : 0;
1473 flags[3] = newBCData_ ? 1 : 0;
1474
1475 CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
1476
1477 newMatrixData_ = globalFlags[0] > 0 ? true : false;
1478 newVectorData_ = globalFlags[1] > 0 ? true : false;
1479 newConstraintData_ = globalFlags[2] > 0 ? true : false;
1480 newBCData_ = globalFlags[3] > 0 ? true : false;
1481
1482 if (newMatrixData_ || newVectorData_ || newConstraintData_) {
1483
1484 CHK_ERR( eqnCommMgr_->exchangeEqns(logStream()) );
1485
1486 needToCallMatrixLoadComplete_ = true;
1487
1488 //so now the remote contributions should be available, let's get them out
1489 //of the eqn comm mgr and put them into our local matrix structure.
1490
1491 debugOutput("# putting remote contributions into linear system...");
1492
1493 CHK_ERR( unpackRemoteContributions(*eqnCommMgr_, ASSEMBLE_SUM) );
1494
1495 eqnCommMgr_->resetCoefs();
1496
1497 newMatrixData_ = false;
1498 newVectorData_ = false;
1499 newConstraintData_ = false;
1500 }
1501
1502 firstRemEqnExchange_ = false;
1503
1504 debugOutput("#LinSysCoreFilter leaving exchangeRemoteEquations");
1505
1506 return(FEI_SUCCESS);
1507}
1508
1509//------------------------------------------------------------------------------
1510int LinSysCoreFilter::unpackRemoteContributions(EqnCommMgr& eqnCommMgr,
1511 int assemblyMode)
1512{
1513 int numRecvEqns = eqnCommMgr.getNumLocalEqns();
1514 std::vector<int>& recvEqnNumbers = eqnCommMgr.localEqnNumbers();
1515 std::vector<fei::CSVec*>& recvEqns = eqnCommMgr.localEqns();
1516 std::vector<std::vector<double>*>& recvRHSs = *(eqnCommMgr.localRHSsPtr());
1517
1518 bool newCoefs = eqnCommMgr.newCoefData();
1519 bool newRHSs = eqnCommMgr.newRHSData();
1520
1521 int i;
1522 double** coefs = new double*[numRecvEqns];
1523
1524 for(i=0; i<numRecvEqns; i++) {
1525 coefs[i] = &(recvEqns[i]->coefs()[0]);
1526 }
1527
1528 for(i=0; i<numRecvEqns; i++) {
1529
1530 int eqn = recvEqnNumbers[i];
1531 if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
1532 fei::console_out() << "LinSysCoreFilter::unpackRemoteContributions: ERROR, recvEqn "
1533 << eqn << " out of range. (localStartRow_: " << reducedStartRow_
1534 << ", localEndRow_: " << reducedEndRow_ << ", localRank_: "
1535 << localRank_ << ")" << FEI_ENDL;
1536 MPI_Abort(comm_, -1);
1537 }
1538
1539 for(size_t ii=0; ii<recvEqns[i]->size(); ii++) {
1540 if (coefs[i][ii] > 1.e+200) {
1541 fei::console_out() << localRank_ << ": LinSysCoreFilter::unpackRemoteContributions: "
1542 << "WARNING, coefs["<<i<<"]["<<ii<<"]: " << coefs[i][ii]
1543 << FEI_ENDL;
1544 MPI_Abort(comm_, -1);
1545 }
1546 }
1547
1548 if (recvEqns[i]->size() > 0 && newCoefs) {
1549 //contribute this equation to the matrix,
1550 CHK_ERR( giveToLocalReducedMatrix(1, &(recvEqnNumbers[i]),
1551 recvEqns[i]->size(),
1552 &(recvEqns[i]->indices()[0]),
1553 &(coefs[i]), assemblyMode ) );
1554 }
1555
1556 //and now the RHS contributions.
1557 if (newRHSs) {
1558 for(int j=0; j<numRHSs_; j++) {
1559 CHK_ERR( giveToLocalReducedRHS(1, &( (*(recvRHSs[i]))[j] ),
1560 &eqn, assemblyMode) );
1561 }
1562 }
1563 }
1564
1565 delete [] coefs;
1566
1567 return(0);
1568}
1569
1570//------------------------------------------------------------------------------
1571int LinSysCoreFilter::exchangeRemoteBCs(std::vector<int>& essEqns,
1572 std::vector<double>& essAlpha,
1573 std::vector<double>& essGamma)
1574{
1575 //we need to make sure that the right thing happens for essential
1576 //boundary conditions that get applied to nodes on elements that touch
1577 //a processor boundary. (Note that for this case, the BC node itself doesn't
1578 //touch the processor boundary.) For an essential boundary condition, the row
1579 //and column of the corresponding equation must be diagonalized. If there is
1580 //a processor boundary on any side of the element that contains the node,
1581 //then there are column contributions to the matrix on the other processor.
1582 //That other processor must be notified and told to make the adjustments
1583 //necessary to enforce the boundary condition.
1584
1585 std::vector<int>* eqns = &essEqns;
1586
1587 if (problemStructure_->numSlaveEquations() > 0) {
1588 int numEqns = essEqns.size();
1589 eqns = new std::vector<int>(numEqns);
1590 int* eqnsPtr = &(*eqns)[0];
1591
1592 for(int ii=0; ii<numEqns; ++ii) {
1593 problemStructure_->translateToReducedEqn(essEqns[ii], eqnsPtr[ii]);
1594 }
1595 }
1596
1597 FEI_OSTREAM* dbgOut = NULL;
1598 if (Filter::logStream() != NULL) {
1599 dbgOut = logStream();
1600 }
1601
1602 eqnCommMgr_->exchangeRemEssBCs(&(*eqns)[0], eqns->size(),
1603 &essAlpha[0], &essGamma[0],
1604 comm_, dbgOut);
1605
1606 int numRemoteEssBCEqns = eqnCommMgr_->getNumRemEssBCEqns();
1607 if (numRemoteEssBCEqns > 0) {
1608 std::vector<int>& remEssBCEqnNumbers = eqnCommMgr_->remEssBCEqnNumbersPtr();
1609 fei::CSVec** remEssBCEqns = &(eqnCommMgr_->remEssBCEqns()[0]);
1610 std::vector<int> remEssBCEqnLengths(remEssBCEqnNumbers.size());
1611
1612 int** indices = new int*[numRemoteEssBCEqns];
1613 double** coefs = new double*[numRemoteEssBCEqns];
1614
1615 for(int i=0; i<numRemoteEssBCEqns; i++) {
1616 coefs[i] = &(remEssBCEqns[i]->coefs()[0]);
1617 indices[i] = &(remEssBCEqns[i]->indices()[0]);
1618 remEssBCEqnLengths[i] = remEssBCEqns[i]->size();
1619 }
1620
1621 CHK_ERR( enforceRemoteEssBCs(numRemoteEssBCEqns,
1622 &remEssBCEqnNumbers[0], indices,
1623 &remEssBCEqnLengths[0], coefs));
1624
1625 delete [] indices;
1626 delete [] coefs;
1627 }
1628
1629 if (problemStructure_->numSlaveEquations() > 0) {
1630 delete eqns;
1631 }
1632
1633 debugOutput("#LinSysCoreFilter exchangeRemoteBCs");
1634
1635 return(FEI_SUCCESS);
1636}
1637
1638//------------------------------------------------------------------------------
1639int LinSysCoreFilter::loadCRMult(int CRID,
1640 int numCRNodes,
1641 const GlobalID* CRNodes,
1642 const int* CRFields,
1643 const double* CRWeights,
1644 double CRValue)
1645{
1646//
1647// Load Lagrange multiplier constraint relation data
1648//
1649// Question: do we really need to pass CRNodes again? Here, I'm going
1650// to ignore it for now (i.e., not store it, but just check it),
1651// as it got passed during the initialization phase, so all we'll
1652// do here is check for errors...
1653//
1654 if (Filter::logStream() != NULL) {
1655 FEI_OSTREAM& os = *logStream();
1656 os<<"FEI: loadCRMult"<<FEI_ENDL;
1657 os<<"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
1658 os<<"#CRNodes:"<<FEI_ENDL;
1659 int i;
1660 for(i=0; i<numCRNodes; ++i) {
1661 GlobalID nodeID = CRNodes[i];
1662 os << static_cast<int>(nodeID) << " ";
1663 }
1664 os << FEI_ENDL << "#fields:"<<FEI_ENDL;
1665 for(i=0; i<numCRNodes; ++i) os << CRFields[i] << " ";
1666 os << FEI_ENDL << "#field-sizes:"<<FEI_ENDL;
1667 for(i=0; i<numCRNodes; ++i) {
1668 int size = problemStructure_->getFieldSize(CRFields[i]);
1669 os << size << " ";
1670 }
1671 os << FEI_ENDL<<"#weights:"<<FEI_ENDL;
1672 int offset = 0;
1673 for(i=0; i<numCRNodes; ++i) {
1674 int size = problemStructure_->getFieldSize(CRFields[i]);
1675 for(int j=0; j<size; ++j) {
1676 os << CRWeights[offset++] << " ";
1677 }
1678 }
1679 os << FEI_ENDL<<"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
1680 }
1681
1682 ConstraintType* multCR = NULL;
1683 CHK_ERR( problemStructure_->getMultConstRecord(CRID, multCR) );
1684
1685 int i;
1686
1687 int lenList = multCR->getMasters().size();
1688 if (lenList < 1) {
1689 fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
1690 <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
1691 ERReturn(-1);
1692 }
1693
1694 // recall the data stored earlier and ensure that the passed data (here, the
1695 // node list) agrees with the initialization data
1696
1697 std::vector<GlobalID>& CRNode_vec = multCR->getMasters();
1698 GlobalID *CRNodePtr = &CRNode_vec[0];
1699
1700 for(i=0; i<lenList; i++) {
1701 if (CRNodePtr[i] != CRNodes[i]) {
1702 fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
1703 << " in initCRMult than it has in loadCRMult."<<FEI_ENDL;
1704 ERReturn(-1);
1705 }
1706 }
1707
1708 std::vector<int>& CRField_vec = multCR->getMasterFieldIDs();
1709 int *CRFieldPtr = &CRField_vec[0];
1710
1711 for (i = 0; i < lenList; i++) {
1712 if (CRFieldPtr[i] != CRFields[i]) {
1713 fei::console_out() <<"ERROR in FEI, constraint with CRID="<<CRID<<" had different field-list"
1714 <<" in initCRMult than it has in loadCRMult."<<FEI_ENDL;
1715 ERReturn(-1);
1716 }
1717 }
1718
1719 newConstraintData_ = true;
1720
1721 if ((int)iworkSpace_.size() < lenList) {
1722 iworkSpace_.resize(lenList);
1723 }
1724
1725 int* fieldSizes = &iworkSpace_[0];
1726
1727 for (i = 0; i < lenList; i++) {
1728 int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
1729 assert(numSolnParams >= 0);
1730 fieldSizes[i] = numSolnParams;
1731 }
1732
1733 std::vector<double>& CRWeightArray = multCR->getMasterWeights();
1734
1735 int offset = 0;
1736
1737 try {
1738
1739 for(i = 0; i < lenList; i++) {
1740 for(int j = 0; j < fieldSizes[i]; j++) {
1741 CRWeightArray.push_back(CRWeights[offset++]);
1742 }
1743 }
1744
1745 }
1746 catch(std::runtime_error& exc) {
1747 fei::console_out() << exc.what() << FEI_ENDL;
1748 ERReturn(-1);
1749 }
1750
1751 multCR->setRHSValue(CRValue);
1752 double* CRWeightsPtr = &CRWeightArray[0];
1753
1754// next, perform assembly of the various terms into the system arrays
1755// (this is a good candidate for a separate function...)
1756
1757 int irow = multCR->getEqnNumber();
1758 double zero = 0.0;
1759 double* zeroPtr = &zero;
1760 CHK_ERR( giveToMatrix(1, &irow, 1, &irow, &zeroPtr, ASSEMBLE_PUT) );
1761
1762 CHK_ERR( giveToRHS(1, &(CRValue), &irow, ASSEMBLE_PUT));
1763
1764 offset = 0;
1765 for(int j = 0; j < lenList; j++) {
1766 int myFieldID = CRFields[j];
1767
1768 const NodeDescriptor& node = Filter::findNodeDescriptor(CRNodePtr[j]);
1769
1770 //first, store the column coeficients for equation irow, the
1771 //constraint's equation.
1772 storeNodalColumnCoefs(irow, node, myFieldID, fieldSizes[j],
1773 &(CRWeightsPtr[offset]));
1774
1775
1776 //next, store store the transpose of the above. i.e., column irow,
1777 //in equations associated with 'node' at 'myFieldID'.
1778
1779 if (node.getOwnerProc() == localRank_) {
1780
1781 storeNodalRowCoefs(node, myFieldID, fieldSizes[j],
1782 &(CRWeightsPtr[offset]), irow);
1783 }
1784 else {
1785
1786 storeNodalSendEqn(node, myFieldID, irow, &(CRWeightsPtr[offset]));
1787 }
1788
1789 offset += fieldSizes[j];
1790 }
1791
1792 return(FEI_SUCCESS);
1793}
1794
1795//------------------------------------------------------------------------------
1796int LinSysCoreFilter::loadCRPen(int CRID,
1797 int numCRNodes,
1798 const GlobalID* CRNodes,
1799 const int* CRFields,
1800 const double* CRWeights,
1801 double CRValue,
1802 double penValue)
1803{
1804 //
1805 // Load penalty constraint relation data
1806 //
1807
1808 debugOutput("FEI: loadCRPen");
1809
1810 ConstraintType* penCR = NULL;
1811 CHK_ERR( problemStructure_->getPenConstRecord(CRID, penCR) );
1812
1813 int i;
1814 int lenList = penCR->getMasters().size();
1815 if (lenList < 1) {
1816 fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
1817 <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
1818 ERReturn(-1);
1819 }
1820
1821 // recall the data stored earlier and ensure that the passed data (here,
1822 // the node list) agrees with the initialization data
1823
1824 std::vector<GlobalID>& CRNode_vec = penCR->getMasters();
1825 GlobalID* CRNodePtr = &CRNode_vec[0];
1826
1827 for(int j = 0; j < lenList; j++) {
1828 if (CRNodePtr[j] != CRNodes[j]) {
1829 fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
1830 << " in initCRPen than it has in loadCRPen."<<FEI_ENDL;
1831 ERReturn(-1);
1832 }
1833 }
1834
1835 newConstraintData_ = true;
1836
1837 // store the weights and rhs-value in the constraint records.
1838
1839 if ((int)iworkSpace_.size() < lenList) {
1840 iworkSpace_.resize(lenList);
1841 }
1842
1843 int* fieldSizes = &iworkSpace_[0];
1844
1845 for (i = 0; i < lenList; i++) {
1846 int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
1847 assert(numSolnParams >= 0);
1848 fieldSizes[i] = numSolnParams;
1849 }
1850
1851 std::vector<double>& CRWeightArray = penCR->getMasterWeights();
1852
1853 try {
1854
1855 int offset = 0;
1856 for (i = 0; i < lenList; i++) {
1857 for (int j = 0; j < fieldSizes[i]; j++) {
1858 CRWeightArray.push_back(CRWeights[offset++]);
1859 }
1860 }
1861
1862 }
1863 catch(std::runtime_error& exc) {
1864 fei::console_out() << exc.what() << FEI_ENDL;
1865 ERReturn(-1);
1866 }
1867
1868 penCR->setRHSValue(CRValue);
1869
1870 double* CRWeightPtr = &CRWeightArray[0];
1871
1872 int ioffset = 0, joffset = 0;
1873 for(i = 0; i < lenList; i++) {
1874 GlobalID iNodeID = CRNodePtr[i];
1875 int iField = CRFields[i];
1876
1877 const NodeDescriptor& iNode = Filter::findNodeDescriptor(iNodeID);
1878 double* iweights = &(CRWeightPtr[ioffset]);
1879 ioffset += fieldSizes[i];
1880
1881 joffset = 0;
1882 for (int j = 0; j < lenList; j++) {
1883 GlobalID jNodeID = CRNodePtr[j];
1884 int jField = CRFields[j];
1885
1886 const NodeDescriptor& jNode = Filter::findNodeDescriptor(jNodeID);
1887 double* jweights = &(CRWeightPtr[joffset]);
1888 joffset += fieldSizes[j];
1889
1890 double rhsValue = CRValue;
1891 if (j < lenList-1) {
1892 rhsValue = 0.0;
1893 }
1894
1895 if (iNode.getOwnerProc() == localRank_) {
1896
1897 storePenNodeData(iNode, iField, fieldSizes[i], iweights,
1898 jNode, jField, fieldSizes[j], jweights,
1899 penValue, rhsValue);
1900 }
1901 else {
1902 storePenNodeSendData(iNode, iField, fieldSizes[i], iweights,
1903 jNode, jField, fieldSizes[j], jweights,
1904 penValue, rhsValue);
1905 }
1906 }
1907 }
1908
1909 return(FEI_SUCCESS);
1910}
1911
1912//------------------------------------------------------------------------------
1913int LinSysCoreFilter::parameters(int numParams, const char *const* paramStrings) {
1914//
1915// this function takes parameters for setting internal things like solver
1916// and preconditioner choice, etc.
1917//
1918 if (numParams == 0 || paramStrings == NULL) {
1919 debugOutput("#LinSysCoreFilter::parameters--- no parameters.");
1920 }
1921 else {
1922 const char* param1 = snl_fei::getParamValue("AZ_matrix_type",
1923 numParams,
1924 paramStrings);
1925
1926 if (param1 != NULL) {
1927 if (!strcmp(param1, "AZ_VBR_MATRIX") ||
1928 !strcmp(param1, "blockMatrix")) {
1929 blockMatrix_ = true;
1930 }
1931 }
1932 else {
1933 param1 = snl_fei::getParamValue("matrixType",
1934 numParams, paramStrings);
1935 if (param1 != NULL) {
1936 if (!strcmp(param1, "AZ_VBR_MATRIX") ||
1937 !strcmp(param1, "blockMatrix")) {
1938 blockMatrix_ = true;
1939 }
1940 }
1941 }
1942
1943 param1 = snl_fei::getParamValue("outputLevel",
1944 numParams,paramStrings);
1945 if ( param1 != NULL){
1946 std::string str(param1);
1947 FEI_ISTRINGSTREAM isstr(str);
1948 isstr >> outputLevel_;
1949 }
1950
1951 param1 = snl_fei::getParam("resolveConflict",numParams,paramStrings);
1952 if ( param1 != NULL){
1953 resolveConflictRequested_ = true;
1954 }
1955
1956 param1 = snl_fei::getParamValue("internalFei", numParams,paramStrings);
1957 if ( param1 != NULL ){
1958 std::string str(param1);
1959 FEI_ISTRINGSTREAM isstr(str);
1960 isstr >> internalFei_;
1961 }
1962
1963 if (Filter::logStream() != NULL) {
1964
1965 (*logStream())<<"#LinSysCoreFilter::parameters"<<FEI_ENDL
1966 <<"# --- numParams: "<< numParams<<FEI_ENDL;
1967 for(int i=0; i<numParams; i++){
1968 (*logStream())<<"#------ paramStrings["<<i<<"]: "
1969 <<paramStrings[i];
1970 if (paramStrings[i][strlen(paramStrings[i])-1] != '\n') {
1971 (*logStream())<<FEI_ENDL;
1972 }
1973 }
1974 }
1975 }
1976
1977 CHK_ERR( Filter::parameters(numParams, paramStrings) );
1978
1979 debugOutput("#LinSysCoreFilter leaving parameters function");
1980
1981 return(FEI_SUCCESS);
1982}
1983
1984//------------------------------------------------------------------------------
1985int LinSysCoreFilter::loadComplete()
1986{
1987 int len = 4;
1988 std::vector<int> flags(len), globalFlags(len);
1989 flags[0] = newMatrixData_ ? 1 : 0;
1990 flags[1] = newVectorData_ ? 1 : 0;
1991 flags[2] = newConstraintData_ ? 1 : 0;
1992 flags[3] = newBCData_ ? 1 : 0;
1993
1994 CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
1995
1996 newMatrixData_ = globalFlags[0] > 0 ? true : false;
1997 newVectorData_ = globalFlags[1] > 0 ? true : false;
1998 newConstraintData_ = globalFlags[2] > 0 ? true : false;
1999 newBCData_ = globalFlags[3] > 0 ? true : false;
2000
2001 bool called_exchange = false;
2002 if (newMatrixData_ || newVectorData_ || newConstraintData_) {
2003 CHK_ERR( exchangeRemoteEquations() );
2004 called_exchange = true;
2005 }
2006
2007 bool called_implbcs = false;
2008 if (newBCData_) {
2009 CHK_ERR( implementAllBCs() );
2010 called_implbcs = true;
2011 }
2012
2013 if (called_exchange || called_implbcs ||needToCallMatrixLoadComplete_) {
2014 debugOutput("#LinSysCoreFilter calling LinSysCore matrixLoadComplete");
2015
2016 CHK_ERR( lsc_->matrixLoadComplete() );
2017 needToCallMatrixLoadComplete_ = false;
2018 }
2019
2020 newMatrixData_ = false;
2021 newVectorData_ = false;
2022 newConstraintData_ = false;
2023 newBCData_ = false;
2024
2025 return(0);
2026}
2027
2028//------------------------------------------------------------------------------
2029int LinSysCoreFilter::residualNorm(int whichNorm, int numFields,
2030 int* fieldIDs, double* norms, double& residTime)
2031{
2032 //
2033 //This function can do 3 kinds of norms: infinity-norm (denoted
2034 //by whichNorm==0), 1-norm and 2-norm.
2035 //
2036 debugOutput("FEI: residualNorm");
2037
2038 if (whichNorm < 0 || whichNorm > 2) return(-1);
2039
2040 CHK_ERR( loadComplete() );
2041
2042 std::vector<double> residValues(numReducedRows_, 0.0);
2043
2044 double start = fei::utils::cpu_time();
2045
2046 CHK_ERR( formResidual(&(residValues[0]), numReducedRows_) );
2047
2048 residTime = fei::utils::cpu_time() - start;
2049
2050 CHK_ERR( Filter::calculateResidualNorms(whichNorm, numFields,
2051 fieldIDs, norms, residValues) );
2052
2053 return(FEI_SUCCESS);
2054}
2055
2056//------------------------------------------------------------------------------
2057int LinSysCoreFilter::formResidual(double* residValues, int numLocalEqns)
2058{
2059 CHK_ERR( lsc_->formResidual(residValues, numLocalEqns))
2060
2061 return(FEI_SUCCESS);
2062}
2063
2064//------------------------------------------------------------------------------
2065int LinSysCoreFilter::solve(int& status, double& sTime) {
2066
2067 debugOutput("FEI: solve");
2068
2069 CHK_ERR( loadComplete() );
2070
2071 debugOutput("#LinSysCoreFilter in solve, calling launchSolver...");
2072
2073 double start = fei::utils::cpu_time();
2074
2075 CHK_ERR( lsc_->launchSolver(status, iterations_) );
2076
2077 sTime = fei::utils::cpu_time() - start;
2078
2079 debugOutput("#LinSysCoreFilter... back from solver");
2080
2081 //now unpack the locally-owned shared entries of the solution vector into
2082 //the eqn-comm-mgr data structures.
2083 CHK_ERR( unpackSolution() );
2084
2085 debugOutput("#LinSysCoreFilter leaving solve");
2086
2087 if (status != 0) return(1);
2088 else return(FEI_SUCCESS);
2089}
2090
2091//------------------------------------------------------------------------------
2092int LinSysCoreFilter::setNumRHSVectors(int numRHSs, int* rhsIDs){
2093
2094 if (numRHSs < 0) {
2095 fei::console_out() << "LinSysCoreFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
2096 ERReturn(-1);
2097 }
2098
2099 numRHSs_ = numRHSs;
2100
2101 rhsIDs_.resize(numRHSs_);
2102 for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
2103
2104 //(we need to set the number of RHSs in the eqn comm manager)
2105 eqnCommMgr_->setNumRHSs(numRHSs_);
2106
2107 return(FEI_SUCCESS);
2108}
2109
2110//------------------------------------------------------------------------------
2111int LinSysCoreFilter::setCurrentRHS(int rhsID)
2112{
2113 std::vector<int>::iterator iter =
2114 std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
2115
2116 if (iter == rhsIDs_.end()) ERReturn(-1)
2117
2118 int index = iter - rhsIDs_.begin();
2119 currentRHS_ = index;
2120
2121 lsc_->setRHSID(rhsID);
2122
2123 return(FEI_SUCCESS);
2124}
2125
2126//------------------------------------------------------------------------------
2127int LinSysCoreFilter::giveToMatrix_symm_noSlaves(int numPtRows,
2128 const int* ptRowNumbers,
2129 const double* const* coefs,
2130 int mode)
2131{
2132 for(int i=0; i<numPtRows; i++) {
2133 int row = ptRowNumbers[i];
2134 const double* valptr = coefs[i];
2135 if (row < localStartRow_ || row > localEndRow_) {
2136 eqnCommMgr_->addRemoteEqn(row, valptr, ptRowNumbers, numPtRows);
2137 continue;
2138 }
2139
2140 if (mode == ASSEMBLE_SUM) {
2141 if (Filter::logStream() != NULL && 0) {
2142 FEI_OSTREAM& os = *logStream();
2143 os << "# calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
2144 << ", columns: ";
2145 for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
2146 os << FEI_ENDL;
2147 }
2148
2149 CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRowNumbers[i]),
2150 numPtRows, ptRowNumbers,
2151 &valptr) );
2152 }
2153 else {
2154 CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
2155 numPtRows, ptRowNumbers,
2156 &valptr) );
2157 }
2158 }
2159
2160 return(0);
2161}
2162
2163//------------------------------------------------------------------------------
2164int LinSysCoreFilter::giveToBlkMatrix_symm_noSlaves(int numPtRows,
2165 const int* ptRowNumbers,
2166 int numBlkRows,
2167 const int* blkRowNumbers,
2168 const int* blkRowSizes,
2169 const double* const* coefs,
2170 int mode)
2171{
2172 int i;
2173 if ((int)dworkSpace2_.size() < numPtRows) {
2174 dworkSpace2_.resize(numPtRows);
2175 }
2176 const double* * valptr = &dworkSpace2_[0];
2177 for(i=0; i<numPtRows; i++) {
2178 int row = ptRowNumbers[i];
2179 valptr[i] = coefs[i];
2180 if (row < localStartRow_ || row > localEndRow_) {
2181 eqnCommMgr_->addRemoteEqn(row, valptr[i], ptRowNumbers, numPtRows);
2182 continue;
2183 }
2184
2185 if (mode == ASSEMBLE_PUT) {
2186 CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
2187 numPtRows, ptRowNumbers,
2188 &(valptr[i])) );
2189 }
2190 }
2191
2192 int offset = 0;
2193 for(i=0; i<numBlkRows; i++) {
2194 int row = ptRowNumbers[offset];
2195 if (row < localStartRow_ || row > localEndRow_) {
2196 offset += blkRowSizes[i];
2197 continue;
2198 }
2199
2200 if (mode == ASSEMBLE_SUM) {
2201 if (Filter::logStream() != NULL && 0) {
2202 FEI_OSTREAM& os = *logStream();
2203 os << "# calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
2204 << ", columns: ";
2205 for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
2206 os << FEI_ENDL;
2207 }
2208
2209 CHK_ERR(lsc_->sumIntoSystemMatrix(blkRowSizes[i],&(ptRowNumbers[offset]),
2210 numPtRows, ptRowNumbers,
2211 1, &(blkRowNumbers[i]),
2212 numBlkRows, blkRowNumbers,
2213 &(valptr[offset])) );
2214 }
2215
2216 offset += blkRowSizes[i];
2217 }
2218
2219 return(0);
2220}
2221
2222//------------------------------------------------------------------------------
2223int LinSysCoreFilter::giveToMatrix(int numPtRows, const int* ptRows,
2224 int numPtCols, const int* ptCols,
2225 const double* const* values, int mode)
2226{
2227 try {
2228
2229 if (problemStructure_->numSlaveEquations() == 0) {
2230 for(int i=0; i<numPtRows; i++) {
2231 if (ptRows[i] < localStartRow_ || ptRows[i] > localEndRow_) {
2232 eqnCommMgr_->addRemoteEqn(ptRows[i], values[i], ptCols, numPtCols);
2233 continue;
2234 }
2235
2236 if (mode == ASSEMBLE_SUM) {
2237 if (Filter::logStream() != NULL && 0) {
2238 FEI_OSTREAM& os = *logStream();
2239 os << "# calling sumIntoSystemMatrix, row: " << ptRows[i]
2240 << ", columns: ";
2241 for(int j=0; j<numPtCols; ++j) os << ptCols[j] << " ";
2242 os << FEI_ENDL;
2243 }
2244
2245 CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]),
2246 numPtCols, ptCols,
2247 &(values[i])) );
2248 }
2249 else {
2250 CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]),
2251 numPtCols, ptCols,
2252 &(values[i])) );
2253 }
2254 }
2255 }
2256 else {
2257 iworkSpace_.resize(numPtCols);
2258 iworkSpace2_.resize(numPtCols);
2259 int* iworkPtr = &iworkSpace_[0];
2260 int* iworkPtr2= &iworkSpace2_[0];
2261 int offset = 0;
2262 for(int ii=0; ii<numPtCols; ii++) {
2263 int reducedEqn = -1;
2264 bool isSlave = problemStructure_->translateToReducedEqn(ptCols[ii],
2265 reducedEqn);
2266 if (isSlave) {
2267 reducedEqn = -1;
2268 iworkPtr[ii] = reducedEqn;
2269 }
2270 else {
2271 iworkPtr[ii] = reducedEqn;
2272 iworkPtr2[offset++] = reducedEqn;
2273 }
2274 }
2275 iworkSpace2_.resize(offset);
2276
2277 for(int i=0; i<numPtRows; i++) {
2278 int row = ptRows[i];
2279
2280 int reducedRow;
2281 bool isSlave = problemStructure_->translateToReducedEqn(row, reducedRow);
2282 if (isSlave) continue;
2283
2284 if (reducedStartRow_ > reducedRow || reducedRow > reducedEndRow_) {
2285
2286 dworkSpace_.resize(0);
2287 for(int j=0; j<numPtCols; j++) {
2288 if (iworkSpace_[j]>=0) {
2289 if (Filter::logStream() != NULL) {
2290 (*logStream())<<"# giveToMatrix remote("<<reducedRow<<","
2291 <<iworkSpace_[j]<<","<<values[i][j]<<")"<<FEI_ENDL;
2292 }
2293
2294 dworkSpace_.push_back(values[i][j]);
2295 }
2296 }
2297
2298 if (mode == ASSEMBLE_SUM) {
2299 if (Filter::logStream() != NULL) {
2300 (*logStream())<<"sum"<<FEI_ENDL;
2301 }
2302
2303 eqnCommMgr_->addRemoteEqn(reducedRow,
2304 &dworkSpace_[0],
2305 &iworkSpace2_[0],
2306 iworkSpace2_.size());
2307 }
2308 else {
2309 if (Filter::logStream() != NULL) {
2310 (*logStream())<<"put"<<FEI_ENDL;
2311 }
2312
2313 eqnCommMgr_put_->addRemoteEqn(reducedRow,
2314 &dworkSpace_[0],
2315 &iworkSpace2_[0],
2316 iworkSpace2_.size());
2317 }
2318
2319 continue;
2320 }
2321
2322 for(int j=0; j<numPtCols; j++) {
2323
2324 int reducedCol = iworkPtr[j];
2325 if (reducedCol<0) continue;
2326
2327 double* tmpCoef = const_cast<double*>(&(values[i][j]));
2328
2329 if (Filter::logStream() != NULL) {
2330 (*logStream())<< "# giveToMatrix local("<<reducedRow
2331 <<","<<reducedCol<<","<<*tmpCoef<<")"<<FEI_ENDL;
2332 }
2333
2334 if (mode == ASSEMBLE_SUM) {
2335 if (Filter::logStream() != NULL && 0) {
2336 FEI_OSTREAM& os = *logStream();
2337 os << "# calling sumIntoSystemMatrix, row: " << reducedRow
2338 << ", columns: " << reducedCol << FEI_ENDL;
2339 }
2340
2341 CHK_ERR( lsc_->sumIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
2342 &tmpCoef ) );
2343 }
2344 else {
2345 CHK_ERR( lsc_->putIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
2346 &tmpCoef ) );
2347 }
2348 }
2349 }
2350 }
2351
2352 }
2353 catch(std::runtime_error& exc) {
2354 fei::console_out() << exc.what() << FEI_ENDL;
2355 ERReturn(-1);
2356 }
2357
2358 return(FEI_SUCCESS);
2359}
2360
2361//------------------------------------------------------------------------------
2362int LinSysCoreFilter::giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
2363 int numPtCols, const int* ptCols,
2364 const double* const* values, int mode)
2365{
2366 bool specialCase = (!firstRemEqnExchange_ && newConstraintData_
2367 && !newMatrixData_) ? true : false;
2368
2369 double fei_min = std::numeric_limits<double>::min();
2370
2371 for(int i=0; i<numPtRows; i++) {
2372
2373 if (mode == ASSEMBLE_SUM) {
2374 const double* values_i = values[i];
2375
2376 for(int j=0; j<numPtCols; ++j) {
2377 if (specialCase && std::abs(values_i[j]) < fei_min) continue;
2378
2379 const double* valPtr = &(values_i[j]);
2380 CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]), 1, &(ptCols[j]),
2381 &valPtr) );
2382 }
2383 }
2384 else {
2385 CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]), numPtCols, ptCols,
2386 &(values[i])) );
2387 }
2388 }
2389
2390 return(FEI_SUCCESS);
2391}
2392
2393//------------------------------------------------------------------------------
2394int LinSysCoreFilter::sumIntoMatrix(fei::CSRMat& mat)
2395{
2396 const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
2397 const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
2398 const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
2399 const std::vector<double>& pckCoefs = mat.getPackedCoefs();
2400
2401 for(size_t i=0; i<rowNumbers.size(); ++i) {
2402 int row = rowNumbers[i];
2403 int offset = rowOffsets[i];
2404 int rowlen = rowOffsets[i+1]-offset;
2405 const int* indices = &pckColInds[offset];
2406 const double* coefs = &pckCoefs[offset];
2407
2408 if (giveToMatrix(1, &row, rowlen, indices, &coefs, ASSEMBLE_SUM) != 0) {
2409 throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
2410 }
2411 }
2412
2413 return(FEI_SUCCESS);
2414}
2415
2416//------------------------------------------------------------------------------
2417int LinSysCoreFilter::getFromMatrix(int numPtRows, const int* ptRows,
2418 const int* rowColOffsets, const int* ptCols,
2419 int numColsPerRow, double** values)
2420{
2421 //This function may be attempting to retrieve matrix rows that are not
2422 //locally owned. If those rows correspond to finite-element nodes that we
2423 //share, AND if the owning processor is also making this function call, then
2424 //we can communicate with that processor and obtain those matrix rows.
2425 //
2426
2427 ProcEqns remoteProcEqns;
2428
2429 //Let's populate this ProcEqns object with the remote equations and the procs
2430 //that we need to receive the remote equations from.
2431 for(int re=0; re<numPtRows; re++) {
2432 int eqn = ptRows[re];
2433 int owner = problemStructure_->getOwnerProcForEqn(eqn);
2434 if (owner == localRank_) continue;
2435
2436 remoteProcEqns.addEqn(eqn, owner);
2437 }
2438
2439 //so now we know which of the requested equations are remotely owned, and we
2440 //know which processors own them.
2441 //Next we're going to need to know which locally-owned equations are needed
2442 //by other processors.
2443 ProcEqns localProcEqns;
2444 CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) )
2445
2446 //ok, now we know which local equations we'll need to send, so let's extract
2447 //those from the matrix
2448 EqnBuffer localEqns;
2449 CHK_ERR( getEqnsFromMatrix(localProcEqns, localEqns) )
2450
2451 //now we can set the lengths in localProcEqns.
2452 std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
2453 fei::CSVec** localEqnsPtr = (localEqns.eqns().size() ? &(localEqns.eqns()[0]) : 0);
2454 std::vector<int> eqnLengths(eqnNumbers.size());
2455 for(size_t i=0; i<eqnNumbers.size(); ++i) {
2456 eqnLengths[i] = localEqnsPtr[i]->size();
2457 }
2458
2459 localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
2460 eqnNumbers.size());
2461
2462 //now mirror those lengths in the remoteProcEqns objects to get ready for the
2463 //all-to-all exchange of equation data.
2464 CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
2465
2466 EqnBuffer remoteEqns;
2467 //we're now ready to do the exchange.
2468 CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
2469 &remoteProcEqns, &remoteEqns, false));
2470
2471 std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
2472 fei::CSVec** remEqnsPtr = (remoteEqns.eqns().size() ? &(remoteEqns.eqns()[0]) : 0);
2473 std::vector<fei::CSVec*>& remEqns = remoteEqns.eqns();
2474
2475 //now we're ready to fill the values array with the remote coefficients.
2476 for(int i=0; i<numPtRows; i++) {
2477 int row = ptRows[i];
2478
2479 int eqnIndex = fei::binarySearch(row, remEqnNumbers);
2480
2481 //if eqnIndex < 0, this is a local equation, so skip to the next loop iter.
2482 if (eqnIndex < 0) continue;
2483
2484 //the equation is remote, so stride across it copying out the coefs.
2485 //if ptCols is NULL, then we're going to copy all coefficients (the whole
2486 //row) into 'values'.
2487 if (ptCols == NULL) {
2488 for(size_t j=0; j<remEqnsPtr[eqnIndex]->size(); j++) {
2489 values[i][j] = remEqns[eqnIndex]->coefs()[j];
2490 }
2491 continue;
2492 }
2493
2494 for(int j=0; j<numColsPerRow; j++) {
2495 int offset = rowColOffsets[i] + j;
2496 int colIndex = fei::binarySearch(ptCols[offset], remEqns[eqnIndex]->indices());
2497 if (colIndex < 0) ERReturn(-1);
2498
2499 values[i][j] = remEqns[eqnIndex]->coefs()[colIndex];
2500 }
2501 }
2502
2503 //and now, get the local stuff out of the matrix.
2504 for(int i=0; i<numPtRows; i++) {
2505 int row = ptRows[i];
2506 if (row < localStartRow_ || localEndRow_ < row) continue;
2507
2508 int rowLen = 0, checkRowLen;
2509 CHK_ERR( lsc_->getMatrixRowLength(row, rowLen) )
2510 if (rowLen <= 0) ERReturn(-1);
2511
2512 //for each local row, establish some temp arrays and get the row from
2513 //the matrix.
2514
2515 std::vector<double> coefs(rowLen);
2516 std::vector<int> indices(rowLen);
2517
2518 CHK_ERR( lsc_->getMatrixRow(row, &coefs[0], &indices[0], rowLen, checkRowLen) );
2519 if (rowLen != checkRowLen) ERReturn(-1);
2520
2521 //now stride across the list of requested column-indices, and find the
2522 //corresponding location in the matrix row. Copy that location into the
2523 //values array.
2524
2525 //again, if ptCols is NULL, then we're going to copy all coefficients
2526 //(the whole row) into 'values'.
2527 if (ptCols == NULL) {
2528 for(int j=0; j<rowLen; j++) {
2529 values[i][j] = coefs[j];
2530 }
2531 continue;
2532 }
2533
2534 for(int j=0; j<numColsPerRow; j++) {
2535 std::vector<int>::iterator iter =
2536 std::find(indices.begin(), indices.end(), ptCols[rowColOffsets[i]+j]);
2537 if (iter == indices.end()) {
2538 ERReturn(-1);
2539 }
2540
2541 int index = iter - indices.begin();
2542 values[i][j] = coefs[index];
2543 }
2544 }
2545
2546 return(FEI_SUCCESS);
2547}
2548
2549//------------------------------------------------------------------------------
2550int LinSysCoreFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
2551{
2552 //Given a ProcEqns object containing lists of equation-numbers, get the data
2553 //for those equations from the local portion of the matrix and store that data
2554 //in the eqnData object.
2555
2556 std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
2557
2558 for(unsigned p=0; p<eqnNumbers.size(); p++) {
2559 for(unsigned i=0; i<eqnNumbers[p]->size(); i++) {
2560 int eqn = (*(eqnNumbers[p]))[i];
2561
2562 if (localStartRow_ > eqn || localEndRow_ < eqn) continue;
2563
2564 //if this equation is already in eqnData, then don't put it in again...
2565 std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
2566 if (fei::binarySearch(eqn, eqnDataEqns) >= 0) continue;
2567
2568 int len = 0;
2569 CHK_ERR( lsc_->getMatrixRowLength(eqn, len) );
2570
2571 if (len <= 0) continue;
2572 std::vector<double> coefs(len);
2573 std::vector<int> indices(len);
2574 int outputLen = 0;
2575
2576 CHK_ERR( lsc_->getMatrixRow(eqn, &coefs[0], &indices[0],
2577 len, outputLen) );
2578 if (outputLen != len) ERReturn(-1);
2579
2580 CHK_ERR( eqnData.addEqn(eqn, &coefs[0], &indices[0], len, false) );
2581 }
2582 }
2583 return(FEI_SUCCESS);
2584}
2585
2586//------------------------------------------------------------------------------
2587int LinSysCoreFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
2588{
2589 //Given a ProcEqns object containing lists of equation-numbers, get the data
2590 //for those equations from the local portion of the RHS vector and store that
2591 // data in the eqnData object. We're only storing rhs coefs in an EqnBuffer
2592 //that was designed for also storing equations with column-indices. So we'll
2593 //put a bogus column-index in with each equation just to make sure the
2594 //EqnBuffer does the right stuff internally...
2595
2596 int numSendProcs = procEqns.getNumProcs();
2597 std::vector<int>& eqnsPerProc = procEqns.eqnsPerProcPtr();
2598 std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
2599
2600 eqnData.setNumRHSs(1);
2601
2602 for(int p=0; p<numSendProcs; p++) {
2603 for(int i=0; i<eqnsPerProc[p]; i++) {
2604 int reducedEqn;
2605 problemStructure_->translateToReducedEqn((*(eqnNumbers[p]))[i], reducedEqn);
2606
2607 if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
2608
2609 //if this equation is already in eqnData, then don't put it in again...
2610 std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
2611 if (fei::binarySearch(reducedEqn, eqnDataEqns) >= 0) continue;
2612
2613 double rhsValue;
2614
2615 CHK_ERR( lsc_->getFromRHSVector(1, &rhsValue, &reducedEqn) );
2616
2617 int bogusIndex = 19;
2618 CHK_ERR( eqnData.addIndices(reducedEqn, &bogusIndex, 1) );
2619 CHK_ERR( eqnData.addRHS(reducedEqn, 0, rhsValue) );
2620 }
2621 }
2622 return(FEI_SUCCESS);
2623}
2624
2625//------------------------------------------------------------------------------
2626int LinSysCoreFilter::giveToRHS(int num, const double* values,
2627 const int* indices, int mode)
2628{
2629 if (problemStructure_->numSlaveEquations() == 0) {
2630 for(int i=0; i<num; i++) {
2631 if (indices[i] < localStartRow_ || indices[i] > localEndRow_) {
2632 if (mode == ASSEMBLE_SUM) {
2633 eqnCommMgr_->addRemoteRHS(indices[i], currentRHS_, values[i]);
2634 }
2635 else {
2636 eqnCommMgr_put_->addRemoteRHS(indices[i], currentRHS_, values[i]);
2637 }
2638
2639 continue;
2640 }
2641
2642 if (mode == ASSEMBLE_SUM) {
2643 CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
2644 }
2645 else {
2646 CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
2647 }
2648 }
2649 }
2650 else {
2651 for(int i=0; i<num; i++) {
2652 int reducedEqn;
2653 bool isSlave = problemStructure_->
2654 translateToReducedEqn(indices[i], reducedEqn);
2655 if (isSlave) continue;
2656
2657 if (reducedEqn < reducedStartRow_ || reducedEqn > reducedEndRow_) {
2658 if (mode == ASSEMBLE_SUM) {
2659 eqnCommMgr_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
2660 }
2661 else {
2662 eqnCommMgr_put_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
2663 }
2664
2665 continue;
2666 }
2667
2668 if (mode == ASSEMBLE_SUM) {
2669 CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &reducedEqn) );
2670 }
2671 else {
2672 CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &reducedEqn) );
2673 }
2674 }
2675 }
2676
2677 return(FEI_SUCCESS);
2678}
2679
2680//------------------------------------------------------------------------------
2681int LinSysCoreFilter::giveToLocalReducedRHS(int num, const double* values,
2682 const int* indices, int mode)
2683{
2684 for(int i=0; i<num; i++) {
2685
2686 if (mode == ASSEMBLE_SUM) {
2687 CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
2688 }
2689 else {
2690 CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
2691 }
2692 }
2693
2694 return(FEI_SUCCESS);
2695}
2696
2697//------------------------------------------------------------------------------
2698int LinSysCoreFilter::sumIntoRHS(fei::CSVec& vec)
2699{
2700 std::vector<int>& indices = vec.indices();
2701 std::vector<double>& coefs = vec.coefs();
2702
2703 CHK_ERR( giveToRHS(indices.size(), &coefs[0], &indices[0], ASSEMBLE_SUM) );
2704
2705 return(FEI_SUCCESS);
2706}
2707
2708//------------------------------------------------------------------------------
2709int LinSysCoreFilter::getFromRHS(int num, double* values, const int* indices)
2710{
2711 //We need to do similar things here as we do in getFromMatrix, with respect to
2712 //communications to obtain values for equations that are remotely owned.
2713
2714 ProcEqns remoteProcEqns;
2715
2716 //Let's populate this ProcEqns object with the remote equations and the procs
2717 //that we need to receive the remote equations from.
2718 for(int re=0; re<num; re++) {
2719 int eqn = indices[re];
2720 int owner = problemStructure_->getOwnerProcForEqn(eqn);
2721 if (owner==localRank_) continue;
2722
2723 remoteProcEqns.addEqn(eqn, owner);
2724 }
2725
2726 //so now we know which of the requested equations are remotely owned, and we
2727 //know which processors own them.
2728 //Next we're going to need to know which locally-owned equations are needed
2729 //by other processors.
2730 ProcEqns localProcEqns;
2731 CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) );
2732
2733 //ok, now we know which equations we'll need to send, so let's extract
2734 //them from the rhs vector.
2735 EqnBuffer localEqns;
2736 CHK_ERR( getEqnsFromRHS(localProcEqns, localEqns) );
2737
2738 //now we can set the lengths in localProcEqns.
2739 std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
2740 fei::CSVec** localEqnsPtr = &(localEqns.eqns()[0]);
2741 std::vector<int> eqnLengths(eqnNumbers.size());
2742 for(size_t i=0; i<eqnNumbers.size(); ++i) {
2743 eqnLengths[i] = localEqnsPtr[i]->size();
2744 }
2745
2746 localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
2747 eqnNumbers.size());
2748
2749 //now mirror those lengths in the remoteProcEqns objects to get ready for the
2750 //all-to-all exchange of equation data.
2751 CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
2752
2753 EqnBuffer remoteEqns;
2754 //we're now ready to do the exchange.
2755 CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
2756 &remoteProcEqns, &remoteEqns, false))
2757
2758 //now we're ready to get the rhs data we've received from other processors.
2759 std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
2760 std::vector<std::vector<double>*>& remRhsCoefs = *(remoteEqns.rhsCoefsPtr());
2761
2762 for(int i=0; i<num; i++) {
2763 int row = indices[i];
2764
2765 int eqnIndex = fei::binarySearch(row, remEqnNumbers);
2766 if (eqnIndex < 0) continue;
2767
2768 values[i] = (*(remRhsCoefs[eqnIndex]))[0];
2769 }
2770
2771 //and now get the local stuff.
2772 for(int i=0; i<num; i++) {
2773 if (indices[i] < localStartRow_ || indices[i] > localEndRow_) continue;
2774
2775 CHK_ERR( lsc_->getFromRHSVector(num, values, indices) );
2776 }
2777
2778 return(FEI_SUCCESS);
2779}
2780
2781//------------------------------------------------------------------------------
2782int LinSysCoreFilter::getEqnSolnEntry(int eqnNumber, double& solnValue)
2783{
2784 //This function's task is to retrieve the solution-value for a global
2785 //equation-number. eqnNumber may or may not be a slave-equation, and may or
2786 //may not be locally owned. If it is not locally owned, it should at least
2787 //be shared.
2788 //return 0 if the solution is successfully retrieved, otherwise return 1.
2789 //
2790
2791 //
2792 //First and probably most common case: there are no slave equations.
2793 //
2794 if (problemStructure_->numSlaveEquations() == 0) {
2795 if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
2796 //Dig into the eqn-comm-mgr for the shared-remote solution value.
2797 CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
2798 }
2799 else {
2800 //It's local, simply get the solution from the assembled linear system.
2801 CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
2802 }
2803 return(0);
2804 }
2805
2806 //
2807 //If we reach this point, there are slave equations to account for.
2808 //So first translate this equation into 'assembled-linear-system'
2809 //equation-numbers.
2810 //
2811 int reducedEqn;
2812 bool isSlave = problemStructure_->translateToReducedEqn(eqnNumber,reducedEqn);
2813
2814 if (isSlave) {
2815 //This is a slave-equation, so construct its solution-value as the linear-
2816 //combination of the master-equations it is defined in terms of.
2817
2818 std::vector<int>* masterEqns = NULL;
2819 std::vector<double>* masterCoefs = NULL;
2820 CHK_ERR( problemStructure_->getMasterEqnNumbers(eqnNumber, masterEqns) );
2821 CHK_ERR( problemStructure_->getMasterEqnCoefs(eqnNumber, masterCoefs) );
2822
2823 int len = masterEqns->size();
2824 solnValue = 0.0;
2825 CHK_ERR( problemStructure_->getMasterEqnRHS(eqnNumber, solnValue) );
2826
2827 double coef = 0.0;
2828 for(int i=0; i<len; i++) {
2829 int mEqn = (*masterEqns)[i];
2830 int mReducedeqn;
2831 problemStructure_->translateToReducedEqn(mEqn, mReducedeqn);
2832
2833 if (reducedStartRow_ > mReducedeqn || mReducedeqn > reducedEndRow_) {
2834 CHK_ERR( getSharedRemoteSolnEntry(mReducedeqn, coef) );
2835 }
2836 else {
2837 CHK_ERR( getReducedSolnEntry(mReducedeqn, coef) );
2838 }
2839 solnValue += coef * (*masterCoefs)[i];
2840 }
2841 }
2842 else {
2843 //This is not a slave-equation, so retrieve the solution from either the
2844 //assembled linear system or the shared-remote data structures.
2845
2846 if (reducedStartRow_ > reducedEqn || reducedEqn > reducedEndRow_) {
2847 CHK_ERR( getSharedRemoteSolnEntry(reducedEqn, solnValue) );
2848 }
2849 else {
2850 CHK_ERR( getReducedSolnEntry(reducedEqn, solnValue) );
2851 }
2852 }
2853
2854 return(0);
2855}
2856
2857//------------------------------------------------------------------------------
2858int LinSysCoreFilter::getSharedRemoteSolnEntry(int eqnNumber, double& solnValue)
2859{
2860 std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
2861 double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
2862
2863 int index = fei::binarySearch(eqnNumber, remoteEqnNumbers);
2864 if (index < 0) {
2865 fei::console_out() << "LinSysCoreFilter::getSharedRemoteSolnEntry: ERROR, eqn "
2866 << eqnNumber << " not found." << FEI_ENDL;
2867 ERReturn(-1);
2868 }
2869 solnValue = remoteSoln[index];
2870 return(0);
2871}
2872
2873//------------------------------------------------------------------------------
2874int LinSysCoreFilter::getReducedSolnEntry(int eqnNumber, double& solnValue)
2875{
2876 //We may safely assume that this function is called with 'eqnNumber' that is
2877 //local in the underlying assembled linear system. i.e., it isn't a slave-
2878 //equation, it isn't remotely owned, etc.
2879 //
2880 CHK_ERR( lsc_->getSolnEntry(eqnNumber, solnValue) );
2881
2882 return(FEI_SUCCESS);
2883}
2884
2885//------------------------------------------------------------------------------
2886int LinSysCoreFilter::unpackSolution()
2887{
2888 //
2889 //This function should be called after the solver has returned,
2890 //and we know that there is a solution in the underlying vector.
2891 //This function ensures that any locally-owned shared solution values are
2892 //available on the sharing processors.
2893 //
2894 if (Filter::logStream() != NULL) {
2895 (*logStream())<< "# entering unpackSolution, outputLevel: "
2896 <<outputLevel_<<FEI_ENDL;
2897 }
2898
2899 //what we need to do is as follows.
2900 //The eqn comm mgr has a list of what it calls 'recv eqns'. These are
2901 //equations that we own, for which we received contributions from other
2902 //processors. The solution values corresponding to these equations need
2903 //to be made available to those remote contributing processors.
2904
2905 int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
2906 std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
2907
2908 for(int i=0; i<numRecvEqns; i++) {
2909 int eqn = recvEqnNumbers[i];
2910
2911 if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
2912 fei::console_out() << "LinSysCoreFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
2913 << ") out of local range." << FEI_ENDL;
2914 MPI_Abort(comm_, -1);
2915 }
2916
2917 double solnValue = 0.0;
2918
2919 CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
2920
2921 eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
2922 }
2923
2924 eqnCommMgr_->exchangeSoln();
2925
2926 debugOutput("#LinSysCoreFilter leaving unpackSolution");
2927 return(FEI_SUCCESS);
2928}
2929
2930//------------------------------------------------------------------------------
2931void LinSysCoreFilter::setEqnCommMgr(EqnCommMgr* eqnCommMgr)
2932{
2933 delete eqnCommMgr_;
2934 eqnCommMgr_ = eqnCommMgr;
2935}
2936
2937//------------------------------------------------------------------------------
2938int LinSysCoreFilter::getBlockNodeSolution(GlobalID elemBlockID,
2939 int numNodes,
2940 const GlobalID *nodeIDs,
2941 int *offsets,
2942 double *results) {
2943
2944 debugOutput("FEI: getBlockNodeSolution");
2945
2946 int numActiveNodes = problemStructure_->getNumActiveNodes();
2947 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2948
2949 if (numActiveNodes <= 0) return(0);
2950
2951 int numSolnParams = 0;
2952
2953 BlockDescriptor* block = NULL;
2954 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2955
2956 //Traverse the node list, checking if nodes are associated with this block.
2957 //If so, put its 'answers' in the results list.
2958
2959 int offset = 0;
2960 for(int i=0; i<numActiveNodes; i++) {
2961 const NodeDescriptor* node_i = NULL;
2962 nodeDB.getNodeAtIndex(i, node_i);
2963
2964 if (offset == numNodes) break;
2965
2966 GlobalID nodeID = nodeIDs[offset];
2967
2968 //first let's set the offset at which this node's solution coefs start.
2969 offsets[offset++] = numSolnParams;
2970
2971 const NodeDescriptor* node = NULL;
2972 int err = 0;
2973 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
2974 //Don't call the getActiveNodeDesc_ID function unless we have to.
2975
2976 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2977 node = node_i;
2978 }
2979 else {
2980 err = nodeDB.getNodeWithID(nodeID, node);
2981 }
2982
2983 //ok. If err is not 0, meaning nodeID is NOT in the
2984 //activeNodes list, then skip to the next loop iteration.
2985
2986 if (err != 0) {
2987 continue;
2988 }
2989
2990 int numFields = node->getNumFields();
2991 const int* fieldIDs = node->getFieldIDList();
2992
2993 for(int j=0; j<numFields; j++) {
2994 if (block->containsField(fieldIDs[j])) {
2995 int size = problemStructure_->getFieldSize(fieldIDs[j]);
2996 assert(size >= 0);
2997
2998 int thisEqn = -1;
2999 node->getFieldEqnNumber(fieldIDs[j], thisEqn);
3000
3001 double answer;
3002 for(int k=0; k<size; k++) {
3003 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
3004 results[numSolnParams++] = answer;
3005 }
3006 }
3007 }//for(j<numFields)loop
3008 }
3009
3010 offsets[numNodes] = numSolnParams;
3011
3012 return(FEI_SUCCESS);
3013}
3014
3015//------------------------------------------------------------------------------
3016int LinSysCoreFilter::getNodalSolution(int numNodes,
3017 const GlobalID *nodeIDs,
3018 int *offsets,
3019 double *results)
3020{
3021 debugOutput("FEI: getNodalSolution");
3022
3023 int numActiveNodes = problemStructure_->getNumActiveNodes();
3024 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3025
3026 if (numActiveNodes <= 0) return(0);
3027
3028 int numSolnParams = 0;
3029
3030 //Traverse the node list, checking if nodes are local.
3031 //If so, put 'answers' in the results list.
3032
3033 int offset = 0;
3034 for(int i=0; i<numActiveNodes; i++) {
3035 const NodeDescriptor* node_i = NULL;
3036 nodeDB.getNodeAtIndex(i, node_i);
3037
3038 if (offset == numNodes) break;
3039
3040 GlobalID nodeID = nodeIDs[offset];
3041
3042 //first let's set the offset at which this node's solution coefs start.
3043 offsets[offset++] = numSolnParams;
3044
3045 const NodeDescriptor* node = NULL;
3046 int err = 0;
3047 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
3048 //Don't call the getNodeWithID function unless we have to.
3049
3050 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
3051 node = node_i;
3052 }
3053 else {
3054 err = nodeDB.getNodeWithID(nodeID, node);
3055 }
3056
3057 //ok. If err is not 0, meaning nodeID is NOT in the
3058 //activeNodes list, then skip to the next loop iteration.
3059
3060 if (err != 0) {
3061 continue;
3062 }
3063
3064 int numFields = node->getNumFields();
3065 const int* fieldIDs = node->getFieldIDList();
3066
3067 for(int j=0; j<numFields; j++) {
3068 int size = problemStructure_->getFieldSize(fieldIDs[j]);
3069 assert(size >= 0);
3070
3071 int thisEqn = -1;
3072 node->getFieldEqnNumber(fieldIDs[j], thisEqn);
3073
3074 double answer;
3075 for(int k=0; k<size; k++) {
3076 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
3077 results[numSolnParams++] = answer;
3078 }
3079 }//for(j<numFields)loop
3080 }
3081
3082 offsets[numNodes] = numSolnParams;
3083
3084 return(FEI_SUCCESS);
3085}
3086
3087//------------------------------------------------------------------------------
3088int LinSysCoreFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
3089 int fieldID,
3090 int numNodes,
3091 const GlobalID *nodeIDs,
3092 double *results)
3093{
3094 //Note: if the user-supplied nodeIDs list containts nodes which are not in
3095 //the specified element-block, then the corresponding positions in the
3096 //results array are simply not referenced. This is dangerous behavior that
3097 //hasn't gotten me into trouble yet.
3098
3099 debugOutput("FEI: getBlockFieldNodeSolution");
3100
3101 int numActiveNodes = problemStructure_->getNumActiveNodes();
3102 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3103
3104 if (numActiveNodes <= 0) return(0);
3105
3106 BlockDescriptor* block = NULL;
3107 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
3108
3109 int fieldSize = problemStructure_->getFieldSize(fieldID);
3110 if (fieldSize <= 0) ERReturn(-1);
3111
3112 if (!block->containsField(fieldID)) {
3113 fei::console_out() << "LinSysCoreFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
3114 << " not contained in element-block " << (int)elemBlockID << FEI_ENDL;
3115 return(1);
3116 }
3117
3118 //Traverse the node list, checking if nodes are associated with this block.
3119 //If so, put the answers in the results list.
3120
3121 for(int i=0; i<numNodes; i++) {
3122 const NodeDescriptor* node_i = NULL;
3123 nodeDB.getNodeAtIndex(i, node_i);
3124
3125 GlobalID nodeID = nodeIDs[i];
3126
3127 const NodeDescriptor* node = NULL;
3128 int err = 0;
3129 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
3130 //Don't call the getNodeWithID function unless we have to. (getNodeWithID
3131 //does a binary-search over all local nodeIDs, while getNodeAtIndex is
3132 //a direct lookup.) Often the user supplies a nodeIDs list that is in the
3133 //"natural" order, so we don't need to call getNodeWithID at all.
3134
3135 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
3136 node = node_i;
3137 }
3138 else {
3139 err = nodeDB.getNodeWithID(nodeID, node);
3140 }
3141
3142 //ok. If err is not 0, meaning nodeID is NOT in the
3143 //activeNodes list, then skip to the next loop iteration.
3144
3145 if (err != 0) {
3146 continue;
3147 }
3148
3149 int eqnNumber = -1;
3150 bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
3151 if (!hasField) continue;
3152
3153 int offset = fieldSize*i;
3154 for(int j=0; j<fieldSize; j++) {
3155 double answer = 0.0;
3156 CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
3157 results[offset+j] = answer;
3158 }
3159 }
3160
3161 return(FEI_SUCCESS);
3162}
3163
3164//------------------------------------------------------------------------------
3165int LinSysCoreFilter::getNodalFieldSolution(int fieldID,
3166 int numNodes,
3167 const GlobalID *nodeIDs,
3168 double *results)
3169{
3170 debugOutput("FEI: getNodalFieldSolution");
3171
3172 int numActiveNodes = problemStructure_->getNumActiveNodes();
3173 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3174
3175 if (numActiveNodes <= 0) return(0);
3176
3177 int fieldSize = problemStructure_->getFieldSize(fieldID);
3178 if (fieldSize <= 0) ERReturn(-1);
3179
3180 //Traverse the node list, checking if nodes have the specified field.
3181 //If so, put the answers in the results list.
3182
3183 for(int i=0; i<numNodes; i++) {
3184 const NodeDescriptor* node_i = NULL;
3185 nodeDB.getNodeAtIndex(i, node_i);
3186
3187 GlobalID nodeID = nodeIDs[i];
3188
3189 const NodeDescriptor* node = NULL;
3190 int err = 0;
3191 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
3192 //Don't call the getNodeWithID function unless we have to.
3193
3194 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
3195 node = node_i;
3196 }
3197 else {
3198 err = nodeDB.getNodeWithID(nodeID, node);
3199 }
3200
3201 //ok. If err is not 0, meaning nodeID is NOT in the
3202 //activeNodes list, then skip to the next loop iteration.
3203
3204 if (err != 0) {
3205 continue;
3206 }
3207
3208 int eqnNumber = -1;
3209 bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
3210
3211 //If this node doesn't have the specified field, then skip to the
3212 //next loop iteration.
3213 if (!hasField) continue;
3214
3215 int offset = fieldSize*i;
3216 for(int j=0; j<fieldSize; j++) {
3217 double answer = 0.0;
3218 CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
3219 results[offset+j] = answer;
3220 }
3221 }
3222
3223 return(FEI_SUCCESS);
3224}
3225
3226//------------------------------------------------------------------------------
3227int LinSysCoreFilter::putBlockNodeSolution(GlobalID elemBlockID,
3228 int numNodes,
3229 const GlobalID *nodeIDs,
3230 const int *offsets,
3231 const double *estimates) {
3232
3233 debugOutput("FEI: putBlockNodeSolution");
3234
3235 int numActiveNodes = problemStructure_->getNumActiveNodes();
3236
3237 if (numActiveNodes <= 0) return(0);
3238
3239 BlockDescriptor* block = NULL;
3240 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
3241
3242 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3243
3244 //traverse the node list, checking for nodes associated with this block
3245 //when an associated node is found, put its 'answers' into the linear system.
3246
3247 unsigned blk_idx = problemStructure_->getIndexOfBlock(elemBlockID);
3248
3249 for(int i=0; i<numNodes; i++) {
3250 const NodeDescriptor* node = NULL;
3251 int err = nodeDB.getNodeWithID(nodeIDs[i], node);
3252
3253 if (err != 0) continue;
3254
3255 if (!node->hasBlockIndex(blk_idx)) continue;
3256
3257 if (node->getOwnerProc() != localRank_) continue;
3258
3259 int numFields = node->getNumFields();
3260 const int* fieldIDs = node->getFieldIDList();
3261 const int* fieldEqnNumbers = node->getFieldEqnNumbers();
3262
3263 if (fieldEqnNumbers[0] < localStartRow_ ||
3264 fieldEqnNumbers[0] > localEndRow_) continue;
3265
3266 int offs = offsets[i];
3267
3268 for(int j=0; j<numFields; j++) {
3269 int size = problemStructure_->getFieldSize(fieldIDs[j]);
3270
3271 if (block->containsField(fieldIDs[j])) {
3272 for(int k=0; k<size; k++) {
3273 int reducedEqn;
3274 problemStructure_->
3275 translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
3276
3277 CHK_ERR( lsc_->putInitialGuess(&reducedEqn,
3278 &estimates[offs+k], 1) );
3279 }
3280 }
3281 offs += size;
3282 }//for(j<numFields)loop
3283 }
3284
3285 return(FEI_SUCCESS);
3286}
3287
3288//------------------------------------------------------------------------------
3289int LinSysCoreFilter::putBlockFieldNodeSolution(GlobalID elemBlockID,
3290 int fieldID,
3291 int numNodes,
3292 const GlobalID *nodeIDs,
3293 const double *estimates)
3294{
3295 int fieldSize = problemStructure_->getFieldSize(fieldID);
3296
3297 if (Filter::logStream() != NULL) {
3298 FEI_OSTREAM& os = *logStream();
3299 os << "FEI: putBlockFieldNodeSolution" << FEI_ENDL;
3300 os << "#blkID" << FEI_ENDL << (int)elemBlockID << FEI_ENDL
3301 << "#fieldID"<<FEI_ENDL << fieldID << FEI_ENDL
3302 << "#fieldSize"<<FEI_ENDL << fieldSize << FEI_ENDL
3303 << "#numNodes"<<FEI_ENDL << numNodes << FEI_ENDL
3304 << "#nodeIDs" << FEI_ENDL;
3305 int i;
3306 for(i=0; i<numNodes; ++i) os << (int)nodeIDs[i] << FEI_ENDL;
3307 os << "#estimates" << FEI_ENDL;
3308 for(i=0; i<numNodes*fieldSize; ++i) os << estimates[i] << FEI_ENDL;
3309 }
3310
3311 BlockDescriptor* block = NULL;
3312 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
3313 if (!block->containsField(fieldID)) return(1);
3314
3315 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3316
3317 //if we have a negative fieldID, we'll need a list of length numNodes,
3318 //in which to put nodeNumbers for passing to the solver...
3319
3320 std::vector<int> numbers(numNodes);
3321
3322 //if we have a fieldID >= 0, then our numbers list will hold equation numbers
3323 //and we'll need fieldSize*numNodes of them.
3324
3325 std::vector<double> data;
3326
3327 if (fieldID >= 0) {
3328 assert(fieldSize >= 0);
3329 numbers.resize(numNodes*fieldSize);
3330 data.resize(numNodes*fieldSize);
3331 }
3332
3333 int count = 0;
3334
3335 for(int i=0; i<numNodes; i++) {
3336 const NodeDescriptor* node = NULL;
3337 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
3338
3339 if (fieldID < 0) numbers[count++] = node->getNodeNumber();
3340 else {
3341 int eqn = -1;
3342 if (node->getFieldEqnNumber(fieldID, eqn)) {
3343 if (eqn >= localStartRow_ && eqn <= localEndRow_) {
3344 for(int j=0; j<fieldSize; j++) {
3345 data[count] = estimates[i*fieldSize + j];
3346 problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
3347 }
3348 }
3349 }
3350 }
3351 }
3352
3353 if (fieldID < 0) {
3354 CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize,
3355 &numbers[0], numNodes, estimates));
3356 }
3357 else {
3358 CHK_ERR(lsc_->putInitialGuess(&numbers[0], &data[0], count));
3359 }
3360
3361 return(FEI_SUCCESS);
3362}
3363
3364//------------------------------------------------------------------------------
3365int LinSysCoreFilter::getBlockElemSolution(GlobalID elemBlockID,
3366 int numElems,
3367 const GlobalID *elemIDs,
3368 int& numElemDOFPerElement,
3369 double *results)
3370{
3371//
3372// return the elemental solution parameters associated with a
3373// particular block of elements
3374//
3375 debugOutput("FEI: getBlockElemSolution");
3376
3377 BlockDescriptor* block = NULL;
3378 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
3379
3380 numElemDOFPerElement = block->getNumElemDOFPerElement();
3381 if (numElemDOFPerElement <= 0) return(0);
3382
3383 ConnectivityTable& ctable =
3384 problemStructure_->getBlockConnectivity(elemBlockID);
3385 std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
3386
3387 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
3388 double answer;
3389
3390 for(int i=0; i<numElems; i++) {
3391 std::map<GlobalID,int>::const_iterator
3392 iter = elemIDList.find(elemIDs[i]);
3393 if (iter == elemIDList.end()) continue;
3394 int index = iter->second;
3395
3396 int offset = i*numElemDOFPerElement;
3397
3398 for(int j=0; j<numElemDOFPerElement; j++) {
3399 int eqn = elemDOFEqnNumbers[index] + j;
3400
3401 CHK_ERR( getEqnSolnEntry(eqn, answer) )
3402
3403 results[offset+j] = answer;
3404 }
3405 }
3406
3407 return(FEI_SUCCESS);
3408}
3409
3410//------------------------------------------------------------------------------
3411int LinSysCoreFilter::putBlockElemSolution(GlobalID elemBlockID,
3412 int numElems,
3413 const GlobalID *elemIDs,
3414 int dofPerElem,
3415 const double *estimates)
3416{
3417 debugOutput("FEI: putBlockElemSolution");
3418
3419 BlockDescriptor* block = NULL;
3420 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
3421
3422 int DOFPerElement = block->getNumElemDOFPerElement();
3423 assert(DOFPerElement == dofPerElem);
3424 if (DOFPerElement <= 0) return(0);
3425
3426 ConnectivityTable& ctable =
3427 problemStructure_->getBlockConnectivity(elemBlockID);
3428 std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
3429
3430 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
3431
3432
3433 for(int i=0; i<numElems; i++) {
3434 std::map<GlobalID,int>::const_iterator
3435 iter = elemIDList.find(elemIDs[i]);
3436 if (iter == elemIDList.end()) continue;
3437
3438 int index = iter->second;
3439
3440 for(int j=0; j<DOFPerElement; j++) {
3441 int reducedEqn;
3442 problemStructure_->
3443 translateToReducedEqn(elemDOFEqnNumbers[index] + j, reducedEqn);
3444 double soln = estimates[i*DOFPerElement + j];
3445
3446 CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
3447 }
3448 }
3449
3450 return(FEI_SUCCESS);
3451}
3452
3453//------------------------------------------------------------------------------
3454int LinSysCoreFilter::getCRMultipliers(int numCRs,
3455 const int* CRIDs,
3456 double* multipliers)
3457{
3458 int multCRsLen = problemStructure_->getNumMultConstRecords();
3459 if (numCRs > multCRsLen) {
3460 return(-1);
3461 }
3462
3463 std::map<GlobalID, ConstraintType*>::const_iterator
3464 cr_iter = problemStructure_->getMultConstRecords().begin(),
3465 cr_end = problemStructure_->getMultConstRecords().end();
3466
3467 int i = 0;
3468 while(cr_iter != cr_end && i < numCRs) {
3469 GlobalID CRID = (*cr_iter).first;
3470 ConstraintType* multCR = (*cr_iter).second;
3471 if (CRID != CRIDs[i]) {
3472 CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[i], multCR) );
3473 }
3474
3475 int eqn = multCR->getEqnNumber();
3476
3477 CHK_ERR( getEqnSolnEntry(eqn, multipliers[i]) );
3478 ++cr_iter;
3479 ++i;
3480 }
3481
3482 return(FEI_SUCCESS);
3483}
3484
3485//------------------------------------------------------------------------------
3486int LinSysCoreFilter::putCRMultipliers(int numMultCRs,
3487 const int* CRIDs,
3488 const double *multEstimates)
3489{
3490 debugOutput("FEI: putCRMultipliers");
3491
3492 for(int j = 0; j < numMultCRs; j++) {
3493 ConstraintType* multCR = NULL;
3494 CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[j], multCR) );
3495
3496 int eqnNumber = multCR->getEqnNumber();
3497 if (eqnNumber < localStartRow_ || eqnNumber > localEndRow_) continue;
3498 CHK_ERR( lsc_->putInitialGuess(&eqnNumber, &(multEstimates[j]), 1));
3499 }
3500
3501 return(FEI_SUCCESS);
3502}
3503
3504//------------------------------------------------------------------------------
3505int LinSysCoreFilter::putNodalFieldData(int fieldID,
3506 int numNodes,
3507 const GlobalID* nodeIDs,
3508 const double* nodeData)
3509{
3510 debugOutput("FEI: putNodalFieldData");
3511
3512 if (fieldID > -1) {
3513 return(putNodalFieldSolution(fieldID, numNodes, nodeIDs, nodeData));
3514 }
3515
3516 int fieldSize = problemStructure_->getFieldSize(fieldID);
3517 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3518
3519 std::vector<int> nodeNumbers(numNodes);
3520
3521 for(int i=0; i<numNodes; i++) {
3522 const NodeDescriptor* node = NULL;
3523 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
3524
3525 int nodeNumber = node->getNodeNumber();
3526 if (nodeNumber < 0) {
3527 fei::console_out() << "LinSysCoreFilter::putNodalFieldData ERROR, node with ID "
3528 << (int)nodeIDs[i] << " doesn't have an associated nodeNumber "
3529 << "assigned. putNodalFieldData shouldn't be called until after the "
3530 << "initComplete method has been called." << FEI_ENDL;
3531 ERReturn(-1);
3532 }
3533
3534 nodeNumbers[i] = nodeNumber;
3535 }
3536
3537 CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize,
3538 &nodeNumbers[0], numNodes, nodeData) );
3539
3540 return(0);
3541}
3542
3543//------------------------------------------------------------------------------
3544int LinSysCoreFilter::putNodalFieldSolution(int fieldID,
3545 int numNodes,
3546 const GlobalID* nodeIDs,
3547 const double* nodeData)
3548{
3549 debugOutput("FEI: putNodalFieldSolution");
3550
3551 if (fieldID < 0) {
3552 return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
3553 }
3554
3555 int fieldSize = problemStructure_->getFieldSize(fieldID);
3556 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
3557
3558 std::vector<int> eqnNumbers(fieldSize);
3559
3560 for(int i=0; i<numNodes; i++) {
3561 const NodeDescriptor* node = NULL;
3562 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
3563
3564 int eqn = -1;
3565 bool hasField = node->getFieldEqnNumber(fieldID, eqn);
3566 if (!hasField) continue;
3567
3568 int reducedEqn = -1;
3569 bool isSlave = problemStructure_->translateToReducedEqn(eqn, reducedEqn);
3570 if (isSlave) continue;
3571
3572 if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
3573
3574 int localLen = fieldSize;
3575 for(int j=0; j<fieldSize; j++) {
3576 int thisEqn = reducedEqn+j;
3577 if (reducedStartRow_ > thisEqn || reducedEndRow_ <thisEqn) {
3578 localLen = j;
3579 }
3580
3581 eqnNumbers[j] = reducedEqn+j;
3582 }
3583
3584 int offset = i*fieldSize;
3585 CHK_ERR( lsc_->putInitialGuess(&eqnNumbers[0],
3586 &nodeData[offset], localLen) );
3587 }
3588
3589 return(0);
3590}
3591
3592//------------------------------------------------------------------------------
3593int LinSysCoreFilter::assembleEqns(int numPtRows,
3594 int numPtCols,
3595 const int* rowNumbers,
3596 const int* colIndices,
3597 const double* const* coefs,
3598 bool structurallySymmetric,
3599 int numBlkEqns, int* blkEqns,
3600 int* blkSizes, bool useBlkEqns,
3601 int mode)
3602{
3603 if (numPtRows == 0) return(FEI_SUCCESS);
3604
3605 bool anySlaves = false;
3606 int numSlaveEqns = problemStructure_->numSlaveEquations();
3607 if (numSlaveEqns > 0) {
3608 rSlave_.resize(numPtRows);
3609 cSlave_.resize(0);
3610 const int* indPtr = colIndices;
3611 for(int r=0; r<numPtRows; r++) {
3612 rSlave_[r] = problemStructure_->isSlaveEqn(rowNumbers[r]) ? 1 : 0;
3613 if (rSlave_[r] == 1) anySlaves = true;
3614
3615 for(int j=0; j<numPtCols; j++) {
3616 int isSlave = problemStructure_->isSlaveEqn(indPtr[j]) ? 1 : 0;
3617 cSlave_.push_back(isSlave);
3618 if (isSlave == 1) anySlaves = true;
3619 }
3620
3621 if (!structurallySymmetric) indPtr += numPtCols;
3622 }
3623 }
3624
3625 if (numSlaveEqns == 0 || !anySlaves) {
3626 if (numSlaveEqns == 0 && structurallySymmetric) {
3627 if (useBlkEqns) {
3628 CHK_ERR( giveToBlkMatrix_symm_noSlaves(numPtRows, rowNumbers,
3629 numBlkEqns, blkEqns, blkSizes,
3630 coefs, mode) );
3631 }
3632 else {
3633 CHK_ERR( giveToMatrix_symm_noSlaves(numPtRows, rowNumbers, coefs, mode) );
3634 }
3635 }
3636 else {
3637 if ((int)dworkSpace2_.size() < numPtRows) {
3638 dworkSpace2_.resize(numPtRows);
3639 }
3640 const double* * coefPtr = &dworkSpace2_[0];
3641 for(int i=0; i<numPtRows; i++) {
3642 coefPtr[i] = coefs[i];
3643 }
3644
3645 if (structurallySymmetric) {
3646 CHK_ERR( giveToMatrix(numPtRows, rowNumbers, numPtRows, rowNumbers,
3647 coefPtr, mode) );
3648 }
3649 else {
3650 const int* indPtr = colIndices;
3651 for(int i=0; i<numPtRows; i++) {
3652 int row = rowNumbers[i];
3653
3654 const double* coefPtr1 = coefs[i];
3655
3656 CHK_ERR(giveToMatrix(1, &row, numPtCols, indPtr, &coefPtr1, mode));
3657 indPtr += numPtCols;
3658 }
3659 }
3660 }
3661 }
3662 else {
3663 int offset = 0;
3664 const int* indicesPtr = colIndices;
3665 for(int i=0; i<numPtRows; i++) {
3666 int row = rowNumbers[i];
3667
3668 const double* coefPtr = coefs[i];
3669 int* colSlave = &cSlave_[offset];
3670 offset += numPtCols;
3671
3672 if (rSlave_[i] == 1) {
3673 //Since this is a slave equation, the non-slave columns of this row go
3674 //into 'Kdi_', and the slave columns go into 'Kdd_'.
3675 for(int jj=0; jj<numPtCols; jj++) {
3676 int col = indicesPtr[jj];
3677 if (colSlave[jj]) {
3678 Kdd_->sumInCoef(row, col, coefPtr[jj]);
3679 }
3680 else {
3681 Kdi_->sumInCoef(row, col, coefPtr[jj]);
3682 }
3683 }
3684
3685 //We also need to put the non-slave rows of column 'row' into 'K_id'.
3686 const int* ii_indicesPtr = colIndices;
3687 for(int ii=0; ii<numPtRows; ii++) {
3688 int rowi = rowNumbers[ii];
3689 if (rSlave_[ii] == 1) continue;
3690
3691 int index = fei::binarySearch(row, ii_indicesPtr, numPtCols);
3692 if (index < 0) continue;
3693
3694 const double* coefs_ii = coefs[ii];
3695
3696 Kid_->sumInCoef(rowi, row, coefs_ii[index]);
3697
3698 if (!structurallySymmetric) ii_indicesPtr += numPtCols;
3699 }
3700
3701 reducedEqnCounter_++;
3702
3703 continue;
3704 }
3705 else {//row is not a slave eqn...
3706
3707 //put all non-slave columns from this row into the assembled matrix.
3708
3709 CHK_ERR( giveToMatrix(1, &row, numPtCols, indicesPtr, &coefPtr, mode) );
3710 }
3711
3712 if (!structurallySymmetric) indicesPtr += numPtCols;
3713 }
3714
3715 if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedEqns() );
3716 }
3717
3718 return(FEI_SUCCESS);
3719}
3720
3721//------------------------------------------------------------------------------
3722int LinSysCoreFilter::assembleReducedEqns()
3723{
3724 fei::FillableMat* D = problemStructure_->getSlaveDependencies();
3725
3726 csrD = *D;
3727 csrKid = *Kid_;
3728 csrKdi = *Kdi_;
3729 csrKdd = *Kdd_;
3730
3731 //form tmpMat1_ = Kid_*D
3732 fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
3733
3734 //form tmpMat2_ = D^T*Kdi_
3735 fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
3736
3737 if (Filter::logStream() != NULL) {
3738 FEI_OSTREAM& os = *Filter::logStream();
3739 os << "# tmpMat1_"<<FEI_ENDL << tmpMat1_ << FEI_ENDL;
3740 os << "# tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
3741 }
3742
3743 //accumulate the above two results into the global system matrix.
3744 CHK_ERR( sumIntoMatrix(tmpMat1_) );
3745 CHK_ERR( sumIntoMatrix(tmpMat2_) );
3746
3747 //form tmpMat1_ = D^T*Kdd_
3748 fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
3749
3750 //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
3751 fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
3752
3753 if (Filter::logStream() != NULL) {
3754 FEI_OSTREAM& os = *Filter::logStream();
3755 os << "# tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
3756 }
3757
3758 //finally, accumulate tmpMat2_ = D^T*Kdd_*D into the global system matrix.
3759 CHK_ERR( sumIntoMatrix(tmpMat2_) );
3760
3761 Kdi_->clear();
3762 Kid_->clear();
3763 Kdd_->clear();
3764 reducedEqnCounter_ = 0;
3765
3766 return(0);
3767}
3768
3769//------------------------------------------------------------------------------
3770int LinSysCoreFilter::assembleRHS(int numValues,
3771 const int* indices,
3772 const double* coefs,
3773 int mode) {
3774//
3775//This function hands the data off to the routine that finally
3776//sticks it into the RHS vector.
3777//
3778
3779 if (problemStructure_->numSlaveEquations() == 0) {
3780 CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
3781 return(FEI_SUCCESS);
3782 }
3783
3784 for(int i = 0; i < numValues; i++) {
3785 int eqn = indices[i];
3786
3787 if (problemStructure_->isSlaveEqn(eqn)) {
3788 fei::add_entry( fd_, eqn, coefs[i]);
3789 reducedRHSCounter_++;
3790 continue;
3791 }
3792
3793 CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
3794 }
3795
3796 if (reducedRHSCounter_ > 300) CHK_ERR( assembleReducedRHS() );
3797
3798 return(FEI_SUCCESS);
3799}
3800
3801//------------------------------------------------------------------------------
3802int LinSysCoreFilter::assembleReducedRHS()
3803{
3804 fei::FillableMat* D = problemStructure_->getSlaveDependencies();
3805
3806 csrD = *D;
3807
3808 //now form tmpVec1_ = D^T*fd_.
3809 fei::multiply_trans_CSRMat_CSVec(csrD, fd_, tmpVec1_);
3810
3811 CHK_ERR( sumIntoRHS(tmpVec1_) );
3812
3813 fd_.clear();
3814 reducedRHSCounter_ = 0;
3815
3816 return(0);
3817}
3818
3819//==============================================================================
3820void LinSysCoreFilter::debugOutput(const char* mesg) {
3821 if (Filter::logStream() != NULL) {
3822 (*logStream())<<mesg<<FEI_ENDL;
3823 }
3824}
void setNumRHSs(int n)
std::vector< int > & eqnNumbers()
int addEqn(int eqnNumber, const double *coefs, const int *indices, int len, bool accumulate, bool create_indices_union=false)
int addRHS(int eqnNumber, int rhsIndex, double value, bool accumulate=true)
std::vector< fei::CSVec * > & eqns()
int addIndices(int eqnNumber, const int *indices, int len)
EqnCommMgr * deepCopy()
int mirrorProcEqns(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
int mirrorProcEqnLengths(ProcEqns &inProcEqns, ProcEqns &outProcEqns)
int getParameters(int &numParams, char **&paramStrings)
virtual int setMultCRComplete()=0
virtual int resetConstraints(double s)=0
virtual int getMatrixRow(int row, double *coefs, int *indices, int len, int &rowLength)=0
virtual int sumIntoSystemMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, int numBlkRows, const int *blkRows, int numBlkCols, const int *blkCols, const double *const *values)=0
virtual int putIntoSystemMatrix(int numPtRows, const int *ptRows, int numPtCols, const int *ptCols, const double *const *values)=0
virtual int getFromRHSVector(int num, double *values, const int *indices)=0
virtual int setLookup(Lookup &lookup)=0
virtual int formResidual(double *values, int len)=0
virtual int sumIntoRHSVector(int num, const double *values, const int *indices)=0
virtual int enforceEssentialBC(int *globalEqn, double *alpha, double *gamma, int len)=0
virtual int setLoadVectors(GlobalID elemBlock, int numElems, const GlobalID *elemIDs, const double *const *load, int numEqnsPerElem, const int *const *eqnIndices)=0
virtual int getMatrixRowLength(int row, int &length)=0
virtual int setConnectivities(GlobalID elemBlock, int numElements, int numNodesPerElem, const GlobalID *elemIDs, const int *const *connNodes)=0
virtual int resetMatrix(double s)=0
virtual int setPenCREqns(int penCRSetID, int numCRs, int numNodesPerCR, int **nodeNumbers, int **eqnNumbers, int *fieldIDs)=0
virtual int matrixLoadComplete()=0
virtual int putIntoRHSVector(int num, const double *values, const int *indices)=0
virtual int resetRHSVector(double s)=0
virtual int setStiffnessMatrices(GlobalID elemBlock, int numElems, const GlobalID *elemIDs, const double *const *const *stiff, int numEqnsPerElem, const int *const *eqnIndices)=0
virtual int launchSolver(int &solveStatus, int &iterations)=0
virtual int setMultCREqns(int multCRSetID, int numCRs, int numNodesPerCR, int **nodeNumbers, int **eqnNumbers, int *fieldIDs, int *multiplierEqnNumbers)=0
virtual int setGlobalOffsets(int len, int *nodeOffsets, int *eqnOffsets, int *blkEqnOffsets)=0
virtual int getSolnEntry(int eqnNumber, double &answer)=0
virtual int putInitialGuess(const int *eqnNumbers, const double *values, int len)=0
virtual int setMatrixStructure(int **ptColIndices, int *ptRrowLengths, int **blkColIndices, int *blkRowLengths, int *ptRowsPerBlkRow)=0
virtual int putNodalFieldData(int fieldID, int fieldSize, int *nodeNumbers, int numNodes, const double *data)=0
virtual int enforceRemoteEssBCs(int numEqns, int *globalEqns, int **colIndices, int *colIndLen, double **coefs)=0
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int getNumNodeDescriptors() const
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
size_t getNumProcs()
void addEqn(int eqnNumber, int proc)
std::vector< int > & eqnsPerProcPtr()
void setProcEqnLengths(int *eqnNumbers, int *eqnLengths, int len)
std::vector< std::vector< int > * > & procEqnNumbersPtr()
int getMasterEqnRHS(int slaveEqn, double &rhsValue)
int getAssociatedNodeNumber(int eqnNumber)
int getFieldSize(int fieldID)
int getBlockDescriptor_index(int index, BlockDescriptor *&block)
int getMasterEqnCoefs(int slaveEqn, std::vector< double > *&masterCoefs)
int getIndexOfBlock(GlobalID blockID) const
int getMasterEqnNumbers(int slaveEqn, std::vector< int > *&masterEqns)
int getOwnerProcForEqn(int eqn)
bool translateToReducedEqn(int eqn, int &reducedEqn)
std::vector< int > rowNumbers
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
std::vector< int > & getMasterFieldIDs()
std::vector< int > & getMasters()
void setRHSValue(double rhs)
std::vector< double > & getMasterWeights()
double cpu_time()
Definition fei_utils.cpp:46
int localProc(MPI_Comm comm)
void multiply_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
void multiply_trans_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
int Allgatherv(MPI_Comm comm, std::vector< T > &sendbuf, std::vector< int > &recvLengths, std::vector< T > &recvbuf)
void multiply_trans_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int numProcs(MPI_Comm comm)