FEI  Version of the Day
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 //------------------------------------------------------------------------------
45 void 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 //------------------------------------------------------------------------------
73 void 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 //------------------------------------------------------------------------------
98 FEDataFilter::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 //------------------------------------------------------------------------------
182 FEDataFilter::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 //------------------------------------------------------------------------------
221 FEDataFilter::~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 //------------------------------------------------------------------------------
236 int 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 //------------------------------------------------------------------------------
284 int 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 //==============================================================================
297 int 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 //------------------------------------------------------------------------------
548 int 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 //------------------------------------------------------------------------------
565 int 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 //------------------------------------------------------------------------------
578 int 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 //------------------------------------------------------------------------------
586 int 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 //------------------------------------------------------------------------------
594 int 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 //------------------------------------------------------------------------------
612 int 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 //------------------------------------------------------------------------------
630 int 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 //------------------------------------------------------------------------------
651 int 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 //------------------------------------------------------------------------------
716 int 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 //------------------------------------------------------------------------------
727 void 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 //------------------------------------------------------------------------------
755 int 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 //------------------------------------------------------------------------------
783 int 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 //------------------------------------------------------------------------------
810 int 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 //------------------------------------------------------------------------------
836 int 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 //------------------------------------------------------------------------------
849 int 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 //------------------------------------------------------------------------------
1006 int 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 //------------------------------------------------------------------------------
1032 int 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 //------------------------------------------------------------------------------
1056 int 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 //------------------------------------------------------------------------------
1081 int 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 //------------------------------------------------------------------------------
1121 int 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 //------------------------------------------------------------------------------
1198 int 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 //------------------------------------------------------------------------------
1278 int 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 //------------------------------------------------------------------------------
1297 int 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 //------------------------------------------------------------------------------
1319 int 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 //------------------------------------------------------------------------------
1354 int 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 //------------------------------------------------------------------------------
1366 int 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 //------------------------------------------------------------------------------
1408 int FEDataFilter::formResidual(double* residValues, int numLocalEqns)
1409 {
1410  //FiniteElementData implementations can't currently do residuals.
1411  return(FEI_SUCCESS);
1412 }
1413 
1414 //------------------------------------------------------------------------------
1415 int 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 //------------------------------------------------------------------------------
1442 int 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 //------------------------------------------------------------------------------
1461 int 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 //------------------------------------------------------------------------------
1474 int 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 //------------------------------------------------------------------------------
1550 int 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 //------------------------------------------------------------------------------
1626 int 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 //------------------------------------------------------------------------------
1635 int FEDataFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
1636 {
1637  ERReturn(-1);
1638 }
1639 
1640 //------------------------------------------------------------------------------
1641 int FEDataFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
1642 {
1643  ERReturn(-1);
1644 }
1645 
1646 //------------------------------------------------------------------------------
1647 int 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 //------------------------------------------------------------------------------
1690 int 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 //------------------------------------------------------------------------------
1726 int FEDataFilter::getFromRHS(int num, double* values, const int* indices)
1727 {
1728  return(-1);
1729 }
1730 
1731 //------------------------------------------------------------------------------
1732 int 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 //------------------------------------------------------------------------------
1753 int 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 //------------------------------------------------------------------------------
1769 int 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 //------------------------------------------------------------------------------
1822 int 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 //------------------------------------------------------------------------------
1867 void FEDataFilter:: setEqnCommMgr(EqnCommMgr* eqnCommMgr)
1868 {
1869  delete eqnCommMgr_;
1870  eqnCommMgr_ = eqnCommMgr;
1871 }
1872 
1873 //------------------------------------------------------------------------------
1874 int 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 //------------------------------------------------------------------------------
1954 int 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 //------------------------------------------------------------------------------
2028 int 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 //------------------------------------------------------------------------------
2097 int 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  NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
2121 
2122  //Traverse the node list, checking if nodes have the specified field.
2123  //If so, put the answers in the results list.
2124 
2125  for(int i=0; i<numNodes; i++) {
2126  NodeDescriptor* node_i = NULL;
2127  nodeDB.getNodeAtIndex(i, node_i);
2128 
2129  GlobalID nodeID = nodeIDs[i];
2130 
2131  NodeDescriptor* node = NULL;
2132  int err = 0;
2133  //Obtain the NodeDescriptor of nodeID in the activeNodes list...
2134  //Don't call the getNodeWithID function unless we have to.
2135 
2136  if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2137  node = node_i;
2138  }
2139  else {
2140  err = nodeDB.getNodeWithID(nodeID, node);
2141  }
2142 
2143  //ok. If err is not 0, meaning nodeID is NOT in the
2144  //activeNodes list, then skip to the next loop iteration.
2145 
2146  if (err != 0) {
2147  continue;
2148  }
2149 
2150  int nodeNumber = node->getNodeNumber();
2151 
2152  int eqnNumber = -1;
2153  bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
2154 
2155  //If this node doesn't have the specified field, then skip to the
2156  //next loop iteration.
2157  if (!hasField) continue;
2158 
2159  std::vector<GlobalID>& shNodeIDs = nodeCommMgr.getSharedNodeIDs();
2160  int shIndex = fei::binarySearch(nodeID, &shNodeIDs[0], shNodeIDs.size());
2161  if (shIndex > -1) {
2162  if (!(problemStructure_->isInLocalElement(nodeNumber))) continue;
2163  }
2164 
2165  int dof_id = problemStructure_->getFieldDofMap().get_dof_id(fieldID, 0);
2166 
2167  int offset = fieldSize*i;
2168 
2169  for(int j=0; j<fieldSize; j++) {
2170  if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
2171  CHK_ERR( getSharedRemoteSolnEntry(eqnNumber+j, results[offset+j]) );
2172  continue;
2173  }
2174 
2175  err = feData_->getSolnEntry(nodeNumber, dof_id+j, results[offset+j]);
2176  if (err != 0) {
2177  fei::console_out() << "FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
2178  << " (nodeID " << nodeID << "), dof_id "<<dof_id
2179  << " couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
2180  ERReturn(-1);
2181  }
2182  }
2183  }
2184 
2185  return(FEI_SUCCESS);
2186 }
2187 
2188 //------------------------------------------------------------------------------
2189 int FEDataFilter::putBlockNodeSolution(GlobalID elemBlockID,
2190  int numNodes,
2191  const GlobalID *nodeIDs,
2192  const int *offsets,
2193  const double *estimates) {
2194 
2195  debugOutput("FEI: putBlockNodeSolution");
2196 
2197  int numActiveNodes = problemStructure_->getNumActiveNodes();
2198 
2199  if (numActiveNodes <= 0) return(0);
2200 
2201  BlockDescriptor* block = NULL;
2202  CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2203 
2204  NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2205 
2206  //traverse the node list, checking for nodes associated with this block
2207  //when an associated node is found, put its 'answers' into the linear system.
2208 
2209  int blk_idx = problemStructure_->getIndexOfBlock(elemBlockID);
2210 
2211  for(int i=0; i<numNodes; i++) {
2212  NodeDescriptor* node = NULL;
2213  int err = nodeDB.getNodeWithID(nodeIDs[i], node);
2214 
2215  if (err != 0) continue;
2216 
2217  if (!node->hasBlockIndex(blk_idx)) continue;
2218 
2219  if (node->getOwnerProc() != localRank_) continue;
2220 
2221  int numFields = node->getNumFields();
2222  const int* fieldIDs = node->getFieldIDList();
2223  const int* fieldEqnNumbers = node->getFieldEqnNumbers();
2224 
2225  if (fieldEqnNumbers[0] < localStartRow_ ||
2226  fieldEqnNumbers[0] > localEndRow_) continue;
2227 
2228  int offs = offsets[i];
2229 
2230  for(int j=0; j<numFields; j++) {
2231  int size = problemStructure_->getFieldSize(fieldIDs[j]);
2232 
2233  if (block->containsField(fieldIDs[j])) {
2234  for(int k=0; k<size; k++) {
2235  int reducedEqn;
2236  problemStructure_->
2237  translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
2238  }
2239  }
2240  offs += size;
2241  }//for(j<numFields)loop
2242  }
2243 
2244  return(FEI_SUCCESS);
2245 }
2246 
2247 //------------------------------------------------------------------------------
2248 int FEDataFilter::putBlockFieldNodeSolution(GlobalID elemBlockID,
2249  int fieldID,
2250  int numNodes,
2251  const GlobalID *nodeIDs,
2252  const double *estimates)
2253 {
2254  debugOutput("FEI: putBlockFieldNodeSolution");
2255 
2256  BlockDescriptor* block = NULL;
2257  CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2258  if (!block->containsField(fieldID)) return(1);
2259 
2260  int fieldSize = problemStructure_->getFieldSize(fieldID);
2261  NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2262 
2263  //if we have a negative fieldID, we'll need a list of length numNodes,
2264  //in which to put nodeNumbers for passing to the solver...
2265 
2266  std::vector<int> numbers(numNodes);
2267 
2268  //if we have a fieldID >= 0, then our numbers list will hold equation numbers
2269  //and we'll need fieldSize*numNodes of them.
2270 
2271  std::vector<double> data;
2272 
2273  if (fieldID >= 0) {
2274  if (fieldSize < 1) {
2275  fei::console_out() << "FEI Warning, putBlockFieldNodeSolution called for field "
2276  << fieldID<<", which has size "<<fieldSize<<FEI_ENDL;
2277  return(0);
2278  }
2279  try {
2280  numbers.resize(numNodes*fieldSize);
2281  data.resize(numNodes*fieldSize);
2282  }
2283  catch(std::runtime_error& exc) {
2284  fei::console_out() << exc.what()<<FEI_ENDL;
2285  ERReturn(-1);
2286  }
2287  }
2288 
2289  int count = 0;
2290 
2291  for(int i=0; i<numNodes; i++) {
2292  NodeDescriptor* node = NULL;
2293  CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2294 
2295  if (fieldID < 0) numbers[count++] = node->getNodeNumber();
2296  else {
2297  int eqn = -1;
2298  if (node->getFieldEqnNumber(fieldID, eqn)) {
2299  if (eqn >= localStartRow_ && eqn <= localEndRow_) {
2300  for(int j=0; j<fieldSize; j++) {
2301  data[count] = estimates[i*fieldSize + j];
2302  problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
2303  }
2304  }
2305  }
2306  }
2307  }
2308 
2309  if (fieldID < 0) {
2310  CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2311  numNodes, &numbers[0],
2312  estimates));
2313  }
2314 
2315  return(FEI_SUCCESS);
2316 }
2317 
2318 //------------------------------------------------------------------------------
2319 int FEDataFilter::getBlockElemSolution(GlobalID elemBlockID,
2320  int numElems,
2321  const GlobalID *elemIDs,
2322  int& numElemDOFPerElement,
2323  double *results)
2324 {
2325 //
2326 // return the elemental solution parameters associated with a
2327 // particular block of elements
2328 //
2329  debugOutput("FEI: getBlockElemSolution");
2330 
2331  BlockDescriptor* block = NULL;
2332  CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2333 
2334  std::map<GlobalID,int>& elemIDList = problemStructure_->
2335  getBlockConnectivity(elemBlockID).elemIDs;
2336 
2337  int len = block->getNumElements();
2338 
2339  //if the user is only asking for a subset of element-solutions, shrink len.
2340  if (len > numElems) len = numElems;
2341 
2342  numElemDOFPerElement = block->getNumElemDOFPerElement();
2343  std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2344  double answer;
2345 
2346 
2347  if (numElemDOFPerElement <= 0) return(0);
2348 
2349  std::map<GlobalID,int>::const_iterator
2350  elemid_end = elemIDList.end(),
2351  elemid_itr = elemIDList.begin();
2352 
2353  for(int i=0; i<len; i++) {
2354  int index = i;
2355 
2356  //if the user-supplied elemIDs are out of order, we need the index of
2357  //the location of this element.
2358  if (elemid_itr->first != elemIDs[i]) {
2359  index = elemid_itr->second;
2360  }
2361 
2362  if (index < 0) continue;
2363 
2364  int offset = i*numElemDOFPerElement;
2365 
2366  for(int j=0; j<numElemDOFPerElement; j++) {
2367  int eqn = elemDOFEqnNumbers[index] + j;
2368 
2369  CHK_ERR( getEqnSolnEntry(eqn, answer) )
2370 
2371  results[offset+j] = answer;
2372  }
2373 
2374  ++elemid_itr;
2375  }
2376 
2377  return(FEI_SUCCESS);
2378 }
2379 
2380 //------------------------------------------------------------------------------
2381 int FEDataFilter::putBlockElemSolution(GlobalID elemBlockID,
2382  int numElems,
2383  const GlobalID *elemIDs,
2384  int dofPerElem,
2385  const double *estimates)
2386 {
2387  debugOutput("FEI: putBlockElemSolution");
2388 
2389  BlockDescriptor* block = NULL;
2390  CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2391 
2392  std::map<GlobalID,int>& elemIDList = problemStructure_->
2393  getBlockConnectivity(elemBlockID).elemIDs;
2394 
2395  int len = block->getNumElements();
2396  if (len > numElems) len = numElems;
2397 
2398  int DOFPerElement = block->getNumElemDOFPerElement();
2399  if (DOFPerElement != dofPerElem) {
2400  fei::console_out() << "FEI ERROR, putBlockElemSolution called with bad 'dofPerElem' ("
2401  <<dofPerElem<<"), block "<<elemBlockID<<" should have dofPerElem=="
2402  <<DOFPerElement<<FEI_ENDL;
2403  ERReturn(-1);
2404  }
2405 
2406  std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2407 
2408  if (DOFPerElement <= 0) return(0);
2409 
2410  std::map<GlobalID,int>::const_iterator
2411  elemid_end = elemIDList.end(),
2412  elemid_itr = elemIDList.begin();
2413 
2414  for(int i=0; i<len; i++) {
2415  int index = i;
2416  if (elemid_itr->first != elemIDs[i]) {
2417  index = elemid_itr->second;
2418  }
2419 
2420  if (index < 0) continue;
2421 
2422  for(int j=0; j<DOFPerElement; j++) {
2423  int reducedEqn;
2424  problemStructure_->
2425  translateToReducedEqn(elemDOFEqnNumbers[i] + j, reducedEqn);
2426 // double soln = estimates[i*DOFPerElement + j];
2427 
2428 // if (useLinSysCore_) {
2429 // CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
2430 // }
2431  }
2432 
2433  ++elemid_itr;
2434  }
2435 
2436  return(FEI_SUCCESS);
2437 }
2438 
2439 //------------------------------------------------------------------------------
2440 int FEDataFilter::getCRMultipliers(int numCRs,
2441  const int* CRIDs,
2442  double* multipliers)
2443 {
2444  for(int i=0; i<numCRs; i++) {
2445  //temporarily, FETI's getMultiplierSoln method isn't implemented.
2446  //CHK_ERR( feData_->getMultiplierSoln(CRIDs[i], multipliers[i]) );
2447  multipliers[i] = -999.99;
2448  }
2449 
2450  return(-1);
2451 }
2452 
2453 //------------------------------------------------------------------------------
2454 int FEDataFilter::putCRMultipliers(int numMultCRs,
2455  const int* CRIDs,
2456  const double *multEstimates)
2457 {
2458  debugOutput("FEI: putCRMultipliers");
2459 
2460  return(FEI_SUCCESS);
2461 }
2462 
2463 //------------------------------------------------------------------------------
2464 int FEDataFilter::putNodalFieldData(int fieldID,
2465  int numNodes,
2466  const GlobalID* nodeIDs,
2467  const double* nodeData)
2468 {
2469  debugOutput("FEI: putNodalFieldData");
2470 
2471  int fieldSize = problemStructure_->getFieldSize(fieldID);
2472  NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2473 
2474  std::vector<int> nodeNumbers(numNodes);
2475 
2476  for(int i=0; i<numNodes; i++) {
2477  NodeDescriptor* node = NULL;
2478  CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2479 
2480  int nodeNumber = node->getNodeNumber();
2481  if (nodeNumber < 0) {
2482  GlobalID nodeID = nodeIDs[i];
2483  fei::console_out() << "FEDataFilter::putNodalFieldData ERROR, node with ID "
2484  << static_cast<int>(nodeID) << " doesn't have an associated nodeNumber "
2485  << "assigned. putNodalFieldData shouldn't be called until after the "
2486  << "initComplete method has been called." << FEI_ENDL;
2487  ERReturn(-1);
2488  }
2489 
2490  nodeNumbers[i] = nodeNumber;
2491  }
2492 
2493  CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2494  numNodes, &nodeNumbers[0],
2495  nodeData) );
2496 
2497  return(0);
2498 }
2499 
2500 //------------------------------------------------------------------------------
2501 int FEDataFilter::putNodalFieldSolution(int fieldID,
2502  int numNodes,
2503  const GlobalID* nodeIDs,
2504  const double* nodeData)
2505 {
2506  debugOutput("FEI: putNodalFieldSolution");
2507 
2508  if (fieldID < 0) {
2509  return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
2510  }
2511 
2512  int fieldSize = problemStructure_->getFieldSize(fieldID);
2513  NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2514 
2515  std::vector<int> eqnNumbers(fieldSize);
2516 
2517  for(int i=0; i<numNodes; i++) {
2518  NodeDescriptor* node = NULL;
2519  CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
2520 
2521  int eqn = -1;
2522  bool hasField = node->getFieldEqnNumber(fieldID, eqn);
2523  if (!hasField) continue;
2524 
2525  }
2526 
2527  return(0);
2528 }
2529 
2530 //------------------------------------------------------------------------------
2531 int FEDataFilter::assembleEqns(int numRows,
2532  int numCols,
2533  const int* rowNumbers,
2534  const int* colIndices,
2535  const double* const* coefs,
2536  bool structurallySymmetric,
2537  int mode)
2538 {
2539  if (numRows == 0) return(FEI_SUCCESS);
2540 
2541  const int* indPtr = colIndices;
2542  for(int i=0; i<numRows; i++) {
2543  int row = rowNumbers[i];
2544 
2545  const double* coefPtr = coefs[i];
2546 
2547  CHK_ERR(giveToMatrix(1, &row, numCols, indPtr, &coefPtr, mode));
2548 
2549  if (!structurallySymmetric) indPtr += numCols;
2550  }
2551 
2552  return(FEI_SUCCESS);
2553 }
2554 
2555 //------------------------------------------------------------------------------
2556 int FEDataFilter::assembleRHS(int numValues,
2557  const int* indices,
2558  const double* coefs,
2559  int mode) {
2560 //
2561 //This function hands the data off to the routine that finally
2562 //sticks it into the RHS vector.
2563 //
2564  if (problemStructure_->numSlaveEquations() == 0) {
2565  CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
2566  return(FEI_SUCCESS);
2567  }
2568 
2569  for(int i = 0; i < numValues; i++) {
2570  int eqn = indices[i];
2571  if (eqn < 0) continue;
2572 
2573  CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
2574  }
2575 
2576  return(FEI_SUCCESS);
2577 }
2578 
2579 //==============================================================================
2580 void FEDataFilter::debugOutput(const char* mesg)
2581 {
2582  if (Filter::logStream() != NULL) {
2583  (*logStream()) << mesg << FEI_ENDL;
2584  }
2585 }
int getParameters(int &numParams, char **&paramStrings)
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
std::vector< int > & getMasterFieldIDs()
virtual void parameters(const fei::ParameterSet &paramset)
Definition: fei_Factory.cpp:38
void getFieldID(int eqnNumber, int &fieldID, int &offset_into_field) const
int getNumNodeDescriptors() const
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int localProc(MPI_Comm comm)
int getNodeWithNumber(int nodeNumber, const NodeDescriptor *&node) const
int getNodeWithEqn(int eqnNumber, const NodeDescriptor *&node) const
std::vector< int > & getMasters()
int numProcs(MPI_Comm comm)