FEI Version of the Day
Loading...
Searching...
No Matches
fei_FEDataFilter.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_macros.hpp>
10
11#include <limits>
12#include <cmath>
13
14#include <fei_defs.h>
15
16#include <fei_CommUtils.hpp>
17#include <fei_TemplateUtils.hpp>
18#include <snl_fei_Constraint.hpp>
20
21#include <fei_LibraryWrapper.hpp>
22#include <SNL_FEI_Structure.hpp>
23#include <fei_FiniteElementData.hpp>
24#include <fei_Lookup.hpp>
25#include <FEI_Implementation.hpp>
26#include <fei_EqnCommMgr.hpp>
27#include <fei_EqnBuffer.hpp>
28#include <fei_NodeDatabase.hpp>
29#include <fei_NodeCommMgr.hpp>
30#include <fei_ProcEqns.hpp>
31#include <fei_BlockDescriptor.hpp>
32#include <fei_ConnectivityTable.hpp>
33#include <snl_fei_Utils.hpp>
34
35#include <fei_FEDataFilter.hpp>
36
37#undef fei_file
38#define fei_file "FEDataFilter.cpp"
39#include <fei_ErrMacros.hpp>
40
41#define ASSEMBLE_PUT 0
42#define ASSEMBLE_SUM 1
43
44//------------------------------------------------------------------------------
45void convert_eqns_to_nodenumbers_and_dof_ids(fei::FieldDofMap<int>& fdmap,
46 const NodeDatabase& nodeDB,
47 int numEqns,
48 const int* eqns,
49 std::vector<int>& nodeNumbers,
50 std::vector<int>& dof_ids)
51{
52 nodeNumbers.resize(numEqns);
53 dof_ids.resize(numEqns);
54
55 for(int i=0; i<numEqns; ++i) {
56 const NodeDescriptor* nodePtr = NULL;
57 int err = nodeDB.getNodeWithEqn(eqns[i], nodePtr);
58 if (err < 0) {
59 nodeNumbers[i] = -1;
60 dof_ids[i] = -1;
61 continue;
62 }
63
64 nodeNumbers[i] = nodePtr->getNodeNumber();
65
66 int fieldID, offset;
67 nodePtr->getFieldID(eqns[i], fieldID, offset);
68 dof_ids[i] = fdmap.get_dof_id(fieldID, offset);
69 }
70}
71
72//------------------------------------------------------------------------------
73void convert_field_and_nodes_to_eqns(const NodeDatabase& nodeDB,
74 int fieldID, int fieldSize,
75 int numNodes, const GlobalID* nodeIDs,
76 std::vector<int>& eqns)
77{
78 eqns.assign(numNodes*fieldSize, -1);
79
80 size_t offset = 0;
81 for(int i=0; i<numNodes; ++i) {
82 const NodeDescriptor* node = NULL;
83 int err = nodeDB.getNodeWithID(nodeIDs[i], node);
84 if (err < 0) {
85 offset += fieldSize;
86 continue;
87 }
88
89 int eqn = 0;
90 node->getFieldEqnNumber(fieldID, eqn);
91 for(int j=0; j<fieldSize; ++j) {
92 eqns[offset++] = eqn+j;
93 }
94 }
95}
96
97//------------------------------------------------------------------------------
98FEDataFilter::FEDataFilter(FEI_Implementation* owner,
99 MPI_Comm comm,
100 SNL_FEI_Structure* probStruct,
101 LibraryWrapper* wrapper,
102 int masterRank)
103 : Filter(probStruct),
104 wrapper_(wrapper),
105 feData_(NULL),
106 useLookup_(true),
107 internalFei_(0),
108 newData_(false),
109 localStartRow_(0),
110 localEndRow_(0),
111 numGlobalEqns_(0),
112 reducedStartRow_(0),
113 reducedEndRow_(0),
114 numReducedRows_(0),
115 iterations_(0),
116 numRHSs_(0),
117 currentRHS_(0),
118 rhsIDs_(),
119 outputLevel_(0),
120 comm_(comm),
121 masterRank_(masterRank),
122 problemStructure_(probStruct),
123 penCRIDs_(),
124 rowIndices_(),
125 rowColOffsets_(0),
126 colIndices_(0),
127 eqnCommMgr_(NULL),
128 eqnCommMgr_put_(NULL),
129 maxElemRows_(0),
130 eStiff_(NULL),
131 eStiff1D_(NULL),
132 eLoad_(NULL),
133 numRegularElems_(0),
134 constraintBlocks_(0, 16),
135 constraintNodeOffsets_(),
136 packedFieldSizes_()
137{
138 localRank_ = fei::localProc(comm_);
139 numProcs_ = fei::numProcs(comm_);
140
141 internalFei_ = 0;
142
143 numRHSs_ = 1;
144 rhsIDs_.resize(numRHSs_);
145 rhsIDs_[0] = 0;
146
147 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
148 createEqnCommMgr_put();
149
150 if (wrapper_->haveFiniteElementData()) {
151 feData_ = wrapper_->getFiniteElementData();
152 }
153 else {
154 fei::console_out() << "FEDataFilter::FEDataFilter ERROR, must be constructed with a "
155 << "FiniteElementData interface. Aborting." << FEI_ENDL;
156#ifndef FEI_SER
157 MPI_Abort(comm_, -1);
158#else
159 abort();
160#endif
161 }
162
163 //We need to get the parameters from the owning FEI_Implementation, if we've
164 //been given a non-NULL FEI_Implementation...
165 if (owner != NULL) {
166 int numParams = -1;
167 char** paramStrings = NULL;
168 int err = owner->getParameters(numParams, paramStrings);
169
170 //Now let's pass them into our own parameter-handling mechanism.
171 err = parameters(numParams, paramStrings);
172 if (err != 0) {
173 fei::console_out() << "FEDataFilter::FEDataFilter ERROR, parameters failed." << FEI_ENDL;
174 MPI_Abort(comm_, -1);
175 }
176 }
177
178 return;
179}
180
181//------------------------------------------------------------------------------
182FEDataFilter::FEDataFilter(const FEDataFilter& src)
183 : Filter(NULL),
184 wrapper_(NULL),
185 feData_(NULL),
186 useLookup_(true),
187 internalFei_(0),
188 newData_(false),
189 localStartRow_(0),
190 localEndRow_(0),
191 numGlobalEqns_(0),
192 reducedStartRow_(0),
193 reducedEndRow_(0),
194 numReducedRows_(0),
195 iterations_(0),
196 numRHSs_(0),
197 currentRHS_(0),
198 rhsIDs_(),
199 outputLevel_(0),
200 comm_(0),
201 masterRank_(0),
202 problemStructure_(NULL),
203 penCRIDs_(),
204 rowIndices_(),
205 rowColOffsets_(0),
206 colIndices_(0),
207 eqnCommMgr_(NULL),
208 eqnCommMgr_put_(NULL),
209 maxElemRows_(0),
210 eStiff_(NULL),
211 eStiff1D_(NULL),
212 eLoad_(NULL),
213 numRegularElems_(0),
214 constraintBlocks_(0, 16),
215 constraintNodeOffsets_(),
216 packedFieldSizes_()
217{
218}
219
220//------------------------------------------------------------------------------
221FEDataFilter::~FEDataFilter() {
222//
223// Destructor function. Free allocated memory, etc.
224//
225 numRHSs_ = 0;
226
227 delete eqnCommMgr_;
228 delete eqnCommMgr_put_;
229
230 delete [] eStiff_;
231 delete [] eStiff1D_;
232 delete [] eLoad_;
233}
234
235//------------------------------------------------------------------------------
236int FEDataFilter::initialize()
237{
238// Determine final sparsity pattern for setting the structure of the
239// underlying sparse matrix.
240//
241 debugOutput("# initialize");
242
243 // now, obtain the global equation info, such as how many equations there
244 // are globally, and what the local starting and ending row-numbers are.
245
246 // let's also get the number of global nodes, and a first-local-node-number.
247 // node-number is a globally 0-based number we are assigning to nodes.
248 // node-numbers are contiguous on a processor -- i.e., a processor owns a
249 // contiguous block of node-numbers. This provides an easier-to-work-with
250 // node numbering than the application-supplied nodeIDs which may not be
251 // assumed to be contiguous or 0-based, or anything else.
252
253 std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
254 localStartRow_ = eqnOffsets[localRank_];
255 localEndRow_ = eqnOffsets[localRank_+1]-1;
256 numGlobalEqns_ = eqnOffsets[numProcs_];
257
258 //--------------------------------------------------------------------------
259 // ----- end active equation calculations -----
260
261 if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
262 eqnCommMgr_ = NULL;
263 if (eqnCommMgr_put_ != NULL) delete eqnCommMgr_put_;
264 eqnCommMgr_put_ = NULL;
265
266 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
267 if (eqnCommMgr_ == NULL) ERReturn(-1);
268
269 int err = createEqnCommMgr_put();
270 if (err != 0) ERReturn(err);
271
272 //(we need to set the number of RHSs in the eqn comm manager)
273 eqnCommMgr_->setNumRHSs(numRHSs_);
274
275 //let's let the underlying linear system know what the global offsets are.
276 //While we're dealing with global offsets, we'll also calculate the starting
277 //and ending 'reduced' rows, etc.
278 CHK_ERR( initLinSysCore() );
279
280 return(FEI_SUCCESS);
281}
282
283//------------------------------------------------------------------------------
284int FEDataFilter::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 FEDataFilter::initLinSysCore()
298{
299 try {
300
301 int err = wrapper_->getFiniteElementData()->setLookup(*problemStructure_);
302
303 if (err != 0) {
304 useLookup_ = false;
305 }
306
307 reducedStartRow_ = localStartRow_;
308 reducedEndRow_ = localEndRow_;
309
310 int numElemBlocks = problemStructure_->getNumElemBlocks();
311 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
312 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
313
314 int numNodes = nodeDB.getNumNodeDescriptors();
315 int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
316 nodeCommMgr.getLocalNodeIDs().size();
317 numNodes -= numRemoteNodes;
318
319 int numSharedNodes = nodeCommMgr.getNumSharedNodes();
320
321 std::vector<int> numElemsPerBlock(numElemBlocks);
322 std::vector<int> numNodesPerElem(numElemBlocks);
323 std::vector<int> elemMatrixSizePerBlock(numElemBlocks);
324
325 for(int blk=0; blk<numElemBlocks; blk++) {
326 BlockDescriptor* block = NULL;
327 CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
328
329 numElemsPerBlock[blk] = block->getNumElements();
330 numNodesPerElem[blk] = block->getNumNodesPerElement();
331
332 int* fieldsPerNode = block->fieldsPerNodePtr();
333 int** fieldIDsTable = block->fieldIDsTablePtr();
334
335 elemMatrixSizePerBlock[blk] = 0;
336
337 for(int nn=0; nn<numNodesPerElem[blk]; nn++) {
338 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
339
340 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
341 elemMatrixSizePerBlock[blk] +=
342 problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
343 }
344 }
345 }
346
347 //Now we need to run the penalty constraint records and figure out how many
348 //extra "element-blocks" to describe. (A penalty constraint will be treated
349 //exactly like an element.) So first, we need to figure out how many different
350 //sizes of constraint connectivities there are, because the constraints with
351 //the same numbers of constrained nodes will be grouped together in blocks.
352
353 if (problemStructure_==NULL) {
354 FEI_COUT << "problemStructrue_ NULL"<<FEI_ENDL;
355 ERReturn(-1);
356 }
357
358 std::map<GlobalID,ConstraintType*>::const_iterator
359 cr_iter = problemStructure_->getPenConstRecords().begin(),
360 cr_end = problemStructure_->getPenConstRecords().end();
361
362 //constraintBlocks will be a sorted list with each "block-id" being the
363 //num-nodes-per-constraint for constraints in that block.
364
365 //numConstraintsPerBlock is the same length as constraintBlocks
366 std::vector<int> numConstraintsPerBlock;
367 std::vector<int> numDofPerConstraint;
368 penCRIDs_.resize(problemStructure_->getNumPenConstRecords());
369
370 int counter = 0;
371 while(cr_iter != cr_end) {
372 penCRIDs_[counter++] = (*cr_iter).first;
373 ConstraintType& cr = *((*cr_iter).second);
374 int nNodes = cr.getMasters().size();
375
376 int insertPoint = -1;
377 int offset = fei::binarySearch(nNodes, constraintBlocks_, insertPoint);
378
379 int nodeOffset = 0;
380 if (offset < 0) {
381 constraintBlocks_.insert(constraintBlocks_.begin()+insertPoint, nNodes);
382 numConstraintsPerBlock.insert(numConstraintsPerBlock.begin()+insertPoint, 1);
383 numDofPerConstraint.insert(numDofPerConstraint.begin()+insertPoint, 0);
384
385 if (insertPoint > 0) {
386 nodeOffset = constraintNodeOffsets_[insertPoint-1] +
387 constraintBlocks_[insertPoint-1];
388 }
389 constraintNodeOffsets_.insert(constraintNodeOffsets_.begin()+insertPoint, nodeOffset);
390 offset = insertPoint;
391 }
392 else {
393 numConstraintsPerBlock[offset]++;
394 ++cr_iter;
395 continue;
396 }
397
398 std::vector<int>& fieldIDsvec = cr.getMasterFieldIDs();
399 int* fieldIDs = &fieldIDsvec[0];
400 for(int k=0; k<nNodes; ++k) {
401 int fieldSize = problemStructure_->getFieldSize(fieldIDs[k]);
402 packedFieldSizes_.insert(packedFieldSizes_.begin()+nodeOffset+k, fieldSize);
403 numDofPerConstraint[offset] += fieldSize;
404 }
405 ++cr_iter;
406 }
407
408 //now combine the elem-block info with the penalty-constraint info.
409 int numBlocksTotal = numElemBlocks + constraintBlocks_.size();
410 for(size_t i=0; i<constraintBlocks_.size(); ++i) {
411 numElemsPerBlock.push_back(numConstraintsPerBlock[i]);
412 numNodesPerElem.push_back(constraintBlocks_[i]);
413 elemMatrixSizePerBlock.push_back(numDofPerConstraint[i]);
414 }
415
416 int numMultCRs = problemStructure_->getNumMultConstRecords();
417
418 CHK_ERR( feData_->describeStructure(numBlocksTotal,
419 &numElemsPerBlock[0],
420 &numNodesPerElem[0],
421 &elemMatrixSizePerBlock[0],
422 numNodes,
423 numSharedNodes,
424 numMultCRs) );
425
426 numRegularElems_ = 0;
427 std::vector<int> numDofPerNode;
428 std::vector<int> dof_ids;
429 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
430
431 for(int i=0; i<numElemBlocks; i++) {
432 BlockDescriptor* block = NULL;
433 CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
434
435 if (block->getNumElements() == 0) continue;
436
437 ConnectivityTable& ctbl =
438 problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
439
440 std::vector<int> cNodeList(block->getNumNodesPerElement());
441
442 int* fieldsPerNode = block->fieldsPerNodePtr();
443 int** fieldIDsTable = block->fieldIDsTablePtr();
444
445 numDofPerNode.resize(0);
446 int total_num_dof = 0;
447 for(int nn=0; nn<numNodesPerElem[i]; nn++) {
448 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
449 numDofPerNode.push_back(0);
450 int indx = numDofPerNode.size()-1;
451
452 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
453 numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
454 }
455 total_num_dof += numDofPerNode[indx];
456 }
457
458 dof_ids.resize(total_num_dof);
459 int doffset = 0;
460 for(int nn=0; nn<numNodesPerElem[i]; ++nn) {
461 for(int nf=0; nf<fieldsPerNode[nn]; ++nf) {
462 int fieldSize = problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
463 for(int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
464 dof_ids[doffset++] = fdmap.get_dof_id(fieldIDsTable[nn][nf], dof_offset);
465 }
466 }
467 }
468
469 int nodesPerElement = block->getNumNodesPerElement();
470 NodeDescriptor** elemConn = &((*ctbl.elem_conn_ptrs)[0]);
471 int offset = 0;
472 int numElems = block->getNumElements();
473 numRegularElems_ += numElems;
474 for(int j=0; j<numElems; j++) {
475
476 for(int k=0; k<nodesPerElement; k++) {
477 NodeDescriptor* node = elemConn[offset++];
478 cNodeList[k] = node->getNodeNumber();
479 }
480
481 CHK_ERR( feData_->setConnectivity(i, ctbl.elemNumbers[j],
482 block->getNumNodesPerElement(),
483 &cNodeList[0],
484 &numDofPerNode[0],
485 &dof_ids[0]) );
486 }
487 }
488
489 std::vector<int> nodeNumbers;
490 cr_iter = problemStructure_->getPenConstRecords().begin();
491 int i = 0;
492 while(cr_iter != cr_end) {
493 ConstraintType& cr = *((*cr_iter).second);
494 std::vector<GlobalID>& nodeIDsvec = cr.getMasters();
495 GlobalID* nodeIDs = &nodeIDsvec[0];
496 int nNodes = cr.getMasters().size();
497 int index = fei::binarySearch(nNodes, constraintBlocks_);
498 if (index < 0) {
499 ERReturn(-1);
500 }
501
502 int total_num_dof = 0;
503 std::vector<int>& masterFieldIDs = cr.getMasterFieldIDs();
504 for(size_t k=0; k<masterFieldIDs.size(); ++k) {
505 total_num_dof += problemStructure_->getFieldSize(masterFieldIDs[k]);
506 }
507
508 dof_ids.resize(total_num_dof);
509 int doffset = 0;
510 for(size_t k=0; k<masterFieldIDs.size(); ++k) {
511 int field_size = problemStructure_->getFieldSize(masterFieldIDs[k]);
512 for(int dof_offset=0; dof_offset<field_size; ++dof_offset) {
513 dof_ids[doffset++] = fdmap.get_dof_id(masterFieldIDs[k], dof_offset);
514 }
515 }
516
517 int blockNum = numElemBlocks + index;
518
519 nodeNumbers.resize(nNodes);
520
521 for(int k=0; k<nNodes; ++k) {
522 const NodeDescriptor* node = Filter::findNode(nodeIDs[k]);
523 if(node == NULL)
524 {
525 nodeNumbers[k] = -1;
526 }
527 else
528 {
529 nodeNumbers[k] = node->getNodeNumber();
530 }
531 }
532
533 int offset = constraintNodeOffsets_[index];
534 CHK_ERR( feData_->setConnectivity(blockNum, numRegularElems_+i++, nNodes, &nodeNumbers[0], &packedFieldSizes_[offset], &dof_ids[0]) );
535 ++cr_iter;
536 }
537
538 }
539 catch(std::runtime_error& exc) {
540 fei::console_out() << exc.what() << FEI_ENDL;
541 ERReturn(-1);
542 }
543
544 return(FEI_SUCCESS);
545}
546
547//------------------------------------------------------------------------------
548int FEDataFilter::resetSystem(double s)
549{
550 //
551 // This puts the value s throughout both the matrix and the vector.
552 //
553 if (Filter::logStream() != NULL) {
554 (*logStream()) << "FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
555 }
556
557 CHK_ERR( feData_->reset() );
558
559 debugOutput("#FEDataFilter leaving resetSystem");
560
561 return(FEI_SUCCESS);
562}
563
564//------------------------------------------------------------------------------
565int FEDataFilter::deleteMultCRs()
566{
567
568 debugOutput("#FEDataFilter::deleteMultCRs");
569
570 int err = feData_->deleteConstraints();
571
572 debugOutput("#FEDataFilter leaving deleteMultCRs");
573
574 return(err);
575}
576
577//------------------------------------------------------------------------------
578int FEDataFilter::resetTheMatrix(double s)
579{
580 //FiniteElementData implementations can't currently reset the matrix without
581 //resetting the rhs vector too.
582 return(FEI_SUCCESS);
583}
584
585//------------------------------------------------------------------------------
586int FEDataFilter::resetTheRHSVector(double s)
587{
588 //FiniteElementData implementations can't currently reset the rhs vector
589 //without resetting the matrix too.
590 return(FEI_SUCCESS);
591}
592
593//------------------------------------------------------------------------------
594int FEDataFilter::resetMatrix(double s)
595{
596 //
597 // This puts the value s throughout both the matrix and the vector.
598 //
599
600 debugOutput("FEI: resetMatrix");
601
602 CHK_ERR( resetTheMatrix(s) );
603
604 eqnCommMgr_->resetCoefs();
605
606 debugOutput("#FEDataFilter leaving resetMatrix");
607
608 return(FEI_SUCCESS);
609}
610
611//------------------------------------------------------------------------------
612int FEDataFilter::resetRHSVector(double s)
613{
614 //
615 // This puts the value s throughout the rhs vector.
616 //
617
618 debugOutput("FEI: resetRHSVector");
619
620 CHK_ERR( resetTheRHSVector(s) );
621
622 eqnCommMgr_->resetCoefs();
623
624 debugOutput("#FEDataFilter leaving resetRHSVector");
625
626 return(FEI_SUCCESS);
627}
628
629//------------------------------------------------------------------------------
630int FEDataFilter::resetInitialGuess(double s)
631{
632 //
633 // This puts the value s throughout the initial guess (solution) vector.
634 //
635 if (Filter::logStream() != NULL) {
636 FEI_OSTREAM& os = *logStream();
637 os << "FEI: resetInitialGuess" << FEI_ENDL;
638 os << "#value to which initial guess is to be set" << FEI_ENDL;
639 os << s << FEI_ENDL;
640 }
641
642 //Actually, the FiniteElementData doesn't currently allow us to alter
643 //values in any initial guess or solution vector.
644
645 debugOutput("#FEDataFilter leaving resetInitialGuess");
646
647 return(FEI_SUCCESS);
648}
649
650//------------------------------------------------------------------------------
651int FEDataFilter::loadNodeBCs(int numNodes,
652 const GlobalID *nodeIDs,
653 int fieldID,
654 const int* offsetsIntoField,
655 const double* prescribedValues)
656{
657 //
658 // load boundary condition information for a given set of nodes
659 //
660 int size = problemStructure_->getFieldSize(fieldID);
661 if (size < 1) {
662 fei::console_out() << "FEI Warning: loadNodeBCs called for fieldID "<<fieldID
663 <<", which was defined with size "<<size<<" (should be positive)."<<FEI_ENDL;
664 return(0);
665 }
666
667 if (Filter::logStream() != NULL) {
668 (*logStream())<<"FEI: loadNodeBCs"<<FEI_ENDL
669 <<"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
670 <<"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
671 <<"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
672 (*logStream())<<"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
673
674 for(int j=0; j<numNodes; j++) {
675 int nodeID = nodeIDs[j];
676 (*logStream())<<nodeID<<" ";
677 (*logStream())<< offsetsIntoField[j]<<" ";
678 (*logStream())<< prescribedValues[j]<<FEI_ENDL;
679 }
680 }
681
682 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
683
684 std::vector<int> essEqns(numNodes);
685 std::vector<double> alpha(numNodes);
686 std::vector<double> gamma(numNodes);
687
688 for(int i=0; i<numNodes; ++i) {
689 NodeDescriptor* node = NULL;
690 nodeDB.getNodeWithID(nodeIDs[i], node);
691 if (node == NULL) {
692 fei::console_out() << "fei_FEDataFilter::loadNodeBCs ERROR, node " << nodeIDs[i]
693 << " not found." << FEI_ENDL;
694 ERReturn(-1);
695 }
696
697 int eqn = -1;
698 if (!node->getFieldEqnNumber(fieldID, eqn)) {
699 ERReturn(-1);
700 }
701
702 essEqns[i] = eqn + offsetsIntoField[i];
703 gamma[i] = prescribedValues[i];
704 alpha[i] = 1.0;
705 }
706
707 if (essEqns.size() > 0) {
708 CHK_ERR( enforceEssentialBCs(&essEqns[0], &alpha[0],
709 &gamma[0], essEqns.size()) );
710 }
711
712 return(FEI_SUCCESS);
713}
714
715//------------------------------------------------------------------------------
716int FEDataFilter::loadElemBCs(int numElems,
717 const GlobalID *elemIDs,
718 int fieldID,
719 const double *const *alpha,
720 const double *const *beta,
721 const double *const *gamma)
722{
723 return(-1);
724}
725
726//------------------------------------------------------------------------------
727void FEDataFilter::allocElemStuff()
728{
729 int nb = problemStructure_->getNumElemBlocks();
730
731 for(int i=0; i<nb; i++) {
732 BlockDescriptor* block = NULL;
733 int err = problemStructure_->getBlockDescriptor_index(i, block);
734 if (err) voidERReturn;
735
736 int numEqns = block->getNumEqnsPerElement();
737 if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
738 }
739
740 eStiff_ = new double*[maxElemRows_];
741 eStiff1D_ = new double[maxElemRows_*maxElemRows_];
742
743 if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
744
745 for(int r=0; r<maxElemRows_; r++) {
746 eStiff_[r] = eStiff1D_ + r*maxElemRows_;
747 }
748
749 eLoad_ = new double[maxElemRows_];
750
751 if (eLoad_ == NULL) voidERReturn
752}
753
754//------------------------------------------------------------------------------
755int FEDataFilter::sumInElem(GlobalID elemBlockID,
756 GlobalID elemID,
757 const GlobalID* elemConn,
758 const double* const* elemStiffness,
759 const double* elemLoad,
760 int elemFormat)
761{
762 if (Filter::logStream() != NULL) {
763 (*logStream()) << "FEI: sumInElem" << FEI_ENDL <<"# elemBlockID " << FEI_ENDL
764 << static_cast<int>(elemBlockID) << FEI_ENDL
765 << "# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
766 BlockDescriptor* block = NULL;
767 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
768 int numNodes = block->getNumNodesPerElement();
769 (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
770 (*logStream()) << "#connected nodes" << FEI_ENDL;
771 for(int i=0; i<numNodes; ++i) {
772 GlobalID nodeID = elemConn[i];
773 (*logStream())<<static_cast<int>(nodeID)<<" ";
774 }
775 (*logStream())<<FEI_ENDL;
776 }
777
778 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
779 elemLoad, elemFormat));
780}
781
782//------------------------------------------------------------------------------
783int FEDataFilter::sumInElemMatrix(GlobalID elemBlockID,
784 GlobalID elemID,
785 const GlobalID* elemConn,
786 const double* const* elemStiffness,
787 int elemFormat)
788{
789 if (Filter::logStream() != NULL) {
790 (*logStream()) << "FEI: sumInElemMatrix"<<FEI_ENDL
791 << "#elemBlockID" << FEI_ENDL << static_cast<int>(elemBlockID)
792 << "# elemID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
793 BlockDescriptor* block = NULL;
794 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
795 int numNodes = block->getNumNodesPerElement();
796 (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
797 (*logStream()) << "#connected nodes" << FEI_ENDL;
798 for(int i=0; i<numNodes; ++i) {
799 GlobalID nodeID = elemConn[i];
800 (*logStream())<<static_cast<int>(nodeID)<<" ";
801 }
802 (*logStream())<<FEI_ENDL;
803 }
804
805 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
806 NULL, elemFormat));
807}
808
809//------------------------------------------------------------------------------
810int FEDataFilter::sumInElemRHS(GlobalID elemBlockID,
811 GlobalID elemID,
812 const GlobalID* elemConn,
813 const double* elemLoad)
814{
815 if (Filter::logStream() != NULL) {
816 (*logStream()) << "FEI: sumInElemRHS"<<FEI_ENDL<<"# elemBlockID " << FEI_ENDL
817 <<static_cast<int>(elemBlockID)
818 << "# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
819 BlockDescriptor* block = NULL;
820 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
821 int numNodes = block->getNumNodesPerElement();
822 (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
823 (*logStream()) << "#connected nodes" << FEI_ENDL;
824 for(int i=0; i<numNodes; ++i) {
825 GlobalID nodeID = elemConn[i];
826 (*logStream())<<static_cast<int>(nodeID)<<" ";
827 }
828 (*logStream())<<FEI_ENDL;
829 }
830
831 return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
832 elemLoad, -1));
833}
834
835//------------------------------------------------------------------------------
836int FEDataFilter::generalElemInput(GlobalID elemBlockID,
837 GlobalID elemID,
838 const GlobalID* elemConn,
839 const double* const* elemStiffness,
840 const double* elemLoad,
841 int elemFormat)
842{
843 (void)elemConn;
844 return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
845 elemFormat) );
846}
847
848//------------------------------------------------------------------------------
849int FEDataFilter::generalElemInput(GlobalID elemBlockID,
850 GlobalID elemID,
851 const double* const* elemStiffness,
852 const double* elemLoad,
853 int elemFormat)
854{
855 //first get the block-descriptor for this elemBlockID...
856
857 BlockDescriptor* block = NULL;
858 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
859
860 //now allocate our local stiffness/load copy if we haven't already.
861
862 if (maxElemRows_ <= 0) allocElemStuff();
863
864 int numElemRows = block->getNumEqnsPerElement();
865
866 //an std::vector.resize operation is free if the size is either shrinking or
867 //staying the same.
868
869 const double* const* stiff = NULL;
870 if (elemStiffness != NULL) stiff = elemStiffness;
871
872 const double* load = NULL;
873 if (elemLoad != NULL) load = elemLoad;
874
875 //we'll make a local dense copy of the element stiffness array
876 //if the stiffness array was passed in using one of the "weird"
877 //element formats, AND if the stiffness array is non-null.
878 if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
879 Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
880 stiff = eStiff_;
881 }
882
883 if (stiff != NULL || load != NULL) newData_ = true;
884
885 if (Filter::logStream() != NULL) {
886 if (stiff != NULL) {
887 (*logStream())
888 << "#numElemRows"<< FEI_ENDL << numElemRows << FEI_ENDL
889 << "#elem-stiff (after being copied into dense-row format)"
890 << FEI_ENDL;
891 for(int i=0; i<numElemRows; i++) {
892 const double* stiff_i = stiff[i];
893 for(int j=0; j<numElemRows; j++) {
894 (*logStream()) << stiff_i[j] << " ";
895 }
896 (*logStream()) << FEI_ENDL;
897 }
898 }
899
900 if (load != NULL) {
901 (*logStream()) << "#elem-load" << FEI_ENDL;
902 for(int i=0; i<numElemRows; i++) {
903 (*logStream()) << load[i] << " ";
904 }
905 (*logStream())<<FEI_ENDL;
906 }
907 }
908
909 //Now we'll proceed to gather the stuff we need to pass the stiffness
910 //data through to the FiniteElementData interface...
911
912 int blockNumber = problemStructure_->getIndexOfBlock(elemBlockID);
913
914 ConnectivityTable& connTable = problemStructure_->
915 getBlockConnectivity(elemBlockID);
916
917 std::map<GlobalID,int>::iterator
918 iter = connTable.elemIDs.find(elemID);
919 if (iter == connTable.elemIDs.end()) {
920 ERReturn(-1);
921 }
922
923 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
924
925 int elemIndex = iter->second;
926
927 int elemNumber = connTable.elemNumbers[elemIndex];
928
929 int numNodes = block->getNumNodesPerElement();
930 int* fieldsPerNode = block->fieldsPerNodePtr();
931 int** fieldIDsTable = block->fieldIDsTablePtr();
932
933 int numDistinctFields = block->getNumDistinctFields();
934 int dof_id = 0;
935 int fieldSize = 0;
936 int total_num_dofs = 0;
937 if (numDistinctFields == 1) {
938 fieldSize = problemStructure_->getFieldSize(fieldIDsTable[0][0]);
939 for(int i=0; i<numNodes; ++i) {
940 total_num_dofs += fieldSize*fieldsPerNode[i];
941 }
942 dof_id = fdmap.get_dof_id(fieldIDsTable[0][0], 0);
943 }
944 else {
945 for(int i=0; i<numNodes; ++i) {
946 for(int nf=0; nf<fieldsPerNode[i]; ++nf) {
947 total_num_dofs += problemStructure_->getFieldSize(fieldIDsTable[i][nf]);
948 }
949 }
950 }
951
952 static std::vector<int> iwork;
953 iwork.resize(2*numNodes+total_num_dofs);
954
955 int* dofsPerNode = &iwork[0];
956 int* nodeNumbers = dofsPerNode+numNodes;
957 int* dof_ids = nodeNumbers+numNodes;
958
959 for(int i=0; i<numNodes; ++i) {
960 dofsPerNode[i] = 0;
961 }
962
963
964 NodeDescriptor** elemNodes =
965 &((*connTable.elem_conn_ptrs)[elemIndex*numNodes]);
966
967 int doffset = 0;
968 for(int nn=0; nn<numNodes; nn++) {
969 NodeDescriptor* node = elemNodes[nn];
970 nodeNumbers[nn] = node->getNodeNumber();
971
972 if (numDistinctFields == 1) {
973 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
974 dofsPerNode[nn] += fieldSize;
975 for(int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
976 dof_ids[doffset++] = dof_id;
977 }
978 }
979 }
980 else {
981 for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
982 int fieldSize = problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
983 int dof_id = fdmap.get_dof_id(fieldIDsTable[nn][nf], 0);
984 dofsPerNode[nn] += fieldSize;
985 for(int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
986 dof_ids[doffset++] = dof_id + dof_offset;
987 }
988 }
989 }
990 }
991
992 if (stiff != NULL) {
993 CHK_ERR( feData_->setElemMatrix(blockNumber, elemNumber, numNodes,
994 nodeNumbers, dofsPerNode, dof_ids, stiff) );
995 }
996
997 if (load != NULL) {
998 CHK_ERR( feData_->setElemVector(blockNumber, elemNumber, numNodes,
999 nodeNumbers, dofsPerNode, dof_ids, load) );
1000 }
1001
1002 return(FEI_SUCCESS);
1003}
1004
1005//------------------------------------------------------------------------------
1006int FEDataFilter::putIntoRHS(int IDType,
1007 int fieldID,
1008 int numIDs,
1009 const GlobalID* IDs,
1010 const double* rhsEntries)
1011{
1012 int fieldSize = problemStructure_->getFieldSize(fieldID);
1013
1014 rowIndices_.resize(fieldSize*numIDs);
1015 int checkNumEqns;
1016
1017 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1018 checkNumEqns,
1019 &rowIndices_[0]));
1020 if (checkNumEqns != numIDs*fieldSize) {
1021 ERReturn(-1);
1022 }
1023
1024 CHK_ERR( exchangeRemoteEquations() );
1025
1026 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
1027
1028 return(0);
1029}
1030
1031//------------------------------------------------------------------------------
1032int FEDataFilter::sumIntoRHS(int IDType,
1033 int fieldID,
1034 int numIDs,
1035 const GlobalID* IDs,
1036 const double* rhsEntries)
1037{
1038 int fieldSize = problemStructure_->getFieldSize(fieldID);
1039
1040 rowIndices_.resize(fieldSize*numIDs);
1041 int checkNumEqns;
1042
1043 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1044 checkNumEqns,
1045 &rowIndices_[0]));
1046 if (checkNumEqns != numIDs*fieldSize) {
1047 ERReturn(-1);
1048 }
1049
1050 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
1051
1052 return(0);
1053}
1054
1055//------------------------------------------------------------------------------
1056int FEDataFilter::sumIntoMatrixDiagonal(int IDType,
1057 int fieldID,
1058 int numIDs,
1059 const GlobalID* IDs,
1060 const double* coefficients)
1061{
1062 int fieldSize = problemStructure_->getFieldSize(fieldID);
1063 const NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1064
1065 std::vector<int> eqns;
1066 convert_field_and_nodes_to_eqns(nodeDB, fieldID, fieldSize, numIDs, IDs, eqns);
1067
1068 std::vector<int> nodeNumbers, dof_ids;
1069 convert_eqns_to_nodenumbers_and_dof_ids(problemStructure_->getFieldDofMap(),
1070 nodeDB, eqns.size(), &eqns[0],
1071 nodeNumbers, dof_ids);
1072
1073 std::vector<int> ones(nodeNumbers.size(), 1);
1074
1075 CHK_ERR( feData_->sumIntoMatrix(nodeNumbers.size(), &nodeNumbers[0], &dof_ids[0],
1076 &ones[0], &nodeNumbers[0], &dof_ids[0], coefficients));
1077 return( 0 );
1078}
1079
1080//------------------------------------------------------------------------------
1081int FEDataFilter::enforceEssentialBCs(const int* eqns,
1082 const double* alpha,
1083 const double* gamma,
1084 int numEqns)
1085{
1086 std::vector<double> values;
1087 std::vector<int> nodeNumbers;
1088 std::vector<int> dof_ids;
1089 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1090
1091 for(int i=0; i<numEqns; i++) {
1092 int reducedEqn = -1;
1093 bool isSlave = problemStructure_->
1094 translateToReducedEqn(eqns[i], reducedEqn);
1095 if (isSlave) continue;
1096
1097 int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqns[i]);
1098
1099 nodeNumbers.push_back(nodeNumber);
1100
1101 const NodeDescriptor* node = NULL;
1102 CHK_ERR( problemStructure_->getNodeDatabase().
1103 getNodeWithNumber(nodeNumber, node));
1104
1105 int fieldID, offset;
1106 node->getFieldID(eqns[i], fieldID, offset);
1107 dof_ids.push_back( fdmap.get_dof_id(fieldID, offset) );
1108
1109 values.push_back(gamma[i]/alpha[i]);
1110 }
1111
1112 CHK_ERR( feData_->setDirichletBCs(nodeNumbers.size(),
1113 &nodeNumbers[0], &dof_ids[0], &values[0]) );
1114
1115 newData_ = true;
1116
1117 return(FEI_SUCCESS);
1118}
1119
1120//------------------------------------------------------------------------------
1121int FEDataFilter::loadFEDataMultCR(int CRID,
1122 int numCRNodes,
1123 const GlobalID* CRNodes,
1124 const int* CRFields,
1125 const double* CRWeights,
1126 double CRValue)
1127{
1128 if (Filter::logStream() != NULL) {
1129 FEI_OSTREAM& os = *logStream();
1130 os<<"FEI: loadCRMult"<<FEI_ENDL;
1131 os<<"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
1132 os<<"#CRNodes:"<<FEI_ENDL;
1133 int i;
1134 for(i=0; i<numCRNodes; ++i) {
1135 GlobalID nodeID = CRNodes[i];
1136 os << static_cast<int>(nodeID) << " ";
1137 }
1138 os << FEI_ENDL << "#fields:"<<FEI_ENDL;
1139 for(i=0; i<numCRNodes; ++i) os << CRFields[i] << " ";
1140 os << FEI_ENDL << "#field-sizes:"<<FEI_ENDL;
1141 for(i=0; i<numCRNodes; ++i) {
1142 int size = problemStructure_->getFieldSize(CRFields[i]);
1143 os << size << " ";
1144 }
1145 os << FEI_ENDL<<"#weights:"<<FEI_ENDL;
1146 int offset = 0;
1147 for(i=0; i<numCRNodes; ++i) {
1148 int size = problemStructure_->getFieldSize(CRFields[i]);
1149 for(int j=0; j<size; ++j) {
1150 os << CRWeights[offset++] << " ";
1151 }
1152 }
1153 os << FEI_ENDL<<"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
1154 }
1155
1156 if (numCRNodes <= 0) return(0);
1157
1158 std::vector<int> nodeNumbers;
1159 std::vector<int> dof_ids;
1160 std::vector<double> weights;
1161
1162 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1163 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1164
1165 double fei_min = std::numeric_limits<double>::min();
1166
1167 int offset = 0;
1168 for(int i=0; i<numCRNodes; i++) {
1169 NodeDescriptor* node = NULL;
1170 CHK_ERR( nodeDB.getNodeWithID(CRNodes[i], node) );
1171
1172 int fieldEqn = -1;
1173 bool hasField = node->getFieldEqnNumber(CRFields[i], fieldEqn);
1174 if (!hasField) ERReturn(-1);
1175
1176 int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
1177 int dof_id = fdmap.get_dof_id(CRFields[i], 0);
1178
1179 for(int f=0; f<fieldSize; f++) {
1180 double weight = CRWeights[offset++];
1181 if (std::abs(weight) > fei_min) {
1182 nodeNumbers.push_back(node->getNodeNumber());
1183 dof_ids.push_back(dof_id+f);
1184 weights.push_back(weight);
1185 }
1186 }
1187 }
1188
1189 CHK_ERR( feData_->setMultiplierCR(CRID, nodeNumbers.size(),
1190 &nodeNumbers[0], &dof_ids[0],
1191 &weights[0], CRValue) );
1192 newData_ = true;
1193
1194 return(0);
1195}
1196
1197//------------------------------------------------------------------------------
1198int FEDataFilter::loadFEDataPenCR(int CRID,
1199 int numCRNodes,
1200 const GlobalID* CRNodes,
1201 const int* CRFields,
1202 const double* CRWeights,
1203 double CRValue,
1204 double penValue)
1205{
1206 if (numCRNodes <= 0) return(0);
1207
1208 std::vector<int> nodeNumbers;
1209 std::vector<int> dofsPerNode;
1210 std::vector<int> dof_ids;
1211 std::vector<double> weights;
1212
1213 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1214 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1215
1216 int offset = 0;
1217 for(int i=0; i<numCRNodes; i++) {
1218 NodeDescriptor* node = NULL;
1219 nodeDB.getNodeWithID(CRNodes[i], node);
1220 if(node == NULL) continue;
1221
1222 int fieldEqn = -1;
1223 bool hasField = node->getFieldEqnNumber(CRFields[i], fieldEqn);
1224 // If a node doesn't have a field, skip it.
1225 if (!hasField) continue;
1226
1227 int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
1228
1229 nodeNumbers.push_back(node->getNodeNumber());
1230 dofsPerNode.push_back(fieldSize);
1231
1232 for(int f=0; f<fieldSize; f++) {
1233 dof_ids.push_back(fdmap.get_dof_id(CRFields[i], f));
1234 double weight = CRWeights[offset++];
1235 weights.push_back(weight);
1236 }
1237 }
1238
1239 std::vector<double*> matrixCoefs(weights.size());
1240 std::vector<double> rhsCoefs(weights.size());
1241 offset = 0;
1242 for(size_t i=0; i<weights.size(); ++i) {
1243 double* coefPtr = new double[weights.size()];
1244 for(size_t j=0; j<weights.size(); ++j) {
1245 coefPtr[j] = weights[i]*weights[j]*penValue;
1246 }
1247 matrixCoefs[i] = coefPtr;
1248 rhsCoefs[i] = weights[i]*penValue*CRValue;
1249 }
1250
1251 int crIndex = fei::binarySearch(CRID, penCRIDs_);
1252
1253 int index = fei::binarySearch(numCRNodes, constraintBlocks_);
1254
1255 int blockNum = problemStructure_->getNumElemBlocks() + index;
1256 int elemNum = numRegularElems_ + crIndex;
1257
1258 CHK_ERR( feData_->setElemMatrix(blockNum, elemNum,
1259 nodeNumbers.size(),
1260 &nodeNumbers[0],
1261 &dofsPerNode[0],
1262 &dof_ids[0],
1263 &matrixCoefs[0]) );
1264
1265 CHK_ERR( feData_->setElemVector(blockNum, elemNum, nodeNumbers.size(),
1266 &nodeNumbers[0], &dofsPerNode[0], &dof_ids[0], &rhsCoefs[0]) );
1267
1268 newData_ = true;
1269
1270 for(size_t i=0; i<weights.size(); ++i) {
1271 delete [] matrixCoefs[i];
1272 }
1273
1274 return(0);
1275}
1276
1277//------------------------------------------------------------------------------
1278int FEDataFilter::loadCRMult(int CRID,
1279 int numCRNodes,
1280 const GlobalID* CRNodes,
1281 const int* CRFields,
1282 const double* CRWeights,
1283 double CRValue)
1284{
1285//
1286// Load Lagrange multiplier constraint relation data
1287//
1288
1289 //Give the constraint data to the underlying solver using this special function...
1290 CHK_ERR( loadFEDataMultCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
1291 CRValue) );
1292
1293 return(FEI_SUCCESS);
1294}
1295
1296//------------------------------------------------------------------------------
1297int FEDataFilter::loadCRPen(int CRID,
1298 int numCRNodes,
1299 const GlobalID* CRNodes,
1300 const int* CRFields,
1301 const double* CRWeights,
1302 double CRValue,
1303 double penValue)
1304{
1305//
1306// Load penalty constraint relation data
1307//
1308
1309 debugOutput("FEI: loadCRPen");
1310
1311 //Give the constraint data to the underlying solver using this special function...
1312 CHK_ERR( loadFEDataPenCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
1313 CRValue, penValue) );
1314
1315 return(FEI_SUCCESS);
1316}
1317
1318//------------------------------------------------------------------------------
1319int FEDataFilter::parameters(int numParams, const char *const* paramStrings)
1320{
1321//
1322// this function takes parameters for setting internal things like solver
1323// and preconditioner choice, etc.
1324//
1325 if (numParams == 0 || paramStrings == NULL) {
1326 debugOutput("#FEDataFilter::parameters --- no parameters.");
1327 }
1328 else {
1329
1330 snl_fei::getIntParamValue("outputLevel",numParams, paramStrings,
1331 outputLevel_);
1332
1333 snl_fei::getIntParamValue("internalFei",numParams,paramStrings,
1334 internalFei_);
1335
1336 if (Filter::logStream() != NULL) {
1337 (*logStream())<<"#FEDataFilter::parameters"<<FEI_ENDL
1338 <<"# --- numParams: "<< numParams<<FEI_ENDL;
1339 for(int i=0; i<numParams; i++){
1340 (*logStream())<<"#------ paramStrings["<<i<<"]: "
1341 <<paramStrings[i]<<FEI_ENDL;
1342 }
1343 }
1344 }
1345
1346 CHK_ERR( Filter::parameters(numParams, paramStrings) );
1347
1348 debugOutput("#FEDataFilter leaving parameters function");
1349
1350 return(FEI_SUCCESS);
1351}
1352
1353//------------------------------------------------------------------------------
1354int FEDataFilter::loadComplete()
1355{
1356 debugOutput("#FEDataFilter calling FEData matrixLoadComplete");
1357
1358 CHK_ERR( feData_->loadComplete() );
1359
1360 newData_ = false;
1361
1362 return(0);
1363}
1364
1365//------------------------------------------------------------------------------
1366int FEDataFilter::residualNorm(int whichNorm, int numFields,
1367 int* fieldIDs, double* norms, double& residTime)
1368{
1369//
1370//This function can do 3 kinds of norms: infinity-norm (denoted
1371//by whichNorm==0), 1-norm and 2-norm.
1372//
1373 debugOutput("FEI: residualNorm");
1374
1375 CHK_ERR( loadComplete() );
1376
1377 //for now, FiniteElementData doesn't do residual calculations.
1378
1379 int fdbNumFields = problemStructure_->getNumFields();
1380 const int* fdbFieldIDs = problemStructure_->getFieldIDsPtr();
1381
1382 int i;
1383
1384 //Since we don't calculate actual residual norms, we'll fill the user's
1385 //array with norm data that is obviously not real norm data.
1386 int offset = 0;
1387 i = 0;
1388 while(offset < numFields && i < fdbNumFields) {
1389 if (fdbFieldIDs[i] >= 0) {
1390 fieldIDs[offset++] = fdbFieldIDs[i];
1391 }
1392 ++i;
1393 }
1394 for(i=0; i<numFields; ++i) {
1395 norms[i] = -99.9;
1396 }
1397
1398 //fill out the end of the array with garbage in case the user-provided
1399 //array is longer than the list of fields we have in fieldDB.
1400 for(i=offset; i<numFields; ++i) {
1401 fieldIDs[i] = -99;
1402 }
1403
1404 return(FEI_SUCCESS);
1405}
1406
1407//------------------------------------------------------------------------------
1408int FEDataFilter::formResidual(double* residValues, int numLocalEqns)
1409{
1410 //FiniteElementData implementations can't currently do residuals.
1411 return(FEI_SUCCESS);
1412}
1413
1414//------------------------------------------------------------------------------
1415int FEDataFilter::solve(int& status, double& sTime) {
1416
1417 debugOutput("FEI: solve");
1418
1419 CHK_ERR( loadComplete() );
1420
1421 debugOutput("#FEDataFilter in solve, calling launchSolver...");
1422
1423 double start = MPI_Wtime();
1424
1425 CHK_ERR( feData_->launchSolver(status, iterations_) );
1426
1427 sTime = MPI_Wtime() - start;
1428
1429 debugOutput("#FEDataFilter... back from solver");
1430
1431 //now unpack the locally-owned shared entries of the solution vector into
1432 //the eqn-comm-mgr data structures.
1433 CHK_ERR( unpackSolution() );
1434
1435 debugOutput("#FEDataFilter leaving solve");
1436
1437 if (status != 0) return(1);
1438 else return(FEI_SUCCESS);
1439}
1440
1441//------------------------------------------------------------------------------
1442int FEDataFilter::setNumRHSVectors(int numRHSs, int* rhsIDs){
1443
1444 if (numRHSs < 0) {
1445 fei::console_out() << "FEDataFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
1446 ERReturn(-1);
1447 }
1448
1449 numRHSs_ = numRHSs;
1450
1451 rhsIDs_.resize(numRHSs_);
1452 for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
1453
1454 //(we need to set the number of RHSs in the eqn comm manager)
1455 eqnCommMgr_->setNumRHSs(numRHSs_);
1456
1457 return(FEI_SUCCESS);
1458}
1459
1460//------------------------------------------------------------------------------
1461int FEDataFilter::setCurrentRHS(int rhsID)
1462{
1463 std::vector<int>::iterator iter =
1464 std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
1465
1466 if (iter == rhsIDs_.end()) ERReturn(-1)
1467
1468 currentRHS_ = iter - rhsIDs_.begin();
1469
1470 return(FEI_SUCCESS);
1471}
1472
1473//------------------------------------------------------------------------------
1474int FEDataFilter::giveToMatrix(int numPtRows, const int* ptRows,
1475 int numPtCols, const int* ptCols,
1476 const double* const* values, int mode)
1477{
1478 //This isn't going to be fast... I need to optimize the whole structure
1479 //of code that's associated with passing data to FiniteElementData.
1480
1481 std::vector<int> rowNodeNumbers, row_dof_ids, colNodeNumbers, col_dof_ids;
1482 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1483 int i;
1484
1485 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1486
1487 //First, we have to get nodeNumbers and dof_ids for each of the
1488 //row-numbers and col-numbers.
1489
1490 for(i=0; i<numPtRows; i++) {
1491 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
1492 if (nodeNumber < 0) ERReturn(-1);
1493 const NodeDescriptor* node = NULL;
1494 CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
1495 int fieldID, offset;
1496 node->getFieldID(ptRows[i], fieldID, offset);
1497
1498 rowNodeNumbers.push_back(nodeNumber);
1499 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1500 }
1501
1502 for(i=0; i<numPtCols; i++) {
1503 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
1504 if (nodeNumber < 0) ERReturn(-1);
1505 const NodeDescriptor* node = NULL;
1506 CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
1507 int fieldID, offset;
1508 node->getFieldID(ptCols[i], fieldID, offset);
1509
1510 colNodeNumbers.push_back(nodeNumber);
1511 col_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1512 }
1513
1514 //now we have to flatten the colNodeNumbers and col_dof_ids out into
1515 //an array of length numPtRows*numPtCols, where the nodeNumbers and
1516 //dof_ids are repeated 'numPtRows' times.
1517
1518 int len = numPtRows*numPtCols;
1519 std::vector<int> allColNodeNumbers(len), all_col_dof_ids(len);
1520 std::vector<double> allValues(len);
1521
1522 int offset = 0;
1523 for(i=0; i<numPtRows; i++) {
1524 for(int j=0; j<numPtCols; j++) {
1525 allColNodeNumbers[offset] = colNodeNumbers[j];
1526 all_col_dof_ids[offset] = col_dof_ids[j];
1527 allValues[offset++] = values[i][j];
1528 }
1529 }
1530
1531 //while we're at it, let's make an array with numPtCols replicated in it
1532 //'numPtRows' times.
1533 std::vector<int> numColsPerRow(numPtRows, numPtCols);
1534
1535 //now we're ready to hand this stuff off to the FiniteElementData
1536 //instantiation.
1537
1538 CHK_ERR( feData_->sumIntoMatrix(numPtRows,
1539 &rowNodeNumbers[0],
1540 &row_dof_ids[0],
1541 &numColsPerRow[0],
1542 &allColNodeNumbers[0],
1543 &all_col_dof_ids[0],
1544 &allValues[0]) );
1545
1546 return(FEI_SUCCESS);
1547}
1548
1549//------------------------------------------------------------------------------
1550int FEDataFilter::giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
1551 int numPtCols, const int* ptCols,
1552 const double* const* values, int mode)
1553{
1554 //This isn't going to be fast... I need to optimize the whole structure
1555 //of code that's associated with passing data to FiniteElementData.
1556
1557 std::vector<int> rowNodeNumbers, row_dof_ids, colNodeNumbers, col_dof_ids;
1558 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1559 int i;
1560
1561 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1562
1563 //First, we have to get nodeNumbers and dof_ids for each of the
1564 //row-numbers and col-numbers.
1565
1566 for(i=0; i<numPtRows; i++) {
1567 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
1568 if (nodeNumber < 0) ERReturn(-1);
1569 const NodeDescriptor* node = NULL;
1570 CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
1571 int fieldID, offset;
1572 node->getFieldID(ptRows[i], fieldID, offset);
1573
1574 rowNodeNumbers.push_back(nodeNumber);
1575 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1576 }
1577
1578 for(i=0; i<numPtCols; i++) {
1579 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
1580 if (nodeNumber < 0) ERReturn(-1);
1581 const NodeDescriptor* node = NULL;
1582 CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
1583 int fieldID, offset;
1584 node->getFieldID(ptCols[i], fieldID, offset);
1585
1586 colNodeNumbers.push_back(nodeNumber);
1587 col_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1588 }
1589
1590 //now we have to flatten the colNodeNumbers and col_dof_ids out into
1591 //an array of length numPtRows*numPtCols, where the nodeNumbers and
1592 //dof_ids are repeated 'numPtRows' times.
1593
1594 int len = numPtRows*numPtCols;
1595 std::vector<int> allColNodeNumbers(len), all_col_dof_ids(len);
1596 std::vector<double> allValues(len);
1597
1598 int offset = 0;
1599 for(i=0; i<numPtRows; i++) {
1600 for(int j=0; j<numPtCols; j++) {
1601 allColNodeNumbers[offset] = colNodeNumbers[j];
1602 all_col_dof_ids[offset] = col_dof_ids[j];
1603 allValues[offset++] = values[i][j];
1604 }
1605 }
1606
1607 //while we're at it, let's make an array with numPtCols replicated in it
1608 //'numPtRows' times.
1609 std::vector<int> numColsPerRow(numPtRows, numPtCols);
1610
1611 //now we're ready to hand this stuff off to the FiniteElementData
1612 //instantiation.
1613
1614 CHK_ERR( feData_->sumIntoMatrix(numPtRows,
1615 &rowNodeNumbers[0],
1616 &row_dof_ids[0],
1617 &numColsPerRow[0],
1618 &allColNodeNumbers[0],
1619 &all_col_dof_ids[0],
1620 &allValues[0]) );
1621
1622 return(FEI_SUCCESS);
1623}
1624
1625//------------------------------------------------------------------------------
1626int FEDataFilter::getFromMatrix(int numPtRows, const int* ptRows,
1627 const int* rowColOffsets, const int* ptCols,
1628 int numColsPerRow, double** values)
1629{
1630 return(-1);
1631
1632}
1633
1634//------------------------------------------------------------------------------
1635int FEDataFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
1636{
1637 ERReturn(-1);
1638}
1639
1640//------------------------------------------------------------------------------
1641int FEDataFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
1642{
1643 ERReturn(-1);
1644}
1645
1646//------------------------------------------------------------------------------
1647int FEDataFilter::giveToRHS(int num, const double* values,
1648 const int* indices, int mode)
1649{
1650 std::vector<int> workspace(num*2);
1651 int* rowNodeNumbers = &workspace[0];
1652 int* row_dof_ids = rowNodeNumbers+num;
1653 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1654 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1655
1656 for(int i=0; i<num; ++i) {
1657 const NodeDescriptor* nodeptr = 0;
1658 int err = nodeDB.getNodeWithEqn(indices[i], nodeptr);
1659 if (err < 0) {
1660 rowNodeNumbers[i] = -1;
1661 row_dof_ids[i] = -1;
1662 continue;
1663 }
1664
1665 rowNodeNumbers[i] = nodeptr->getNodeNumber();
1666
1667 int fieldID, offset;
1668 nodeptr->getFieldID(indices[i], fieldID, offset);
1669
1670 row_dof_ids[i] = fdmap.get_dof_id(fieldID, offset);
1671 }
1672
1673 if (mode == ASSEMBLE_SUM) {
1674 CHK_ERR( feData_->sumIntoRHSVector(num,
1675 rowNodeNumbers,
1676 row_dof_ids,
1677 values) );
1678 }
1679 else {
1680 CHK_ERR( feData_->putIntoRHSVector(num,
1681 rowNodeNumbers,
1682 row_dof_ids,
1683 values) );
1684 }
1685
1686 return(FEI_SUCCESS);
1687}
1688
1689//------------------------------------------------------------------------------
1690int FEDataFilter::giveToLocalReducedRHS(int num, const double* values,
1691 const int* indices, int mode)
1692{
1693 std::vector<int> rowNodeNumbers, row_dof_ids;
1694 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1695 fei::FieldDofMap<int>& fdmap = problemStructure_->getFieldDofMap();
1696
1697 for(int i=0; i<num; i++) {
1698 int nodeNumber = problemStructure_->getAssociatedNodeNumber(indices[i]);
1699 if (nodeNumber < 0) ERReturn(-1);
1700
1701 const NodeDescriptor* node = NULL;
1702 CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
1703
1704 int fieldID, offset;
1705 node->getFieldID(indices[i], fieldID, offset);
1706
1707 rowNodeNumbers.push_back(nodeNumber);
1708 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1709 }
1710
1711 if (mode == ASSEMBLE_SUM) {
1712 CHK_ERR( feData_->sumIntoRHSVector(rowNodeNumbers.size(),
1713 &rowNodeNumbers[0],
1714 &row_dof_ids[0], values) );
1715 }
1716 else {
1717 CHK_ERR( feData_->putIntoRHSVector(rowNodeNumbers.size(),
1718 &rowNodeNumbers[0],
1719 &row_dof_ids[0], values) );
1720 }
1721
1722 return(FEI_SUCCESS);
1723}
1724
1725//------------------------------------------------------------------------------
1726int FEDataFilter::getFromRHS(int num, double* values, const int* indices)
1727{
1728 return(-1);
1729}
1730
1731//------------------------------------------------------------------------------
1732int FEDataFilter::getEqnSolnEntry(int eqnNumber, double& solnValue)
1733{
1734 //This function's task is to retrieve the solution-value for a global
1735 //equation-number. eqnNumber may or may not be a slave-equation, and may or
1736 //may not be locally owned. If it is not locally owned, it should at least
1737 //be shared.
1738 //return 0 if the solution is successfully retrieved, otherwise return 1.
1739 //
1740
1741 if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
1742 //Dig into the eqn-comm-mgr for the shared-remote solution value.
1743 CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
1744 }
1745 else {
1746 //It's local, simply get the solution from the assembled linear system.
1747 CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
1748 }
1749 return(0);
1750}
1751
1752//------------------------------------------------------------------------------
1753int FEDataFilter::getSharedRemoteSolnEntry(int eqnNumber, double& solnValue)
1754{
1755 std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
1756 double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
1757
1758 int index = fei::binarySearch(eqnNumber, remoteEqnNumbers);
1759 if (index < 0) {
1760 fei::console_out() << "FEDataFilter::getSharedRemoteSolnEntry: ERROR, eqn "
1761 << eqnNumber << " not found." << FEI_ENDL;
1762 ERReturn(-1);
1763 }
1764 solnValue = remoteSoln[index];
1765 return(0);
1766}
1767
1768//------------------------------------------------------------------------------
1769int FEDataFilter::getReducedSolnEntry(int eqnNumber, double& solnValue)
1770{
1771 //We may safely assume that this function is called with 'eqnNumber' that is
1772 //local in the underlying assembled linear system. i.e., it isn't a slave-
1773 //equation, it isn't remotely owned, etc.
1774 //
1775
1776 int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqnNumber);
1777
1778 //if nodeNumber < 0, it probably means we're trying to look up the
1779 //node for a lagrange-multiplier (which doesn't exist). In that
1780 //case, we're just going to ignore the request and return for now...
1781 if (nodeNumber < 0) {solnValue = -999.99; return(FEI_SUCCESS);}
1782
1783 const NodeDescriptor* node = NULL;
1784 problemStructure_->getNodeDatabase().getNodeWithNumber(nodeNumber, node);
1785 if(node == NULL) {
1786 // KHP: If a node doesn't exist, we still need to
1787 // return a solution value....Zero seems like a logical
1788 // choice however, FEI_SUCCESS seems wrong however I don't
1789 // want to trip any asserts or other error conditions.
1790 solnValue = 0.0;
1791 return FEI_SUCCESS;
1792 }
1793
1794 int eqn = problemStructure_->translateFromReducedEqn(eqnNumber);
1795 int fieldID, offset;
1796 node->getFieldID(eqn, fieldID, offset);
1797 int dof_id = problemStructure_->getFieldDofMap().get_dof_id(fieldID, offset);
1798
1799 bool fetiHasNode = true;
1800 GlobalID nodeID = node->getGlobalNodeID();
1801 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
1802 std::vector<GlobalID>& shNodeIDs = nodeCommMgr.getSharedNodeIDs();
1803 int shIndex = fei::binarySearch(nodeID, shNodeIDs);
1804 if (shIndex >= 0) {
1805 if (!(problemStructure_->isInLocalElement(nodeNumber)) ) fetiHasNode = false;
1806 }
1807
1808 if (fetiHasNode) {
1809 int err = feData_->getSolnEntry(nodeNumber, dof_id, solnValue);
1810 if (err != 0) {
1811 fei::console_out() << "FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
1812 << " (nodeID " << node->getGlobalNodeID() << "), dof_id "<<dof_id
1813 << " couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
1814 ERReturn(-1);
1815 }
1816 }
1817
1818 return(FEI_SUCCESS);
1819}
1820
1821//------------------------------------------------------------------------------
1822int FEDataFilter::unpackSolution()
1823{
1824 //
1825 //This function should be called after the solver has returned,
1826 //and we know that there is a solution in the underlying vector.
1827 //This function ensures that any locally-owned shared solution values are
1828 //available on the sharing processors.
1829 //
1830 if (Filter::logStream() != NULL) {
1831 (*logStream())<< "# entering unpackSolution, outputLevel: "
1832 <<outputLevel_<<FEI_ENDL;
1833 }
1834
1835 //what we need to do is as follows.
1836 //The eqn comm mgr has a list of what it calls 'recv eqns'. These are
1837 //equations that we own, for which we received contributions from other
1838 //processors. The solution values corresponding to these equations need
1839 //to be made available to those remote contributing processors.
1840
1841 int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
1842 std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
1843
1844 for(int i=0; i<numRecvEqns; i++) {
1845 int eqn = recvEqnNumbers[i];
1846
1847 if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
1848 fei::console_out() << "FEDataFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
1849 << ") out of local range." << FEI_ENDL;
1850 MPI_Abort(comm_, -1);
1851 }
1852
1853 double solnValue = 0.0;
1854
1855 CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
1856
1857 eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
1858 }
1859
1860 eqnCommMgr_->exchangeSoln();
1861
1862 debugOutput("#FEDataFilter leaving unpackSolution");
1863 return(FEI_SUCCESS);
1864}
1865
1866//------------------------------------------------------------------------------
1867void FEDataFilter:: setEqnCommMgr(EqnCommMgr* eqnCommMgr)
1868{
1869 delete eqnCommMgr_;
1870 eqnCommMgr_ = eqnCommMgr;
1871}
1872
1873//------------------------------------------------------------------------------
1874int FEDataFilter::getBlockNodeSolution(GlobalID elemBlockID,
1875 int numNodes,
1876 const GlobalID *nodeIDs,
1877 int *offsets,
1878 double *results)
1879{
1880 debugOutput("FEI: getBlockNodeSolution");
1881
1882 int numActiveNodes = problemStructure_->getNumActiveNodes();
1883 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1884
1885 if (numActiveNodes <= 0) return(0);
1886
1887 int numSolnParams = 0;
1888
1889 BlockDescriptor* block = NULL;
1890 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1891
1892 //Traverse the node list, checking if nodes are associated with this block.
1893 //If so, put its 'answers' in the results list.
1894
1895 int offset = 0;
1896 for(int i=0; i<numActiveNodes; i++) {
1897 NodeDescriptor* node_i = NULL;
1898 nodeDB.getNodeAtIndex(i, node_i);
1899
1900 if (offset == numNodes) break;
1901
1902 GlobalID nodeID = nodeIDs[offset];
1903
1904 //first let's set the offset at which this node's solution coefs start.
1905 offsets[offset++] = numSolnParams;
1906
1907 NodeDescriptor* node = NULL;
1908 int err = 0;
1909 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
1910 //Don't call the getActiveNodeDesc_ID function unless we have to.
1911
1912 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
1913 node = node_i;
1914 }
1915 else {
1916 err = nodeDB.getNodeWithID(nodeID, node);
1917 }
1918
1919 //ok. If err is not 0, meaning nodeID is NOT in the
1920 //activeNodes list, then skip to the next loop iteration.
1921
1922 if (err != 0) {
1923 continue;
1924 }
1925
1926 int numFields = node->getNumFields();
1927 const int* fieldIDs = node->getFieldIDList();
1928
1929 for(int j=0; j<numFields; j++) {
1930 if (block->containsField(fieldIDs[j])) {
1931 int size = problemStructure_->getFieldSize(fieldIDs[j]);
1932 if (size < 1) {
1933 continue;
1934 }
1935
1936 int thisEqn = -1;
1937 node->getFieldEqnNumber(fieldIDs[j], thisEqn);
1938
1939 double answer;
1940 for(int k=0; k<size; k++) {
1941 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
1942 results[numSolnParams++] = answer;
1943 }
1944 }
1945 }//for(j<numFields)loop
1946 }
1947
1948 offsets[numNodes] = numSolnParams;
1949
1950 return(FEI_SUCCESS);
1951}
1952
1953//------------------------------------------------------------------------------
1954int FEDataFilter::getNodalSolution(int numNodes,
1955 const GlobalID *nodeIDs,
1956 int *offsets,
1957 double *results)
1958{
1959 debugOutput("FEI: getNodalSolution");
1960
1961 int numActiveNodes = problemStructure_->getNumActiveNodes();
1962 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1963
1964 if (numActiveNodes <= 0) return(0);
1965
1966 int numSolnParams = 0;
1967
1968 //Traverse the node list, checking if nodes are local.
1969 //If so, put 'answers' in the results list.
1970
1971 int offset = 0;
1972 for(int i=0; i<numActiveNodes; i++) {
1973 NodeDescriptor* node_i = NULL;
1974 nodeDB.getNodeAtIndex(i, node_i);
1975
1976 if (offset == numNodes) break;
1977
1978 GlobalID nodeID = nodeIDs[offset];
1979
1980 //first let's set the offset at which this node's solution coefs start.
1981 offsets[offset++] = numSolnParams;
1982
1983 NodeDescriptor* node = NULL;
1984 int err = 0;
1985 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
1986 //Don't call the getNodeWithID function unless we have to.
1987
1988 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
1989 node = node_i;
1990 }
1991 else {
1992 err = nodeDB.getNodeWithID(nodeID, node);
1993 }
1994
1995 //ok. If err is not 0, meaning nodeID is NOT in the
1996 //activeNodes list, then skip to the next loop iteration.
1997
1998 if (err != 0) {
1999 continue;
2000 }
2001
2002 int numFields = node->getNumFields();
2003 const int* fieldIDs = node->getFieldIDList();
2004
2005 for(int j=0; j<numFields; j++) {
2006 int size = problemStructure_->getFieldSize(fieldIDs[j]);
2007 if (size < 1) {
2008 continue;
2009 }
2010
2011 int thisEqn = -1;
2012 node->getFieldEqnNumber(fieldIDs[j], thisEqn);
2013
2014 double answer;
2015 for(int k=0; k<size; k++) {
2016 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
2017 results[numSolnParams++] = answer;
2018 }
2019 }//for(j<numFields)loop
2020 }
2021
2022 offsets[numNodes] = numSolnParams;
2023
2024 return(FEI_SUCCESS);
2025}
2026
2027//------------------------------------------------------------------------------
2028int FEDataFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
2029 int fieldID,
2030 int numNodes,
2031 const GlobalID *nodeIDs,
2032 double *results)
2033{
2034 debugOutput("FEI: getBlockFieldNodeSolution");
2035
2036 int numActiveNodes = problemStructure_->getNumActiveNodes();
2037 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2038
2039 if (numActiveNodes <= 0) return(0);
2040
2041 BlockDescriptor* block = NULL;
2042 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2043
2044 int fieldSize = problemStructure_->getFieldSize(fieldID);
2045 if (fieldSize <= 0) ERReturn(-1);
2046
2047 if (!block->containsField(fieldID)) {
2048 fei::console_out() << "FEDataFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
2049 << " not contained in element-block " << static_cast<int>(elemBlockID) << FEI_ENDL;
2050 return(1);
2051 }
2052
2053 //Traverse the node list, checking if nodes are associated with this block.
2054 //If so, put the answers in the results list.
2055
2056 for(int i=0; i<numNodes; i++) {
2057 NodeDescriptor* node_i = NULL;
2058 nodeDB.getNodeAtIndex(i, node_i);
2059
2060 GlobalID nodeID = nodeIDs[i];
2061
2062 NodeDescriptor* node = NULL;
2063 int err = 0;
2064 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
2065 //Don't call the getActiveNodeDesc_ID function unless we have to.
2066
2067 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2068 node = node_i;
2069 }
2070 else {
2071 err = nodeDB.getNodeWithID(nodeID, node);
2072 }
2073
2074 //ok. If err is not 0, meaning nodeID is NOT in the
2075 //activeNodes list, then skip to the next loop iteration.
2076
2077 if (err != 0) {
2078 continue;
2079 }
2080
2081 int eqnNumber = -1;
2082 bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
2083 if (!hasField) continue;
2084
2085 int offset = fieldSize*i;
2086 for(int j=0; j<fieldSize; j++) {
2087 double answer = 0.0;
2088 CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
2089 results[offset+j] = answer;
2090 }
2091 }
2092
2093 return(FEI_SUCCESS);
2094}
2095
2096//------------------------------------------------------------------------------
2097int FEDataFilter::getNodalFieldSolution(int fieldID,
2098 int numNodes,
2099 const GlobalID *nodeIDs,
2100 double *results)
2101{
2102 debugOutput("FEI: getNodalFieldSolution");
2103
2104 int numActiveNodes = problemStructure_->getNumActiveNodes();
2105 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2106
2107 if (numActiveNodes <= 0) return(0);
2108
2109 if (problemStructure_->numSlaveEquations() != 0) {
2110 fei::console_out() << "FEDataFilter::getEqnSolnEntry ERROR FETI-support is not currently"
2111 << " compatible with the FEI's constraint reduction." << FEI_ENDL;
2112 ERReturn(-1);
2113 }
2114
2115 int fieldSize = problemStructure_->getFieldSize(fieldID);
2116 if (fieldSize <= 0) {
2117 ERReturn(-1);
2118 }
2119
2120 //Traverse the node list, checking if nodes have the specified field.
2121 //If so, put the answers in the results list.
2122
2123 for(int i=0; i<numNodes; i++) {
2124 NodeDescriptor* node_i = NULL;
2125 nodeDB.getNodeAtIndex(i, node_i);
2126
2127 GlobalID nodeID = nodeIDs[i];
2128
2129 NodeDescriptor* node = NULL;
2130 int err = 0;
2131 //Obtain the NodeDescriptor of nodeID in the activeNodes list...
2132 //Don't call the getNodeWithID function unless we have to.
2133
2134 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2135 node = node_i;
2136 }
2137 else {
2138 err = nodeDB.getNodeWithID(nodeID, node);
2139 }
2140
2141 //ok. If err is not 0, meaning nodeID is NOT in the
2142 //activeNodes list, then skip to the next loop iteration.
2143
2144 if (err != 0) {
2145 continue;
2146 }
2147
2148 int nodeNumber = node->getNodeNumber();
2149
2150 int eqnNumber = -1;
2151 bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
2152
2153 //If this node doesn't have the specified field, then skip to the
2154 //next loop iteration.
2155 if (!hasField) continue;
2156
2157 int dof_id = problemStructure_->getFieldDofMap().get_dof_id(fieldID, 0);
2158
2159 int offset = fieldSize*i;
2160
2161 for(int j=0; j<fieldSize; j++) {
2162 if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
2163 CHK_ERR( getSharedRemoteSolnEntry(eqnNumber+j, results[offset+j]) );
2164 continue;
2165 }
2166
2167 err = feData_->getSolnEntry(nodeNumber, dof_id+j, results[offset+j]);
2168 if (err != 0) {
2169 fei::console_out() << "FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
2170 << " (nodeID " << nodeID << "), dof_id "<<dof_id
2171 << " couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
2172 ERReturn(-1);
2173 }
2174 }
2175 }
2176
2177 return(FEI_SUCCESS);
2178}
2179
2180//------------------------------------------------------------------------------
2181int FEDataFilter::putBlockNodeSolution(GlobalID elemBlockID,
2182 int numNodes,
2183 const GlobalID *nodeIDs,
2184 const int *offsets,
2185 const double *estimates) {
2186
2187 debugOutput("FEI: putBlockNodeSolution");
2188
2189 int numActiveNodes = problemStructure_->getNumActiveNodes();
2190
2191 if (numActiveNodes <= 0) return(0);
2192
2193 BlockDescriptor* block = NULL;
2194 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2195
2196 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2197
2198 //traverse the node list, checking for nodes associated with this block
2199 //when an associated node is found, put its 'answers' into the linear system.
2200
2201 int blk_idx = problemStructure_->getIndexOfBlock(elemBlockID);
2202
2203 for(int i=0; i<numNodes; i++) {
2204 NodeDescriptor* node = NULL;
2205 int err = nodeDB.getNodeWithID(nodeIDs[i], node);
2206
2207 if (err != 0) continue;
2208
2209 if (!node->hasBlockIndex(blk_idx)) continue;
2210
2211 if (node->getOwnerProc() != localRank_) continue;
2212
2213 int numFields = node->getNumFields();
2214 const int* fieldIDs = node->getFieldIDList();
2215 const int* fieldEqnNumbers = node->getFieldEqnNumbers();
2216
2217 if (fieldEqnNumbers[0] < localStartRow_ ||
2218 fieldEqnNumbers[0] > localEndRow_) continue;
2219
2220 int offs = offsets[i];
2221
2222 for(int j=0; j<numFields; j++) {
2223 int size = problemStructure_->getFieldSize(fieldIDs[j]);
2224
2225 if (block->containsField(fieldIDs[j])) {
2226 for(int k=0; k<size; k++) {
2227 int reducedEqn;
2228 problemStructure_->
2229 translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
2230 }
2231 }
2232 offs += size;
2233 }//for(j<numFields)loop
2234 }
2235
2236 return(FEI_SUCCESS);
2237}
2238
2239//------------------------------------------------------------------------------
2240int FEDataFilter::putBlockFieldNodeSolution(GlobalID elemBlockID,
2241 int fieldID,
2242 int numNodes,
2243 const GlobalID *nodeIDs,
2244 const double *estimates)
2245{
2246 debugOutput("FEI: putBlockFieldNodeSolution");
2247
2248 BlockDescriptor* block = NULL;
2249 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2250 if (!block->containsField(fieldID)) return(1);
2251
2252 int fieldSize = problemStructure_->getFieldSize(fieldID);
2253 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2254
2255 //if we have a negative fieldID, we'll need a list of length numNodes,
2256 //in which to put nodeNumbers for passing to the solver...
2257
2258 std::vector<int> numbers(numNodes);
2259
2260 //if we have a fieldID >= 0, then our numbers list will hold equation numbers
2261 //and we'll need fieldSize*numNodes of them.
2262
2263 std::vector<double> data;
2264
2265 if (fieldID >= 0) {
2266 if (fieldSize < 1) {
2267 fei::console_out() << "FEI Warning, putBlockFieldNodeSolution called for field "
2268 << fieldID<<", which has size "<<fieldSize<<FEI_ENDL;
2269 return(0);
2270 }
2271 try {
2272 numbers.resize(numNodes*fieldSize);
2273 data.resize(numNodes*fieldSize);
2274 }
2275 catch(std::runtime_error& exc) {
2276 fei::console_out() << exc.what()<<FEI_ENDL;
2277 ERReturn(-1);
2278 }
2279 }
2280
2281 int count = 0;
2282
2283 for(int i=0; i<numNodes; i++) {
2284 NodeDescriptor* node = NULL;
2285 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2286
2287 if (fieldID < 0) numbers[count++] = node->getNodeNumber();
2288 else {
2289 int eqn = -1;
2290 if (node->getFieldEqnNumber(fieldID, eqn)) {
2291 if (eqn >= localStartRow_ && eqn <= localEndRow_) {
2292 for(int j=0; j<fieldSize; j++) {
2293 data[count] = estimates[i*fieldSize + j];
2294 problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
2295 }
2296 }
2297 }
2298 }
2299 }
2300
2301 if (fieldID < 0) {
2302 CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2303 numNodes, &numbers[0],
2304 estimates));
2305 }
2306
2307 return(FEI_SUCCESS);
2308}
2309
2310//------------------------------------------------------------------------------
2311int FEDataFilter::getBlockElemSolution(GlobalID elemBlockID,
2312 int numElems,
2313 const GlobalID *elemIDs,
2314 int& numElemDOFPerElement,
2315 double *results)
2316{
2317//
2318// return the elemental solution parameters associated with a
2319// particular block of elements
2320//
2321 debugOutput("FEI: getBlockElemSolution");
2322
2323 BlockDescriptor* block = NULL;
2324 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2325
2326 std::map<GlobalID,int>& elemIDList = problemStructure_->
2327 getBlockConnectivity(elemBlockID).elemIDs;
2328
2329 int len = block->getNumElements();
2330
2331 //if the user is only asking for a subset of element-solutions, shrink len.
2332 if (len > numElems) len = numElems;
2333
2334 numElemDOFPerElement = block->getNumElemDOFPerElement();
2335 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2336 double answer;
2337
2338
2339 if (numElemDOFPerElement <= 0) return(0);
2340
2341 std::map<GlobalID,int>::const_iterator
2342 elemid_itr = elemIDList.begin();
2343
2344 for(int i=0; i<len; i++) {
2345 int index = i;
2346
2347 //if the user-supplied elemIDs are out of order, we need the index of
2348 //the location of this element.
2349 if (elemid_itr->first != elemIDs[i]) {
2350 index = elemid_itr->second;
2351 }
2352
2353 if (index < 0) continue;
2354
2355 int offset = i*numElemDOFPerElement;
2356
2357 for(int j=0; j<numElemDOFPerElement; j++) {
2358 int eqn = elemDOFEqnNumbers[index] + j;
2359
2360 CHK_ERR( getEqnSolnEntry(eqn, answer) )
2361
2362 results[offset+j] = answer;
2363 }
2364
2365 ++elemid_itr;
2366 }
2367
2368 return(FEI_SUCCESS);
2369}
2370
2371//------------------------------------------------------------------------------
2372int FEDataFilter::putBlockElemSolution(GlobalID elemBlockID,
2373 int numElems,
2374 const GlobalID *elemIDs,
2375 int dofPerElem,
2376 const double *estimates)
2377{
2378 debugOutput("FEI: putBlockElemSolution");
2379
2380 BlockDescriptor* block = NULL;
2381 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2382
2383 std::map<GlobalID,int>& elemIDList = problemStructure_->
2384 getBlockConnectivity(elemBlockID).elemIDs;
2385
2386 int len = block->getNumElements();
2387 if (len > numElems) len = numElems;
2388
2389 int DOFPerElement = block->getNumElemDOFPerElement();
2390 if (DOFPerElement != dofPerElem) {
2391 fei::console_out() << "FEI ERROR, putBlockElemSolution called with bad 'dofPerElem' ("
2392 <<dofPerElem<<"), block "<<elemBlockID<<" should have dofPerElem=="
2393 <<DOFPerElement<<FEI_ENDL;
2394 ERReturn(-1);
2395 }
2396
2397 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2398
2399 if (DOFPerElement <= 0) return(0);
2400
2401 std::map<GlobalID,int>::const_iterator
2402 elemid_itr = elemIDList.begin();
2403
2404 for(int i=0; i<len; i++) {
2405 int index = i;
2406 if (elemid_itr->first != elemIDs[i]) {
2407 index = elemid_itr->second;
2408 }
2409
2410 if (index < 0) continue;
2411
2412 for(int j=0; j<DOFPerElement; j++) {
2413 int reducedEqn;
2414 problemStructure_->
2415 translateToReducedEqn(elemDOFEqnNumbers[i] + j, reducedEqn);
2416// double soln = estimates[i*DOFPerElement + j];
2417
2418// if (useLinSysCore_) {
2419// CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
2420// }
2421 }
2422
2423 ++elemid_itr;
2424 }
2425
2426 return(FEI_SUCCESS);
2427}
2428
2429//------------------------------------------------------------------------------
2430int FEDataFilter::getCRMultipliers(int numCRs,
2431 const int* CRIDs,
2432 double* multipliers)
2433{
2434 for(int i=0; i<numCRs; i++) {
2435 //temporarily, FETI's getMultiplierSoln method isn't implemented.
2436 //CHK_ERR( feData_->getMultiplierSoln(CRIDs[i], multipliers[i]) );
2437 multipliers[i] = -999.99;
2438 }
2439
2440 return(-1);
2441}
2442
2443//------------------------------------------------------------------------------
2444int FEDataFilter::putCRMultipliers(int numMultCRs,
2445 const int* CRIDs,
2446 const double *multEstimates)
2447{
2448 debugOutput("FEI: putCRMultipliers");
2449
2450 return(FEI_SUCCESS);
2451}
2452
2453//------------------------------------------------------------------------------
2454int FEDataFilter::putNodalFieldData(int fieldID,
2455 int numNodes,
2456 const GlobalID* nodeIDs,
2457 const double* nodeData)
2458{
2459 debugOutput("FEI: putNodalFieldData");
2460
2461 int fieldSize = problemStructure_->getFieldSize(fieldID);
2462 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2463
2464 std::vector<int> nodeNumbers(numNodes);
2465
2466 for(int i=0; i<numNodes; i++) {
2467 NodeDescriptor* node = NULL;
2468 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2469
2470 int nodeNumber = node->getNodeNumber();
2471 if (nodeNumber < 0) {
2472 GlobalID nodeID = nodeIDs[i];
2473 fei::console_out() << "FEDataFilter::putNodalFieldData ERROR, node with ID "
2474 << static_cast<int>(nodeID) << " doesn't have an associated nodeNumber "
2475 << "assigned. putNodalFieldData shouldn't be called until after the "
2476 << "initComplete method has been called." << FEI_ENDL;
2477 ERReturn(-1);
2478 }
2479
2480 nodeNumbers[i] = nodeNumber;
2481 }
2482
2483 CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2484 numNodes, &nodeNumbers[0],
2485 nodeData) );
2486
2487 return(0);
2488}
2489
2490//------------------------------------------------------------------------------
2491int FEDataFilter::putNodalFieldSolution(int fieldID,
2492 int numNodes,
2493 const GlobalID* nodeIDs,
2494 const double* nodeData)
2495{
2496 debugOutput("FEI: putNodalFieldSolution");
2497
2498 if (fieldID < 0) {
2499 return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
2500 }
2501
2502 int fieldSize = problemStructure_->getFieldSize(fieldID);
2503 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2504
2505 std::vector<int> eqnNumbers(fieldSize);
2506
2507 for(int i=0; i<numNodes; i++) {
2508 NodeDescriptor* node = NULL;
2509 CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2510
2511 int eqn = -1;
2512 bool hasField = node->getFieldEqnNumber(fieldID, eqn);
2513 if (!hasField) continue;
2514
2515 }
2516
2517 return(0);
2518}
2519
2520//------------------------------------------------------------------------------
2521int FEDataFilter::assembleEqns(int numRows,
2522 int numCols,
2523 const int* rowNumbers,
2524 const int* colIndices,
2525 const double* const* coefs,
2526 bool structurallySymmetric,
2527 int mode)
2528{
2529 if (numRows == 0) return(FEI_SUCCESS);
2530
2531 const int* indPtr = colIndices;
2532 for(int i=0; i<numRows; i++) {
2533 int row = rowNumbers[i];
2534
2535 const double* coefPtr = coefs[i];
2536
2537 CHK_ERR(giveToMatrix(1, &row, numCols, indPtr, &coefPtr, mode));
2538
2539 if (!structurallySymmetric) indPtr += numCols;
2540 }
2541
2542 return(FEI_SUCCESS);
2543}
2544
2545//------------------------------------------------------------------------------
2546int FEDataFilter::assembleRHS(int numValues,
2547 const int* indices,
2548 const double* coefs,
2549 int mode) {
2550//
2551//This function hands the data off to the routine that finally
2552//sticks it into the RHS vector.
2553//
2554 if (problemStructure_->numSlaveEquations() == 0) {
2555 CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
2556 return(FEI_SUCCESS);
2557 }
2558
2559 for(int i = 0; i < numValues; i++) {
2560 int eqn = indices[i];
2561 if (eqn < 0) continue;
2562
2563 CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
2564 }
2565
2566 return(FEI_SUCCESS);
2567}
2568
2569//==============================================================================
2570void FEDataFilter::debugOutput(const char* mesg)
2571{
2572 if (Filter::logStream() != NULL) {
2573 (*logStream()) << mesg << FEI_ENDL;
2574 }
2575}
EqnCommMgr * deepCopy()
int getParameters(int &numParams, char **&paramStrings)
virtual int setDirichletBCs(int numBCs, const int *nodeNumbers, const int *dof_ids, const double *values)=0
virtual int setMultiplierCR(int CRID, int numNodes, const int *nodeNumbers, const int *dof_ids, const double *coefWeights, double rhsValue)=0
virtual int setElemMatrix(int elemBlockID, int elemID, int numNodes, const int *nodeNumbers, const int *numDofPerNode, const int *dof_ids, const double *const *coefs)=0
virtual int sumIntoMatrix(int numRowNodes, const int *rowNodeNumbers, const int *row_dof_ids, const int *numColNodesPerRow, const int *colNodeNumbers, const int *col_dof_ids, const double *coefs)=0
virtual int describeStructure(int numElemBlocks, const int *numElemsPerBlock, const int *numNodesPerElem, const int *elemMatrixSizePerBlock, int totalNumNodes, int numSharedNodes, int numMultCRs)=0
virtual int setConnectivity(int elemBlockID, int elemID, int numNodes, const int *nodeNumbers, const int *numDofPerNode, const int *dof_ids)=0
virtual int reset()=0
virtual int deleteConstraints()=0
virtual int setElemVector(int elemBlockID, int elemID, int numNodes, const int *nodeNumbers, const int *numDofPerNode, const int *dof_ids, const double *coefs)=0
virtual int loadComplete()=0
virtual int launchSolver(int &solveStatus, int &iterations)=0
virtual int getSolnEntry(int nodeNumber, int dof_id, double &value)=0
virtual int putNodalFieldData(int fieldID, int fieldSize, int numNodes, const int *nodeNumbers, const double *coefs)=0
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
int getNodeWithNumber(int nodeNumber, const NodeDescriptor *&node) const
int getNodeWithEqn(int eqnNumber, const NodeDescriptor *&node) const
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int getNumNodeDescriptors() const
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
void getFieldID(int eqnNumber, int &fieldID, int &offset_into_field) const
int getAssociatedNodeNumber(int eqnNumber)
int getFieldSize(int fieldID)
int getBlockDescriptor_index(int index, BlockDescriptor *&block)
int translateFromReducedEqn(int reducedEqn)
const int * getFieldIDsPtr()
int getIndexOfBlock(GlobalID blockID) const
bool isInLocalElement(int nodeNumber)
bool translateToReducedEqn(int eqn, int &reducedEqn)
std::vector< int > & getMasterFieldIDs()
std::vector< int > & getMasters()
int localProc(MPI_Comm comm)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int numProcs(MPI_Comm comm)