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