FEI Version of the Day
FEI_Implementation.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_CommUtils.hpp>
00010 #include <fei_iostream.hpp>
00011 #include <fei_fstream.hpp>
00012 #include <fei_sstream.hpp>
00013 
00014 #include <fei_utils.hpp>
00015 
00016 #include <fei_Data.hpp>
00017 #include <fei_LinearSystemCore.hpp>
00018 #include <fei_FiniteElementData.hpp>
00019 
00020 #include <fei_LibraryWrapper.hpp>
00021 #include <fei_Filter.hpp>
00022 
00023 #include <fei_LinSysCoreFilter.hpp>
00024 #include <fei_FEDataFilter.hpp>
00025 
00026 #include <SNL_FEI_Structure.hpp>
00027 #include <fei_BlockDescriptor.hpp>
00028 #include <fei_NodeDatabase.hpp>
00029 #include <fei_ConnectivityTable.hpp>
00030 #include <FEI_Implementation.hpp>
00031 #include <snl_fei_Utils.hpp>
00032 
00033 #undef fei_file
00034 #define fei_file "FEI_Implementation.cpp"
00035 
00036 #include <fei_ErrMacros.hpp>
00037 
00038 //------------------------------------------------------------------------------
00039 FEI_Implementation::FEI_Implementation(fei::SharedPtr<LibraryWrapper> libWrapper,
00040                                        MPI_Comm comm, int masterRank)
00041  : wrapper_(libWrapper),
00042    linSysCore_(NULL),
00043    lscArray_(),
00044    haveLinSysCore_(false),
00045    haveFEData_(false),
00046    problemStructure_(NULL),
00047    filter_(NULL),
00048    numInternalFEIs_(0),
00049    internalFEIsAllocated_(false),
00050    matrixIDs_(),
00051    numRHSIDs_(),
00052    rhsIDs_(),
00053    IDsAllocated_(false),
00054    matScalars_(),
00055    matScalarsSet_(false),
00056    rhsScalars_(),
00057    rhsScalarsSet_(false),
00058    index_soln_filter_(0),
00059    index_current_filter_(0),
00060    index_current_rhs_row_(0),
00061    solveType_(-1),
00062    setSolveTypeCalled_(false),
00063    initPhaseIsComplete_(false),
00064    aggregateSystemFormed_(false),
00065    newMatrixDataLoaded_(0),
00066    soln_fei_matrix_(NULL),
00067    soln_fei_vector_(NULL),
00068    comm_(comm),
00069    masterRank_(0),
00070    localRank_(0),
00071    numProcs_(1),
00072    outputLevel_(0),
00073    solveCounter_(1),
00074    debugOutput_(0),
00075    dbgOStreamPtr_(NULL),
00076    dbgFileOpened_(false),
00077    dbgFStreamPtr_(NULL),
00078    initTime_(0.0),
00079    loadTime_(0.0),
00080    solveTime_(0.0),
00081    solnReturnTime_(0.0),
00082    numParams_(0),
00083    paramStrings_(NULL)
00084 {
00085   localRank_ = fei::localProc(comm);
00086   numProcs_ = fei::numProcs(comm);
00087 
00088    problemStructure_ = new SNL_FEI_Structure(comm_);
00089 
00090    if (problemStructure_ == NULL) {
00091      messageAbort("problemStructure allocation failed");
00092    }
00093 
00094    //If we have a FiniteElementData instance as the underlying
00095    //data receptacle and solver, then we'll set the shared-node-ownership rule
00096    //to make sure shared nodes are owned by a proc which contains them in
00097    //local elements.
00098    haveFEData_ = wrapper_->haveFiniteElementData();
00099    if (haveFEData_) {
00100      haveFEData_ = true;
00101      NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
00102      nodeCommMgr.setSharedOwnershipRule(NodeCommMgr::PROC_WITH_LOCAL_ELEM);
00103    }
00104 
00105    haveLinSysCore_ = wrapper_->haveLinearSystemCore();
00106    if (haveLinSysCore_) {
00107      linSysCore_ = wrapper_->getLinearSystemCore();
00108      lscArray_.push_back(linSysCore_);
00109    }
00110 
00111    numInternalFEIs_ = 1;
00112    matrixIDs_.resize(1);
00113    matrixIDs_[0] = 0;
00114    numRHSIDs_.resize(1);
00115    numRHSIDs_[0] = 1;
00116    rhsIDs_.resize(1);
00117    rhsIDs_[0] = new int[1];
00118    rhsIDs_[0][0] = 0;
00119    rhsScalars_.resize(numInternalFEIs_);
00120    for(int ii=0; ii<numInternalFEIs_; ii++) rhsScalars_[ii] = NULL;
00121 
00122    //  and the time spent in the constructor is...
00123 
00124    return;
00125 }
00126 
00127 //------------------------------------------------------------------------------
00128 FEI_Implementation::~FEI_Implementation()
00129 {
00130   //
00131   //  Destructor function. Free allocated memory, etc.
00132   //
00133 
00134   if (debugOutput_) {
00135     (*dbgOStreamPtr_) << "FEI: destructor" << FEI_ENDL;
00136   }
00137 
00138   if (soln_fei_matrix_) {
00139       linSysCore_->destroyMatrixData(*soln_fei_matrix_);
00140       delete soln_fei_matrix_;
00141       soln_fei_matrix_ = NULL;
00142    }
00143 
00144    if (soln_fei_vector_) {
00145       linSysCore_->destroyVectorData(*soln_fei_vector_);
00146       delete soln_fei_vector_;
00147       soln_fei_vector_ = NULL;
00148    }
00149 
00150    deleteIDs();
00151 
00152    if (internalFEIsAllocated_) {
00153       for(int j=0; j<numInternalFEIs_; j++){
00154          delete filter_[j];
00155       }
00156       delete [] filter_;
00157    }
00158 
00159    deleteRHSScalars();
00160 
00161    internalFEIsAllocated_ = false;
00162    numInternalFEIs_ = 0;
00163 
00164    delete problemStructure_;
00165 
00166    for(int k=0; k<numParams_; k++) delete [] paramStrings_[k];
00167    delete [] paramStrings_;
00168 
00169    if (dbgFileOpened_ == true) {
00170      dbgFStreamPtr_->close(); delete dbgFStreamPtr_;
00171    }
00172    else delete dbgOStreamPtr_;
00173 }
00174 
00175 
00176 //------------------------------------------------------------------------------
00177 void FEI_Implementation::deleteIDs()
00178 {
00179   matrixIDs_.resize(0);
00180   for(size_t i=0; i<rhsIDs_.size(); i++) {
00181     delete [] rhsIDs_[i];
00182   }
00183   rhsIDs_.resize(0);
00184   numRHSIDs_.resize(0);
00185 }
00186 
00187 //------------------------------------------------------------------------------
00188 void FEI_Implementation::deleteRHSScalars()
00189 {
00190   for(size_t i=0; i<rhsScalars_.size(); i++) {
00191     delete [] rhsScalars_[i];
00192   }
00193   rhsScalars_.resize(0);
00194 }
00195 
00196 //------------------------------------------------------------------------------
00197 int FEI_Implementation::setCurrentMatrix(int matID)
00198 {
00199    if (debugOutput_) {
00200      (*dbgOStreamPtr_) << "FEI: setCurrentMatrix" << FEI_ENDL << "#matrix-id"
00201                        << FEI_ENDL<<matID<<FEI_ENDL;
00202    }
00203 
00204    index_current_filter_ = -1;
00205 
00206    for(int i=0; i<numInternalFEIs_; i++){
00207       if (matrixIDs_[i] == matID) index_current_filter_ = i;
00208    }
00209 
00210    if (debugOutput_) {
00211      (*dbgOStreamPtr_) << "#--- ID: " << matID
00212                        << ", ind: "<<index_current_filter_<<FEI_ENDL;
00213    }
00214 
00215    //if matID wasn't found, return non-zero (error)
00216    if (index_current_filter_ == -1) {
00217       fei::console_out() << "FEI_Implementation::setCurrentMatrix: ERROR, invalid matrix ID "
00218            << "supplied" << FEI_ENDL;
00219       return(-1);
00220    }
00221 
00222    debugOut("#FEI_Implementation leaving setCurrentMatrix");
00223 
00224    return(0);
00225 }
00226 
00227 //------------------------------------------------------------------------------
00228 int FEI_Implementation::getParameters(int& numParams, char**& paramStrings)
00229 {
00230   numParams = numParams_;
00231   paramStrings = paramStrings_;
00232   return(0);
00233 }
00234 
00235 //------------------------------------------------------------------------------
00236 int FEI_Implementation::setCurrentRHS(int rhsID)
00237 {
00238   if (debugOutput_) {
00239     (*dbgOStreamPtr_) << "FEI: setCurrentRHS" << FEI_ENDL << "#rhs-id"
00240                       << FEI_ENDL<<rhsID<<FEI_ENDL;
00241   }
00242 
00243   bool found = false;
00244 
00245   for(int j=0; j<numInternalFEIs_; j++){
00246     int index = fei::searchList(rhsID, rhsIDs_[j], numRHSIDs_[j]);
00247     if (index >= 0) {
00248       index_current_rhs_row_ = j;
00249       CHK_ERR( filter_[index_current_rhs_row_]->setCurrentRHS(rhsID) )
00250         found = true;
00251       break;
00252     }
00253   }
00254 
00255   if (!found) {
00256     fei::console_out() << "FEI_Implementation::setCurrentRHS: ERROR, invalid RHS ID" 
00257          << FEI_ENDL;
00258     ERReturn(-1);
00259   }
00260 
00261   debugOut("#FEI_Implementation leaving setCurrentRHS");
00262 
00263   return(0);
00264 }
00265 
00266 //------------------------------------------------------------------------------
00267 int FEI_Implementation::setSolveType(int solveType)
00268 {
00269   if (debugOutput_) {
00270     (*dbgOStreamPtr_)<<"FEI: setSolveType"<<FEI_ENDL;
00271     (*dbgOStreamPtr_)<<solveType<<FEI_ENDL;
00272   }
00273 
00274   solveType_ = solveType;
00275 
00276   if (solveType_ == FEI_SINGLE_SYSTEM) {
00277     if (matrixIDs_.size() > 1) {
00278       messageAbort("setSolveType: solve-type is FEI_SINGLE_SYSTEM, but setIDLists() has been called with numMatrices > 1.");
00279     }
00280   }
00281   else if (solveType_ == FEI_EIGEN_SOLVE) {
00282   }
00283   else if (solveType_ == FEI_AGGREGATE_SUM) {
00284     //solving a linear-combination of separately
00285     //assembled matrices and rhs vectors
00286   }
00287   else if (solveType_ == FEI_AGGREGATE_PRODUCT) {
00288     //solving a product of separately assembled
00289     //matrices -- i.e., (C^T*M*C)x = rhs
00290   }
00291   else if (solveType_ == 4) {
00292     //4 means we'll be doing a multi-level solution
00293   }
00294 
00295   return(0);
00296 }
00297 
00298 //------------------------------------------------------------------------------
00299 int FEI_Implementation::setIDLists(int numMatrices, const int* matrixIDs,
00300                                    int numRHSs, const int* rhsIDs)
00301 {
00302    if (debugOutput_) {
00303      FEI_OSTREAM& os = *dbgOStreamPtr_;
00304      os << "FEI: setIDLists" << FEI_ENDL
00305         << "#num-matrices" << FEI_ENDL << numMatrices << FEI_ENDL
00306         << "#matrixIDs" << FEI_ENDL;
00307      int i;
00308      for(i=0; i<numMatrices; ++i) os << matrixIDs[i] << " ";
00309      os << FEI_ENDL << "#num-rhs's" << FEI_ENDL;
00310      for(i=0; i<numRHSs; ++i) os << rhsIDs[i] << " ";
00311      os << FEI_ENDL;
00312    }
00313 
00314    deleteIDs();
00315 
00316    // We will try to assign the rhs's evenly over the matrices. i.e., give
00317    // roughly equal numbers of rhs's to each matrix.
00318 
00319    //first, let's make sure we have at least 1 matrixID to which we can assign
00320    //rhs's...
00321    int myNumMatrices = numMatrices;
00322    if (myNumMatrices == 0) myNumMatrices = 1;
00323 
00324    matrixIDs_.resize(myNumMatrices);
00325 
00326    if (rhsScalars_.size() != 0) deleteRHSScalars();
00327 
00328    numInternalFEIs_ = myNumMatrices;
00329 
00330    if (numMatrices == 0) {
00331       matrixIDs_[0] = 0;
00332    }
00333    else {
00334       for(int i=0; i<numMatrices; i++) matrixIDs_[i] = matrixIDs[i];
00335    }
00336 
00337    int quotient = numRHSs/myNumMatrices;
00338    int rem = numRHSs%numMatrices;
00339 
00340    //the allocateInternalFEIs function which will be called later from within
00341    //initComplete(), takes a list of matrixIDs, and a list
00342    //of numRHSsPerMatrix, and then a table of rhsIDs, where the table has a row
00343    //for each matrixID. Each of those rows is a list of the rhsIDs assigned to
00344    //the corresponding matrix. Is that clear???
00345 
00346    numRHSIDs_.resize(myNumMatrices);
00347    rhsIDs_.resize(myNumMatrices);
00348 
00349    int offset = 0;
00350    for(int i=0; i<myNumMatrices; i++) {
00351       numRHSIDs_[i] = quotient;
00352       if (i < rem) numRHSIDs_[i]++;
00353 
00354       rhsIDs_[i] = numRHSIDs_[i] > 0 ? new int[numRHSIDs_[i]] : NULL ;
00355 
00356       for(int j=0; j<numRHSIDs_[i]; j++) {
00357         rhsIDs_[i][j] = rhsIDs[offset+j];
00358       }
00359 
00360       offset += numRHSIDs_[i];
00361    }
00362 
00363    return(0);
00364 }
00365 
00366 //------------------------------------------------------------------------------
00367 int FEI_Implementation::initFields(int numFields,
00368                                    const int *fieldSizes,
00369                                    const int *fieldIDs,
00370                                    const int *fieldTypes)
00371 {
00372     CHK_ERR( problemStructure_->initFields(numFields, fieldSizes, fieldIDs, fieldTypes) );
00373 
00374     return(0);
00375 }
00376  
00377 //------------------------------------------------------------------------------
00378 int FEI_Implementation::initElemBlock(GlobalID elemBlockID,
00379                                       int numElements,
00380                                       int numNodesPerElement,
00381                                       const int* numFieldsPerNode,
00382                                       const int* const* nodalFieldIDs,
00383                                       int numElemDofFieldsPerElement,
00384                                       const int* elemDOFFieldIDs,
00385                                       int interleaveStrategy)
00386 {
00387    CHK_ERR( problemStructure_->initElemBlock(elemBlockID,
00388                                              numElements,
00389                                              numNodesPerElement,
00390                                              numFieldsPerNode,
00391                                              nodalFieldIDs,
00392                                              numElemDofFieldsPerElement,
00393                                              elemDOFFieldIDs,
00394                                              interleaveStrategy) );
00395 
00396    return(0);
00397 }
00398 
00399 //------------------------------------------------------------------------------
00400 int FEI_Implementation::initElem(GlobalID elemBlockID,
00401                                  GlobalID elemID,
00402                                  const GlobalID* elemConn)
00403 {
00404    CHK_ERR( problemStructure_->initElem(elemBlockID, elemID, elemConn) );
00405 
00406    return(0);
00407 }
00408 
00409 //------------------------------------------------------------------------------
00410 int FEI_Implementation::initSlaveVariable(GlobalID slaveNodeID, 
00411                                           int slaveFieldID,
00412                                           int offsetIntoSlaveField,
00413                                           int numMasterNodes,
00414                                           const GlobalID* masterNodeIDs,
00415                                           const int* masterFieldIDs,
00416                                           const double* weights,
00417                                           double rhsValue)
00418 {
00419    CHK_ERR( problemStructure_->initSlaveVariable(slaveNodeID, slaveFieldID,
00420                                             offsetIntoSlaveField,
00421                                                    numMasterNodes, masterNodeIDs,
00422                                             masterFieldIDs, weights, rhsValue));
00423 
00424    return(0);
00425 }
00426 
00427 //------------------------------------------------------------------------------
00428 int FEI_Implementation::deleteMultCRs()
00429 {
00430   debugOut("FEI: deleteMultCRs");
00431 
00432   CHK_ERR( problemStructure_->deleteMultCRs() );
00433 
00434   int err = -1;
00435   if (internalFEIsAllocated_) {
00436     err = filter_[index_current_filter_]->deleteMultCRs();
00437   }
00438 
00439   return(err);
00440 }
00441 
00442 //------------------------------------------------------------------------------
00443 int FEI_Implementation::initSharedNodes(int numSharedNodes,
00444                                         const GlobalID *sharedNodeIDs,
00445                                         const int* numProcsPerNode,
00446                                         const int *const *sharingProcIDs)
00447 {
00448   //
00449   //  In this function we simply accumulate the incoming data into
00450   // internal arrays in the problemStructure_ object.
00451   //
00452   CHK_ERR( problemStructure_->initSharedNodes(numSharedNodes,
00453                                               sharedNodeIDs,
00454                                               numProcsPerNode,
00455                                               sharingProcIDs));
00456 
00457   return(0);
00458 }
00459 
00460 //------------------------------------------------------------------------------
00461 int FEI_Implementation::initCRMult(int numCRNodes,
00462                                    const GlobalID* CRNodes,
00463                                    const int *CRFields,
00464                                    int& CRID)
00465 {
00466 //
00467 // Store Lagrange Multiplier constraint data into internal structures, and
00468 // return an identifier 'CRID' by which this constraint may be referred to
00469 // later.
00470 //
00471 
00472    CHK_ERR( problemStructure_->initCRMult(numCRNodes,
00473                                           CRNodes,
00474                                           CRFields,
00475                                           CRID));
00476 
00477    return(0);
00478 }
00479 
00480 //------------------------------------------------------------------------------
00481 int FEI_Implementation::initCRPen(int numCRNodes,
00482                                   const GlobalID* CRNodes,
00483                                   const int *CRFields,
00484                                   int& CRID)
00485 {
00486 //
00487 // Store penalty constraint data and return an identifier 'CRID' by which the
00488 // constraint may be referred to later.
00489 //
00490 
00491    CHK_ERR( problemStructure_->initCRPen(numCRNodes,
00492                                          CRNodes,
00493                                          CRFields,
00494                                          CRID));
00495 
00496    return(0);
00497 }
00498 
00499 //------------------------------------------------------------------------------
00500 int FEI_Implementation::initComplete()
00501 {
00502     bool generateGraph = !haveFEData_;
00503 
00504     CHK_ERR( problemStructure_->initComplete(generateGraph) );
00505 
00506     //now allocate one or more internal instances of Filter, depending on
00507     //whether the user has indicated that they're doing an aggregate solve
00508     //etc., via the functions setSolveType() and setIDLists().
00509 
00510     CHK_ERR( allocateInternalFEIs() );
00511 
00512     for(int i=0; i<numInternalFEIs_; ++i) {
00513       CHK_ERR( filter_[i]->initialize() );
00514     }
00515 
00516     problemStructure_->destroyMatIndices();
00517 
00518     initPhaseIsComplete_ = true;
00519     return(0);
00520 }
00521 
00522 //------------------------------------------------------------------------------
00523 int FEI_Implementation::resetSystem(double s)
00524 {
00525 //
00526 //  This puts the value s throughout both the matrix and the vector.
00527 //
00528    if (!internalFEIsAllocated_) return(0);
00529 
00530    CHK_ERR( filter_[index_current_filter_]->resetSystem(s))
00531  
00532    return(0);
00533 }
00534 
00535 //------------------------------------------------------------------------------
00536 int FEI_Implementation::resetMatrix(double s)
00537 {
00538    if (!internalFEIsAllocated_) return(0);
00539 
00540    CHK_ERR( filter_[index_current_filter_]->resetMatrix(s))
00541 
00542    return(0);
00543 }
00544 
00545 //------------------------------------------------------------------------------
00546 int FEI_Implementation::resetRHSVector(double s)
00547 {
00548    if (!internalFEIsAllocated_) return(0);
00549 
00550    CHK_ERR( filter_[index_current_filter_]->resetRHSVector(s))
00551 
00552    return(0);
00553 }
00554 
00555 //------------------------------------------------------------------------------
00556 int FEI_Implementation::resetInitialGuess(double s)
00557 {
00558    if (!internalFEIsAllocated_) return(0);
00559 
00560    CHK_ERR( filter_[index_current_filter_]->resetInitialGuess(s))
00561 
00562    return(0);
00563 }
00564 
00565 //------------------------------------------------------------------------------
00566 int FEI_Implementation::loadNodeBCs(int numNodes,
00567                                     const GlobalID *nodeIDs,
00568                                     int fieldID,
00569                                     const int* offsetsIntoField,
00570                                     const double* prescribedValues)
00571 {
00572    if (!internalFEIsAllocated_)
00573       notAllocatedAbort("FEI_Implementation::loadNodeBCs");
00574 
00575    int index = index_current_filter_;
00576    if (solveType_ == 2) index = index_soln_filter_;
00577 
00578    CHK_ERR( filter_[index]->loadNodeBCs(numNodes,
00579                                nodeIDs, fieldID,
00580                                offsetsIntoField, prescribedValues));
00581 
00582    return(0);
00583 }
00584 
00585 //------------------------------------------------------------------------------
00586 int FEI_Implementation::loadElemBCs(int numElems,
00587                                     const GlobalID *elemIDs,
00588                                     int fieldID,
00589                                     const double *const *alpha,
00590                                     const double *const *beta,
00591                                     const double *const *gamma)
00592 {
00593    if (!internalFEIsAllocated_)
00594       notAllocatedAbort("FEI_Implementation::loadElemBCs");
00595 
00596    int index = index_current_filter_;
00597    if (solveType_ == 2) index = index_soln_filter_;
00598 
00599    CHK_ERR( filter_[index]->loadElemBCs(numElems,
00600                                elemIDs, fieldID,
00601                                alpha, beta, gamma))
00602 
00603    return(0);
00604 }
00605 
00606 //------------------------------------------------------------------------------
00607 int FEI_Implementation::sumInElem(GlobalID elemBlockID,
00608                         GlobalID elemID,
00609                         const GlobalID* elemConn,
00610                         const double* const* elemStiffness,
00611                         const double* elemLoad,
00612                         int elemFormat)
00613 {
00614   if (!internalFEIsAllocated_) {
00615     notAllocatedAbort("FEI_Implementation::sumInElem");
00616   }
00617 
00618   CHK_ERR( filter_[index_current_filter_]->sumInElem(elemBlockID, elemID,
00619                                                      elemConn, elemStiffness,
00620                                                      elemLoad, elemFormat));
00621 
00622   newMatrixDataLoaded_ = 1;
00623 
00624   return(0);
00625 }
00626 
00627 //------------------------------------------------------------------------------
00628 int FEI_Implementation::sumInElemMatrix(GlobalID elemBlockID,
00629                         GlobalID elemID,
00630                         const GlobalID* elemConn,
00631                         const double* const* elemStiffness,
00632                         int elemFormat)
00633 {
00634   if (!internalFEIsAllocated_)
00635     notAllocatedAbort("FEI_Implementation::sumInElemMatrix");
00636 
00637   CHK_ERR( filter_[index_current_filter_]->sumInElemMatrix(elemBlockID,
00638                                           elemID, elemConn,
00639                                           elemStiffness, elemFormat))
00640 
00641   newMatrixDataLoaded_ = 1;
00642 
00643 
00644   return(0);
00645 }
00646 
00647 //------------------------------------------------------------------------------
00648 int FEI_Implementation::sumInElemRHS(GlobalID elemBlockID,
00649                         GlobalID elemID,
00650                         const GlobalID* elemConn,
00651                         const double* elemLoad)
00652 {
00653    if (!internalFEIsAllocated_)
00654       notAllocatedAbort("FEI_Implementation::sumInElemRHS");
00655 
00656    CHK_ERR( filter_[index_current_rhs_row_]->sumInElemRHS(elemBlockID,
00657                                           elemID, elemConn, elemLoad))
00658 
00659    newMatrixDataLoaded_ = 1;
00660 
00661    return(0);
00662 }
00663 
00664 //------------------------------------------------------------------------------
00665 int FEI_Implementation::loadCRMult(int CRID,
00666                                    int numCRNodes,
00667                                    const GlobalID* CRNodes,
00668                                    const int* CRFields,
00669                                    const double* CRWeights,
00670                                    double CRValue)
00671 {
00672    if (!internalFEIsAllocated_)
00673       notAllocatedAbort("FEI_Implementation::loadCRMult");
00674 
00675    newMatrixDataLoaded_ = 1;
00676 
00677    CHK_ERR( filter_[index_current_filter_]->loadCRMult(CRID,
00678                                            numCRNodes, CRNodes,
00679                                            CRFields, CRWeights, CRValue));
00680 
00681    return(0);
00682 }
00683 
00684 //------------------------------------------------------------------------------
00685 int FEI_Implementation::loadCRPen(int CRID,
00686                                   int numCRNodes,
00687                                   const GlobalID* CRNodes,
00688                                   const int* CRFields,
00689                                   const double* CRWeights,
00690                                   double CRValue,
00691                                   double penValue)
00692 {
00693    if (!internalFEIsAllocated_)
00694       notAllocatedAbort("FEI_Implementation::loadCRPen");
00695 
00696    CHK_ERR( filter_[index_current_filter_]->loadCRPen(CRID,
00697                                           numCRNodes, CRNodes,
00698                                           CRFields, CRWeights,
00699                                           CRValue, penValue))
00700 
00701    newMatrixDataLoaded_ = 1;
00702 
00703    return(0);
00704 }
00705 
00706 //------------------------------------------------------------------------------
00707 int FEI_Implementation::sumIntoRHS(int IDType,
00708                                    int fieldID,
00709                                    int numIDs,
00710                                    const GlobalID* IDs,
00711                                    const double* rhsEntries)
00712 {
00713   if (!internalFEIsAllocated_)
00714     notAllocatedAbort("FEI_Implementation::sumIntoRHS");
00715 
00716   CHK_ERR( filter_[index_current_rhs_row_]->sumIntoRHS(IDType, fieldID,
00717                                                        numIDs, IDs, rhsEntries) );
00718   newMatrixDataLoaded_ = 1;
00719 
00720   return(0);
00721 }
00722 
00723 //------------------------------------------------------------------------------
00724 int FEI_Implementation::sumIntoMatrixDiagonal(int IDType,
00725                              int fieldID,
00726                              int numIDs,
00727                              const GlobalID* IDs,
00728                              const double* coefficients)
00729 {
00730   if (!internalFEIsAllocated_)
00731     notAllocatedAbort("FEI_Implementation::sumIntoMatrixDiagonal");
00732 
00733   CHK_ERR( filter_[index_current_filter_]->sumIntoMatrixDiagonal(IDType, fieldID,
00734                                                        numIDs, IDs, coefficients) );
00735   newMatrixDataLoaded_ = 1;
00736 
00737   return(0);
00738 }
00739 
00740 //------------------------------------------------------------------------------
00741 int FEI_Implementation::putIntoRHS(int IDType,
00742                                    int fieldID,
00743                                    int numIDs,
00744                                    const GlobalID* IDs,
00745                                    const double* rhsEntries)
00746 {
00747   if (!internalFEIsAllocated_)
00748     notAllocatedAbort("FEI_Implementation::putIntoRHS");
00749 
00750   CHK_ERR( filter_[index_current_rhs_row_]->putIntoRHS(IDType, fieldID,
00751                                                        numIDs, IDs, rhsEntries) );
00752   newMatrixDataLoaded_ = 1;
00753 
00754   return(0);
00755 }
00756 
00757 //------------------------------------------------------------------------------
00758 int FEI_Implementation::setMatScalars(int numScalars,
00759                                       const int* IDs, const double* scalars)
00760 {
00761    for(int i=0; i<numScalars; i++){
00762       std::vector<int>::iterator iter =
00763           std::find(matrixIDs_.begin(), matrixIDs_.end(), IDs[i]);
00764       if (iter != matrixIDs_.end()) {
00765          int index = iter - matrixIDs_.begin();
00766          matScalars_[index] = scalars[i];
00767       }
00768       else {
00769          fei::console_out() << "FEI_Implementation::setMatScalars: ERROR, invalid ID supplied"
00770               << FEI_ENDL;
00771          return(1);
00772       }
00773    }
00774 
00775    aggregateSystemFormed_ = false;
00776    matScalarsSet_ = true;
00777 
00778    return(0);
00779 }
00780 
00781 //------------------------------------------------------------------------------
00782 int FEI_Implementation::setRHSScalars(int numScalars,
00783                                       const int* IDs, const double* scalars)
00784 {
00785   for(int i=0; i<numScalars; i++){
00786      bool found = false;
00787 
00788      for(int j=0; j<numInternalFEIs_; j++){
00789          int index = fei::searchList(IDs[i], rhsIDs_[j], numRHSIDs_[j]);
00790          if (index>=0) {
00791              rhsScalars_[j][index] = scalars[i];
00792              found = true;
00793              break;
00794          }
00795       }
00796 
00797       if (!found) {
00798          fei::console_out() << "FEI_Implementation::setRHSScalars: ERROR, invalid RHS ID supplied"
00799              << FEI_ENDL;
00800          return(1);
00801       }
00802    }
00803 
00804    aggregateSystemFormed_ = false;
00805    rhsScalarsSet_ = true;
00806 
00807    return(0);
00808 }
00809 
00810 //------------------------------------------------------------------------------
00811 int FEI_Implementation::parameters(int numParams, const char *const* paramStrings)
00812 {
00813   // this function takes parameters and passes them to the internal
00814   // fei objects.
00815 
00816   if (numParams == 0 || paramStrings == NULL) {
00817     debugOut("#--- no parameters");
00818     return(0);
00819   }
00820 
00821   // merge these parameters with any others we may have, for later use.
00822   snl_fei::mergeStringLists(paramStrings_, numParams_,
00823                                    paramStrings, numParams);
00824 
00825   snl_fei::getIntParamValue("numMatrices", numParams,paramStrings, numInternalFEIs_);
00826 
00827   snl_fei::getIntParamValue("outputLevel", numParams,paramStrings, outputLevel_);
00828 
00829   const char* param = snl_fei::getParamValue("debugOutput",numParams,paramStrings);
00830   if (param != NULL) {
00831     setDebugOutput(param,"FEI_log");
00832   }
00833 
00834   if (debugOutput_) {
00835     (*dbgOStreamPtr_)<<"FEI: parameters"<<FEI_ENDL;
00836     (*dbgOStreamPtr_)<<"#FEI_Implementation, num-params "<<FEI_ENDL
00837                      <<numParams<<FEI_ENDL;
00838     (*dbgOStreamPtr_)<<"# "<<numParams<<" parameter lines follow:"<<FEI_ENDL;
00839     for(int i=0; i<numParams; i++){
00840       (*dbgOStreamPtr_)<<paramStrings[i]<<FEI_ENDL;
00841     }
00842   }
00843 
00844   if (haveLinSysCore_) {
00845     linSysCore_->parameters(numParams, (char**)paramStrings);
00846   }
00847   if (haveFEData_) {
00848     wrapper_->getFiniteElementData()->parameters(numParams, (char**)paramStrings);
00849   }
00850 
00851   problemStructure_->parameters(numParams, paramStrings);
00852 
00853   if (internalFEIsAllocated_){
00854     for(int i=0; i<numInternalFEIs_; i++){
00855       CHK_ERR( filter_[i]->parameters(numParams, paramStrings) );
00856     }
00857   }
00858 
00859   debugOut("#FEI_Implementation leaving parameters method");
00860 
00861   return(0);
00862 }
00863 
00864 //------------------------------------------------------------------------------
00865 void FEI_Implementation::setDebugOutput(const char* path, const char* name)
00866 {
00867   //
00868   //This function turns on debug output, and opens a file to put it in.
00869   //
00870   if (dbgFileOpened_) {
00871     dbgFStreamPtr_->close();
00872   }
00873 
00874   dbgFileOpened_ = false;
00875   delete dbgOStreamPtr_;
00876 
00877   FEI_OSTRINGSTREAM osstr;
00878   if (path != NULL) {
00879     osstr << path << "/";
00880   }
00881   osstr << name << "." << numProcs_ << "." << localRank_;
00882 
00883   debugOutput_ = 1;
00884   dbgFStreamPtr_ = new FEI_OFSTREAM(osstr.str().c_str(), IOS_APP);
00885   if (!dbgFStreamPtr_ || dbgFStreamPtr_->bad()){
00886     fei::console_out() << "couldn't open debug output file: " << osstr.str() << FEI_ENDL;
00887     debugOutput_ = 0;
00888   }
00889 
00890   if (debugOutput_) {
00891     const char* version_str = NULL;
00892     version(version_str);
00893 
00894     (*dbgFStreamPtr_) << version_str << FEI_ENDL;
00895 
00896     problemStructure_->setDbgOut(*dbgFStreamPtr_, path, "_0");
00897     dbgOStreamPtr_ = dbgFStreamPtr_;
00898     dbgOStreamPtr_->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00899     dbgFileOpened_ = true;
00900 
00901     if (internalFEIsAllocated_) {
00902       for(int i=0; i<numInternalFEIs_; ++i) {
00903         filter_[i]->setLogStream(dbgOStreamPtr_);
00904       }
00905     }
00906   }
00907 }
00908 
00909 //------------------------------------------------------------------------------
00910 int FEI_Implementation::loadComplete(bool applyBCs,
00911                                      bool globalAssemble)
00912 {
00913   (void)applyBCs;
00914   (void)globalAssemble;
00915 
00916   buildLinearSystem();
00917 
00918   return(0);
00919 }
00920 
00921 //------------------------------------------------------------------------------
00922 int FEI_Implementation::residualNorm(int whichNorm, int numFields,
00923                                      int* fieldIDs, double* norms)
00924 {
00925    buildLinearSystem();
00926 
00927    double residTime = 0.0;
00928 
00929    int err = filter_[index_soln_filter_]->residualNorm(whichNorm, numFields,
00930                                                  fieldIDs, norms, residTime);
00931 
00932    solveTime_ += residTime;
00933 
00934    return(err);
00935 }
00936 
00937 //------------------------------------------------------------------------------
00938 int FEI_Implementation::solve(int& status)
00939 {
00940    buildLinearSystem();
00941 
00942    double sTime = 0.0;
00943 
00944    int err = filter_[index_soln_filter_]->solve(status, sTime);
00945 
00946    solveTime_ += sTime;
00947 
00948    return(err);
00949 }
00950 
00951 //------------------------------------------------------------------------------
00952 int FEI_Implementation::iterations(int& itersTaken) const {
00953   itersTaken = filter_[index_soln_filter_]->iterations();
00954   return(0);
00955 }
00956 
00957 //------------------------------------------------------------------------------
00958 int FEI_Implementation::version(const char*& versionString)
00959 {
00960   versionString = fei::utils::version();
00961   return(0);
00962 }
00963 
00964 //------------------------------------------------------------------------------
00965 int FEI_Implementation::cumulative_cpu_times(double& initTime,
00966                                              double& loadTime,
00967                                              double& solveTime,
00968                                              double& solnReturnTime)
00969 {
00970    initTime = initTime_;
00971    loadTime = loadTime_;
00972    solveTime = solveTime_;
00973    solnReturnTime = solnReturnTime_;
00974 
00975    return(0);
00976 }
00977 
00978 //------------------------------------------------------------------------------
00979 int FEI_Implementation::getBlockNodeSolution(GlobalID elemBlockID,
00980                                              int numNodes,
00981                                              const GlobalID *nodeIDs,
00982                                              int *offsets,
00983                                              double *results)
00984 {
00985    CHK_ERR(filter_[index_soln_filter_]->getBlockNodeSolution(elemBlockID,
00986                                                numNodes,
00987                                                nodeIDs,
00988                                                offsets,
00989                                                results))
00990    return(0);
00991 }
00992 
00993 //------------------------------------------------------------------------------
00994 int FEI_Implementation::getNodalSolution(int numNodes,
00995                                          const GlobalID *nodeIDs,
00996                                          int *offsets,
00997                                          double *results)
00998 {
00999    CHK_ERR(filter_[index_soln_filter_]->getNodalSolution(numNodes,
01000                                                          nodeIDs,
01001                                                          offsets,
01002                                                          results))
01003 
01004    return(0);
01005 }
01006 
01007 //------------------------------------------------------------------------------
01008 int FEI_Implementation::getBlockFieldNodeSolution(GlobalID elemBlockID,
01009                                         int fieldID,
01010                                         int numNodes, 
01011                                         const GlobalID *nodeIDs, 
01012                                         double *results)
01013 {
01014    CHK_ERR( filter_[index_soln_filter_]->getBlockFieldNodeSolution(elemBlockID,
01015                                                      fieldID, numNodes,
01016                                                      nodeIDs, results))
01017 
01018    return(0);
01019 }
01020 
01021 //------------------------------------------------------------------------------
01022 int FEI_Implementation::putBlockNodeSolution(GlobalID elemBlockID,
01023                                              int numNodes,
01024                                              const GlobalID *nodeIDs,
01025                                              const int *offsets,
01026                                              const double *estimates)
01027 {
01028    CHK_ERR( filter_[index_soln_filter_]->putBlockNodeSolution(elemBlockID,
01029                                                 numNodes, nodeIDs,
01030                                                 offsets, estimates))
01031    return(0);
01032 }
01033 
01034 //------------------------------------------------------------------------------
01035 int FEI_Implementation::putBlockFieldNodeSolution(GlobalID elemBlockID,
01036                                         int fieldID,
01037                                         int numNodes,
01038                                         const GlobalID *nodeIDs,
01039                                         const double *estimates)
01040 {
01041    int err = filter_[index_soln_filter_]->putBlockFieldNodeSolution(elemBlockID,
01042                                                      fieldID, numNodes,
01043                                                      nodeIDs, estimates);
01044    return(err);
01045 }
01046 
01047 //------------------------------------------------------------------------------
01048 int FEI_Implementation::getBlockElemSolution(GlobalID elemBlockID,  
01049                                    int numElems, 
01050                                    const GlobalID *elemIDs,
01051                                    int& numElemDOFPerElement,
01052                                    double *results)
01053 {
01054    CHK_ERR( filter_[index_soln_filter_]->getBlockElemSolution(elemBlockID,
01055                                                        numElems, elemIDs,
01056                                                        numElemDOFPerElement,
01057                                                        results))
01058     return(0);
01059 } 
01060       
01061 //------------------------------------------------------------------------------
01062 int FEI_Implementation::putBlockElemSolution(GlobalID elemBlockID,
01063                                    int numElems,
01064                                    const GlobalID *elemIDs,
01065                                    int dofPerElem,
01066                                    const double *estimates)
01067 {
01068    CHK_ERR( filter_[index_soln_filter_]->putBlockElemSolution(elemBlockID,
01069                                                        numElems, elemIDs,
01070                                                        dofPerElem, estimates))
01071    return(0);
01072 }
01073 
01074 //------------------------------------------------------------------------------
01075 int FEI_Implementation::getNumCRMultipliers(int& numMultCRs)
01076 {
01077   numMultCRs = problemStructure_->getNumMultConstRecords();
01078   return(0);
01079 }
01080 
01081 //------------------------------------------------------------------------------
01082 int FEI_Implementation::getCRMultIDList(int numMultCRs,
01083                                         int* multIDs)
01084 {
01085    if (numMultCRs > problemStructure_->getNumMultConstRecords()) return(-1);
01086 
01087    std::map<GlobalID,snl_fei::Constraint<GlobalID>*>::const_iterator
01088      cr_iter = problemStructure_->getMultConstRecords().begin(),
01089      cr_end = problemStructure_->getMultConstRecords().end();
01090    int i = 0;
01091    while(cr_iter != cr_end) {
01092       multIDs[i++] = (*cr_iter).first;
01093       ++cr_iter;
01094    }
01095 
01096    return(0);
01097 }
01098 
01099 //------------------------------------------------------------------------------
01100 int FEI_Implementation::getCRMultipliers(int numMultCRs,
01101                                          const int* CRIDs,
01102                                          double* multipliers)
01103 {
01104    CHK_ERR( filter_[index_soln_filter_]->getCRMultipliers(numMultCRs,
01105                                                      CRIDs, multipliers))
01106    return(0);
01107 }
01108 
01109 //------------------------------------------------------------------------------
01110 int FEI_Implementation::putCRMultipliers(int numMultCRs,
01111                                          const int* CRIDs,
01112                                          const double *multEstimates)
01113 {
01114     return(
01115     filter_[index_soln_filter_]->putCRMultipliers(numMultCRs,
01116                                             CRIDs,
01117                                             multEstimates)
01118     );
01119 }
01120 
01121 //-----------------------------------------------------------------------------
01122 //  some utility functions to aid in using the "put" functions for passing
01123 //  an initial guess to the solver
01124 //-----------------------------------------------------------------------------
01125 
01126 //------------------------------------------------------------------------------
01127 int FEI_Implementation::getBlockElemIDList(GlobalID elemBlockID,
01128                                            int numElems,
01129                                            GlobalID* elemIDs)
01130 {
01131   //
01132   //  return the list of element IDs for a given block... the length parameter
01133   //  lenElemIDList should be used to check memory allocation for the calling
01134   //  method, as the calling method should have gotten a copy of this param
01135   //  from a call to getNumBlockElements before allocating memory for elemIDList
01136   //
01137   ConnectivityTable& connTable = problemStructure_->
01138     getBlockConnectivity(elemBlockID);
01139 
01140   std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
01141   int len = elemIDList.size();
01142   if (len > numElems) len = numElems;
01143 
01144   fei::copyKeysToArray(elemIDList, len, elemIDs);
01145 
01146   return(FEI_SUCCESS);
01147 }
01148 
01149 //------------------------------------------------------------------------------
01150 int FEI_Implementation::getBlockNodeIDList(GlobalID elemBlockID,
01151                                            int numNodes,
01152                                            GlobalID *nodeIDs)
01153 {
01154   //
01155   //  similar comments as for getBlockElemIDList(), except for returning the
01156   //  active node list
01157   //
01158   int numActiveNodes = problemStructure_->getNumActiveNodes();
01159   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01160 
01161   int offset = 0;
01162   for(int i=0; i<numActiveNodes; i++) {
01163     const NodeDescriptor* node = NULL;
01164     nodeDB.getNodeAtIndex(i, node);
01165     if (node==NULL) continue;
01166     if (node->containedInBlock(elemBlockID))
01167       nodeIDs[offset++] = node->getGlobalNodeID();
01168     if (offset == numNodes) break;
01169   }
01170 
01171   return(FEI_SUCCESS);
01172 }
01173 
01174 //------------------------------------------------------------------------------
01175 int FEI_Implementation::getNumNodesPerElement(GlobalID blockID,
01176                                               int& nodesPerElem) const
01177 {
01178   //
01179   //  return the number of nodes associated with elements of a given block ID
01180   //
01181   BlockDescriptor* block = NULL;
01182   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01183 
01184   nodesPerElem = block->getNumNodesPerElement();
01185   return(FEI_SUCCESS);
01186 }
01187  
01188 //------------------------------------------------------------------------------
01189 int FEI_Implementation::getNumEqnsPerElement(GlobalID blockID, int& numEqns)
01190 const
01191 {
01192   //
01193   //  return the number of eqns associated with elements of a given block ID
01194   //
01195   BlockDescriptor* block = NULL;
01196   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01197 
01198   numEqns = block->getNumEqnsPerElement();
01199   return(FEI_SUCCESS);
01200 }
01201 
01202 //------------------------------------------------------------------------------
01203 int FEI_Implementation::getNumSolnParams(GlobalID nodeID, int& numSolnParams)
01204 const {
01205   //
01206   //  return the number of solution parameters at a given node
01207   //
01208   const NodeDescriptor* node = NULL;
01209   int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
01210 
01211   if (err != 0) {
01212     ERReturn(-1);
01213   }
01214 
01215   numSolnParams = node->getNumNodalDOF();
01216   return(0);
01217 }
01218  
01219 //------------------------------------------------------------------------------
01220 int FEI_Implementation::getNumElemBlocks(int& numElemBlocks) const
01221 {
01222   //
01223   //  the number of element blocks
01224   //
01225   numElemBlocks = problemStructure_->getNumElemBlocks();
01226   return( 0 );
01227 }
01228 
01229 //------------------------------------------------------------------------------
01230 int FEI_Implementation::getNumBlockActNodes(GlobalID blockID, int& numNodes)
01231 const {
01232   //
01233   //  return the number of active nodes associated with a given element block ID
01234   //
01235   BlockDescriptor* block = NULL;
01236   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01237 
01238   numNodes = block->getNumActiveNodes();
01239   return(FEI_SUCCESS);
01240 }
01241 
01242 //------------------------------------------------------------------------------
01243 int FEI_Implementation::getNumBlockActEqns(GlobalID blockID, int& numEqns)
01244 const {
01245 //
01246 // return the number of active equations associated with a given element
01247 // block ID
01248 //
01249   BlockDescriptor* block = NULL;
01250   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01251 
01252   numEqns = block->getTotalNumEqns();
01253   return(FEI_SUCCESS);
01254 }
01255 
01256 //------------------------------------------------------------------------------
01257 int FEI_Implementation::getNumBlockElements(GlobalID blockID, int& numElems) const {
01258 //
01259 //  return the number of elements associated with a given elem blockID
01260 //
01261   BlockDescriptor* block = NULL;
01262   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01263 
01264   numElems = block->getNumElements();
01265   return(FEI_SUCCESS);
01266 }
01267 
01268 //------------------------------------------------------------------------------
01269 int FEI_Implementation::getNumBlockElemDOF(GlobalID blockID,
01270                                            int& DOFPerElem) const
01271 {
01272   //
01273   //  return the number of elem equations associated with a given elem blockID
01274   //
01275   BlockDescriptor* block = NULL;
01276   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01277 
01278   DOFPerElem = block->getNumElemDOFPerElement();
01279 
01280   return(FEI_SUCCESS);
01281 }
01282 
01283 //------------------------------------------------------------------------------
01284 int FEI_Implementation::getFieldSize(int fieldID,
01285                                      int& numScalars)
01286 {
01287   //
01288   //  return the number of scalars associated with a given fieldID
01289   //
01290 
01291   numScalars = problemStructure_->getFieldSize(fieldID);
01292   return(0);  
01293 }
01294 
01295 //------------------------------------------------------------------------------
01296 int FEI_Implementation::getEqnNumbers(GlobalID ID, int idType,
01297                                       int fieldID, int& numEqns,
01298                                       int* eqnNumbers)
01299 {
01300   //
01301   // Translate from an ID/fieldID pair to a list of equation-numbers
01302   //
01303 
01304   return( problemStructure_->getEqnNumbers(ID, idType, fieldID,
01305                                            numEqns, eqnNumbers) );
01306 }
01307 
01308 //------------------------------------------------------------------------------
01309 int FEI_Implementation::getNodalFieldSolution(int fieldID,
01310                                               int numNodes,
01311                                               const GlobalID* nodeIDs,
01312                                               double* results)
01313 {
01314   return( filter_[index_soln_filter_]->getNodalFieldSolution(fieldID, numNodes,
01315                                                        nodeIDs, results) );
01316 }
01317 
01318 //------------------------------------------------------------------------------
01319 int FEI_Implementation::getNumLocalNodes(int& numNodes)
01320 {
01321   numNodes = problemStructure_->getNodeDatabase().getNodeIDs().size();
01322   return( 0 );
01323 }
01324 
01325 //------------------------------------------------------------------------------
01326 int FEI_Implementation::getLocalNodeIDList(int& numNodes,
01327                                            GlobalID* nodeIDs,
01328                                            int lenNodeIDs)
01329 {
01330   std::map<GlobalID,int>& nodes =
01331     problemStructure_->getNodeDatabase().getNodeIDs();
01332   numNodes = nodes.size();
01333   int len = numNodes;
01334   if (lenNodeIDs < len) len = lenNodeIDs;
01335 
01336   fei::copyKeysToArray(nodes, len, nodeIDs);
01337 
01338   return( 0 );
01339 }
01340 
01341 //------------------------------------------------------------------------------
01342 int FEI_Implementation::putNodalFieldData(int fieldID,
01343                                           int numNodes,
01344                                           const GlobalID* nodeIDs,
01345                                           const double* nodeData)
01346 {
01347   return( filter_[index_soln_filter_]->putNodalFieldData(fieldID, numNodes,
01348                                                    nodeIDs, nodeData) );
01349 }
01350 
01351 //------------------------------------------------------------------------------
01352 void FEI_Implementation::buildLinearSystem()
01353 {
01354   //
01355   //Private function.
01356   //
01357   //At the point when this function is called, all matrix assembly has
01358   //already taken place, with the data having been directed into the
01359   //appropriate Filter instance in the filter_ list. Now it's
01360   //time to finalize the linear system (get A,x,b ready to give to a
01361   //solver), performing this list of last-minute tasks:
01362   //
01363   //1. Have each Filter instance exchangeRemoteEquations.
01364   //2. Aggregate the system (form a linear combination of LHS's and
01365   //   RHS's) if solveType_==2.
01366   //3. call loadComplete to have the 'master' Filter instance
01367   //   (filter_[index_soln_filter_]) enforce any boundary conditions
01368   //   that have been loaded.
01369   //
01370   debugOut("#   buildLinearSystem");
01371 
01372   //
01373   //figure out if new matrix-data was loaded on any processor.
01374   int anyDataLoaded = newMatrixDataLoaded_;
01375 #ifndef FEI_SER
01376   if (numProcs_ > 1) {
01377     if (MPI_Allreduce(&newMatrixDataLoaded_, &anyDataLoaded, 1, MPI_INT,
01378                       MPI_SUM, comm_) != MPI_SUCCESS) voidERReturn
01379   }
01380 #endif
01381 
01382   if (anyDataLoaded) {
01383 #ifndef FEI_SER
01384     for(int i=0; i<numInternalFEIs_; i++) {
01385       filter_[i]->exchangeRemoteEquations();
01386     }
01387 #endif
01388     newMatrixDataLoaded_ = 0;
01389   }
01390 
01391   if (solveType_ == 2){
01392     //solveType_ == 2 means this is a linear-combination solve --
01393     //i.e., we're solving an aggregate system which is the sum of
01394     //several individual matrices and rhs's.
01395 
01396     if (!aggregateSystemFormed_) {
01397       if (!matScalarsSet_ || !rhsScalarsSet_) {
01398         FEI_COUT << "FEI_Implementation: WARNING: solveType_==2, aggregating system, but setMatScalars and/or setRHSScalars not yet called." << FEI_ENDL;
01399       }
01400 
01401       aggregateSystem();
01402     }
01403   }
01404 
01405   filter_[index_soln_filter_]->loadComplete();
01406 
01407   debugOut("#   leaving buildLinearSystem");
01408 }
01409 
01410 //------------------------------------------------------------------------------
01411 int FEI_Implementation::aggregateSystem()
01412 {
01413   debugOut("#   aggregateSystem");
01414    if (!haveLinSysCore_) ERReturn(-1);
01415 
01416    if (soln_fei_matrix_ == NULL) {
01417       soln_fei_matrix_ = new Data();
01418 
01419       CHK_ERR( lscArray_[index_soln_filter_]->
01420                copyOutMatrix(1.0, *soln_fei_matrix_) );
01421    }
01422 
01423    if (soln_fei_vector_ == NULL) {
01424       soln_fei_vector_ = new Data();
01425 
01426       CHK_ERR( lscArray_[index_soln_filter_]->
01427                      setRHSID(rhsIDs_[index_soln_filter_][0]) );
01428 
01429       CHK_ERR( lscArray_[index_soln_filter_]->
01430                     copyOutRHSVector(1.0, *soln_fei_vector_) );
01431    }
01432 
01433    Data tmp;
01434    Data tmpv;
01435 
01436    for(int i=0; i<numInternalFEIs_; i++){
01437 
01438       if (i == index_soln_filter_) {
01439          tmp.setTypeName(soln_fei_matrix_->getTypeName());
01440          tmp.setDataPtr(soln_fei_matrix_->getDataPtr());
01441          CHK_ERR( lscArray_[index_soln_filter_]->
01442                         copyInMatrix(matScalars_[i], tmp) );
01443       }
01444       else {
01445          CHK_ERR( lscArray_[i]->getMatrixPtr(tmp) );
01446 
01447          CHK_ERR( lscArray_[index_soln_filter_]->
01448                         sumInMatrix(matScalars_[i], tmp) );
01449       }
01450 
01451       for(int j=0; j<numRHSIDs_[i]; j++){
01452          if ((i == index_soln_filter_) && (j == 0)) {
01453             tmpv.setTypeName(soln_fei_vector_->getTypeName());
01454             tmpv.setDataPtr(soln_fei_vector_->getDataPtr());
01455          }
01456          else {
01457             CHK_ERR( lscArray_[i]->setRHSID(rhsIDs_[i][j]) );
01458             CHK_ERR( lscArray_[i]->getRHSVectorPtr(tmpv) );
01459          }
01460 
01461          if (i == index_soln_filter_) {
01462             CHK_ERR( lscArray_[index_soln_filter_]->
01463                          copyInRHSVector(rhsScalars_[i][j], tmpv) );
01464          }
01465          else {
01466             CHK_ERR( lscArray_[index_soln_filter_]->
01467                          sumInRHSVector(rhsScalars_[i][j], tmpv) );
01468          }
01469       }
01470    }
01471 
01472    aggregateSystemFormed_ = true;
01473 
01474    debugOut("#   leaving aggregateSystem");
01475 
01476    return(0);
01477 }
01478 
01479 //==============================================================================
01480 int FEI_Implementation::allocateInternalFEIs(){
01481 //
01482 //This is a private FEI_Implementation function, to be called from within
01483 //setSolveType or the other overloading of allocateInternalFEIs.
01484 //Assumes that numInternalFEIs_ has already been set.
01485 //It is also safe to assume that problemStructure_->initComplete() has already
01486 //been called.
01487 //
01488 
01489    if (internalFEIsAllocated_) return(0);
01490 
01491    matScalars_.resize(numInternalFEIs_);
01492 
01493    if (rhsScalars_.size() != 0) deleteRHSScalars();
01494 
01495    rhsScalars_.resize(numInternalFEIs_);
01496 
01497    for(int i=0; i<numInternalFEIs_; i++){
01498       matScalars_[i] = 1.0;
01499 
01500       rhsScalars_[i] = new double[numRHSIDs_[i]];
01501 
01502       for(int j=0; j<numRHSIDs_[i]; j++){
01503          rhsScalars_[i][j] = 1.0;
01504       }
01505    }
01506 
01507    IDsAllocated_ = true;
01508 
01509    if (numInternalFEIs_ > 0) {
01510       index_soln_filter_ = 0;
01511       index_current_filter_ = 0;
01512       filter_ = new Filter*[numInternalFEIs_];
01513 
01514       if (haveLinSysCore_) {
01515         if (numRHSIDs_[0] == 0) {
01516           int dummyID = -1;
01517           linSysCore_->setNumRHSVectors(1, &dummyID);
01518         }
01519         else {
01520           linSysCore_->setNumRHSVectors(numRHSIDs_[0], rhsIDs_[0]);
01521         }
01522 
01523         for(int i=1; i<numInternalFEIs_; i++) {
01524           fei::SharedPtr<LinearSystemCore> lsc(linSysCore_->clone());
01525           lsc->parameters(numParams_, paramStrings_);
01526 
01527           if (numRHSIDs_[i] == 0) {
01528             int dummyID = -1;
01529             lsc->setNumRHSVectors(1, &dummyID);
01530           }
01531           else {
01532             lsc->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
01533           }
01534 
01535           lscArray_.push_back(lsc);
01536         }
01537       }
01538 
01539       for(int i=0; i<numInternalFEIs_; i++){
01540 
01541         if (haveLinSysCore_) {
01542           filter_[i] = new LinSysCoreFilter(this, comm_, problemStructure_,
01543                                             lscArray_[i].get(), masterRank_);
01544         }
01545         else if (haveFEData_) {
01546           filter_[i] = new FEDataFilter(this, comm_, problemStructure_,
01547                                         wrapper_.get(), masterRank_);
01548         }
01549         else {
01550           fei::console_out() << "FEI_Implementation: ERROR, don't have LinearSystemCore"
01551                << " or FiniteElementData implementation..." << FEI_ENDL;
01552           ERReturn(-1);
01553         }
01554 
01555         filter_[i]->setLogStream(dbgOStreamPtr_);
01556 
01557         FEI_OSTRINGSTREAM osstr;
01558         osstr<<"internalFei "<< i;
01559         std::string osstr_str = osstr.str();
01560         const char* param = osstr_str.c_str();
01561         filter_[i]->parameters(1, &param);
01562 
01563         if (debugOutput_) {
01564           (*dbgOStreamPtr_)<<"#-- fei["<<i<<"]->setNumRHSVectors "
01565                            <<numRHSIDs_[i]<<FEI_ENDL;
01566         }
01567 
01568         if (numRHSIDs_[i] == 0) {
01569           int dummyID = -1;
01570           filter_[i]->setNumRHSVectors(1, &dummyID);
01571         }
01572         else {
01573           filter_[i]->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
01574         }
01575       }
01576 
01577       internalFEIsAllocated_ = true;
01578    }
01579    else {
01580      needParametersAbort("FEI_Implementation::allocateInternalFEIs");
01581    }
01582 
01583    return(0);
01584 }
01585 
01586 //==============================================================================
01587 void FEI_Implementation::debugOut(const char* msg) {
01588   if (debugOutput_) { (*dbgOStreamPtr_)<<msg<<FEI_ENDL; }
01589 }
01590 
01591 //==============================================================================
01592 void FEI_Implementation::debugOut(const char* msg, int whichFEI) {
01593    if (debugOutput_) {
01594      (*dbgOStreamPtr_)<<msg<<", -> fei["<<whichFEI<<"]"<<FEI_ENDL;
01595    }
01596 }
01597 
01598 //==============================================================================
01599 void FEI_Implementation::messageAbort(const char* msg){
01600 
01601     fei::console_out() << "FEI_Implementation: ERROR " << msg << " Aborting." << FEI_ENDL;
01602     MPI_Abort(comm_, -1);
01603 }
01604 
01605 //==============================================================================
01606 void FEI_Implementation::notAllocatedAbort(const char* name){
01607 
01608     fei::console_out() << name
01609          << FEI_ENDL << "ERROR, internal data structures not allocated."
01610          << FEI_ENDL << "'setIDLists' and/or 'setSolveType' must be called"
01611          << FEI_ENDL << "first to identify solveType and number of matrices"
01612          << FEI_ENDL << "to be assembled." << FEI_ENDL;
01613     MPI_Abort(comm_, -1);
01614 }
01615 
01616 //==============================================================================
01617 void FEI_Implementation::needParametersAbort(const char* name){
01618 
01619    fei::console_out() << name
01620      << FEI_ENDL << "FEI_Implementation: ERROR, numMatrices has not been specified."
01621      << FEI_ENDL << "FEI_Implementation: 'parameters' must be called up front with"
01622      << FEI_ENDL << "FEI_Implementation: the string 'numMatrices n' to specify that"
01623      << FEI_ENDL << "FEI_Implementation: n matrices will be assembled." << FEI_ENDL;
01624     MPI_Abort(comm_, -1);
01625 }
01626 
01627 //==============================================================================
01628 void FEI_Implementation::badParametersAbort(const char* name){
01629 
01630    fei::console_out() << name
01631         << FEI_ENDL << "FEI_Implementation: ERROR, inconsistent 'solveType' and"
01632         << FEI_ENDL << "FEI_Implementation: 'numMatrices' parameters specified."
01633         << FEI_ENDL << "FEI_Implementation: Aborting."
01634         << FEI_ENDL;
01635     MPI_Abort(comm_, -1);
01636 }
01637 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends