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_CERR << "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_CERR << "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 {
00371     CHK_ERR( problemStructure_->initFields(numFields, fieldSizes, fieldIDs) );
00372 
00373     return(0);
00374 }
00375  
00376 //------------------------------------------------------------------------------
00377 int FEI_Implementation::initElemBlock(GlobalID elemBlockID,
00378                                       int numElements,
00379                                       int numNodesPerElement,
00380                                       const int* numFieldsPerNode,
00381                                       const int* const* nodalFieldIDs,
00382                                       int numElemDofFieldsPerElement,
00383                                       const int* elemDOFFieldIDs,
00384                                       int interleaveStrategy)
00385 {
00386    CHK_ERR( problemStructure_->initElemBlock(elemBlockID,
00387                                              numElements,
00388                                              numNodesPerElement,
00389                                              numFieldsPerNode,
00390                                              nodalFieldIDs,
00391                                              numElemDofFieldsPerElement,
00392                                              elemDOFFieldIDs,
00393                                              interleaveStrategy) );
00394 
00395    return(0);
00396 }
00397 
00398 //------------------------------------------------------------------------------
00399 int FEI_Implementation::initElem(GlobalID elemBlockID,
00400                                  GlobalID elemID,
00401                                  const GlobalID* elemConn)
00402 {
00403    CHK_ERR( problemStructure_->initElem(elemBlockID, elemID, elemConn) );
00404 
00405    return(0);
00406 }
00407 
00408 //------------------------------------------------------------------------------
00409 int FEI_Implementation::initSlaveVariable(GlobalID slaveNodeID, 
00410                                           int slaveFieldID,
00411                                           int offsetIntoSlaveField,
00412                                           int numMasterNodes,
00413                                           const GlobalID* masterNodeIDs,
00414                                           const int* masterFieldIDs,
00415                                           const double* weights,
00416                                           double rhsValue)
00417 {
00418    CHK_ERR( problemStructure_->initSlaveVariable(slaveNodeID, slaveFieldID,
00419                                             offsetIntoSlaveField,
00420                                                    numMasterNodes, masterNodeIDs,
00421                                             masterFieldIDs, weights, rhsValue));
00422 
00423    return(0);
00424 }
00425 
00426 //------------------------------------------------------------------------------
00427 int FEI_Implementation::deleteMultCRs()
00428 {
00429   debugOut("FEI: deleteMultCRs");
00430 
00431   CHK_ERR( problemStructure_->deleteMultCRs() );
00432 
00433   int err = -1;
00434   if (internalFEIsAllocated_) {
00435     err = filter_[index_current_filter_]->deleteMultCRs();
00436   }
00437 
00438   return(err);
00439 }
00440 
00441 //------------------------------------------------------------------------------
00442 int FEI_Implementation::initSharedNodes(int numSharedNodes,
00443                                         const GlobalID *sharedNodeIDs,
00444                                         const int* numProcsPerNode,
00445                                         const int *const *sharingProcIDs)
00446 {
00447   //
00448   //  In this function we simply accumulate the incoming data into
00449   // internal arrays in the problemStructure_ object.
00450   //
00451   CHK_ERR( problemStructure_->initSharedNodes(numSharedNodes,
00452                                               sharedNodeIDs,
00453                                               numProcsPerNode,
00454                                               sharingProcIDs));
00455 
00456   return(0);
00457 }
00458 
00459 //------------------------------------------------------------------------------
00460 int FEI_Implementation::initCRMult(int numCRNodes,
00461                                    const GlobalID* CRNodes,
00462                                    const int *CRFields,
00463                                    int& CRID)
00464 {
00465 //
00466 // Store Lagrange Multiplier constraint data into internal structures, and
00467 // return an identifier 'CRID' by which this constraint may be referred to
00468 // later.
00469 //
00470 
00471    CHK_ERR( problemStructure_->initCRMult(numCRNodes,
00472                                           CRNodes,
00473                                           CRFields,
00474                                           CRID));
00475 
00476    return(0);
00477 }
00478 
00479 //------------------------------------------------------------------------------
00480 int FEI_Implementation::initCRPen(int numCRNodes,
00481                                   const GlobalID* CRNodes,
00482                                   const int *CRFields,
00483                                   int& CRID)
00484 {
00485 //
00486 // Store penalty constraint data and return an identifier 'CRID' by which the
00487 // constraint may be referred to later.
00488 //
00489 
00490    CHK_ERR( problemStructure_->initCRPen(numCRNodes,
00491                                          CRNodes,
00492                                          CRFields,
00493                                          CRID));
00494 
00495    return(0);
00496 }
00497 
00498 //------------------------------------------------------------------------------
00499 int FEI_Implementation::initComplete()
00500 {
00501     bool generateGraph = !haveFEData_;
00502 
00503     CHK_ERR( problemStructure_->initComplete(generateGraph) );
00504 
00505     //now allocate one or more internal instances of Filter, depending on
00506     //whether the user has indicated that they're doing an aggregate solve
00507     //etc., via the functions setSolveType() and setIDLists().
00508 
00509     CHK_ERR( allocateInternalFEIs() );
00510 
00511     for(int i=0; i<numInternalFEIs_; ++i) {
00512       CHK_ERR( filter_[i]->initialize() );
00513     }
00514 
00515     problemStructure_->destroyMatIndices();
00516 
00517     initPhaseIsComplete_ = true;
00518     return(0);
00519 }
00520 
00521 //------------------------------------------------------------------------------
00522 int FEI_Implementation::resetSystem(double s)
00523 {
00524 //
00525 //  This puts the value s throughout both the matrix and the vector.
00526 //
00527    if (!internalFEIsAllocated_) return(0);
00528 
00529    CHK_ERR( filter_[index_current_filter_]->resetSystem(s))
00530  
00531    return(0);
00532 }
00533 
00534 //------------------------------------------------------------------------------
00535 int FEI_Implementation::resetMatrix(double s)
00536 {
00537    if (!internalFEIsAllocated_) return(0);
00538 
00539    CHK_ERR( filter_[index_current_filter_]->resetMatrix(s))
00540 
00541    return(0);
00542 }
00543 
00544 //------------------------------------------------------------------------------
00545 int FEI_Implementation::resetRHSVector(double s)
00546 {
00547    if (!internalFEIsAllocated_) return(0);
00548 
00549    CHK_ERR( filter_[index_current_filter_]->resetRHSVector(s))
00550 
00551    return(0);
00552 }
00553 
00554 //------------------------------------------------------------------------------
00555 int FEI_Implementation::resetInitialGuess(double s)
00556 {
00557    if (!internalFEIsAllocated_) return(0);
00558 
00559    CHK_ERR( filter_[index_current_filter_]->resetInitialGuess(s))
00560 
00561    return(0);
00562 }
00563 
00564 //------------------------------------------------------------------------------
00565 int FEI_Implementation::loadNodeBCs(int numNodes,
00566                                     const GlobalID *nodeIDs,
00567                                     int fieldID,
00568                                     const int* offsetsIntoField,
00569                                     const double* prescribedValues)
00570 {
00571    if (!internalFEIsAllocated_)
00572       notAllocatedAbort("FEI_Implementation::loadNodeBCs");
00573 
00574    int index = index_current_filter_;
00575    if (solveType_ == 2) index = index_soln_filter_;
00576 
00577    CHK_ERR( filter_[index]->loadNodeBCs(numNodes,
00578                                nodeIDs, fieldID,
00579                                offsetsIntoField, prescribedValues));
00580 
00581    return(0);
00582 }
00583 
00584 //------------------------------------------------------------------------------
00585 int FEI_Implementation::loadElemBCs(int numElems,
00586                                     const GlobalID *elemIDs,
00587                                     int fieldID,
00588                                     const double *const *alpha,
00589                                     const double *const *beta,
00590                                     const double *const *gamma)
00591 {
00592    if (!internalFEIsAllocated_)
00593       notAllocatedAbort("FEI_Implementation::loadElemBCs");
00594 
00595    int index = index_current_filter_;
00596    if (solveType_ == 2) index = index_soln_filter_;
00597 
00598    CHK_ERR( filter_[index]->loadElemBCs(numElems,
00599                                elemIDs, fieldID,
00600                                alpha, beta, gamma))
00601 
00602    return(0);
00603 }
00604 
00605 //------------------------------------------------------------------------------
00606 int FEI_Implementation::sumInElem(GlobalID elemBlockID,
00607                         GlobalID elemID,
00608                         const GlobalID* elemConn,
00609                         const double* const* elemStiffness,
00610                         const double* elemLoad,
00611                         int elemFormat)
00612 {
00613   if (!internalFEIsAllocated_) {
00614     notAllocatedAbort("FEI_Implementation::sumInElem");
00615   }
00616 
00617   CHK_ERR( filter_[index_current_filter_]->sumInElem(elemBlockID, elemID,
00618                                                      elemConn, elemStiffness,
00619                                                      elemLoad, elemFormat));
00620 
00621   newMatrixDataLoaded_ = 1;
00622 
00623   return(0);
00624 }
00625 
00626 //------------------------------------------------------------------------------
00627 int FEI_Implementation::sumInElemMatrix(GlobalID elemBlockID,
00628                         GlobalID elemID,
00629                         const GlobalID* elemConn,
00630                         const double* const* elemStiffness,
00631                         int elemFormat)
00632 {
00633   if (!internalFEIsAllocated_)
00634     notAllocatedAbort("FEI_Implementation::sumInElemMatrix");
00635 
00636   CHK_ERR( filter_[index_current_filter_]->sumInElemMatrix(elemBlockID,
00637                                           elemID, elemConn,
00638                                           elemStiffness, elemFormat))
00639 
00640   newMatrixDataLoaded_ = 1;
00641 
00642 
00643   return(0);
00644 }
00645 
00646 //------------------------------------------------------------------------------
00647 int FEI_Implementation::sumInElemRHS(GlobalID elemBlockID,
00648                         GlobalID elemID,
00649                         const GlobalID* elemConn,
00650                         const double* elemLoad)
00651 {
00652    if (!internalFEIsAllocated_)
00653       notAllocatedAbort("FEI_Implementation::sumInElemRHS");
00654 
00655    CHK_ERR( filter_[index_current_rhs_row_]->sumInElemRHS(elemBlockID,
00656                                           elemID, elemConn, elemLoad))
00657 
00658    newMatrixDataLoaded_ = 1;
00659 
00660    return(0);
00661 }
00662 
00663 //------------------------------------------------------------------------------
00664 int FEI_Implementation::loadCRMult(int CRID,
00665                                    int numCRNodes,
00666                                    const GlobalID* CRNodes,
00667                                    const int* CRFields,
00668                                    const double* CRWeights,
00669                                    double CRValue)
00670 {
00671    if (!internalFEIsAllocated_)
00672       notAllocatedAbort("FEI_Implementation::loadCRMult");
00673 
00674    newMatrixDataLoaded_ = 1;
00675 
00676    CHK_ERR( filter_[index_current_filter_]->loadCRMult(CRID,
00677                                            numCRNodes, CRNodes,
00678                                            CRFields, CRWeights, CRValue));
00679 
00680    return(0);
00681 }
00682 
00683 //------------------------------------------------------------------------------
00684 int FEI_Implementation::loadCRPen(int CRID,
00685                                   int numCRNodes,
00686                                   const GlobalID* CRNodes,
00687                                   const int* CRFields,
00688                                   const double* CRWeights,
00689                                   double CRValue,
00690                                   double penValue)
00691 {
00692    if (!internalFEIsAllocated_)
00693       notAllocatedAbort("FEI_Implementation::loadCRPen");
00694 
00695    CHK_ERR( filter_[index_current_filter_]->loadCRPen(CRID,
00696                                           numCRNodes, CRNodes,
00697                                           CRFields, CRWeights,
00698                                           CRValue, penValue))
00699 
00700    newMatrixDataLoaded_ = 1;
00701 
00702    return(0);
00703 }
00704 
00705 //------------------------------------------------------------------------------
00706 int FEI_Implementation::sumIntoRHS(int IDType,
00707                                    int fieldID,
00708                                    int numIDs,
00709                                    const GlobalID* IDs,
00710                                    const double* rhsEntries)
00711 {
00712   if (!internalFEIsAllocated_)
00713     notAllocatedAbort("FEI_Implementation::sumIntoRHS");
00714 
00715   CHK_ERR( filter_[index_current_rhs_row_]->sumIntoRHS(IDType, fieldID,
00716                                                        numIDs, IDs, rhsEntries) );
00717   newMatrixDataLoaded_ = 1;
00718 
00719   return(0);
00720 }
00721 
00722 //------------------------------------------------------------------------------
00723 int FEI_Implementation::sumIntoMatrixDiagonal(int IDType,
00724                              int fieldID,
00725                              int numIDs,
00726                              const GlobalID* IDs,
00727                              const double* coefficients)
00728 {
00729   if (!internalFEIsAllocated_)
00730     notAllocatedAbort("FEI_Implementation::sumIntoMatrixDiagonal");
00731 
00732   CHK_ERR( filter_[index_current_filter_]->sumIntoMatrixDiagonal(IDType, fieldID,
00733                                                        numIDs, IDs, coefficients) );
00734   newMatrixDataLoaded_ = 1;
00735 
00736   return(0);
00737 }
00738 
00739 //------------------------------------------------------------------------------
00740 int FEI_Implementation::putIntoRHS(int IDType,
00741                                    int fieldID,
00742                                    int numIDs,
00743                                    const GlobalID* IDs,
00744                                    const double* rhsEntries)
00745 {
00746   if (!internalFEIsAllocated_)
00747     notAllocatedAbort("FEI_Implementation::putIntoRHS");
00748 
00749   CHK_ERR( filter_[index_current_rhs_row_]->putIntoRHS(IDType, fieldID,
00750                                                        numIDs, IDs, rhsEntries) );
00751   newMatrixDataLoaded_ = 1;
00752 
00753   return(0);
00754 }
00755 
00756 //------------------------------------------------------------------------------
00757 int FEI_Implementation::setMatScalars(int numScalars,
00758                                       const int* IDs, const double* scalars)
00759 {
00760    for(int i=0; i<numScalars; i++){
00761       std::vector<int>::iterator iter =
00762           std::find(matrixIDs_.begin(), matrixIDs_.end(), IDs[i]);
00763       if (iter != matrixIDs_.end()) {
00764          int index = iter - matrixIDs_.begin();
00765          matScalars_[index] = scalars[i];
00766       }
00767       else {
00768          FEI_CERR << "FEI_Implementation::setMatScalars: ERROR, invalid ID supplied"
00769               << FEI_ENDL;
00770          return(1);
00771       }
00772    }
00773 
00774    aggregateSystemFormed_ = false;
00775    matScalarsSet_ = true;
00776 
00777    return(0);
00778 }
00779 
00780 //------------------------------------------------------------------------------
00781 int FEI_Implementation::setRHSScalars(int numScalars,
00782                                       const int* IDs, const double* scalars)
00783 {
00784   for(int i=0; i<numScalars; i++){
00785      bool found = false;
00786 
00787      for(int j=0; j<numInternalFEIs_; j++){
00788          int index = fei::searchList(IDs[i], rhsIDs_[j], numRHSIDs_[j]);
00789          if (index>=0) {
00790              rhsScalars_[j][index] = scalars[i];
00791              found = true;
00792              break;
00793          }
00794       }
00795 
00796       if (!found) {
00797          FEI_CERR << "FEI_Implementation::setRHSScalars: ERROR, invalid RHS ID supplied"
00798              << FEI_ENDL;
00799          return(1);
00800       }
00801    }
00802 
00803    aggregateSystemFormed_ = false;
00804    rhsScalarsSet_ = true;
00805 
00806    return(0);
00807 }
00808 
00809 //------------------------------------------------------------------------------
00810 int FEI_Implementation::parameters(int numParams, const char *const* paramStrings)
00811 {
00812   // this function takes parameters and passes them to the internal
00813   // fei objects.
00814 
00815   if (numParams == 0 || paramStrings == NULL) {
00816     debugOut("#--- no parameters");
00817     return(0);
00818   }
00819 
00820   // merge these parameters with any others we may have, for later use.
00821   snl_fei::mergeStringLists(paramStrings_, numParams_,
00822                                    paramStrings, numParams);
00823 
00824   snl_fei::getIntParamValue("numMatrices", numParams,paramStrings, numInternalFEIs_);
00825 
00826   snl_fei::getIntParamValue("outputLevel", numParams,paramStrings, outputLevel_);
00827 
00828   const char* param = snl_fei::getParamValue("debugOutput",numParams,paramStrings);
00829   if (param != NULL) {
00830     setDebugOutput(param,"FEI_log");
00831   }
00832 
00833   if (debugOutput_) {
00834     (*dbgOStreamPtr_)<<"FEI: parameters"<<FEI_ENDL;
00835     (*dbgOStreamPtr_)<<"#FEI_Implementation, num-params "<<FEI_ENDL
00836                      <<numParams<<FEI_ENDL;
00837     (*dbgOStreamPtr_)<<"# "<<numParams<<" parameter lines follow:"<<FEI_ENDL;
00838     for(int i=0; i<numParams; i++){
00839       (*dbgOStreamPtr_)<<paramStrings[i]<<FEI_ENDL;
00840     }
00841   }
00842 
00843   if (haveLinSysCore_) {
00844     linSysCore_->parameters(numParams, (char**)paramStrings);
00845   }
00846   if (haveFEData_) {
00847     wrapper_->getFiniteElementData()->parameters(numParams, (char**)paramStrings);
00848   }
00849 
00850   problemStructure_->parameters(numParams, paramStrings);
00851 
00852   if (internalFEIsAllocated_){
00853     for(int i=0; i<numInternalFEIs_; i++){
00854       CHK_ERR( filter_[i]->parameters(numParams, paramStrings) );
00855     }
00856   }
00857 
00858   debugOut("#FEI_Implementation leaving parameters method");
00859 
00860   return(0);
00861 }
00862 
00863 //------------------------------------------------------------------------------
00864 void FEI_Implementation::setDebugOutput(const char* path, const char* name)
00865 {
00866   //
00867   //This function turns on debug output, and opens a file to put it in.
00868   //
00869   if (dbgFileOpened_) {
00870     dbgFStreamPtr_->close();
00871   }
00872 
00873   dbgFileOpened_ = false;
00874   delete dbgOStreamPtr_;
00875 
00876   FEI_OSTRINGSTREAM osstr;
00877   if (path != NULL) {
00878     osstr << path << "/";
00879   }
00880   osstr << name << "." << numProcs_ << "." << localRank_;
00881 
00882   debugOutput_ = 1;
00883   dbgFStreamPtr_ = new FEI_OFSTREAM(osstr.str().c_str(), IOS_APP);
00884   if (!dbgFStreamPtr_ || dbgFStreamPtr_->bad()){
00885     FEI_CERR << "couldn't open debug output file: " << osstr.str() << FEI_ENDL;
00886     debugOutput_ = 0;
00887   }
00888 
00889   if (debugOutput_) {
00890     const char* version_str = NULL;
00891     version(version_str);
00892 
00893     (*dbgFStreamPtr_) << version_str << FEI_ENDL;
00894 
00895     problemStructure_->setDbgOut(*dbgFStreamPtr_, path, "_0");
00896     dbgOStreamPtr_ = dbgFStreamPtr_;
00897     dbgOStreamPtr_->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00898     dbgFileOpened_ = true;
00899 
00900     if (internalFEIsAllocated_) {
00901       for(int i=0; i<numInternalFEIs_; ++i) {
00902         filter_[i]->setLogStream(dbgOStreamPtr_);
00903       }
00904     }
00905   }
00906 }
00907 
00908 //------------------------------------------------------------------------------
00909 int FEI_Implementation::loadComplete(bool applyBCs,
00910                                      bool globalAssemble)
00911 {
00912   (void)applyBCs;
00913   (void)globalAssemble;
00914 
00915   buildLinearSystem();
00916 
00917   return(0);
00918 }
00919 
00920 //------------------------------------------------------------------------------
00921 int FEI_Implementation::residualNorm(int whichNorm, int numFields,
00922                                      int* fieldIDs, double* norms)
00923 {
00924    buildLinearSystem();
00925 
00926    double residTime = 0.0;
00927 
00928    int err = filter_[index_soln_filter_]->residualNorm(whichNorm, numFields,
00929                                                  fieldIDs, norms, residTime);
00930 
00931    solveTime_ += residTime;
00932 
00933    return(err);
00934 }
00935 
00936 //------------------------------------------------------------------------------
00937 int FEI_Implementation::solve(int& status)
00938 {
00939    buildLinearSystem();
00940 
00941    double sTime = 0.0;
00942 
00943    int err = filter_[index_soln_filter_]->solve(status, sTime);
00944 
00945    solveTime_ += sTime;
00946 
00947    return(err);
00948 }
00949 
00950 //------------------------------------------------------------------------------
00951 int FEI_Implementation::iterations(int& itersTaken) const {
00952   itersTaken = filter_[index_soln_filter_]->iterations();
00953   return(0);
00954 }
00955 
00956 //------------------------------------------------------------------------------
00957 int FEI_Implementation::version(const char*& versionString)
00958 {
00959   versionString = fei::utils::version();
00960   return(0);
00961 }
00962 
00963 //------------------------------------------------------------------------------
00964 int FEI_Implementation::cumulative_cpu_times(double& initTime,
00965                                              double& loadTime,
00966                                              double& solveTime,
00967                                              double& solnReturnTime)
00968 {
00969    initTime = initTime_;
00970    loadTime = loadTime_;
00971    solveTime = solveTime_;
00972    solnReturnTime = solnReturnTime_;
00973 
00974    return(0);
00975 }
00976 
00977 //------------------------------------------------------------------------------
00978 int FEI_Implementation::getBlockNodeSolution(GlobalID elemBlockID,
00979                                              int numNodes,
00980                                              const GlobalID *nodeIDs,
00981                                              int *offsets,
00982                                              double *results)
00983 {
00984    CHK_ERR(filter_[index_soln_filter_]->getBlockNodeSolution(elemBlockID,
00985                                                numNodes,
00986                                                nodeIDs,
00987                                                offsets,
00988                                                results))
00989    return(0);
00990 }
00991 
00992 //------------------------------------------------------------------------------
00993 int FEI_Implementation::getNodalSolution(int numNodes,
00994                                          const GlobalID *nodeIDs,
00995                                          int *offsets,
00996                                          double *results)
00997 {
00998    CHK_ERR(filter_[index_soln_filter_]->getNodalSolution(numNodes,
00999                                                          nodeIDs,
01000                                                          offsets,
01001                                                          results))
01002 
01003    return(0);
01004 }
01005 
01006 //------------------------------------------------------------------------------
01007 int FEI_Implementation::getBlockFieldNodeSolution(GlobalID elemBlockID,
01008                                         int fieldID,
01009                                         int numNodes, 
01010                                         const GlobalID *nodeIDs, 
01011                                         double *results)
01012 {
01013    CHK_ERR( filter_[index_soln_filter_]->getBlockFieldNodeSolution(elemBlockID,
01014                                                      fieldID, numNodes,
01015                                                      nodeIDs, results))
01016 
01017    return(0);
01018 }
01019 
01020 //------------------------------------------------------------------------------
01021 int FEI_Implementation::putBlockNodeSolution(GlobalID elemBlockID,
01022                                              int numNodes,
01023                                              const GlobalID *nodeIDs,
01024                                              const int *offsets,
01025                                              const double *estimates)
01026 {
01027    CHK_ERR( filter_[index_soln_filter_]->putBlockNodeSolution(elemBlockID,
01028                                                 numNodes, nodeIDs,
01029                                                 offsets, estimates))
01030    return(0);
01031 }
01032 
01033 //------------------------------------------------------------------------------
01034 int FEI_Implementation::putBlockFieldNodeSolution(GlobalID elemBlockID,
01035                                         int fieldID,
01036                                         int numNodes,
01037                                         const GlobalID *nodeIDs,
01038                                         const double *estimates)
01039 {
01040    int err = filter_[index_soln_filter_]->putBlockFieldNodeSolution(elemBlockID,
01041                                                      fieldID, numNodes,
01042                                                      nodeIDs, estimates);
01043    return(err);
01044 }
01045 
01046 //------------------------------------------------------------------------------
01047 int FEI_Implementation::getBlockElemSolution(GlobalID elemBlockID,  
01048                                    int numElems, 
01049                                    const GlobalID *elemIDs,
01050                                    int& numElemDOFPerElement,
01051                                    double *results)
01052 {
01053    CHK_ERR( filter_[index_soln_filter_]->getBlockElemSolution(elemBlockID,
01054                                                        numElems, elemIDs,
01055                                                        numElemDOFPerElement,
01056                                                        results))
01057     return(0);
01058 } 
01059       
01060 //------------------------------------------------------------------------------
01061 int FEI_Implementation::putBlockElemSolution(GlobalID elemBlockID,
01062                                    int numElems,
01063                                    const GlobalID *elemIDs,
01064                                    int dofPerElem,
01065                                    const double *estimates)
01066 {
01067    CHK_ERR( filter_[index_soln_filter_]->putBlockElemSolution(elemBlockID,
01068                                                        numElems, elemIDs,
01069                                                        dofPerElem, estimates))
01070    return(0);
01071 }
01072 
01073 //------------------------------------------------------------------------------
01074 int FEI_Implementation::getNumCRMultipliers(int& numMultCRs)
01075 {
01076   numMultCRs = problemStructure_->getNumMultConstRecords();
01077   return(0);
01078 }
01079 
01080 //------------------------------------------------------------------------------
01081 int FEI_Implementation::getCRMultIDList(int numMultCRs,
01082                                         int* multIDs)
01083 {
01084    if (numMultCRs > problemStructure_->getNumMultConstRecords()) return(-1);
01085 
01086    std::map<GlobalID,snl_fei::Constraint<GlobalID>*>::const_iterator
01087      cr_iter = problemStructure_->getMultConstRecords().begin(),
01088      cr_end = problemStructure_->getMultConstRecords().end();
01089    int i = 0;
01090    while(cr_iter != cr_end) {
01091       multIDs[i++] = (*cr_iter).first;
01092       ++cr_iter;
01093    }
01094 
01095    return(0);
01096 }
01097 
01098 //------------------------------------------------------------------------------
01099 int FEI_Implementation::getCRMultipliers(int numMultCRs,
01100                                          const int* CRIDs,
01101                                          double* multipliers)
01102 {
01103    CHK_ERR( filter_[index_soln_filter_]->getCRMultipliers(numMultCRs,
01104                                                      CRIDs, multipliers))
01105    return(0);
01106 }
01107 
01108 //------------------------------------------------------------------------------
01109 int FEI_Implementation::putCRMultipliers(int numMultCRs,
01110                                          const int* CRIDs,
01111                                          const double *multEstimates)
01112 {
01113     return(
01114     filter_[index_soln_filter_]->putCRMultipliers(numMultCRs,
01115                                             CRIDs,
01116                                             multEstimates)
01117     );
01118 }
01119 
01120 //-----------------------------------------------------------------------------
01121 //  some utility functions to aid in using the "put" functions for passing
01122 //  an initial guess to the solver
01123 //-----------------------------------------------------------------------------
01124 
01125 //------------------------------------------------------------------------------
01126 int FEI_Implementation::getBlockElemIDList(GlobalID elemBlockID,
01127                                            int numElems,
01128                                            GlobalID* elemIDs)
01129 {
01130   //
01131   //  return the list of element IDs for a given block... the length parameter
01132   //  lenElemIDList should be used to check memory allocation for the calling
01133   //  method, as the calling method should have gotten a copy of this param
01134   //  from a call to getNumBlockElements before allocating memory for elemIDList
01135   //
01136   ConnectivityTable& connTable = problemStructure_->
01137     getBlockConnectivity(elemBlockID);
01138 
01139   std::map<GlobalID,int>& elemIDList = connTable.elemIDs;
01140   int len = elemIDList.size();
01141   if (len > numElems) len = numElems;
01142 
01143   fei::copyKeysToArray(elemIDList, len, elemIDs);
01144 
01145   return(FEI_SUCCESS);
01146 }
01147 
01148 //------------------------------------------------------------------------------
01149 int FEI_Implementation::getBlockNodeIDList(GlobalID elemBlockID,
01150                                            int numNodes,
01151                                            GlobalID *nodeIDs)
01152 {
01153   //
01154   //  similar comments as for getBlockElemIDList(), except for returning the
01155   //  active node list
01156   //
01157   int numActiveNodes = problemStructure_->getNumActiveNodes();
01158   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01159 
01160   int offset = 0;
01161   for(int i=0; i<numActiveNodes; i++) {
01162     const NodeDescriptor* node = NULL;
01163     CHK_ERR( nodeDB.getNodeAtIndex(i, node) );
01164     if (node->containedInBlock(elemBlockID))
01165       nodeIDs[offset++] = node->getGlobalNodeID();
01166     if (offset == numNodes) break;
01167   }
01168 
01169   return(FEI_SUCCESS);
01170 }
01171 
01172 //------------------------------------------------------------------------------
01173 int FEI_Implementation::getNumNodesPerElement(GlobalID blockID,
01174                                               int& nodesPerElem) const
01175 {
01176   //
01177   //  return the number of nodes associated with elements of a given block ID
01178   //
01179   BlockDescriptor* block = NULL;
01180   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01181 
01182   nodesPerElem = block->getNumNodesPerElement();
01183   return(FEI_SUCCESS);
01184 }
01185  
01186 //------------------------------------------------------------------------------
01187 int FEI_Implementation::getNumEqnsPerElement(GlobalID blockID, int& numEqns)
01188 const
01189 {
01190   //
01191   //  return the number of eqns associated with elements of a given block ID
01192   //
01193   BlockDescriptor* block = NULL;
01194   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01195 
01196   numEqns = block->getNumEqnsPerElement();
01197   return(FEI_SUCCESS);
01198 }
01199 
01200 //------------------------------------------------------------------------------
01201 int FEI_Implementation::getNumSolnParams(GlobalID nodeID, int& numSolnParams)
01202 const {
01203   //
01204   //  return the number of solution parameters at a given node
01205   //
01206   const NodeDescriptor* node = NULL;
01207   int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
01208 
01209   if (err != 0) {
01210     ERReturn(-1);
01211   }
01212 
01213   numSolnParams = node->getNumNodalDOF();
01214   return(0);
01215 }
01216  
01217 //------------------------------------------------------------------------------
01218 int FEI_Implementation::getNumElemBlocks(int& numElemBlocks) const
01219 {
01220   //
01221   //  the number of element blocks
01222   //
01223   numElemBlocks = problemStructure_->getNumElemBlocks();
01224   return( 0 );
01225 }
01226 
01227 //------------------------------------------------------------------------------
01228 int FEI_Implementation::getNumBlockActNodes(GlobalID blockID, int& numNodes)
01229 const {
01230   //
01231   //  return the number of active nodes associated with a given element block ID
01232   //
01233   BlockDescriptor* block = NULL;
01234   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01235 
01236   numNodes = block->getNumActiveNodes();
01237   return(FEI_SUCCESS);
01238 }
01239 
01240 //------------------------------------------------------------------------------
01241 int FEI_Implementation::getNumBlockActEqns(GlobalID blockID, int& numEqns)
01242 const {
01243 //
01244 // return the number of active equations associated with a given element
01245 // block ID
01246 //
01247   BlockDescriptor* block = NULL;
01248   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01249 
01250   numEqns = block->getTotalNumEqns();
01251   return(FEI_SUCCESS);
01252 }
01253 
01254 //------------------------------------------------------------------------------
01255 int FEI_Implementation::getNumBlockElements(GlobalID blockID, int& numElems) const {
01256 //
01257 //  return the number of elements associated with a given elem blockID
01258 //
01259   BlockDescriptor* block = NULL;
01260   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01261 
01262   numElems = block->getNumElements();
01263   return(FEI_SUCCESS);
01264 }
01265 
01266 //------------------------------------------------------------------------------
01267 int FEI_Implementation::getNumBlockElemDOF(GlobalID blockID,
01268                                            int& DOFPerElem) const
01269 {
01270   //
01271   //  return the number of elem equations associated with a given elem blockID
01272   //
01273   BlockDescriptor* block = NULL;
01274   CHK_ERR( problemStructure_->getBlockDescriptor(blockID, block) );
01275 
01276   DOFPerElem = block->getNumElemDOFPerElement();
01277 
01278   return(FEI_SUCCESS);
01279 }
01280 
01281 //------------------------------------------------------------------------------
01282 int FEI_Implementation::getFieldSize(int fieldID,
01283                                      int& numScalars)
01284 {
01285   //
01286   //  return the number of scalars associated with a given fieldID
01287   //
01288 
01289   numScalars = problemStructure_->getFieldSize(fieldID);
01290   return(0);  
01291 }
01292 
01293 //------------------------------------------------------------------------------
01294 int FEI_Implementation::getEqnNumbers(GlobalID ID, int idType,
01295                                       int fieldID, int& numEqns,
01296                                       int* eqnNumbers)
01297 {
01298   //
01299   // Translate from an ID/fieldID pair to a list of equation-numbers
01300   //
01301 
01302   return( problemStructure_->getEqnNumbers(ID, idType, fieldID,
01303                                            numEqns, eqnNumbers) );
01304 }
01305 
01306 //------------------------------------------------------------------------------
01307 int FEI_Implementation::getNodalFieldSolution(int fieldID,
01308                                               int numNodes,
01309                                               const GlobalID* nodeIDs,
01310                                               double* results)
01311 {
01312   return( filter_[index_soln_filter_]->getNodalFieldSolution(fieldID, numNodes,
01313                                                        nodeIDs, results) );
01314 }
01315 
01316 //------------------------------------------------------------------------------
01317 int FEI_Implementation::getNumLocalNodes(int& numNodes)
01318 {
01319   numNodes = problemStructure_->getNodeDatabase().getNodeIDs().size();
01320   return( 0 );
01321 }
01322 
01323 //------------------------------------------------------------------------------
01324 int FEI_Implementation::getLocalNodeIDList(int& numNodes,
01325                                            GlobalID* nodeIDs,
01326                                            int lenNodeIDs)
01327 {
01328   std::map<GlobalID,int>& nodes =
01329     problemStructure_->getNodeDatabase().getNodeIDs();
01330   numNodes = nodes.size();
01331   int len = numNodes;
01332   if (lenNodeIDs < len) len = lenNodeIDs;
01333 
01334   fei::copyKeysToArray(nodes, len, nodeIDs);
01335 
01336   return( 0 );
01337 }
01338 
01339 //------------------------------------------------------------------------------
01340 int FEI_Implementation::putNodalFieldData(int fieldID,
01341                                           int numNodes,
01342                                           const GlobalID* nodeIDs,
01343                                           const double* nodeData)
01344 {
01345   return( filter_[index_soln_filter_]->putNodalFieldData(fieldID, numNodes,
01346                                                    nodeIDs, nodeData) );
01347 }
01348 
01349 //------------------------------------------------------------------------------
01350 void FEI_Implementation::buildLinearSystem()
01351 {
01352   //
01353   //Private function.
01354   //
01355   //At the point when this function is called, all matrix assembly has
01356   //already taken place, with the data having been directed into the
01357   //appropriate Filter instance in the filter_ list. Now it's
01358   //time to finalize the linear system (get A,x,b ready to give to a
01359   //solver), performing this list of last-minute tasks:
01360   //
01361   //1. Have each Filter instance exchangeRemoteEquations.
01362   //2. Aggregate the system (form a linear combination of LHS's and
01363   //   RHS's) if solveType_==2.
01364   //3. call loadComplete to have the 'master' Filter instance
01365   //   (filter_[index_soln_filter_]) enforce any boundary conditions
01366   //   that have been loaded.
01367   //
01368   debugOut("#   buildLinearSystem");
01369 
01370   //
01371   //figure out if new matrix-data was loaded on any processor.
01372   int anyDataLoaded = newMatrixDataLoaded_;
01373 #ifndef FEI_SER
01374   if (numProcs_ > 1) {
01375     if (MPI_Allreduce(&newMatrixDataLoaded_, &anyDataLoaded, 1, MPI_INT,
01376                       MPI_SUM, comm_) != MPI_SUCCESS) voidERReturn
01377   }
01378 #endif
01379 
01380   if (anyDataLoaded) {
01381 #ifndef FEI_SER
01382     for(int i=0; i<numInternalFEIs_; i++) {
01383       filter_[i]->exchangeRemoteEquations();
01384     }
01385 #endif
01386     newMatrixDataLoaded_ = 0;
01387   }
01388 
01389   if (solveType_ == 2){
01390     //solveType_ == 2 means this is a linear-combination solve --
01391     //i.e., we're solving an aggregate system which is the sum of
01392     //several individual matrices and rhs's.
01393 
01394     if (!aggregateSystemFormed_) {
01395       if (!matScalarsSet_ || !rhsScalarsSet_) {
01396         FEI_COUT << "FEI_Implementation: WARNING: solveType_==2, aggregating system, but setMatScalars and/or setRHSScalars not yet called." << FEI_ENDL;
01397       }
01398 
01399       aggregateSystem();
01400     }
01401   }
01402 
01403   filter_[index_soln_filter_]->loadComplete();
01404 
01405   debugOut("#   leaving buildLinearSystem");
01406 }
01407 
01408 //------------------------------------------------------------------------------
01409 int FEI_Implementation::aggregateSystem()
01410 {
01411   debugOut("#   aggregateSystem");
01412    if (!haveLinSysCore_) ERReturn(-1);
01413 
01414    if (soln_fei_matrix_ == NULL) {
01415       soln_fei_matrix_ = new Data();
01416 
01417       CHK_ERR( lscArray_[index_soln_filter_]->
01418                copyOutMatrix(1.0, *soln_fei_matrix_) );
01419    }
01420 
01421    if (soln_fei_vector_ == NULL) {
01422       soln_fei_vector_ = new Data();
01423 
01424       CHK_ERR( lscArray_[index_soln_filter_]->
01425                      setRHSID(rhsIDs_[index_soln_filter_][0]) );
01426 
01427       CHK_ERR( lscArray_[index_soln_filter_]->
01428                     copyOutRHSVector(1.0, *soln_fei_vector_) );
01429    }
01430 
01431    Data tmp;
01432    Data tmpv;
01433 
01434    for(int i=0; i<numInternalFEIs_; i++){
01435 
01436       if (i == index_soln_filter_) {
01437          tmp.setTypeName(soln_fei_matrix_->getTypeName());
01438          tmp.setDataPtr(soln_fei_matrix_->getDataPtr());
01439          CHK_ERR( lscArray_[index_soln_filter_]->
01440                         copyInMatrix(matScalars_[i], tmp) );
01441       }
01442       else {
01443          CHK_ERR( lscArray_[i]->getMatrixPtr(tmp) );
01444 
01445          CHK_ERR( lscArray_[index_soln_filter_]->
01446                         sumInMatrix(matScalars_[i], tmp) );
01447       }
01448 
01449       for(int j=0; j<numRHSIDs_[i]; j++){
01450          if ((i == index_soln_filter_) && (j == 0)) {
01451             tmpv.setTypeName(soln_fei_vector_->getTypeName());
01452             tmpv.setDataPtr(soln_fei_vector_->getDataPtr());
01453          }
01454          else {
01455             CHK_ERR( lscArray_[i]->setRHSID(rhsIDs_[i][j]) );
01456             CHK_ERR( lscArray_[i]->getRHSVectorPtr(tmpv) );
01457          }
01458 
01459          if (i == index_soln_filter_) {
01460             CHK_ERR( lscArray_[index_soln_filter_]->
01461                          copyInRHSVector(rhsScalars_[i][j], tmpv) );
01462          }
01463          else {
01464             CHK_ERR( lscArray_[index_soln_filter_]->
01465                          sumInRHSVector(rhsScalars_[i][j], tmpv) );
01466          }
01467       }
01468    }
01469 
01470    aggregateSystemFormed_ = true;
01471 
01472    debugOut("#   leaving aggregateSystem");
01473 
01474    return(0);
01475 }
01476 
01477 //==============================================================================
01478 int FEI_Implementation::allocateInternalFEIs(){
01479 //
01480 //This is a private FEI_Implementation function, to be called from within
01481 //setSolveType or the other overloading of allocateInternalFEIs.
01482 //Assumes that numInternalFEIs_ has already been set.
01483 //It is also safe to assume that problemStructure_->initComplete() has already
01484 //been called.
01485 //
01486 
01487    if (internalFEIsAllocated_) return(0);
01488 
01489    matScalars_.resize(numInternalFEIs_);
01490 
01491    if (rhsScalars_.size() != 0) deleteRHSScalars();
01492 
01493    rhsScalars_.resize(numInternalFEIs_);
01494 
01495    for(int i=0; i<numInternalFEIs_; i++){
01496       matScalars_[i] = 1.0;
01497 
01498       rhsScalars_[i] = new double[numRHSIDs_[i]];
01499 
01500       for(int j=0; j<numRHSIDs_[i]; j++){
01501          rhsScalars_[i][j] = 1.0;
01502       }
01503    }
01504 
01505    IDsAllocated_ = true;
01506 
01507    if (numInternalFEIs_ > 0) {
01508       index_soln_filter_ = 0;
01509       index_current_filter_ = 0;
01510       filter_ = new Filter*[numInternalFEIs_];
01511 
01512       if (haveLinSysCore_) {
01513         if (numRHSIDs_[0] == 0) {
01514           int dummyID = -1;
01515           linSysCore_->setNumRHSVectors(1, &dummyID);
01516         }
01517         else {
01518           linSysCore_->setNumRHSVectors(numRHSIDs_[0], rhsIDs_[0]);
01519         }
01520 
01521         for(int i=1; i<numInternalFEIs_; i++) {
01522           fei::SharedPtr<LinearSystemCore> lsc(linSysCore_->clone());
01523           lsc->parameters(numParams_, paramStrings_);
01524 
01525           if (numRHSIDs_[i] == 0) {
01526             int dummyID = -1;
01527             lsc->setNumRHSVectors(1, &dummyID);
01528           }
01529           else {
01530             lsc->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
01531           }
01532 
01533           lscArray_.push_back(lsc);
01534         }
01535       }
01536 
01537       for(int i=0; i<numInternalFEIs_; i++){
01538 
01539         if (haveLinSysCore_) {
01540           filter_[i] = new LinSysCoreFilter(this, comm_, problemStructure_,
01541                                             lscArray_[i].get(), masterRank_);
01542         }
01543         else if (haveFEData_) {
01544           filter_[i] = new FEDataFilter(this, comm_, problemStructure_,
01545                                         wrapper_.get(), masterRank_);
01546         }
01547         else {
01548           FEI_CERR << "FEI_Implementation: ERROR, don't have LinearSystemCore"
01549                << " or FiniteElementData implementation..." << FEI_ENDL;
01550           ERReturn(-1);
01551         }
01552 
01553         filter_[i]->setLogStream(dbgOStreamPtr_);
01554 
01555         FEI_OSTRINGSTREAM osstr;
01556         osstr<<"internalFei "<< i;
01557         std::string osstr_str = osstr.str();
01558         const char* param = osstr_str.c_str();
01559         filter_[i]->parameters(1, &param);
01560 
01561         if (debugOutput_) {
01562           (*dbgOStreamPtr_)<<"#-- fei["<<i<<"]->setNumRHSVectors "
01563                            <<numRHSIDs_[i]<<FEI_ENDL;
01564         }
01565 
01566         if (numRHSIDs_[i] == 0) {
01567           int dummyID = -1;
01568           filter_[i]->setNumRHSVectors(1, &dummyID);
01569         }
01570         else {
01571           filter_[i]->setNumRHSVectors(numRHSIDs_[i], rhsIDs_[i]);
01572         }
01573       }
01574 
01575       internalFEIsAllocated_ = true;
01576    }
01577    else {
01578      needParametersAbort("FEI_Implementation::allocateInternalFEIs");
01579    }
01580 
01581    return(0);
01582 }
01583 
01584 //==============================================================================
01585 void FEI_Implementation::debugOut(const char* msg) {
01586   if (debugOutput_) { (*dbgOStreamPtr_)<<msg<<FEI_ENDL; }
01587 }
01588 
01589 //==============================================================================
01590 void FEI_Implementation::debugOut(const char* msg, int whichFEI) {
01591    if (debugOutput_) {
01592      (*dbgOStreamPtr_)<<msg<<", -> fei["<<whichFEI<<"]"<<FEI_ENDL;
01593    }
01594 }
01595 
01596 //==============================================================================
01597 void FEI_Implementation::messageAbort(const char* msg){
01598 
01599     FEI_CERR << "FEI_Implementation: ERROR " << msg << " Aborting." << FEI_ENDL;
01600     MPI_Abort(comm_, -1);
01601 }
01602 
01603 //==============================================================================
01604 void FEI_Implementation::notAllocatedAbort(const char* name){
01605 
01606     FEI_CERR << name
01607          << FEI_ENDL << "ERROR, internal data structures not allocated."
01608          << FEI_ENDL << "'setIDLists' and/or 'setSolveType' must be called"
01609          << FEI_ENDL << "first to identify solveType and number of matrices"
01610          << FEI_ENDL << "to be assembled." << FEI_ENDL;
01611     MPI_Abort(comm_, -1);
01612 }
01613 
01614 //==============================================================================
01615 void FEI_Implementation::needParametersAbort(const char* name){
01616 
01617    FEI_CERR << name
01618      << FEI_ENDL << "FEI_Implementation: ERROR, numMatrices has not been specified."
01619      << FEI_ENDL << "FEI_Implementation: 'parameters' must be called up front with"
01620      << FEI_ENDL << "FEI_Implementation: the string 'numMatrices n' to specify that"
01621      << FEI_ENDL << "FEI_Implementation: n matrices will be assembled." << FEI_ENDL;
01622     MPI_Abort(comm_, -1);
01623 }
01624 
01625 //==============================================================================
01626 void FEI_Implementation::badParametersAbort(const char* name){
01627 
01628    FEI_CERR << name
01629         << FEI_ENDL << "FEI_Implementation: ERROR, inconsistent 'solveType' and"
01630         << FEI_ENDL << "FEI_Implementation: 'numMatrices' parameters specified."
01631         << FEI_ENDL << "FEI_Implementation: Aborting."
01632         << FEI_ENDL;
01633     MPI_Abort(comm_, -1);
01634 }
01635 

Generated on Tue Jul 13 09:27:45 2010 for FEI by  doxygen 1.4.7