9 #include <fei_CommUtils.hpp>
10 #include <fei_iostream.hpp>
11 #include <fei_fstream.hpp>
12 #include <fei_sstream.hpp>
14 #include <fei_utils.hpp>
16 #include <fei_Data.hpp>
17 #include <fei_LinearSystemCore.hpp>
18 #include <fei_FiniteElementData.hpp>
20 #include <fei_LibraryWrapper.hpp>
21 #include <fei_Filter.hpp>
23 #include <fei_LinSysCoreFilter.hpp>
24 #include <fei_FEDataFilter.hpp>
26 #include <SNL_FEI_Structure.hpp>
27 #include <fei_BlockDescriptor.hpp>
28 #include <fei_NodeDatabase.hpp>
29 #include <fei_ConnectivityTable.hpp>
30 #include <FEI_Implementation.hpp>
31 #include <snl_fei_Utils.hpp>
34 #define fei_file "FEI_Implementation.cpp"
36 #include <fei_ErrMacros.hpp>
40 MPI_Comm comm,
int masterRank)
41 : wrapper_(libWrapper),
44 haveLinSysCore_(false),
46 problemStructure_(NULL),
49 internalFEIsAllocated_(false),
55 matScalarsSet_(false),
57 rhsScalarsSet_(false),
58 index_soln_filter_(0),
59 index_current_filter_(0),
60 index_current_rhs_row_(0),
62 setSolveTypeCalled_(false),
63 initPhaseIsComplete_(false),
64 aggregateSystemFormed_(false),
65 newMatrixDataLoaded_(0),
66 soln_fei_matrix_(NULL),
67 soln_fei_vector_(NULL),
76 dbgFileOpened_(false),
90 if (problemStructure_ == NULL) {
91 messageAbort(
"problemStructure allocation failed");
98 haveFEData_ = wrapper_->haveFiniteElementData();
101 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
102 nodeCommMgr.setSharedOwnershipRule(NodeCommMgr::PROC_WITH_LOCAL_ELEM);
105 haveLinSysCore_ = wrapper_->haveLinearSystemCore();
106 if (haveLinSysCore_) {
107 linSysCore_ = wrapper_->getLinearSystemCore();
108 lscArray_.push_back(linSysCore_);
111 numInternalFEIs_ = 1;
112 matrixIDs_.resize(1);
114 numRHSIDs_.resize(1);
117 rhsIDs_[0] =
new int[1];
119 rhsScalars_.resize(numInternalFEIs_);
120 for(
int ii=0; ii<numInternalFEIs_; ii++) rhsScalars_[ii] = NULL;
135 (*dbgOStreamPtr_) <<
"FEI: destructor" << FEI_ENDL;
138 if (soln_fei_matrix_) {
140 delete soln_fei_matrix_;
141 soln_fei_matrix_ = NULL;
144 if (soln_fei_vector_) {
146 delete soln_fei_vector_;
147 soln_fei_vector_ = NULL;
152 if (internalFEIsAllocated_) {
153 for(
int j=0; j<numInternalFEIs_; j++){
161 internalFEIsAllocated_ =
false;
162 numInternalFEIs_ = 0;
164 delete problemStructure_;
166 for(
int k=0; k<numParams_; k++)
delete [] paramStrings_[k];
167 delete [] paramStrings_;
169 if (dbgFileOpened_ ==
true) {
170 dbgFStreamPtr_->close();
delete dbgFStreamPtr_;
172 else delete dbgOStreamPtr_;
177 void FEI_Implementation::deleteIDs()
179 matrixIDs_.resize(0);
180 for(
size_t i=0; i<rhsIDs_.size(); i++) {
181 delete [] rhsIDs_[i];
184 numRHSIDs_.resize(0);
188 void FEI_Implementation::deleteRHSScalars()
190 for(
size_t i=0; i<rhsScalars_.size(); i++) {
191 delete [] rhsScalars_[i];
193 rhsScalars_.resize(0);
200 (*dbgOStreamPtr_) <<
"FEI: setCurrentMatrix" << FEI_ENDL <<
"#matrix-id"
201 << FEI_ENDL<<matID<<FEI_ENDL;
204 index_current_filter_ = -1;
206 for(
int i=0; i<numInternalFEIs_; i++){
207 if (matrixIDs_[i] == matID) index_current_filter_ = i;
211 (*dbgOStreamPtr_) <<
"#--- ID: " << matID
212 <<
", ind: "<<index_current_filter_<<FEI_ENDL;
216 if (index_current_filter_ == -1) {
217 fei::console_out() <<
"FEI_Implementation::setCurrentMatrix: ERROR, invalid matrix ID "
218 <<
"supplied" << FEI_ENDL;
222 debugOut(
"#FEI_Implementation leaving setCurrentMatrix");
230 numParams = numParams_;
231 paramStrings = paramStrings_;
239 (*dbgOStreamPtr_) <<
"FEI: setCurrentRHS" << FEI_ENDL <<
"#rhs-id"
240 << FEI_ENDL<<rhsID<<FEI_ENDL;
245 for(
int j=0; j<numInternalFEIs_; j++){
248 index_current_rhs_row_ = j;
249 CHK_ERR( filter_[index_current_rhs_row_]->
setCurrentRHS(rhsID) )
256 fei::console_out() <<
"FEI_Implementation::setCurrentRHS: ERROR, invalid RHS ID"
261 debugOut(
"#FEI_Implementation leaving setCurrentRHS");
270 (*dbgOStreamPtr_)<<
"FEI: setSolveType"<<FEI_ENDL;
271 (*dbgOStreamPtr_)<<solveType<<FEI_ENDL;
274 solveType_ = solveType;
276 if (solveType_ == FEI_SINGLE_SYSTEM) {
277 if (matrixIDs_.size() > 1) {
278 messageAbort(
"setSolveType: solve-type is FEI_SINGLE_SYSTEM, but setIDLists() has been called with numMatrices > 1.");
281 else if (solveType_ == FEI_EIGEN_SOLVE) {
283 else if (solveType_ == FEI_AGGREGATE_SUM) {
287 else if (solveType_ == FEI_AGGREGATE_PRODUCT) {
291 else if (solveType_ == 4) {
300 int numRHSs,
const int* rhsIDs)
303 FEI_OSTREAM& os = *dbgOStreamPtr_;
304 os <<
"FEI: setIDLists" << FEI_ENDL
305 <<
"#num-matrices" << FEI_ENDL << numMatrices << FEI_ENDL
306 <<
"#matrixIDs" << FEI_ENDL;
308 for(i=0; i<numMatrices; ++i) os << matrixIDs[i] <<
" ";
309 os << FEI_ENDL <<
"#num-rhs's" << FEI_ENDL;
310 for(i=0; i<numRHSs; ++i) os << rhsIDs[i] <<
" ";
321 int myNumMatrices = numMatrices;
322 if (myNumMatrices == 0) myNumMatrices = 1;
324 matrixIDs_.resize(myNumMatrices);
326 if (rhsScalars_.size() != 0) deleteRHSScalars();
328 numInternalFEIs_ = myNumMatrices;
330 if (numMatrices == 0) {
334 for(
int i=0; i<numMatrices; i++) matrixIDs_[i] = matrixIDs[i];
337 int quotient = numRHSs/myNumMatrices;
338 int rem = numRHSs%numMatrices;
346 numRHSIDs_.resize(myNumMatrices);
347 rhsIDs_.resize(myNumMatrices);
350 for(
int i=0; i<myNumMatrices; i++) {
351 numRHSIDs_[i] = quotient;
352 if (i < rem) numRHSIDs_[i]++;
354 rhsIDs_[i] = numRHSIDs_[i] > 0 ?
new int[numRHSIDs_[i]] : NULL ;
356 for(
int j=0; j<numRHSIDs_[i]; j++) {
357 rhsIDs_[i][j] = rhsIDs[offset+j];
360 offset += numRHSIDs_[i];
368 const int *fieldSizes,
370 const int *fieldTypes)
372 CHK_ERR( problemStructure_->initFields(numFields, fieldSizes, fieldIDs, fieldTypes) );
380 int numNodesPerElement,
381 const int* numFieldsPerNode,
382 const int*
const* nodalFieldIDs,
383 int numElemDofFieldsPerElement,
384 const int* elemDOFFieldIDs,
385 int interleaveStrategy)
387 CHK_ERR( problemStructure_->initElemBlock(elemBlockID,
392 numElemDofFieldsPerElement,
394 interleaveStrategy) );
402 const GlobalID* elemConn)
404 CHK_ERR( problemStructure_->initElem(elemBlockID, elemID, elemConn) );
412 int offsetIntoSlaveField,
414 const GlobalID* masterNodeIDs,
415 const int* masterFieldIDs,
416 const double* weights,
419 CHK_ERR( problemStructure_->initSlaveVariable(slaveNodeID, slaveFieldID,
420 offsetIntoSlaveField,
421 numMasterNodes, masterNodeIDs,
422 masterFieldIDs, weights, rhsValue));
430 debugOut(
"FEI: deleteMultCRs");
432 CHK_ERR( problemStructure_->deleteMultCRs() );
435 if (internalFEIsAllocated_) {
436 err = filter_[index_current_filter_]->deleteMultCRs();
444 const GlobalID *sharedNodeIDs,
445 const int* numProcsPerNode,
446 const int *
const *sharingProcIDs)
452 CHK_ERR( problemStructure_->initSharedNodes(numSharedNodes,
462 const GlobalID* CRNodes,
472 CHK_ERR( problemStructure_->initCRMult(numCRNodes,
482 const GlobalID* CRNodes,
491 CHK_ERR( problemStructure_->initCRPen(numCRNodes,
502 bool generateGraph = !haveFEData_;
504 CHK_ERR( problemStructure_->initComplete(generateGraph) );
510 CHK_ERR( allocateInternalFEIs() );
512 for(
int i=0; i<numInternalFEIs_; ++i) {
513 CHK_ERR( filter_[i]->initialize() );
516 problemStructure_->destroyMatIndices();
518 initPhaseIsComplete_ =
true;
528 if (!internalFEIsAllocated_)
return(0);
530 CHK_ERR( filter_[index_current_filter_]->
resetSystem(s))
538 if (!internalFEIsAllocated_)
return(0);
540 CHK_ERR( filter_[index_current_filter_]->
resetMatrix(s))
548 if (!internalFEIsAllocated_)
return(0);
558 if (!internalFEIsAllocated_)
return(0);
567 const GlobalID *nodeIDs,
569 const int* offsetsIntoField,
570 const double* prescribedValues)
572 if (!internalFEIsAllocated_)
573 notAllocatedAbort(
"FEI_Implementation::loadNodeBCs");
575 int index = index_current_filter_;
576 if (solveType_ == 2) index = index_soln_filter_;
580 offsetsIntoField, prescribedValues));
587 const GlobalID *elemIDs,
589 const double *
const *alpha,
590 const double *
const *beta,
591 const double *
const *gamma)
593 if (!internalFEIsAllocated_)
594 notAllocatedAbort(
"FEI_Implementation::loadElemBCs");
596 int index = index_current_filter_;
597 if (solveType_ == 2) index = index_soln_filter_;
609 const GlobalID* elemConn,
610 const double*
const* elemStiffness,
611 const double* elemLoad,
614 if (!internalFEIsAllocated_) {
615 notAllocatedAbort(
"FEI_Implementation::sumInElem");
618 CHK_ERR( filter_[index_current_filter_]->
sumInElem(elemBlockID, elemID,
619 elemConn, elemStiffness,
620 elemLoad, elemFormat));
622 newMatrixDataLoaded_ = 1;
630 const GlobalID* elemConn,
631 const double*
const* elemStiffness,
634 if (!internalFEIsAllocated_)
635 notAllocatedAbort(
"FEI_Implementation::sumInElemMatrix");
639 elemStiffness, elemFormat))
641 newMatrixDataLoaded_ = 1;
650 const GlobalID* elemConn,
651 const double* elemLoad)
653 if (!internalFEIsAllocated_)
654 notAllocatedAbort(
"FEI_Implementation::sumInElemRHS");
656 CHK_ERR( filter_[index_current_rhs_row_]->
sumInElemRHS(elemBlockID,
657 elemID, elemConn, elemLoad))
659 newMatrixDataLoaded_ = 1;
667 const GlobalID* CRNodes,
669 const double* CRWeights,
672 if (!internalFEIsAllocated_)
673 notAllocatedAbort(
"FEI_Implementation::loadCRMult");
675 newMatrixDataLoaded_ = 1;
677 CHK_ERR( filter_[index_current_filter_]->
loadCRMult(CRID,
679 CRFields, CRWeights, CRValue));
687 const GlobalID* CRNodes,
689 const double* CRWeights,
693 if (!internalFEIsAllocated_)
694 notAllocatedAbort(
"FEI_Implementation::loadCRPen");
696 CHK_ERR( filter_[index_current_filter_]->
loadCRPen(CRID,
701 newMatrixDataLoaded_ = 1;
711 const double* rhsEntries)
713 if (!internalFEIsAllocated_)
714 notAllocatedAbort(
"FEI_Implementation::sumIntoRHS");
716 CHK_ERR( filter_[index_current_rhs_row_]->
sumIntoRHS(IDType, fieldID,
717 numIDs, IDs, rhsEntries) );
718 newMatrixDataLoaded_ = 1;
724 int FEI_Implementation::sumIntoMatrixDiagonal(
int IDType,
728 const double* coefficients)
730 if (!internalFEIsAllocated_)
731 notAllocatedAbort(
"FEI_Implementation::sumIntoMatrixDiagonal");
733 CHK_ERR( filter_[index_current_filter_]->sumIntoMatrixDiagonal(IDType, fieldID,
734 numIDs, IDs, coefficients) );
735 newMatrixDataLoaded_ = 1;
745 const double* rhsEntries)
747 if (!internalFEIsAllocated_)
748 notAllocatedAbort(
"FEI_Implementation::putIntoRHS");
750 CHK_ERR( filter_[index_current_rhs_row_]->
putIntoRHS(IDType, fieldID,
751 numIDs, IDs, rhsEntries) );
752 newMatrixDataLoaded_ = 1;
759 const int* IDs,
const double* scalars)
761 for(
int i=0; i<numScalars; i++){
762 std::vector<int>::iterator iter =
763 std::find(matrixIDs_.begin(), matrixIDs_.end(), IDs[i]);
764 if (iter != matrixIDs_.end()) {
765 int index = iter - matrixIDs_.begin();
766 matScalars_[index] = scalars[i];
769 fei::console_out() <<
"FEI_Implementation::setMatScalars: ERROR, invalid ID supplied"
775 aggregateSystemFormed_ =
false;
776 matScalarsSet_ =
true;
783 const int* IDs,
const double* scalars)
785 for(
int i=0; i<numScalars; i++){
788 for(
int j=0; j<numInternalFEIs_; j++){
791 rhsScalars_[j][index] = scalars[i];
798 fei::console_out() <<
"FEI_Implementation::setRHSScalars: ERROR, invalid RHS ID supplied"
804 aggregateSystemFormed_ =
false;
805 rhsScalarsSet_ =
true;
816 if (numParams == 0 || paramStrings == NULL) {
817 debugOut(
"#--- no parameters");
822 snl_fei::mergeStringLists(paramStrings_, numParams_,
823 paramStrings, numParams);
825 snl_fei::getIntParamValue(
"numMatrices", numParams,paramStrings, numInternalFEIs_);
827 snl_fei::getIntParamValue(
"outputLevel", numParams,paramStrings, outputLevel_);
829 const char* param = snl_fei::getParamValue(
"debugOutput",numParams,paramStrings);
831 setDebugOutput(param,
"FEI_log");
835 (*dbgOStreamPtr_)<<
"FEI: parameters"<<FEI_ENDL;
836 (*dbgOStreamPtr_)<<
"#FEI_Implementation, num-params "<<FEI_ENDL
837 <<numParams<<FEI_ENDL;
838 (*dbgOStreamPtr_)<<
"# "<<numParams<<
" parameter lines follow:"<<FEI_ENDL;
839 for(
int i=0; i<numParams; i++){
840 (*dbgOStreamPtr_)<<paramStrings[i]<<FEI_ENDL;
844 if (haveLinSysCore_) {
845 linSysCore_->
parameters(numParams, (
char**)paramStrings);
848 wrapper_->getFiniteElementData()->parameters(numParams, (
char**)paramStrings);
851 problemStructure_->
parameters(numParams, paramStrings);
853 if (internalFEIsAllocated_){
854 for(
int i=0; i<numInternalFEIs_; i++){
855 CHK_ERR( filter_[i]->
parameters(numParams, paramStrings) );
859 debugOut(
"#FEI_Implementation leaving parameters method");
865 void FEI_Implementation::setDebugOutput(
const char* path,
const char* name)
870 if (dbgFileOpened_) {
871 dbgFStreamPtr_->close();
874 dbgFileOpened_ =
false;
875 delete dbgOStreamPtr_;
877 FEI_OSTRINGSTREAM osstr;
879 osstr << path <<
"/";
881 osstr << name <<
"." << numProcs_ <<
"." << localRank_;
884 dbgFStreamPtr_ =
new FEI_OFSTREAM(osstr.str().c_str(), IOS_APP);
885 if (!dbgFStreamPtr_ || dbgFStreamPtr_->bad()){
886 fei::console_out() <<
"couldn't open debug output file: " << osstr.str() << FEI_ENDL;
891 const char* version_str = NULL;
894 (*dbgFStreamPtr_) << version_str << FEI_ENDL;
896 problemStructure_->setDbgOut(*dbgFStreamPtr_, path,
"_0");
897 dbgOStreamPtr_ = dbgFStreamPtr_;
898 dbgOStreamPtr_->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
899 dbgFileOpened_ =
true;
901 if (internalFEIsAllocated_) {
902 for(
int i=0; i<numInternalFEIs_; ++i) {
903 filter_[i]->setLogStream(dbgOStreamPtr_);
914 (void)globalAssemble;
923 int* fieldIDs,
double* norms)
927 double residTime = 0.0;
929 int err = filter_[index_soln_filter_]->residualNorm(whichNorm, numFields,
930 fieldIDs, norms, residTime);
932 solveTime_ += residTime;
944 int err = filter_[index_soln_filter_]->solve(status, sTime);
953 itersTaken = filter_[index_soln_filter_]->iterations();
968 double& solnReturnTime)
970 initTime = initTime_;
971 loadTime = loadTime_;
972 solveTime = solveTime_;
973 solnReturnTime = solnReturnTime_;
981 const GlobalID *nodeIDs,
995 const GlobalID *nodeIDs,
1011 const GlobalID *nodeIDs,
1024 const GlobalID *nodeIDs,
1026 const double *estimates)
1030 offsets, estimates))
1038 const GlobalID *nodeIDs,
1039 const double *estimates)
1041 int err = filter_[index_soln_filter_]->putBlockFieldNodeSolution(elemBlockID,
1043 nodeIDs, estimates);
1050 const GlobalID *elemIDs,
1051 int& numElemDOFPerElement,
1056 numElemDOFPerElement,
1064 const GlobalID *elemIDs,
1066 const double *estimates)
1070 dofPerElem, estimates))
1077 numMultCRs = problemStructure_->getNumMultConstRecords();
1085 if (numMultCRs > problemStructure_->getNumMultConstRecords())
return(-1);
1087 std::map<GlobalID,snl_fei::Constraint<GlobalID>*>::const_iterator
1088 cr_iter = problemStructure_->getMultConstRecords().begin(),
1089 cr_end = problemStructure_->getMultConstRecords().end();
1091 while(cr_iter != cr_end) {
1092 multIDs[i++] = (*cr_iter).first;
1102 double* multipliers)
1105 CRIDs, multipliers))
1112 const double *multEstimates)
1138 getBlockConnectivity(elemBlockID);
1140 std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
1141 int len = elemIDList.size();
1142 if (len > numElems) len = numElems;
1146 return(FEI_SUCCESS);
1158 int numActiveNodes = problemStructure_->getNumActiveNodes();
1159 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1163 for(
int i=0; i<numActiveNodes; i++) {
1166 if (node==NULL)
continue;
1167 if (node->hasBlockIndex(blk_idx))
1168 nodeIDs[offset++] = node->getGlobalNodeID();
1169 if (offset == numNodes)
break;
1172 return(FEI_SUCCESS);
1177 int& nodesPerElem)
const
1183 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1185 nodesPerElem = block->getNumNodesPerElement();
1186 return(FEI_SUCCESS);
1197 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1199 numEqns = block->getNumEqnsPerElement();
1200 return(FEI_SUCCESS);
1210 int err = problemStructure_->getNodeDatabase().
getNodeWithID(nodeID, node);
1216 numSolnParams = node->getNumNodalDOF();
1237 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1239 numNodes = block->getNumActiveNodes();
1240 return(FEI_SUCCESS);
1251 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1253 numEqns = block->getTotalNumEqns();
1254 return(FEI_SUCCESS);
1263 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1265 numElems = block->getNumElements();
1266 return(FEI_SUCCESS);
1271 int& DOFPerElem)
const
1277 CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
1279 DOFPerElem = block->getNumElemDOFPerElement();
1281 return(FEI_SUCCESS);
1298 int fieldID,
int& numEqns,
1305 return( problemStructure_->getEqnNumbers(ID, idType, fieldID,
1306 numEqns, eqnNumbers) );
1312 const GlobalID* nodeIDs,
1316 nodeIDs, results) );
1322 numNodes = problemStructure_->getNodeDatabase().
getNodeIDs().size();
1331 std::map<GlobalID,int>& nodes =
1332 problemStructure_->getNodeDatabase().
getNodeIDs();
1333 numNodes = nodes.size();
1335 if (lenNodeIDs < len) len = lenNodeIDs;
1345 const GlobalID* nodeIDs,
1346 const double* nodeData)
1349 nodeIDs, nodeData) );
1353 void FEI_Implementation::buildLinearSystem()
1371 debugOut(
"# buildLinearSystem");
1375 int anyDataLoaded = newMatrixDataLoaded_;
1377 if (numProcs_ > 1) {
1378 if (MPI_Allreduce(&newMatrixDataLoaded_, &anyDataLoaded, 1, MPI_INT,
1379 MPI_SUM, comm_) != MPI_SUCCESS) voidERReturn
1383 if (anyDataLoaded) {
1385 for(
int i=0; i<numInternalFEIs_; i++) {
1386 filter_[i]->exchangeRemoteEquations();
1389 newMatrixDataLoaded_ = 0;
1392 if (solveType_ == 2){
1397 if (!aggregateSystemFormed_) {
1398 if (!matScalarsSet_ || !rhsScalarsSet_) {
1399 FEI_COUT <<
"FEI_Implementation: WARNING: solveType_==2, aggregating system, but setMatScalars and/or setRHSScalars not yet called." << FEI_ENDL;
1406 filter_[index_soln_filter_]->loadComplete();
1408 debugOut(
"# leaving buildLinearSystem");
1412 int FEI_Implementation::aggregateSystem()
1414 debugOut(
"# aggregateSystem");
1415 if (!haveLinSysCore_) ERReturn(-1);
1417 if (soln_fei_matrix_ == NULL) {
1418 soln_fei_matrix_ =
new Data();
1420 CHK_ERR( lscArray_[index_soln_filter_]->
1421 copyOutMatrix(1.0, *soln_fei_matrix_) );
1424 if (soln_fei_vector_ == NULL) {
1425 soln_fei_vector_ =
new Data();
1427 CHK_ERR( lscArray_[index_soln_filter_]->
1428 setRHSID(rhsIDs_[index_soln_filter_][0]) );
1430 CHK_ERR( lscArray_[index_soln_filter_]->
1431 copyOutRHSVector(1.0, *soln_fei_vector_) );
1437 for(
int i=0; i<numInternalFEIs_; i++){
1439 if (i == index_soln_filter_) {
1442 CHK_ERR( lscArray_[index_soln_filter_]->
1443 copyInMatrix(matScalars_[i], tmp) );
1446 CHK_ERR( lscArray_[i]->getMatrixPtr(tmp) );
1448 CHK_ERR( lscArray_[index_soln_filter_]->
1449 sumInMatrix(matScalars_[i], tmp) );
1452 for(
int j=0; j<numRHSIDs_[i]; j++){
1453 if ((i == index_soln_filter_) && (j == 0)) {
1458 CHK_ERR( lscArray_[i]->setRHSID(rhsIDs_[i][j]) );
1459 CHK_ERR( lscArray_[i]->getRHSVectorPtr(tmpv) );
1462 if (i == index_soln_filter_) {
1463 CHK_ERR( lscArray_[index_soln_filter_]->
1464 copyInRHSVector(rhsScalars_[i][j], tmpv) );
1467 CHK_ERR( lscArray_[index_soln_filter_]->
1468 sumInRHSVector(rhsScalars_[i][j], tmpv) );
1473 aggregateSystemFormed_ =
true;
1475 debugOut(
"# leaving aggregateSystem");
1481 int FEI_Implementation::allocateInternalFEIs(){
1490 if (internalFEIsAllocated_)
return(0);
1492 matScalars_.resize(numInternalFEIs_);
1494 if (rhsScalars_.size() != 0) deleteRHSScalars();
1496 rhsScalars_.resize(numInternalFEIs_);
1498 for(
int i=0; i<numInternalFEIs_; i++){
1499 matScalars_[i] = 1.0;
1501 rhsScalars_[i] =
new double[numRHSIDs_[i]];
1503 for(
int j=0; j<numRHSIDs_[i]; j++){
1504 rhsScalars_[i][j] = 1.0;
1508 IDsAllocated_ =
true;
1510 if (numInternalFEIs_ > 0) {
1511 index_soln_filter_ = 0;
1512 index_current_filter_ = 0;
1513 filter_ =
new Filter*[numInternalFEIs_];
1515 if (haveLinSysCore_) {
1516 if (numRHSIDs_[0] == 0) {
1524 for(
int i=1; i<numInternalFEIs_; i++) {
1528 if (numRHSIDs_[i] == 0) {
1530 lsc->setNumRHSVectors(1, &dummyID);
1533 lsc->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
1536 lscArray_.push_back(lsc);
1540 for(
int i=0; i<numInternalFEIs_; i++){
1542 if (haveLinSysCore_) {
1544 lscArray_[i].
get(), masterRank_);
1546 else if (haveFEData_) {
1547 filter_[i] =
new FEDataFilter(
this, comm_, problemStructure_,
1548 wrapper_.
get(), masterRank_);
1551 fei::console_out() <<
"FEI_Implementation: ERROR, don't have LinearSystemCore"
1552 <<
" or FiniteElementData implementation..." << FEI_ENDL;
1556 filter_[i]->setLogStream(dbgOStreamPtr_);
1558 FEI_OSTRINGSTREAM osstr;
1559 osstr<<
"internalFei "<< i;
1560 std::string osstr_str = osstr.str();
1561 const char* param = osstr_str.c_str();
1562 filter_[i]->parameters(1, ¶m);
1565 (*dbgOStreamPtr_)<<
"#-- fei["<<i<<
"]->setNumRHSVectors "
1566 <<numRHSIDs_[i]<<FEI_ENDL;
1569 if (numRHSIDs_[i] == 0) {
1571 filter_[i]->setNumRHSVectors(1, &dummyID);
1574 filter_[i]->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
1578 internalFEIsAllocated_ =
true;
1581 needParametersAbort(
"FEI_Implementation::allocateInternalFEIs");
1588 void FEI_Implementation::debugOut(
const char* msg) {
1589 if (debugOutput_) { (*dbgOStreamPtr_)<<msg<<FEI_ENDL; }
1593 void FEI_Implementation::debugOut(
const char* msg,
int whichFEI) {
1595 (*dbgOStreamPtr_)<<msg<<
", -> fei["<<whichFEI<<
"]"<<FEI_ENDL;
1600 void FEI_Implementation::messageAbort(
const char* msg){
1602 fei::console_out() <<
"FEI_Implementation: ERROR " << msg <<
" Aborting." << FEI_ENDL;
1603 MPI_Abort(comm_, -1);
1607 void FEI_Implementation::notAllocatedAbort(
const char* name){
1610 << FEI_ENDL <<
"ERROR, internal data structures not allocated."
1611 << FEI_ENDL <<
"'setIDLists' and/or 'setSolveType' must be called"
1612 << FEI_ENDL <<
"first to identify solveType and number of matrices"
1613 << FEI_ENDL <<
"to be assembled." << FEI_ENDL;
1614 MPI_Abort(comm_, -1);
1618 void FEI_Implementation::needParametersAbort(
const char* name){
1621 << FEI_ENDL <<
"FEI_Implementation: ERROR, numMatrices has not been specified."
1622 << FEI_ENDL <<
"FEI_Implementation: 'parameters' must be called up front with"
1623 << FEI_ENDL <<
"FEI_Implementation: the string 'numMatrices n' to specify that"
1624 << FEI_ENDL <<
"FEI_Implementation: n matrices will be assembled." << FEI_ENDL;
1625 MPI_Abort(comm_, -1);
1629 void FEI_Implementation::badParametersAbort(
const char* name){
1632 << FEI_ENDL <<
"FEI_Implementation: ERROR, inconsistent 'solveType' and"
1633 << FEI_ENDL <<
"FEI_Implementation: 'numMatrices' parameters specified."
1634 << FEI_ENDL <<
"FEI_Implementation: Aborting."
1636 MPI_Abort(comm_, -1);
int getBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
int resetSystem(double s=0.0)
int putBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, const double *estimates)
int getParameters(int &numParams, char **¶mStrings)
int getFieldSize(int fieldID)
int setRHSScalars(int numScalars, const int *IDs, const double *scalars)
virtual ~FEI_Implementation()
std::map< GlobalID, int > & getNodeIDs()
int getNodalFieldSolution(int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
int putBlockNodeSolution(GlobalID elemBlockID, int numNodes, const GlobalID *nodeIDs, const int *offsets, const double *estimates)
int sumInElem(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, const double *elemLoad, int elemFormat)
int loadElemBCs(int numElems, const GlobalID *elemIDs, int fieldID, const double *const *alpha, const double *const *beta, const double *const *gamma)
FEI_Implementation(fei::SharedPtr< LibraryWrapper > libWrapper, MPI_Comm comm, int masterRank=0)
int loadNodeBCs(int numNodes, const GlobalID *nodeIDs, int fieldID, const int *offsetsIntoField, const double *prescribedValues)
int getNodalSolution(int numNodes, const GlobalID *nodeIDs, int *offsets, double *results)
int initSlaveVariable(GlobalID slaveNodeID, int slaveFieldID, int offsetIntoSlaveField, int numMasterNodes, const GlobalID *masterNodeIDs, const int *masterFieldIDs, const double *weights, double rhsValue)
int iterations(int &itersTaken) const
void setTypeName(const char *name)
int setIDLists(int numMatrices, const int *matrixIDs, int numRHSs, const int *rhsIDs)
int getBlockFieldNodeSolution(GlobalID elemBlockID, int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)
int getNumElemBlocks(int &numElemBlocks) const
int sumInElemMatrix(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *const *elemStiffness, int elemFormat)
int initCRPen(int numCRNodes, const GlobalID *CRNodes, const int *CRFields, int &CRID)
int version(const char *&versionString)
int initCRMult(int numCRNodes, const GlobalID *CRNodes, const int *CRFields, int &CRID)
int getEqnNumbers(GlobalID ID, int idType, int fieldID, int &numEqns, int *eqnNumbers)
int resetInitialGuess(double s=0.0)
int sumInElemRHS(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn, const double *elemLoad)
void copyKeysToArray(const MAP_TYPE &map_obj, unsigned lenList, int *list)
int sumIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *rhsEntries)
int initElemBlock(GlobalID elemBlockID, int numElements, int numNodesPerElement, const int *numFieldsPerNode, const int *const *nodalFieldIDs, int numElemDofFieldsPerElement, const int *elemDOFFieldIDs, int interleaveStrategy)
int setSolveType(int solveType)
int getNumBlockActEqns(GlobalID blockID, int &numEqns) const
int loadComplete(bool applyBCs=true, bool globalAssemble=true)
int getBlockNodeIDList(GlobalID elemBlockID, int numNodes, GlobalID *nodeIDs)
int getBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int &numElemDOFPerElement, double *results)
int initElem(GlobalID elemBlockID, GlobalID elemID, const GlobalID *elemConn)
int getCRMultIDList(int numMultCRs, int *multIDs)
virtual LinearSystemCore * clone()=0
int setCurrentRHS(int rhsID)
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int getNumEqnsPerElement(GlobalID blockID, int &numEqns) const
int getBlockElemIDList(GlobalID elemBlockID, int numElems, GlobalID *elemIDs)
int cumulative_cpu_times(double &initTime, double &loadTime, double &solveTime, double &solnReturnTime)
int putBlockElemSolution(GlobalID elemBlockID, int numElems, const GlobalID *elemIDs, int dofPerElem, const double *estimates)
int putNodalFieldData(int fieldID, int numNodes, const GlobalID *nodeIDs, const double *nodeData)
int parameters(int numParams, const char *const *paramStrings)
int initFields(int numFields, const int *fieldSizes, const int *fieldIDs, const int *fieldTypes=NULL)
virtual int setNumRHSVectors(int numRHSs, const int *rhsIDs)=0
void * getDataPtr() const
std::ostream & console_out()
void setDataPtr(void *ptr)
virtual int destroyVectorData(Data &data)=0
int getNumCRMultipliers(int &numMultCRs)
int getNumNodesPerElement(GlobalID blockID, int &nodesPerElem) const
const char * getTypeName() const
int putIntoRHS(int IDType, int fieldID, int numIDs, const GlobalID *IDs, const double *rhsEntries)
int localProc(MPI_Comm comm)
int loadCRPen(int CRID, int numCRNodes, const GlobalID *CRNodes, const int *CRFields, const double *CRWeights, double CRValue, double penValue)
int getCRMultipliers(int numCRs, const int *CRIDs, double *multipliers)
virtual int parameters(int numParams, const char *const *params)=0
int getNumBlockElements(GlobalID blockID, int &numElems) const
int resetMatrix(double s=0.0)
int initSharedNodes(int numSharedNodes, const GlobalID *sharedNodeIDs, const int *numProcsPerNode, const int *const *sharingProcIDs)
int getNumSolnParams(GlobalID nodeID, int &numSolnParams) const
int getNumLocalNodes(int &numNodes)
int loadCRMult(int CRID, int numCRNodes, const GlobalID *CRNodes, const int *CRFields, const double *CRWeights, double CRValue)
int getLocalNodeIDList(int &numNodes, GlobalID *nodeIDs, int lenNodeIDs)
int putCRMultipliers(int numMultCRs, const int *CRIDs, const double *multEstimates)
int setMatScalars(int numScalars, const int *IDs, const double *scalars)
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
int parameters(int numParams, const char *const *paramStrings)
int resetRHSVector(double s=0.0)
int getFieldSize(int fieldID, int &numScalars)
int residualNorm(int whichNorm, int numFields, int *fieldIDs, double *norms)
int getNumBlockElemDOF(GlobalID blockID, int &DOFPerElem) const
int searchList(const T &item, const T *list, int len)
virtual int destroyMatrixData(Data &data)=0
int getNumBlockActNodes(GlobalID blockID, int &numNodes) const
int numProcs(MPI_Comm comm)
int getIndexOfBlock(GlobalID blockID) const
int setCurrentMatrix(int matID)