FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_FEI_Impl.cpp
Go to the documentation of this file.
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 <fei_utils.hpp>
12 
13 #include <fei_FEI_Impl.hpp>
14 #include <fei_Record.hpp>
15 #include <fei_TemplateUtils.hpp>
16 #include <fei_ParameterSet.hpp>
17 #include <fei_base.hpp>
18 
19 #include <fei_Pattern.hpp>
20 #include <fei_LibraryWrapper.hpp>
21 #include <fei_Data.hpp>
22 #include <fei_defs.h>
23 
24 #include <stdexcept>
25 #include <cmath>
26 
27 #undef fei_file
28 #define fei_file "fei_FEI_Impl.cpp"
29 
30 #include <fei_ErrMacros.hpp>
31 
33  MPI_Comm comm,
34  int masterRank)
35  : wrapper_(1),
36  nodeIDType_(0),
37  elemIDType_(1),
38  constraintIDType_(2),
39  factory_(1),
40  rowSpace_(NULL),
41  matGraph_(),
42  x_(),
43  b_(),
44  A_(),
45  linSys_(),
46  newData_(false),
47  soln_fei_matrix_(NULL),
48  soln_fei_vector_(NULL),
49  comm_(comm),
50  masterRank_(masterRank),
51  localProc_(0),
52  numProcs_(1),
53  numParams_(0),
54  paramStrings_(NULL),
55  matrixIDs_(),
56  rhsIDs_(),
57  matScalars_(),
58  matScalarsSet_(false),
59  rhsScalars_(),
60  rhsScalarsSet_(false),
61  constraintID_(0),
62  index_soln_(0),
63  index_current_(0),
64  index_current_rhs_row_(0),
65  solveType_(0),
66  iterations_(0),
67  setSolveTypeCalled_(false),
68  initPhaseIsComplete_(false),
69  aggregateSystemFormed_(false),
70  newMatrixDataLoaded_(0),
71  solveCounter_(0),
72  initTime_(0.0),
73  loadTime_(0.0),
74  solveTime_(0.0),
75  solnReturnTime_(0.0),
76  iwork_(),
77  nodeset_(),
78  nodeset_filled_(false),
79  block_dof_per_elem_(),
80  any_blocks_have_elem_dof_(false)
81 {
82  wrapper_[0] = wrapper;
83 
84  factory_[0].reset(new snl_fei::Factory(comm, wrapper_[0]));
85  createdFactory_ = true;
86 
88 }
89 
91  MPI_Comm comm,
92  int masterRank)
93  : wrapper_(1),
94  nodeIDType_(0),
95  elemIDType_(1),
96  constraintIDType_(2),
97  factory_(1),
98  createdFactory_(false),
99  rowSpace_(NULL),
100  matGraph_(),
101  x_(),
102  b_(),
103  A_(),
104  linSys_(),
105  newData_(false),
106  soln_fei_matrix_(NULL),
107  soln_fei_vector_(NULL),
108  comm_(comm),
109  masterRank_(masterRank),
110  localProc_(0),
111  numProcs_(1),
112  numParams_(0),
113  paramStrings_(NULL),
114  matrixIDs_(),
115  rhsIDs_(),
116  matScalars_(),
117  matScalarsSet_(false),
118  rhsScalars_(),
119  rhsScalarsSet_(false),
120  constraintID_(0),
121  index_soln_(0),
122  index_current_(0),
123  index_current_rhs_row_(0),
124  solveType_(0),
125  iterations_(0),
126  setSolveTypeCalled_(false),
127  initPhaseIsComplete_(false),
128  aggregateSystemFormed_(false),
129  newMatrixDataLoaded_(0),
130  solveCounter_(0),
131  initTime_(0.0),
132  loadTime_(0.0),
133  solveTime_(0.0),
134  solnReturnTime_(0.0),
135  iwork_(),
136  nodeset_(),
137  nodeset_filled_(false),
138  block_dof_per_elem_(),
139  any_blocks_have_elem_dof_(false)
140 {
141  wrapper_[0].reset(0);
142 
143  const snl_fei::Factory* snlfactory
144  = dynamic_cast<const snl_fei::Factory*>(factory);
145  if (snlfactory != NULL) {
146  wrapper_[0] = snlfactory->get_LibraryWrapper();
147  }
148 
149  factory_[0] = factory->clone();
150 
152 }
153 
155 {
156  for(int k=0; k<numParams_; k++) {
157  delete [] paramStrings_[k];
158  }
159  delete [] paramStrings_;
160 
161  if (soln_fei_matrix_ != NULL && wrapper_[0].get() != NULL) {
162  fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
163  if (lsc.get() != NULL) {
164  lsc->destroyMatrixData(*soln_fei_matrix_);
165  delete soln_fei_matrix_;
166  soln_fei_matrix_ = NULL;
167  }
168  }
169 
170  if (soln_fei_vector_ != NULL && wrapper_[0].get() != NULL) {
171  fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
172  if (lsc.get() != NULL) {
173  lsc->destroyVectorData(*soln_fei_vector_);
174  delete soln_fei_vector_;
175  soln_fei_vector_ = NULL;
176  }
177  }
178 }
179 
181 {
182  localProc_ = fei::localProc(comm_);
183  numProcs_ = fei::numProcs(comm_);
184 
185  constraintID_ = localProc_*100000;
186 
187  matrixIDs_.resize(1);
188  matrixIDs_[0] = 0;
189  A_.resize(1);
190  rhsIDs_.resize(1);
191  rhsIDs_[0] = 0;
192  b_.resize(1);
193 
194  rowSpace_ = factory_[0]->createVectorSpace(comm_, (const char*)NULL);
195 
196  rowSpace_->defineIDTypes(1, &nodeIDType_);
197  rowSpace_->defineIDTypes(1, &elemIDType_);
198  rowSpace_->defineIDTypes(1, &constraintIDType_);
199 
200  matGraph_ = factory_[0]->createMatrixGraph(rowSpace_, rowSpace_,(const char*)NULL);
201  if (matGraph_.get() == NULL) {
202  voidERReturn;
203  }
204 }
205 
206 int fei::FEI_Impl::setIDLists(int numMatrices,
207  const int* matrixIDs,
208  int numRHSs,
209  const int* rhsIDs)
210 {
211  if (numMatrices < 1 || numRHSs < 1) {
212  fei::console_out() << "fei::FEI_Impl::setIDLists ERROR, numMatrices and numRHSs "
213  << "must both be greater than 0."<<FEI_ENDL;
214  ERReturn(-1);
215  }
216 
217  matrixIDs_.resize(0);
218  A_.resize(numMatrices);
219  for(int i=0; i<numMatrices; ++i) {
220  fei::sortedListInsert(matrixIDs[i], matrixIDs_);
221  }
222  if ((int)matrixIDs_.size() != numMatrices) {
223  fei::console_out() << "fei::FEI_Impl::setIDLists ERROR creating matrixIDs_ list."<<FEI_ENDL;
224  ERReturn(-1);
225  }
226 
227  rhsIDs_.resize(0);
228  b_.resize(numRHSs);
229  for(int i=0; i<numRHSs; ++i) {
230  fei::sortedListInsert(rhsIDs[i], rhsIDs_);
231  }
232  if ((int)rhsIDs_.size() != numRHSs) {
233  fei::console_out() << "fei::FEI_Impl::setIDLists ERROR creating rhsIDs_ list."<<FEI_ENDL;
234  ERReturn(-1);
235  }
236 
237  if (wrapper_[0].get() != NULL) {
238  fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
239  if (linSysCore.get() != NULL) {
240  linSysCore->setNumRHSVectors(rhsIDs_.size(), &rhsIDs_[0]);
241  }
242  }
243 
244  return(0);
245 }
246 
248 {
249  return( linSys_ );
250 }
251 
252 int fei::FEI_Impl::parameters(int numParams,
253  const char *const* paramStrings)
254 {
255  // merge these parameters with any others we may have, for later use.
256  snl_fei::mergeStringLists(paramStrings_, numParams_,
257  paramStrings, numParams);
258 
259  if (wrapper_[0].get() != NULL) {
260  if (wrapper_[0]->haveLinearSystemCore()) {
261  CHK_ERR( wrapper_[0]->getLinearSystemCore()->parameters(numParams, (char**)paramStrings) );
262  }
263  if (wrapper_[0]->haveFiniteElementData()) {
264  CHK_ERR( wrapper_[0]->getFiniteElementData()->parameters(numParams, (char**)paramStrings) );
265  }
266  }
267 
268  std::vector<std::string> stdstrings;
269  fei::utils::char_ptrs_to_strings(numParams, paramStrings, stdstrings);
270  fei::ParameterSet paramset;
271  fei::utils::parse_strings(stdstrings, " ", paramset);
272  factory_[0]->parameters(paramset);
273 
274  return(0);
275 }
276 
277 int fei::FEI_Impl::setSolveType(int solveType)
278 {
279  solveType_ = solveType;
280  setSolveTypeCalled_ = true;
281 
282  return(0);
283 }
284 
285 int fei::FEI_Impl::initFields(int numFields,
286  const int *fieldSizes,
287  const int *fieldIDs,
288  const int *fieldTypes)
289 {
290  rowSpace_->defineFields(numFields, fieldIDs, fieldSizes, fieldTypes);
291 
292  return(0);
293 }
294 
296  int numElements,
297  int numNodesPerElement,
298  const int* numFieldsPerNode,
299  const int* const* nodalFieldIDs,
300  int numElemDofFieldsPerElement,
301  const int* elemDOFFieldIDs,
302  int interleaveStrategy)
303 {
304  //define pattern that describes the layout of fields for elements in
305  //this element-block
306 
307  int numIDs = numNodesPerElement;
308  if (numElemDofFieldsPerElement > 0) ++numIDs;
309  std::vector<int> idTypes;
310  std::vector<int> numFieldsPerID;
311  std::vector<int> fieldIDs;
312 
313  int i, j;
314  for(i=0; i<numNodesPerElement; ++i) {
315  idTypes.push_back(nodeIDType_);
316  numFieldsPerID.push_back(numFieldsPerNode[i]);
317 
318  for(j=0; j<numFieldsPerNode[i]; ++j) {
319  fieldIDs.push_back(nodalFieldIDs[i][j]);
320  }
321  }
322 
323  if (numElemDofFieldsPerElement>0) {
324  idTypes.push_back(elemIDType_);
325  numFieldsPerID.push_back(numElemDofFieldsPerElement);
326  for(i=0; i<numElemDofFieldsPerElement; ++i) {
327  fieldIDs.push_back(elemDOFFieldIDs[i]);
328  }
329 
330  block_dof_per_elem_.insert(std::pair<int,int>(elemBlockID, numElemDofFieldsPerElement));
331  any_blocks_have_elem_dof_ = true;
332  }
333 
334  int pattID = matGraph_->definePattern(numIDs,
335  &idTypes[0],
336  &numFieldsPerID[0],
337  &fieldIDs[0]);
338 
339  //initialize connectivity-block
340  CHK_ERR( matGraph_->initConnectivityBlock(elemBlockID, numElements,
341  pattID) );
342 
343  return(0);
344 }
345 
347  GlobalID elemID,
348  const GlobalID* elemConn)
349 {
350  bool elemdof = false;
351 
352  if (any_blocks_have_elem_dof_) {
353  std::map<int,int>::const_iterator
354  b_iter = block_dof_per_elem_.find(elemBlockID);
355  if (b_iter != block_dof_per_elem_.end()) {
356  int numIDs = matGraph_->getNumIDsPerConnectivityList(elemBlockID);
357  iwork_.resize(numIDs);
358  int* iworkPtr = &iwork_[0];
359  for(int i=0; i<numIDs-1; ++i) {
360  iworkPtr[i] = elemConn[i];
361  }
362  iworkPtr[numIDs-1] = elemID;
363 
364  CHK_ERR( matGraph_->initConnectivity(elemBlockID, elemID, iworkPtr) );
365  elemdof = true;
366  }
367  }
368 
369  if (!elemdof) {
370  CHK_ERR( matGraph_->initConnectivity(elemBlockID, elemID, elemConn) );
371  }
372 
373  nodeset_filled_ = false;
374 
375  return(0);
376 }
377 
379  int slaveFieldID,
380  int offsetIntoSlaveField,
381  int numMasterNodes,
382  const GlobalID* masterNodeIDs,
383  const int* masterFieldIDs,
384  const double* weights,
385  double rhsValue)
386 {
387  throw std::runtime_error("FEI_Impl::initSlaveVariable not implemented.");
388  return(0);
389 }
390 
392 {
393  throw std::runtime_error("FEI_Impl::deleteMultCRs not implemented.");
394  return(0);
395 }
396 
397 int fei::FEI_Impl::initSharedNodes(int numSharedNodes,
398  const GlobalID *sharedNodeIDs,
399  const int* numProcsPerNode,
400  const int *const *sharingProcIDs)
401 {
402  CHK_ERR( rowSpace_->initSharedIDs(numSharedNodes,
403  nodeIDType_,
404  sharedNodeIDs,
405  numProcsPerNode,
406  sharingProcIDs) );
407 
408  return(0);
409 }
410 
411 int fei::FEI_Impl::initCRMult(int numCRNodes,
412  const GlobalID* CRNodeIDs,
413  const int *CRFieldIDs,
414  int& CRID)
415 {
416  iwork_.assign(numCRNodes,nodeIDType_);
417 
418  CRID = constraintID_++;
419 
420  CHK_ERR( matGraph_->initLagrangeConstraint(CRID,
421  constraintIDType_,
422  numCRNodes,
423  &iwork_[0],
424  CRNodeIDs,
425  CRFieldIDs) );
426 
427  return(0);
428 }
429 
430 int fei::FEI_Impl::initCRPen(int numCRNodes,
431  const GlobalID* CRNodeIDs,
432  const int *CRFieldIDs,
433  int& CRID)
434 {
435  iwork_.assign(numCRNodes, nodeIDType_);
436 
437  CRID = constraintID_++;
438 
439  CHK_ERR( matGraph_->initPenaltyConstraint(CRID,
440  constraintIDType_,
441  numCRNodes,
442  &iwork_[0],
443  CRNodeIDs,
444  CRFieldIDs) );
445 
446  return(0);
447 }
448 
450 {
451  CHK_ERR( matGraph_->initComplete() );
452 
453  if (matrixIDs_.size() < 1 || factory_.size() < 1 ||
454  A_.size() < 1 || b_.size() < 1) {
455  ERReturn(-1);
456  }
457 
458  A_[0] = factory_[0]->createMatrix(matGraph_);
459 
460  std::vector<std::string> stdstrings;
461  fei::utils::char_ptrs_to_strings(numParams_, paramStrings_, stdstrings);
462  fei::ParameterSet params;
463  fei::utils::parse_strings(stdstrings, " ", params);
464 
465  CHK_ERR( A_[0]->parameters(params) );
466 
467  b_[0] = factory_[0]->createVector(matGraph_);
468 
469  if (matrixIDs_.size() > 1) {
470  bool multiple_factories = false;
471 
472  if (wrapper_[0].get() != NULL) {
473  multiple_factories = true;
474 
475  fei::SharedPtr<LinearSystemCore> linsyscore = wrapper_[0]->getLinearSystemCore();
476  if (linsyscore.get() == NULL) {
477  fei::console_out() << "fei::FEI_Impl ERROR, multiple matrix/rhs assembly not supported "
478  << "non-null LibraryWrapper holds null LinearSystemCore."<<FEI_ENDL;
479  ERReturn(-1);
480  }
481 
482  wrapper_.resize(matrixIDs_.size());
483  factory_.resize(matrixIDs_.size());
484  for(unsigned i=1; i<matrixIDs_.size(); ++i) {
485  fei::SharedPtr<LinearSystemCore> lscclone(linsyscore->clone());
486  wrapper_[i].reset(new LibraryWrapper(lscclone));
487  factory_[i].reset(new snl_fei::Factory(comm_, wrapper_[i]));
488  }
489  }
490 
492  for(unsigned i=1; i<matrixIDs_.size(); ++i) {
493  factory = multiple_factories ? factory_[i] : factory_[0];
494 
495  A_[i] = factory->createMatrix(matGraph_);
496  CHK_ERR( A_[i]->parameters(params) );
497  }
498 
499  for(unsigned i=1; i<rhsIDs_.size(); ++i) {
500  factory = multiple_factories ? factory_[i] : factory_[0];
501 
502  b_[i] = factory->createVector(matGraph_);
503  }
504 
505  if (wrapper_[0].get() != NULL) {
507  = wrapper_[0]->getLinearSystemCore();
508 
509  linsyscore->setNumRHSVectors(1, &(rhsIDs_[0]));
510 
511  unsigned num = rhsIDs_.size();
512  if (matrixIDs_.size() < num) num = matrixIDs_.size();
513  for(unsigned i=1; i<num; ++i) {
514  fei::SharedPtr<LinearSystemCore> lsc = wrapper_[i]->getLinearSystemCore();
515 
516  if (i==num-1 && rhsIDs_.size() > num) {
517  int numRHSs = rhsIDs_.size() - matrixIDs_.size() + 1;
518  lsc->setNumRHSVectors(numRHSs, &(rhsIDs_[i]));
519  }
520  else {
521  lsc->setNumRHSVectors(1, &(rhsIDs_[i]));
522  }
523  }
524 
525  if (rhsIDs_.size() < matrixIDs_.size()) {
526  int dummyID = -1;
527  for(unsigned i=rhsIDs_.size(); i<matrixIDs_.size(); ++i) {
528  wrapper_[i]->getLinearSystemCore()->setNumRHSVectors(1, &dummyID);
529  }
530  }
531  }
532 
533  for(unsigned i=1; i<matrixIDs_.size(); ++i) {
534  factory = multiple_factories ? factory_[i] : factory_[0];
535 
536  A_[i] = factory->createMatrix(matGraph_);
537  CHK_ERR( A_[i]->parameters(params) );
538  }
539 
540  for(unsigned i=1; i<rhsIDs_.size(); ++i) {
541  factory = multiple_factories ? factory_[i] : factory_[0];
542 
543  b_[i] = factory->createVector(matGraph_);
544  }
545  }
546 
547  x_ = factory_[0]->createVector(matGraph_, true);
548 
549  linSys_ = factory_[0]->createLinearSystem(matGraph_);
550 
551  CHK_ERR( linSys_->parameters(numParams_, paramStrings_) );
552 
553  linSys_->setMatrix(A_[0]);
554  linSys_->setRHS(b_[0]);
555  linSys_->setSolutionVector(x_);
556 
557  return(0);
558 }
559 
561 {
562  std::vector<int>::const_iterator
563  iter = std::lower_bound(matrixIDs_.begin(), matrixIDs_.end(), matID);
564  if (iter == matrixIDs_.end() || *iter != matID) {
565  fei::console_out() << "fei::FEI_Impl::setCurrentMatrix: matID ("<<matID
566  <<") not found." <<FEI_ENDL;
567  return(-1);
568  }
569 
570  index_current_ = iter-matrixIDs_.begin();
571 
572  return(0);
573 }
574 
576 {
577  std::vector<int>::const_iterator
578  iter = std::lower_bound(rhsIDs_.begin(), rhsIDs_.end(), rhsID);
579  if (iter == rhsIDs_.end() || *iter != rhsID) {
580  fei::console_out() << "fei::FEI_Impl::setCurrentRHS: rhsID ("<<rhsID<<") not found."
581  << FEI_ENDL;
582  return(-1);
583  }
584 
585  index_current_rhs_row_ = iter - rhsIDs_.begin();
586 
587  return(0);
588 }
589 
591 {
592  int err = A_[index_current_]->putScalar(s);
593  err += x_->putScalar(s);
594  err += b_[index_current_rhs_row_]->putScalar(s);
595 
596  return(err);
597 }
598 
600 {
601  return( A_[index_current_]->putScalar(s) );
602 }
603 
605 {
606  return( b_[index_current_rhs_row_]->putScalar(s) );
607 }
608 
610 {
611  return( x_->putScalar(s) );
612 }
613 
614 int fei::FEI_Impl::loadNodeBCs(int numNodes,
615  const GlobalID *nodeIDs,
616  int fieldID,
617  const int* offsetsIntoField,
618  const double* prescribedValues)
619 {
620  CHK_ERR( linSys_->loadEssentialBCs(numNodes, nodeIDs,
621  nodeIDType_, fieldID,
622  offsetsIntoField, prescribedValues) );
623 
624  newData_ = true;
625 
626  return(0);
627 }
628 
629 int fei::FEI_Impl::loadElemBCs(int numElems,
630  const GlobalID* elemIDs,
631  int fieldID,
632  const double *const *alpha,
633  const double *const *beta,
634  const double *const *gamma)
635 {
636  throw std::runtime_error("FEI_Impl::loadElemBCs not implemented.");
637  return(0);
638 }
639 
641  GlobalID elemID,
642  const GlobalID* elemConn,
643  const double* const* elemStiffness,
644  const double* elemLoad,
645  int elemFormat)
646 {
647  CHK_ERR( A_[index_current_]->sumIn(elemBlockID, elemID, elemStiffness, elemFormat) );
648 
649  int num = matGraph_->getConnectivityNumIndices(elemBlockID);
650  std::vector<int> indices(num);
651  CHK_ERR( matGraph_->getConnectivityIndices(elemBlockID, elemID, num,
652  &indices[0], num) );
653  CHK_ERR( b_[index_current_rhs_row_]->sumIn(num, &indices[0], elemLoad, 0) );
654 
655  newData_ = true;
656 
657  return(0);
658 }
659 
661  GlobalID elemID,
662  const GlobalID* elemConn,
663  const double* const* elemStiffness,
664  int elemFormat)
665 {
666  CHK_ERR( A_[index_current_]->sumIn(elemBlockID, elemID, elemStiffness, elemFormat) );
667 
668  newData_ = true;
669 
670  return(0);
671 }
672 
674  GlobalID elemID,
675  const GlobalID* elemConn,
676  const double* elemLoad)
677 {
678  int num = matGraph_->getConnectivityNumIndices(elemBlockID);
679  std::vector<int> indices(num);
680  CHK_ERR( matGraph_->getConnectivityIndices(elemBlockID, elemID, num,
681  &indices[0], num) );
682  CHK_ERR( b_[index_current_rhs_row_]->sumIn(num, &indices[0], elemLoad, 0) );
683 
684  newData_ = true;
685 
686  return(0);
687 }
688 
689 int fei::FEI_Impl::loadCRMult(int CRMultID,
690  int numCRNodes,
691  const GlobalID* CRNodeIDs,
692  const int* CRFieldIDs,
693  const double* CRWeights,
694  double CRValue)
695 {
696  newData_ = true;
697 
698  CHK_ERR( linSys_->loadLagrangeConstraint(CRMultID, CRWeights, CRValue) );
699 
700  return(0);
701 }
702 
703 int fei::FEI_Impl::loadCRPen(int CRPenID,
704  int numCRNodes,
705  const GlobalID* CRNodeIDs,
706  const int* CRFieldIDs,
707  const double* CRWeights,
708  double CRValue,
709  double penValue)
710 {
711  newData_ = true;
712 
713  CHK_ERR( linSys_->loadPenaltyConstraint(CRPenID, CRWeights, penValue, CRValue) );
714 
715  return(0);
716 }
717 
719  int fieldID,
720  int numIDs,
721  const GlobalID* IDs,
722  const double* coefficients)
723 {
724  CHK_ERR( inputRHS(IDType, fieldID, numIDs, IDs, coefficients, false) );
725 
726  newData_ = true;
727 
728  return(0);
729 }
730 
732  int fieldID,
733  int numIDs,
734  const GlobalID* IDs,
735  const double* coefficients)
736 {
737  CHK_ERR( inputRHS(IDType, fieldID, numIDs, IDs, coefficients, true) );
738 
739  newData_ = true;
740 
741  return(0);
742 }
743 
744 int fei::FEI_Impl::inputRHS(int IDType,
745  int fieldID,
746  int numIDs,
747  const GlobalID* IDs,
748  const double* coefficients,
749  bool sumInto)
750 {
751  int fieldSize = rowSpace_->getFieldSize(fieldID);
752 
753  int offset = 0;
754  for(int i=0; i<numIDs; ++i) {
755  int globalIndex = 0;
756  CHK_ERR( rowSpace_->getGlobalIndex(IDType, IDs[i], globalIndex) );
757 
758  for(int j=0; j<fieldSize; ++j) {
759  int eqn = globalIndex+j;
760  if (sumInto) {
761  CHK_ERR( b_[index_current_rhs_row_]->sumIn(1, &eqn, &(coefficients[offset++])) );
762  }
763  else {
764  CHK_ERR( b_[index_current_rhs_row_]->copyIn(1, &eqn, &(coefficients[offset++])) );
765  }
766  }
767  }
768 
769  return(0);
770 }
771 
772 int fei::FEI_Impl::setMatScalars(int numScalars,
773  const int* IDs,
774  const double* scalars)
775 {
776  matScalars_.resize(matrixIDs_.size());
777 
778  for(int i=0; i<numScalars; ++i) {
779  std::vector<int>::const_iterator
780  iter = std::lower_bound(matrixIDs_.begin(), matrixIDs_.end(), IDs[i]);
781  if (iter == matrixIDs_.end() || *iter != IDs[i]) {
782  continue;
783  }
784 
785  unsigned index = iter - matrixIDs_.begin();
786  matScalars_[index] = scalars[i];
787  }
788 
789  matScalarsSet_ = true;
790 
791  return(0);
792 }
793 
794 int fei::FEI_Impl::setRHSScalars(int numScalars,
795  const int* IDs,
796  const double* scalars)
797 {
798  rhsScalars_.resize(rhsIDs_.size());
799 
800  for(int i=0; i<numScalars; ++i) {
801  std::vector<int>::const_iterator
802  iter = std::lower_bound(rhsIDs_.begin(), rhsIDs_.end(), IDs[i]);
803  if (iter == rhsIDs_.end() || *iter != IDs[i]) {
804  continue;
805  }
806 
807  unsigned index = iter - rhsIDs_.begin();
808  rhsScalars_[index] = scalars[i];
809  }
810 
811  rhsScalarsSet_ = true;
812 
813  return(0);
814 }
815 
816 int fei::FEI_Impl::loadComplete(bool applyBCs,
817  bool globalAssemble)
818 {
819  if (!newData_) {
820  return(0);
821  }
822 
823  if (linSys_.get() == NULL) {
824  fei::console_out() << "fei::FEI_Impl::loadComplete: loadComplete can not be called"
825  << " until after initComplete has been called."<<FEI_ENDL;
826  return(-1);
827  }
828 
829  if (solveType_ == FEI_AGGREGATE_SUM) {
830  for(unsigned i=0; i<A_.size(); ++i) {
831  CHK_ERR( A_[i]->gatherFromOverlap() );
832  }
833 
834  for(unsigned j=0; j<b_.size(); ++j) {
835  CHK_ERR( b_[j]->gatherFromOverlap() );
836  }
837 
838  CHK_ERR( aggregateSystem() );
839  }
840 
841  CHK_ERR( linSys_->loadComplete(applyBCs, globalAssemble) );
842 
843  if (b_.size() > 1) {
844  int rhs_counter = 0;
845  for(unsigned i=0; i<b_.size(); ++i) {
846  FEI_OSTRINGSTREAM osstr;
847  osstr << "rhs_" << rhs_counter++;
848  CHK_ERR( linSys_->putAttribute(osstr.str().c_str(), b_[i].get()) );
849  }
850  }
851 
852  newData_ = false;
853 
854  return(0);
855 }
856 
858 {
859  if (wrapper_[0].get() == NULL) {
860  ERReturn(-1);
861  }
862 
863  if (wrapper_[0].get() != NULL) {
864  CHK_ERR( aggregateSystem_LinSysCore() );
865  }
866 
867  return(0);
868 }
869 
871 {
872  fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
873  if (lsc.get() == NULL) return(-1);
874 
875  if (soln_fei_matrix_ == NULL) {
876  soln_fei_matrix_ = new Data();
877 
878  CHK_ERR( lsc->copyOutMatrix(1.0, *soln_fei_matrix_) );
879  }
880 
881  if (soln_fei_vector_ == NULL) {
882  soln_fei_vector_ = new Data();
883 
884  CHK_ERR( lsc->setRHSID(rhsIDs_[0]) );
885  CHK_ERR( lsc->copyOutRHSVector(1.0, *soln_fei_vector_) );
886  }
887 
888  Data tmp, tmpv;
889  for(unsigned i=0; i<matrixIDs_.size(); ++i) {
890  fei::SharedPtr<LinearSystemCore> lsci = wrapper_[i]->getLinearSystemCore();
891  if (lsci.get() == NULL) return(-1);
892 
893  if (i==0) {
894  CHK_ERR( lsci->copyInMatrix(matScalars_[i], *soln_fei_matrix_) );
895  }
896  else {
897  CHK_ERR( lsci->getMatrixPtr(tmp) );
898  CHK_ERR( lsc->sumInMatrix(matScalars_[i], tmp) );
899  }
900  }
901 
902  int last_mat = matrixIDs_.size() - 1;
903  for(unsigned i=0; i<rhsIDs_.size(); ++i) {
904  fei::SharedPtr<LinearSystemCore> lsci = (int)i<last_mat ?
905  wrapper_[i]->getLinearSystemCore() : wrapper_[last_mat]->getLinearSystemCore();
906 
907  if (i==0) {
908  CHK_ERR( lsci->copyInRHSVector(rhsScalars_[i], *soln_fei_vector_) );
909  }
910  else {
911  CHK_ERR( lsci->setRHSID(rhsIDs_[i]) );
912  CHK_ERR( lsci->getRHSVectorPtr(tmpv) );
913  CHK_ERR( lsc->sumInRHSVector(rhsScalars_[i], tmpv) );
914  }
915  }
916 
917  return(0);
918 }
919 
920 int fei::FEI_Impl::residualNorm(int whichNorm,
921  int numFields,
922  int* fieldIDs,
923  double* norms)
924 {
925  CHK_ERR( loadComplete() );
926 
928  matGraph_->getRowSpace();
929  int numLocalEqns = rowspace->getNumIndices_Owned();
930 
931  std::vector<double> residValues(numLocalEqns);
932  double* residValuesPtr = &residValues[0];
933 
934  std::vector<int> globalEqnOffsets;
935  rowspace->getGlobalIndexOffsets(globalEqnOffsets);
936  int firstLocalOffset = globalEqnOffsets[localProc_];
937 
938  if (wrapper_[0].get() == NULL) {
939  fei::SharedPtr<fei::Vector> r = factory_[0]->createVector(rowspace);
940 
941  //form r = A*x
942  CHK_ERR( A_[index_current_]->multiply(x_.get(), r.get()) );
943 
944  //form r = b - A*x
945  CHK_ERR( r->update(1.0, b_[index_current_rhs_row_].get(), -1.0) );
946 
947  //now put the values from r into the residValues array.
948  for(int ii=0; ii<numLocalEqns; ++ii) {
949  int index = firstLocalOffset+ii;
950  CHK_ERR( r->copyOut(1, &index, &(residValuesPtr[ii]) ) );
951  }
952  }
953  else {
954  fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
955  if (linSysCore.get() != NULL) {
956  CHK_ERR( linSysCore->formResidual(residValuesPtr, numLocalEqns) );
957  }
958  else {
959  fei::console_out() << "FEI_Impl::residualNorm: warning: residualNorm not implemented"
960  << " for FiniteElementData."<<FEI_ENDL;
961  int offset = 0;
962  for(int ii=0; ii<numFields; ++ii) {
963  int fieldSize = rowspace->getFieldSize(fieldIDs[ii]);
964  for(int jj=0; jj<fieldSize; ++jj) {
965  norms[offset++] = -99.9;
966  }
967  }
968  return(0);
969  }
970  }
971 
972  int numLocalNodes = rowspace->getNumOwnedIDs(nodeIDType_);
973  std::vector<int> nodeIDs(numLocalNodes);
974  int* nodeIDsPtr = &nodeIDs[0];
975  int check;
976  CHK_ERR( rowspace->getOwnedIDs(nodeIDType_, numLocalNodes,
977  nodeIDsPtr, check) );
978  std::vector<int> indices(numLocalEqns);
979  int* indicesPtr = &indices[0];
980 
981  std::vector<double> tmpNorms(numFields, 0.0);
982 
983  for(int i=0; i<numFields; ++i) {
984  int fieldSize = rowspace->getFieldSize(fieldIDs[i]);
985 
986  CHK_ERR( rowspace->getGlobalIndices(numLocalNodes, nodeIDsPtr,
987  nodeIDType_, fieldIDs[i],
988  indicesPtr) );
989 
990  double tmp = 0.0;
991  for(int j=0; j<fieldSize*numLocalNodes; ++j) {
992  if (indicesPtr[j] < 0) {
993  continue;
994  }
995 
996  double val = std::abs(residValuesPtr[indicesPtr[j]-firstLocalOffset]);
997  switch(whichNorm) {
998  case 0:
999  if (val > tmp) tmp = val;
1000  break;
1001  case 1:
1002  tmp += val;
1003  break;
1004  case 2:
1005  tmp += val*val;
1006  break;
1007  default:
1008  fei::console_out() << "FEI_Impl::residualNorm: whichNorm="<<whichNorm<<" not recognized"
1009  << FEI_ENDL;
1010  return(-1);
1011  }
1012  }
1013  tmpNorms[i] = tmp;
1014  }
1015 
1016  std::vector<double> normsArray(numFields);
1017  for(int i=0; i<numFields; ++i) {
1018  normsArray[i] = norms[i];
1019  }
1020 
1021  switch(whichNorm) {
1022  case 0:
1023  CHK_ERR( fei::GlobalMax(comm_, tmpNorms, normsArray) );
1024  break;
1025  default:
1026  CHK_ERR( fei::GlobalSum(comm_, tmpNorms, normsArray) );
1027  }
1028 
1029  for(int i=0; i<numFields; ++i) {
1030  norms[i] = normsArray[i];
1031  }
1032 
1033  if (whichNorm == 2) {
1034  for(int ii=0; ii<numFields; ++ii) {
1035  norms[ii] = std::sqrt(norms[ii]);
1036  }
1037  }
1038 
1039  return(0);
1040 }
1041 
1042 int fei::FEI_Impl::solve(int& status)
1043 {
1044  CHK_ERR( loadComplete() );
1045 
1046  fei::SharedPtr<fei::Solver> solver = factory_[0]->createSolver();
1047 
1048  std::vector<std::string> stdstrings;
1049  fei::utils::char_ptrs_to_strings(numParams_, paramStrings_, stdstrings);
1050  fei::ParameterSet params;
1051  fei::utils::parse_strings(stdstrings, " ", params);
1052 
1053  int err = solver->solve(linSys_.get(),
1054  NULL,//preconditioningMatrix
1055  params, iterations_, status);
1056 
1057  CHK_ERR( x_->scatterToOverlap() );
1058 
1059  return(err);
1060 }
1061 
1062 int fei::FEI_Impl::iterations(int& itersTaken) const
1063 {
1064  itersTaken = iterations_;
1065 
1066  return(0);
1067 }
1068 
1069 int fei::FEI_Impl::version(const char*& versionString)
1070 {
1071  versionString = fei::utils::version();
1072 
1073  return(0);
1074 }
1075 
1077  double& loadTime,
1078  double& solveTime,
1079  double& solnReturnTime)
1080 {
1081  initTime = initTime_;
1082  loadTime = loadTime_;
1083  solveTime = solveTime_;
1084  solnReturnTime = solnReturnTime_;
1085 
1086  return(0);
1087 }
1088 
1090  int numNodes,
1091  const GlobalID *nodeIDs,
1092  int *offsets,
1093  double *results)
1094 {
1095  (void)elemBlockID;
1096  return ( getNodalSolution(numNodes, nodeIDs, offsets, results) );
1097 }
1098 
1100  const GlobalID* nodeIDs,
1101  int* offsets,
1102  double* results)
1103 {
1104  int j;
1105  int offset = 0;
1106  std::vector<int> fieldIDs;
1107  for(int i=0; i<numNodes; ++i) {
1108  offsets[i] = offset;
1109 
1110  GlobalID nodeID = nodeIDs[i];
1111 
1112  int numFields = rowSpace_->getNumFields(nodeIDType_, nodeID);
1113  iwork_.resize( numFields*2 );
1114  int* fieldSizes = &iwork_[0];
1115  rowSpace_->getFields(nodeIDType_, nodeID, fieldIDs);
1116 
1117  int numDOF = 0;
1118  for(j=0; j<numFields; ++j) {
1119  fieldSizes[j] = rowSpace_->getFieldSize(fieldIDs[j]);
1120  numDOF += fieldSizes[j];
1121  }
1122 
1123  if (!rowSpace_->isLocal(nodeIDType_, nodeID)) {
1124  offset += numDOF;
1125  continue;
1126  }
1127 
1128  int results_offset = offset;
1129  for(j=0; j<numFields; ++j) {
1130  CHK_ERR( x_->copyOutFieldData(fieldIDs[j], nodeIDType_,
1131  1, &nodeID,
1132  &(results[results_offset])));
1133 
1134  results_offset += fieldSizes[j];
1135  }
1136 
1137  offset += numDOF;
1138  }
1139 
1140  offsets[numNodes] = offset;
1141 
1142  return(0);
1143 }
1144 
1146  int fieldID,
1147  int numNodes,
1148  const GlobalID *nodeIDs,
1149  double *results)
1150 {
1151  throw std::runtime_error("FEI_Impl::getBlockFieldNodeSolution not implemented.");
1152  return(0);
1153 }
1154 
1156  int numElems,
1157  const GlobalID *elemIDs,
1158  int& numElemDOFPerElement,
1159  double *results)
1160 {
1161  std::map<int,int>::const_iterator b_iter = block_dof_per_elem_.find(elemBlockID);
1162  if (b_iter == block_dof_per_elem_.end()) return(-1);
1163  numElemDOFPerElement = (*b_iter).second;
1164 
1165  fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(elemBlockID);
1166  if (block==NULL) {
1167  FEI_OSTRINGSTREAM osstr;
1168  osstr<< "fei::FEI_Impl::getBlockElemSolution ERROR, elemBlockID "
1169  << elemBlockID << " not valid.";
1170  throw std::runtime_error(osstr.str());
1171  }
1172 
1173  fei::Pattern* pattern = block->getRowPattern();
1174 
1175  int numIDs = pattern->getNumIDs();
1176  const int* idTypes = pattern->getIDTypes();
1177  const int* fIDs= pattern->getFieldIDs();
1178  const int* numFieldsPerID = pattern->getNumFieldsPerID();
1179 
1180  std::vector<int> fieldIDs;
1181  int foffset = 0;
1182  int i;
1183  for(i=0; i<numIDs; ++i) {
1184  if (idTypes[i] == elemIDType_) {
1185  for(int j=0; j<numFieldsPerID[i]; ++j) {
1186  fieldIDs.push_back(fIDs[foffset++]);
1187  }
1188  break;
1189  }
1190  foffset += numFieldsPerID[i];
1191  }
1192 
1193  if (fieldIDs.size() < 1) {
1194  ERReturn(-1);
1195  }
1196 
1197  int offset = 0;
1198  for(i=0; i<numElems; ++i) {
1199  foffset = offset;
1200  for(size_t j=0; j<fieldIDs.size(); ++j) {
1201  int fieldSize;
1202  getFieldSize(fieldIDs[j], fieldSize);
1203 
1204  CHK_ERR( x_->copyOutFieldData(fieldIDs[j], elemIDType_, 1, &(elemIDs[i]),
1205  &(results[foffset])) );
1206  foffset += fieldSize;
1207  }
1208 
1209  offset += numElemDOFPerElement;
1210  }
1211 
1212  return(0);
1213 }
1214 
1216 {
1217  numMultCRs = rowSpace_->getNumOwnedAndSharedIDs(constraintIDType_);
1218  return(0);
1219 }
1220 
1221 int fei::FEI_Impl::getCRMultIDList(int numMultCRs, int* multIDs)
1222 {
1223  int checkNum;
1224  CHK_ERR( rowSpace_->getOwnedAndSharedIDs(constraintIDType_, numMultCRs,
1225  multIDs, checkNum) );
1226  return(0);
1227 }
1228 
1230  const int* CRIDs,
1231  double *multipliers)
1232 {
1233  iwork_.resize(numCRs);
1234 
1235  for(int i=0; i<numCRs; ++i) {
1236  CHK_ERR( rowSpace_->getGlobalIndex(constraintIDType_, CRIDs[i], iwork_[i]));
1237  }
1238 
1239  CHK_ERR( x_->copyOut(numCRs, &iwork_[0], multipliers) );
1240 
1241  return(0);
1242 }
1243 
1245  int numNodes,
1246  const GlobalID *nodeIDs,
1247  const int *offsets,
1248  const double *estimates)
1249 {
1250  throw std::runtime_error("FEI_Impl::putBlockNodeSolution not implemented.");
1251  return(0);
1252 }
1253 
1255  int fieldID,
1256  int numNodes,
1257  const GlobalID *nodeIDs,
1258  const double *estimates)
1259 {
1260  throw std::runtime_error("FEI_Impl::putBlockFieldNodeSolution not implemented.");
1261  return(0);
1262 }
1263 
1265  int numElems,
1266  const GlobalID *elemIDs,
1267  int dofPerElem,
1268  const double *estimates)
1269 {
1270  throw std::runtime_error("FEI_Impl::putBlockElemSolution not implemented.");
1271  return(0);
1272 }
1273 
1275  const int* CRIDs,
1276  const double* multEstimates)
1277 {
1278  throw std::runtime_error("FEI_Impl::putCRMultipliers not implemented.");
1279  return(0);
1280 }
1281 
1283  int numNodes,
1284  GlobalID *nodeIDs)
1285 {
1286  if (!nodeset_filled_ || elemBlockID != nodeset_blockid_) {
1287  CHK_ERR( fillNodeset(elemBlockID) );
1288  }
1289 
1290  fei::copySetToArray(nodeset_, numNodes, nodeIDs);
1291  return(0);
1292 }
1293 
1294 int fei::FEI_Impl::fillNodeset(int blockID) const
1295 {
1296  if (nodeset_filled_ && blockID == nodeset_blockid_) {
1297  return(0);
1298  }
1299 
1300  fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(blockID);
1301  if (block==NULL) {
1302  FEI_OSTRINGSTREAM osstr;
1303  osstr<< "fei::FEI_Impl::fillNodeset ERROR, blockID "
1304  << blockID << " not valid.";
1305  throw std::runtime_error(osstr.str());
1306  }
1307 
1308  int nodeType = 0;
1309  snl_fei::RecordCollection* nodeRecords = NULL;
1310  matGraph_->getRowSpace()->getRecordCollection(nodeType, nodeRecords);
1311  std::vector<int>& nodes = block->getRowConnectivities();
1312 
1313  nodeset_.clear();
1314 
1315  for(unsigned i=0; i<nodes.size(); ++i) {
1316  nodeset_.insert(nodeRecords->getRecordWithLocalID(nodes[i])->getID());
1317  }
1318 
1319  nodeset_filled_ = true;
1320  nodeset_blockid_ = blockID;
1321 
1322  return(0);
1323 }
1324 
1326  int numElems,
1327  GlobalID* elemIDs)
1328 {
1329  fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(elemBlockID);
1330  if (block==NULL) {
1331  FEI_OSTRINGSTREAM osstr;
1332  osstr<< "fei::FEI_Impl::getBlockElemIDList ERROR, elemBlockID "
1333  << elemBlockID << " not valid.";
1334  throw std::runtime_error(osstr.str());
1335  }
1336 
1337  std::map<int,int>& elemIDSet = block->getConnectivityIDs();
1338 
1339  fei::copyKeysToArray(elemIDSet, numElems, elemIDs);
1340 
1341  return(0);
1342 }
1343 
1344 int fei::FEI_Impl::getNumSolnParams(GlobalID nodeID, int& numSolnParams) const
1345 {
1346  numSolnParams = rowSpace_->getNumDegreesOfFreedom( nodeIDType_, nodeID);
1347  return(0);
1348 }
1349 
1350 int fei::FEI_Impl::getNumElemBlocks(int& numElemBlocks) const
1351 {
1352  numElemBlocks = matGraph_->getNumConnectivityBlocks();
1353  return(0);
1354 }
1355 
1356 int fei::FEI_Impl::getNumBlockActNodes(GlobalID blockID, int& numNodes) const
1357 {
1358  if (!nodeset_filled_ || blockID != nodeset_blockid_) {
1359  CHK_ERR( fillNodeset(blockID) );
1360  }
1361 
1362  numNodes = nodeset_.size();
1363  return(0);
1364 }
1365 
1366 int fei::FEI_Impl::getNumBlockActEqns(GlobalID blockID, int& numEqns) const
1367 {
1368  throw std::runtime_error("fei::FEI_Impl::getNumBlockActEqns not implemented.");
1369  return(0);
1370 }
1371 
1372 int fei::FEI_Impl::getNumNodesPerElement(GlobalID blockID, int& nodesPerElem) const
1373 {
1374  nodesPerElem = matGraph_->getNumIDsPerConnectivityList(blockID);
1375  return(0);
1376 }
1377 
1378 int fei::FEI_Impl::getNumEqnsPerElement(GlobalID blockID, int& numEqns) const
1379 {
1380  numEqns = matGraph_->getConnectivityNumIndices(blockID);
1381  return(0);
1382 }
1383 
1384 int fei::FEI_Impl::getNumBlockElements(GlobalID blockID, int& numElems) const
1385 {
1386  const fei::ConnectivityBlock* cblock =
1387  matGraph_->getConnectivityBlock(blockID);
1388  if (cblock != NULL) {
1389  numElems = cblock->getConnectivityIDs().size();
1390  }
1391  else {
1392  numElems = 0;
1393  }
1394  return(0);
1395 }
1396 
1397 int fei::FEI_Impl::getNumBlockElemDOF(GlobalID blockID, int& DOFPerElem) const
1398 {
1399  std::map<int,int>::const_iterator b_iter = block_dof_per_elem_.find(blockID);
1400  if (b_iter == block_dof_per_elem_.end()) DOFPerElem = 0;
1401  else DOFPerElem = (*b_iter).second;
1402 
1403  return(0);
1404 }
1405 
1406 int fei::FEI_Impl::getParameters(int& numParams, char**& paramStrings)
1407 {
1408  numParams = numParams_;
1409  paramStrings = paramStrings_;
1410  return(0);
1411 }
1412 
1413 int fei::FEI_Impl::getFieldSize(int fieldID, int& numScalars)
1414 {
1415  numScalars = rowSpace_->getFieldSize(fieldID);
1416  return(0);
1417 }
1418 
1420  int idType,
1421  int fieldID,
1422  int& numEqns,
1423  int* eqnNumbers)
1424 {
1425  numEqns = rowSpace_->getFieldSize(fieldID);
1426  CHK_ERR( rowSpace_->getGlobalIndices(1, &ID, idType, fieldID, eqnNumbers) );
1427  return(0);
1428 }
1429 
1431  int numNodes,
1432  const GlobalID* nodeIDs,
1433  double* results)
1434 {
1435  CHK_ERR( x_->copyOutFieldData(fieldID, nodeIDType_, numNodes,
1436  nodeIDs, results) );
1437  return(0);
1438 }
1439 
1441 {
1442  numNodes = rowSpace_->getNumOwnedAndSharedIDs(nodeIDType_);
1443  return(0);
1444 }
1445 
1447  GlobalID* nodeIDs,
1448  int lenNodeIDs)
1449 {
1450  CHK_ERR( rowSpace_->getOwnedAndSharedIDs(nodeIDType_, lenNodeIDs, nodeIDs, numNodes) );
1451  return(0);
1452 }
1453 
1455  int numNodes,
1456  const GlobalID* nodeIDs,
1457  const double* nodeData)
1458 {
1459  int err = 0;
1460  if (fieldID < 0) {
1461  bool data_passed = false;
1462  if (wrapper_[0].get() != NULL) {
1463  std::vector<int> numbers(numNodes);
1464  for(int i=0; i<numNodes; ++i) {
1465  err = rowSpace_->getGlobalBlkIndex(nodeIDType_, nodeIDs[i], numbers[i]);
1466  if (err != 0) {
1467  fei::console_out() << "fei::FEI_Impl::putNodalFieldData ERROR, nodeID "
1468  << nodeIDs[i] << " not found."<<FEI_ENDL;
1469  ERReturn(-1);
1470  }
1471  }
1472 
1473  int fieldSize = 0;
1474  try {
1475  fieldSize = rowSpace_->getFieldSize(fieldID);
1476  }
1477  catch (std::runtime_error& exc) {
1478  fei::console_out() << "fei::FEI_Impl::putNodalFieldData ERROR: " <<exc.what()<<FEI_ENDL;
1479  ERReturn(-1);
1480  }
1481 
1482  fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
1483  if (linSysCore.get() != NULL) {
1484  linSysCore->putNodalFieldData(fieldID, fieldSize,
1485  &numbers[0], numNodes, nodeData);
1486  data_passed = true;
1487  }
1488  else {
1489  //If we enter this block, we're probably dealing with a FiniteElementData
1490  //instance.
1491  fei::SharedPtr<FiniteElementData> fedata = wrapper_[0]->getFiniteElementData();
1492  if (fedata.get() != NULL) {
1493  fedata->putNodalFieldData(fieldID, fieldSize, numNodes,
1494  &numbers[0], nodeData);
1495  data_passed = true;
1496  }
1497  }
1498  }
1499 
1500  if (!data_passed) {
1501  //If we get to here and data_passed is false, wrapper_[0] is probably NULL.
1502  if (wrapper_[0].get() == NULL) {
1503  fei::SharedPtr<fei::Vector> dataVector =factory_[0]->createVector(matGraph_);
1504 
1505  CHK_ERR( dataVector->copyInFieldData(fieldID, nodeIDType_,
1506  numNodes, nodeIDs, nodeData) );
1507  if (fieldID == -3) {
1508  CHK_ERR( linSys_->putAttribute("coordinates", dataVector.get()) );
1509  }
1510  else {
1511  FEI_OSTRINGSTREAM osstr;
1512  osstr << "fieldID:" << fieldID;
1513  CHK_ERR( linSys_->putAttribute(osstr.str().c_str(), dataVector.get()) );
1514  }
1515  }
1516  else {
1517  fei::console_out() << "fei::FEI_Impl::putNodalFieldData ERROR, non-null LibraryWrapper"
1518  << " contains neither LinearSystemCore or FiniteElementData. " <<FEI_ENDL;
1519  ERReturn(-1);
1520  }
1521  }
1522  return(0);
1523  }
1524 
1525  CHK_ERR( x_->copyInFieldData(fieldID, nodeIDType_,
1526  numNodes, nodeIDs, nodeData));
1527  return(0);
1528 }
virtual int setRHSID(int rhsID)=0
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int setSolveType(int solveType)
int sortedListInsert(const T &item, std::vector< T > &list)
int getNumLocalNodes(int &numNodes)
int putNodalFieldData(int fieldID, int numNodes, const GlobalID *nodeIDs, const double *nodeData)
void char_ptrs_to_strings(int numStrings, const char *const *charstrings, std::vector< std::string > &stdstrings)
Definition: fei_utils.cpp:164
int sumIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *coefficients)
int getNumCRMultipliers(int &numMultCRs)
const int * getIDTypes() const
Definition: fei_Pattern.hpp:67
const std::map< int, int > & getConnectivityIDs() const
int iterations(int &itersTaken) const
int version(const char *&versionString)
#define FEI_AGGREGATE_SUM
Definition: fei_defs.h:67
int fillNodeset(int blockID) const
const int * getNumFieldsPerID() const
Definition: fei_Pattern.hpp:73
int GlobalID
Definition: fei_defs.h:60
int setIDLists(int numMatrices, const int *matrixIDs, int numRHSs, const int *rhsIDs)
fei::SharedPtr< fei::LinearSystem > getLinearSystem()
virtual int copyOutRHSVector(double scalar, Data &data)=0
int getBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int &numElemDOFPerElement, double *results)
int getParameters(int &numParams, char **&paramStrings)
int getOwnedIDs(int idtype, int lenList, int *IDs, int &numLocalIDs)
int initCRPen(int numCRNodes, const GlobalID *CRNodeIDs, const int *CRFieldIDs, int &CRID)
int aggregateSystem_LinSysCore()
virtual int putNodalFieldData(int fieldID, int fieldSize, int *nodeNumbers, int numNodes, const double *data)=0
int initFields(int numFields, const int *fieldSizes, const int *fieldIDs, const int *fieldTypes=NULL)
virtual ~FEI_Impl()
const int * getFieldIDs() const
Definition: fei_Pattern.hpp:76
virtual int sumInRHSVector(double scalar, const Data &data)=0
void copySetToArray(const SET_TYPE &set_obj, int lenList, int *list)
int getFieldSize(int fieldID, int &numScalars)
int getNumSolnParams(GlobalID nodeID, int &numSolnParams) const
int getBlockNodeIDList(GlobalID elemBlockID, int numNodes, GlobalID *nodeIDs)
fei::SharedPtr< LibraryWrapper > get_LibraryWrapper() const
static unsigned getFieldSize(int fieldID, fei::VectorSpace *space1, fei::VectorSpace *space2)
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
int putBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, const double *estimates)
int setCurrentMatrix(int matrixID)
int getCRMultipliers(int numCRs, const int *CRIDs, double *multipliers)
int cumulative_cpu_times(double &initTime, double &loadTime, double &solveTime, double &solnReturnTime)
virtual int getMatrixPtr(Data &data)=0
int mergeStringLists(char **&strings, int &numStrings, const char *const *stringsToMerge, int numStringsToMerge)
void reset(T *p=0)
int setCurrentRHS(int rhsID)
void copyKeysToArray(const MAP_TYPE &map_obj, unsigned lenList, int *list)
virtual int update(double a, const fei::Vector *x, double b)=0
#define MPI_Comm
Definition: fei_mpi.h:56
int putBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, const int *offsets, const double *estimates)
std::vector< int > & getRowConnectivities()
int parameters(int numParams, const char *const *paramStrings)
int resetSystem(double s=0.0)
int getGlobalIndices(int numIDs, const int *IDs, int idType, int fieldID, int *globalIndices)
int loadComplete(bool applyBCs=true, bool globalAssemble=true)
virtual int copyOutMatrix(double scalar, Data &data)=0
virtual int putNodalFieldData(int fieldID, int fieldSize, int numNodes, const int *nodeNumbers, const double *coefs)=0
virtual LinearSystemCore * clone()=0
int getBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
int getEqnNumbers(GlobalID ID, int idType, int fieldID, int &numEqns, int *eqnNumbers)
int getNodalSolution(int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
int initElemBlock(GlobalID elemBlockID, int numElements, int numNodesPerElement, const int *numFieldsPerNode, const int *const *nodalFieldIDs, int numElemDofFieldsPerElement, const int *elemDOFFieldIDs, int interleaveStrategy)
int putCRMultipliers(int numMultCRs, const int *CRIDs, const double *multEstimates)
int getCRMultIDList(int numMultCRs, int *multIDs)
int getNumNodesPerElement(GlobalID blockID, int &nodesPerElem) const
int getNumOwnedIDs(int idType)
int sumInElem(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, const double *elemLoad, int elemFormat)
int initSharedNodes(int numSharedNodes, const GlobalID *sharedNodeIDs, const int *numProcsPerNode, const int *const *sharingProcIDs)
std::vector< fei::SharedPtr< LibraryWrapper > > wrapper_
#define ERReturn(a)
T * get() const
int getNumBlockElemDOF(GlobalID blockID, int &DOFPerElem) const
int sumInElemMatrix(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, int elemFormat)
int putBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int dofPerElem, const double *estimates)
void basic_initializations()
int resetInitialGuess(double s=0.0)
int resetRHSVector(double s=0.0)
virtual int copyInFieldData(int fieldID, int idType, int numIDs, const int *IDs, const double *data, int vectorIndex=0)=0
int setMatScalars(int numScalars, const int *IDs, const double *scalars)
#define voidERReturn
virtual int getRHSVectorPtr(Data &data)=0
virtual fei::SharedPtr< fei::Vector > createVector(fei::SharedPtr< fei::VectorSpace > vecSpace, int numVectors=1)=0
virtual int copyInMatrix(double scalar, const Data &data)=0
virtual fei::SharedPtr< Factory > clone() const =0
int initSlaveVariable(GlobalID slaveNodeID, int slaveFieldID, int offsetIntoSlaveField, int numMasterNodes, const GlobalID *masterNodeIDs, const int *masterFieldIDs, const double *weights, double rhsValue)
int getNumElemBlocks(int &numElemBlocks) const
#define FEI_ENDL
virtual int setNumRHSVectors(int numRHSs, const int *rhsIDs)=0
std::ostream & console_out()
fei::Record< int > * getRecordWithLocalID(int lid)
int getNumBlockActNodes(GlobalID blockID, int &numNodes) const
void getGlobalIndexOffsets(std::vector< int > &globalOffsets) const
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
virtual int destroyVectorData(Data &data)=0
int getNumIDs() const
Definition: fei_Pattern.hpp:64
int residualNorm(int whichNorm, int numFields, int *fieldIDs, double *norms)
int solve(int &status)
int inputRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *rhsEntries, bool sumInto)
int putIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *coefficients)
int localProc(MPI_Comm comm)
int resetMatrix(double s=0.0)
const char * version()
Definition: fei_utils.hpp:53
std::vector< fei::SharedPtr< fei::Factory > > factory_
int loadCRPen(int CRPenID, int numCRNodes, const GlobalID *CRNodeIDs, const int *CRFieldIDs, const double *CRWeights, double CRValue, double penValue)
int getNumBlockActEqns(GlobalID blockID, int &numEqns) const
virtual int solve(fei::LinearSystem *linearSystem, fei::Matrix *preconditioningMatrix, const fei::ParameterSet &parameterSet, int &iterationsTaken, int &status)
Definition: fei_Solver.cpp:65
int initCRMult(int numCRNodes, const GlobalID *CRNodeIDs, const int *CRFieldIDs, int &CRID)
virtual int sumInMatrix(double scalar, const Data &data)=0
int loadCRMult(int CRMultID, int numCRNodes, const GlobalID *CRNodeIDs, const int *CRFieldIDs, const double *CRWeights, double CRValue)
unsigned getFieldSize(int fieldID)
int getNumEqnsPerElement(GlobalID blockID, int &numEqns) const
int getNumBlockElements(GlobalID blockID, int &numElems) const
int initElem(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn)
int getBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
int loadNodeBCs(int numNodes, const GlobalID *nodeIDs, int fieldID, const int *offsetsIntoField, const double *prescribedValues)
int setRHSScalars(int numScalars, const int *IDs, const double *scalars)
int getLocalNodeIDList(int &numNodes, GlobalID *nodeIDs, int lenNodeIDs)
int loadElemBCs(int numElems, const GlobalID *elemIDs, int fieldID, const double *const *alpha, const double *const *beta, const double *const *gamma)
virtual fei::SharedPtr< fei::Matrix > createMatrix(fei::SharedPtr< fei::MatrixGraph > matrixGraph)=0
virtual int formResidual(double *values, int len)=0
#define CHK_ERR(a)
int getBlockElemIDList(GlobalID elemBlockID, int numElems, GlobalID *elemIDs)
int getNodalFieldSolution(int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
#define FEI_OSTRINGSTREAM
Definition: fei_sstream.hpp:32
int sumInElemRHS(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *elemLoad)
virtual int copyOut(int numValues, const int *indices, double *values, int vectorIndex=0) const =0
const fei::Pattern * getRowPattern() const
virtual int destroyMatrixData(Data &data)=0
virtual int copyInRHSVector(double scalar, const Data &data)=0
int numProcs(MPI_Comm comm)
int getNumIndices_Owned() const
FEI_Impl(fei::SharedPtr< LibraryWrapper > wrapper, MPI_Comm comm, int masterRank=0)