fei_FEI_Impl.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_macros.hpp>
00010 
00011 #include <fei_utils.hpp>
00012 
00013 #include <fei_FEI_Impl.hpp>
00014 #include <fei_Record.hpp>
00015 #include <fei_TemplateUtils.hpp>
00016 #include <fei_ParameterSet.hpp>
00017 #include <fei_base.hpp>
00018 
00019 #include <fei_Pattern.hpp>
00020 #include <fei_LibraryWrapper.hpp>
00021 #include <fei_Data.hpp>
00022 #include <fei_defs.h>
00023 
00024 #include <stdexcept>
00025 #include <cmath>
00026 
00027 #undef fei_file
00028 #define fei_file "fei_FEI_Impl.cpp"
00029 
00030 #include <fei_ErrMacros.hpp>
00031 
00032 fei::FEI_Impl::FEI_Impl(fei::SharedPtr<LibraryWrapper> wrapper,
00033           MPI_Comm comm,
00034           int masterRank)
00035   : wrapper_(1),
00036     nodeIDType_(0),
00037     elemIDType_(1),
00038     constraintIDType_(2),
00039     factory_(1),
00040     rowSpace_(NULL),
00041     matGraph_(),
00042     x_(),
00043     b_(),
00044     A_(),
00045     linSys_(),
00046     newData_(false),
00047     soln_fei_matrix_(NULL),
00048     soln_fei_vector_(NULL),
00049     comm_(comm),
00050     masterRank_(masterRank),
00051     localProc_(0),
00052     numProcs_(1),
00053     numParams_(0),
00054     paramStrings_(NULL),
00055     matrixIDs_(),
00056     rhsIDs_(),
00057     matScalars_(),
00058     matScalarsSet_(false),
00059     rhsScalars_(),
00060     rhsScalarsSet_(false),
00061     constraintID_(0),
00062     index_soln_(0),
00063     index_current_(0),
00064     index_current_rhs_row_(0),
00065     solveType_(0),
00066     iterations_(0),
00067     setSolveTypeCalled_(false),
00068     initPhaseIsComplete_(false),
00069     aggregateSystemFormed_(false),
00070     newMatrixDataLoaded_(0),
00071     solveCounter_(0),
00072     initTime_(0.0),
00073     loadTime_(0.0),
00074     solveTime_(0.0),
00075     solnReturnTime_(0.0),
00076     iwork_(),
00077     nodeset_(),
00078     nodeset_filled_(false),
00079     block_dof_per_elem_(),
00080     any_blocks_have_elem_dof_(false)
00081 {
00082   wrapper_[0] = wrapper;
00083 
00084   factory_[0].reset(new snl_fei::Factory(comm, wrapper_[0]));
00085   createdFactory_ = true;
00086 
00087   basic_initializations();
00088 }
00089 
00090 fei::FEI_Impl::FEI_Impl(const fei::Factory* factory,
00091           MPI_Comm comm,
00092           int masterRank)
00093   : wrapper_(1),
00094     nodeIDType_(0),
00095     elemIDType_(1),
00096     constraintIDType_(2),
00097     factory_(1),
00098     createdFactory_(false),
00099     rowSpace_(NULL),
00100     matGraph_(),
00101     x_(),
00102     b_(),
00103     A_(),
00104     linSys_(),
00105     newData_(false),
00106     soln_fei_matrix_(NULL),
00107     soln_fei_vector_(NULL),
00108     comm_(comm),
00109     masterRank_(masterRank),
00110     localProc_(0),
00111     numProcs_(1),
00112     numParams_(0),
00113     paramStrings_(NULL),
00114     matrixIDs_(),
00115     rhsIDs_(),
00116     matScalars_(),
00117     matScalarsSet_(false),
00118     rhsScalars_(),
00119     rhsScalarsSet_(false),
00120     constraintID_(0),
00121     index_soln_(0),
00122     index_current_(0),
00123     index_current_rhs_row_(0),
00124     solveType_(0),
00125     iterations_(0),
00126     setSolveTypeCalled_(false),
00127     initPhaseIsComplete_(false),
00128     aggregateSystemFormed_(false),
00129     newMatrixDataLoaded_(0),
00130     solveCounter_(0),
00131     initTime_(0.0),
00132     loadTime_(0.0),
00133     solveTime_(0.0),
00134     solnReturnTime_(0.0),
00135     iwork_(),
00136     nodeset_(),
00137     nodeset_filled_(false),
00138     block_dof_per_elem_(),
00139     any_blocks_have_elem_dof_(false)
00140 {
00141   wrapper_[0].reset(0);
00142 
00143   const snl_fei::Factory* snlfactory
00144     = dynamic_cast<const snl_fei::Factory*>(factory);
00145   if (snlfactory != NULL) {
00146     wrapper_[0] = snlfactory->get_LibraryWrapper();
00147   }
00148 
00149   factory_[0] = factory->clone();
00150 
00151   basic_initializations();
00152 }
00153 
00154 fei::FEI_Impl::~FEI_Impl()
00155 {
00156   for(int k=0; k<numParams_; k++) {
00157     delete [] paramStrings_[k];
00158   }
00159   delete [] paramStrings_;
00160 
00161   if (soln_fei_matrix_ != NULL && wrapper_[0].get() != NULL) {
00162     fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
00163     if (lsc.get() != NULL) {
00164       lsc->destroyMatrixData(*soln_fei_matrix_);
00165       delete soln_fei_matrix_;
00166       soln_fei_matrix_ = NULL;
00167     }
00168   }
00169 
00170   if (soln_fei_vector_ != NULL && wrapper_[0].get() != NULL) {
00171     fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
00172     if (lsc.get() != NULL) {
00173       lsc->destroyVectorData(*soln_fei_vector_);
00174       delete soln_fei_vector_;
00175       soln_fei_vector_ = NULL;
00176     }
00177   }
00178 }
00179 
00180 void fei::FEI_Impl::basic_initializations()
00181 {
00182   localProc_ = fei::localProc(comm_);
00183   numProcs_ = fei::numProcs(comm_);
00184 
00185   constraintID_ = localProc_*100000;
00186 
00187   matrixIDs_.resize(1);
00188   matrixIDs_[0] = 0;
00189   A_.resize(1);
00190   rhsIDs_.resize(1);
00191   rhsIDs_[0] = 0;
00192   b_.resize(1);
00193 
00194   rowSpace_ = factory_[0]->createVectorSpace(comm_, (const char*)NULL);
00195 
00196   rowSpace_->defineIDTypes(1, &nodeIDType_);
00197   rowSpace_->defineIDTypes(1, &elemIDType_);
00198   rowSpace_->defineIDTypes(1, &constraintIDType_);
00199 
00200   matGraph_ = factory_[0]->createMatrixGraph(rowSpace_, rowSpace_,(const char*)NULL);
00201   if (matGraph_.get() == NULL) {
00202     voidERReturn;
00203   }
00204 }
00205 
00206 int fei::FEI_Impl::setIDLists(int numMatrices,
00207              const int* matrixIDs,
00208              int numRHSs,
00209              const int* rhsIDs)
00210 {
00211   if (numMatrices < 1 || numRHSs < 1) {
00212     FEI_CERR << "fei::FEI_Impl::setIDLists ERROR, numMatrices and numRHSs "
00213    << "must both be greater than 0."<<FEI_ENDL;
00214     ERReturn(-1);
00215   }
00216 
00217   matrixIDs_.resize(0);
00218   A_.resize(numMatrices);
00219   for(int i=0; i<numMatrices; ++i) {
00220     fei::sortedListInsert(matrixIDs[i], matrixIDs_);
00221   }
00222   if ((int)matrixIDs_.size() != numMatrices) {
00223     FEI_CERR << "fei::FEI_Impl::setIDLists ERROR creating matrixIDs_ list."<<FEI_ENDL;
00224     ERReturn(-1);
00225   }
00226 
00227   rhsIDs_.resize(0);
00228   b_.resize(numRHSs);
00229   for(int i=0; i<numRHSs; ++i) {
00230     fei::sortedListInsert(rhsIDs[i], rhsIDs_);
00231   }
00232   if ((int)rhsIDs_.size() != numRHSs) {
00233     FEI_CERR << "fei::FEI_Impl::setIDLists ERROR creating rhsIDs_ list."<<FEI_ENDL;
00234     ERReturn(-1);
00235   }
00236 
00237   if (wrapper_[0].get() != NULL) {
00238     fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
00239     if (linSysCore.get() != NULL) {
00240       linSysCore->setNumRHSVectors(rhsIDs_.size(), &rhsIDs_[0]);
00241     }
00242   }
00243 
00244   return(0);
00245 }
00246 
00247 fei::SharedPtr<fei::LinearSystem> fei::FEI_Impl::getLinearSystem()
00248 {
00249   return( linSys_ );
00250 }
00251 
00252 int fei::FEI_Impl::parameters(int numParams, 
00253              const char *const* paramStrings)
00254 {
00255   // merge these parameters with any others we may have, for later use.
00256   snl_fei::mergeStringLists(paramStrings_, numParams_,
00257            paramStrings, numParams);
00258 
00259   if (wrapper_[0].get() != NULL) {
00260     if (wrapper_[0]->haveLinearSystemCore()) {
00261       CHK_ERR( wrapper_[0]->getLinearSystemCore()->parameters(numParams, (char**)paramStrings) );
00262     }
00263     if (wrapper_[0]->haveFiniteElementData()) {
00264       CHK_ERR( wrapper_[0]->getFiniteElementData()->parameters(numParams, (char**)paramStrings) );
00265     }
00266   }
00267 
00268   std::vector<std::string> stdstrings;
00269   fei::utils::char_ptrs_to_strings(numParams, paramStrings, stdstrings);
00270   fei::ParameterSet paramset;
00271   fei::utils::parse_strings(stdstrings, " ", paramset);
00272   factory_[0]->parameters(paramset);
00273 
00274   return(0);
00275 }
00276 
00277 int fei::FEI_Impl::setSolveType(int solveType)
00278 {
00279   solveType_ = solveType;
00280   setSolveTypeCalled_ = true;
00281 
00282   return(0);
00283 }
00284 
00285 int fei::FEI_Impl::initFields(int numFields, 
00286              const int *fieldSizes, 
00287              const int *fieldIDs)
00288 {
00289   rowSpace_->defineFields(numFields, fieldIDs, fieldSizes);
00290 
00291   return(0);
00292 }
00293 
00294 int fei::FEI_Impl::initElemBlock(GlobalID elemBlockID,
00295           int numElements,
00296           int numNodesPerElement,
00297           const int* numFieldsPerNode,
00298           const int* const* nodalFieldIDs,
00299           int numElemDofFieldsPerElement,
00300           const int* elemDOFFieldIDs,
00301           int interleaveStrategy)
00302 {
00303   //define pattern that describes the layout of fields for elements in
00304   //this element-block
00305 
00306   int numIDs = numNodesPerElement;
00307   if (numElemDofFieldsPerElement > 0) ++numIDs;
00308   std::vector<int> idTypes;
00309   std::vector<int> numFieldsPerID;
00310   std::vector<int> fieldIDs;
00311 
00312   int i, j;
00313   for(i=0; i<numNodesPerElement; ++i) {
00314     idTypes.push_back(nodeIDType_);
00315     numFieldsPerID.push_back(numFieldsPerNode[i]);
00316 
00317     for(j=0; j<numFieldsPerNode[i]; ++j) {
00318       fieldIDs.push_back(nodalFieldIDs[i][j]);
00319     }
00320   }
00321 
00322   if (numElemDofFieldsPerElement>0) {
00323     idTypes.push_back(elemIDType_);
00324     numFieldsPerID.push_back(numElemDofFieldsPerElement);
00325     for(i=0; i<numElemDofFieldsPerElement; ++i) {
00326       fieldIDs.push_back(elemDOFFieldIDs[i]);
00327     }
00328 
00329     block_dof_per_elem_.insert(std::pair<int,int>(elemBlockID, numElemDofFieldsPerElement));
00330     any_blocks_have_elem_dof_ = true;
00331   }
00332 
00333   int pattID = matGraph_->definePattern(numIDs,
00334          &idTypes[0],
00335          &numFieldsPerID[0],
00336          &fieldIDs[0]);
00337 
00338   //initialize connectivity-block
00339   CHK_ERR( matGraph_->initConnectivityBlock(elemBlockID, numElements,
00340               pattID) );
00341 
00342   return(0);
00343 }
00344 
00345 int fei::FEI_Impl::initElem(GlobalID elemBlockID,
00346            GlobalID elemID,
00347            const GlobalID* elemConn)
00348 {
00349   bool elemdof = false;
00350 
00351   if (any_blocks_have_elem_dof_) {
00352     std::map<int,int>::const_iterator
00353       b_iter = block_dof_per_elem_.find(elemBlockID);
00354     if (b_iter != block_dof_per_elem_.end()) {
00355       int numIDs = matGraph_->getNumIDsPerConnectivityList(elemBlockID);
00356       iwork_.resize(numIDs);
00357       int* iworkPtr = &iwork_[0];
00358       for(int i=0; i<numIDs-1; ++i) {
00359   iworkPtr[i] = elemConn[i];
00360       }
00361       iworkPtr[numIDs-1] = elemID;
00362 
00363       CHK_ERR( matGraph_->initConnectivity(elemBlockID, elemID, iworkPtr) );
00364       elemdof = true;
00365     }
00366   }
00367 
00368   if (!elemdof) {
00369     CHK_ERR( matGraph_->initConnectivity(elemBlockID, elemID, elemConn) );
00370   }
00371 
00372   nodeset_filled_ = false;
00373 
00374   return(0);
00375 }
00376 
00377 int fei::FEI_Impl::initSlaveVariable(GlobalID slaveNodeID, 
00378               int slaveFieldID,
00379               int offsetIntoSlaveField,
00380               int numMasterNodes,
00381               const GlobalID* masterNodeIDs,
00382               const int* masterFieldIDs,
00383               const double* weights,
00384               double rhsValue)
00385 {
00386   throw std::runtime_error("FEI_Impl::initSlaveVariable not implemented.");
00387   return(0);
00388 }
00389 
00390 int fei::FEI_Impl::deleteMultCRs()
00391 {
00392   throw std::runtime_error("FEI_Impl::deleteMultCRs not implemented.");
00393   return(0);
00394 }
00395 
00396 int fei::FEI_Impl::initSharedNodes(int numSharedNodes,
00397             const GlobalID *sharedNodeIDs,  
00398             const int* numProcsPerNode, 
00399             const int *const *sharingProcIDs)
00400 {
00401   CHK_ERR( rowSpace_->initSharedIDs(numSharedNodes,
00402             nodeIDType_,
00403             sharedNodeIDs,
00404             numProcsPerNode,
00405             sharingProcIDs) );
00406 
00407   return(0);
00408 }
00409 
00410 int fei::FEI_Impl::initCRMult(int numCRNodes,
00411              const GlobalID* CRNodeIDs,
00412              const int *CRFieldIDs,
00413              int& CRID)
00414 {
00415   iwork_.assign(numCRNodes,nodeIDType_);
00416 
00417   CRID = constraintID_++;
00418 
00419   CHK_ERR( matGraph_->initLagrangeConstraint(CRID,
00420                constraintIDType_,
00421                numCRNodes,
00422                &iwork_[0],
00423                CRNodeIDs,
00424                CRFieldIDs) );
00425 
00426   return(0);
00427 }
00428 
00429 int fei::FEI_Impl::initCRPen(int numCRNodes,
00430             const GlobalID* CRNodeIDs, 
00431             const int *CRFieldIDs,
00432             int& CRID)
00433 {
00434   iwork_.assign(numCRNodes, nodeIDType_);
00435 
00436   CRID = constraintID_++;
00437 
00438   CHK_ERR( matGraph_->initPenaltyConstraint(CRID,
00439               constraintIDType_,
00440               numCRNodes,
00441               &iwork_[0],
00442               CRNodeIDs,
00443               CRFieldIDs) );
00444 
00445   return(0);
00446 }
00447 
00448 int fei::FEI_Impl::initComplete()
00449 {
00450   CHK_ERR( matGraph_->initComplete() );
00451 
00452   if (matrixIDs_.size() < 1 || factory_.size() < 1 ||
00453       A_.size() < 1 || b_.size() < 1) {
00454     ERReturn(-1);
00455   }
00456 
00457   A_[0] = factory_[0]->createMatrix(matGraph_);
00458 
00459   std::vector<std::string> stdstrings;
00460   fei::utils::char_ptrs_to_strings(numParams_, paramStrings_, stdstrings);
00461   fei::ParameterSet params;
00462   fei::utils::parse_strings(stdstrings, " ", params);
00463 
00464   CHK_ERR( A_[0]->parameters(params) );
00465 
00466   b_[0] = factory_[0]->createVector(matGraph_);
00467 
00468   if (matrixIDs_.size() > 1) {
00469     bool multiple_factories = false;
00470 
00471     if (wrapper_[0].get() != NULL) {
00472       multiple_factories = true;
00473 
00474       fei::SharedPtr<LinearSystemCore> linsyscore = wrapper_[0]->getLinearSystemCore();
00475       if (linsyscore.get() == NULL) {
00476   FEI_CERR << "fei::FEI_Impl ERROR, multiple matrix/rhs assembly not supported "
00477        << "non-null LibraryWrapper holds null LinearSystemCore."<<FEI_ENDL;
00478   ERReturn(-1);
00479       }
00480 
00481       wrapper_.resize(matrixIDs_.size());
00482       factory_.resize(matrixIDs_.size());
00483       for(unsigned i=1; i<matrixIDs_.size(); ++i) {
00484   fei::SharedPtr<LinearSystemCore> lscclone(linsyscore->clone());
00485   wrapper_[i].reset(new LibraryWrapper(lscclone));
00486   factory_[i].reset(new snl_fei::Factory(comm_, wrapper_[i]));
00487       }
00488     }
00489 
00490     fei::SharedPtr<fei::Factory> factory;
00491     for(unsigned i=1; i<matrixIDs_.size(); ++i) {
00492       factory = multiple_factories ? factory_[i] : factory_[0];
00493 
00494       A_[i] = factory->createMatrix(matGraph_);
00495       CHK_ERR( A_[i]->parameters(params) );
00496     }
00497 
00498     for(unsigned i=1; i<rhsIDs_.size(); ++i) {
00499       factory = multiple_factories ? factory_[i] : factory_[0];
00500 
00501       b_[i] = factory->createVector(matGraph_);
00502     }
00503 
00504     if (wrapper_[0].get() != NULL) {
00505       fei::SharedPtr<LinearSystemCore> linsyscore
00506         = wrapper_[0]->getLinearSystemCore();
00507 
00508       linsyscore->setNumRHSVectors(1, &(rhsIDs_[0]));
00509 
00510       unsigned num = rhsIDs_.size();
00511       if (matrixIDs_.size() < num) num = matrixIDs_.size();
00512       for(unsigned i=1; i<num; ++i) {
00513   fei::SharedPtr<LinearSystemCore> lsc = wrapper_[i]->getLinearSystemCore();
00514 
00515   if (i==num-1 && rhsIDs_.size() > num) {
00516     int numRHSs = rhsIDs_.size() - matrixIDs_.size() + 1;
00517     lsc->setNumRHSVectors(numRHSs, &(rhsIDs_[i]));
00518   }
00519   else {
00520     lsc->setNumRHSVectors(1, &(rhsIDs_[i]));
00521   }
00522       }
00523 
00524       if (rhsIDs_.size() < matrixIDs_.size()) {
00525   int dummyID = -1;
00526   for(unsigned i=rhsIDs_.size(); i<matrixIDs_.size(); ++i) {
00527     wrapper_[i]->getLinearSystemCore()->setNumRHSVectors(1, &dummyID);
00528   }
00529       }
00530     }
00531 
00532     for(unsigned i=1; i<matrixIDs_.size(); ++i) {
00533       factory = multiple_factories ? factory_[i] : factory_[0];
00534 
00535       A_[i] = factory->createMatrix(matGraph_);
00536       CHK_ERR( A_[i]->parameters(params) );
00537     }
00538 
00539     for(unsigned i=1; i<rhsIDs_.size(); ++i) {
00540       factory = multiple_factories ? factory_[i] : factory_[0];
00541 
00542       b_[i] = factory->createVector(matGraph_);
00543     }
00544   }
00545 
00546   x_ = factory_[0]->createVector(matGraph_, true);
00547 
00548   linSys_ = factory_[0]->createLinearSystem(matGraph_);
00549 
00550   CHK_ERR( linSys_->parameters(numParams_, paramStrings_) );
00551 
00552   linSys_->setMatrix(A_[0]);
00553   linSys_->setRHS(b_[0]);
00554   linSys_->setSolutionVector(x_);
00555 
00556   return(0);
00557 }
00558 
00559 int fei::FEI_Impl::setCurrentMatrix(int matID)
00560 {
00561   std::vector<int>::const_iterator
00562     iter = std::lower_bound(matrixIDs_.begin(), matrixIDs_.end(), matID);
00563   if (iter == matrixIDs_.end() || *iter != matID) {
00564     FEI_CERR << "fei::FEI_Impl::setCurrentMatrix: matID ("<<matID
00565        <<") not found." <<FEI_ENDL;
00566     return(-1);
00567   }
00568 
00569   index_current_ = iter-matrixIDs_.begin();
00570 
00571   return(0);
00572 }
00573 
00574 int fei::FEI_Impl::setCurrentRHS(int rhsID)
00575 {
00576   std::vector<int>::const_iterator
00577     iter = std::lower_bound(rhsIDs_.begin(), rhsIDs_.end(), rhsID);
00578   if (iter == rhsIDs_.end() || *iter != rhsID) {
00579     FEI_CERR << "fei::FEI_Impl::setCurrentRHS: rhsID ("<<rhsID<<") not found."
00580    << FEI_ENDL;
00581     return(-1);
00582   }
00583 
00584   index_current_rhs_row_ = iter - rhsIDs_.begin();
00585 
00586   return(0);
00587 }
00588 
00589 int fei::FEI_Impl::resetSystem(double s)
00590 {
00591   int err = A_[index_current_]->putScalar(s);
00592   err += x_->putScalar(s);
00593   err += b_[index_current_rhs_row_]->putScalar(s);
00594 
00595   return(err);
00596 }
00597 
00598 int fei::FEI_Impl::resetMatrix(double s)
00599 {
00600   return( A_[index_current_]->putScalar(s) );
00601 }
00602 
00603 int fei::FEI_Impl::resetRHSVector(double s)
00604 {
00605   return( b_[index_current_rhs_row_]->putScalar(s) );
00606 }
00607 
00608 int fei::FEI_Impl::resetInitialGuess(double s)
00609 {
00610   return( x_->putScalar(s) );
00611 }
00612 
00613 int fei::FEI_Impl::loadNodeBCs(int numNodes,
00614                                 const GlobalID *nodeIDs,
00615                                 int fieldID,
00616                                 const int* offsetsIntoField,
00617                                 const double* prescribedValues)
00618 {
00619   CHK_ERR( linSys_->loadEssentialBCs(numNodes, nodeIDs,
00620                                      nodeIDType_, fieldID,
00621                                      offsetsIntoField, prescribedValues) );
00622 
00623   newData_ = true;
00624 
00625   return(0);
00626 }
00627 
00628 int fei::FEI_Impl::loadElemBCs(int numElems,
00629                     const GlobalID* elemIDs,  
00630                     int fieldID,
00631                     const double *const *alpha,  
00632                     const double *const *beta,  
00633                     const double *const *gamma)
00634 {
00635   throw std::runtime_error("FEI_Impl::loadElemBCs not implemented.");
00636   return(0);
00637 }
00638 
00639 int fei::FEI_Impl::sumInElem(GlobalID elemBlockID,
00640                  GlobalID elemID,
00641                  const GlobalID* elemConn,
00642                  const double* const* elemStiffness,
00643                  const double* elemLoad,
00644                  int elemFormat)
00645 {
00646   CHK_ERR( A_[index_current_]->sumIn(elemBlockID, elemID, elemStiffness, elemFormat) );
00647 
00648   int num = matGraph_->getConnectivityNumIndices(elemBlockID);
00649   std::vector<int> indices(num);
00650   CHK_ERR( matGraph_->getConnectivityIndices(elemBlockID, elemID, num,
00651                                              &indices[0], num) );
00652   CHK_ERR( b_[index_current_rhs_row_]->sumIn(num, &indices[0], elemLoad, 0) );
00653 
00654   newData_ = true;
00655 
00656   return(0);
00657 }
00658 
00659 int fei::FEI_Impl::sumInElemMatrix(GlobalID elemBlockID,
00660             GlobalID elemID,
00661             const GlobalID* elemConn,
00662             const double* const* elemStiffness,
00663             int elemFormat)
00664 {
00665   CHK_ERR( A_[index_current_]->sumIn(elemBlockID, elemID, elemStiffness, elemFormat) );
00666 
00667   newData_ = true;
00668 
00669   return(0);
00670 }
00671 
00672 int fei::FEI_Impl::sumInElemRHS(GlobalID elemBlockID,
00673          GlobalID elemID,
00674          const GlobalID* elemConn,
00675          const double* elemLoad)
00676 {
00677   int num = matGraph_->getConnectivityNumIndices(elemBlockID);
00678   std::vector<int> indices(num);
00679   CHK_ERR( matGraph_->getConnectivityIndices(elemBlockID, elemID, num,
00680                                              &indices[0], num) );
00681   CHK_ERR( b_[index_current_rhs_row_]->sumIn(num, &indices[0], elemLoad, 0) );
00682 
00683   newData_ = true;
00684 
00685   return(0);
00686 }
00687 
00688 int fei::FEI_Impl::loadCRMult(int CRMultID,
00689              int numCRNodes,
00690              const GlobalID* CRNodeIDs,
00691              const int* CRFieldIDs,
00692              const double* CRWeights,
00693              double CRValue)
00694 {
00695   newData_ = true;
00696 
00697   CHK_ERR( linSys_->loadLagrangeConstraint(CRMultID, CRWeights, CRValue) );
00698 
00699   return(0);
00700 }
00701 
00702 int fei::FEI_Impl::loadCRPen(int CRPenID,
00703                  int numCRNodes,
00704                  const GlobalID* CRNodeIDs,
00705                  const int* CRFieldIDs,
00706                  const double* CRWeights,
00707                  double CRValue,
00708                  double penValue)
00709 {
00710   newData_ = true;
00711 
00712   CHK_ERR( linSys_->loadPenaltyConstraint(CRPenID, CRWeights, penValue, CRValue) );
00713 
00714   return(0);
00715 }
00716 
00717 int fei::FEI_Impl::putIntoRHS(int IDType,
00718              int fieldID,
00719              int numIDs,
00720              const GlobalID* IDs,
00721              const double* coefficients)
00722 {
00723   CHK_ERR( inputRHS(IDType, fieldID, numIDs, IDs, coefficients, false) );
00724 
00725   newData_ = true;
00726 
00727   return(0);
00728 }
00729 
00730 int fei::FEI_Impl::sumIntoRHS(int IDType,
00731              int fieldID,
00732              int numIDs,
00733              const GlobalID* IDs,
00734              const double* coefficients)
00735 {
00736   CHK_ERR( inputRHS(IDType, fieldID, numIDs, IDs, coefficients, true) );
00737 
00738   newData_ = true;
00739 
00740   return(0);
00741 }
00742 
00743 int fei::FEI_Impl::inputRHS(int IDType,
00744            int fieldID,
00745            int numIDs,
00746            const GlobalID* IDs,
00747            const double* coefficients,
00748            bool sumInto)
00749 {
00750   int fieldSize = rowSpace_->getFieldSize(fieldID);
00751 
00752   int offset = 0;
00753   for(int i=0; i<numIDs; ++i) {
00754     int globalIndex = 0;
00755     CHK_ERR( rowSpace_->getGlobalIndex(IDType, IDs[i], globalIndex) );
00756 
00757     for(int j=0; j<fieldSize; ++j) {
00758       int eqn = globalIndex+j;
00759       if (sumInto) {
00760   CHK_ERR( b_[index_current_rhs_row_]->sumIn(1, &eqn, &(coefficients[offset++])) );
00761       }
00762       else {
00763   CHK_ERR( b_[index_current_rhs_row_]->copyIn(1, &eqn, &(coefficients[offset++])) );
00764       }
00765     }
00766   }
00767 
00768   return(0);
00769 }
00770 
00771 int fei::FEI_Impl::setMatScalars(int numScalars,
00772           const int* IDs, 
00773           const double* scalars)
00774 {
00775   matScalars_.resize(matrixIDs_.size());
00776 
00777   for(int i=0; i<numScalars; ++i) {
00778     std::vector<int>::const_iterator
00779       iter = std::lower_bound(matrixIDs_.begin(), matrixIDs_.end(), IDs[i]);
00780     if (iter == matrixIDs_.end() || *iter != IDs[i]) {
00781       continue;
00782     }
00783 
00784     unsigned index = iter - matrixIDs_.begin();
00785     matScalars_[index] = scalars[i];
00786   }
00787 
00788   matScalarsSet_ = true;
00789 
00790   return(0);
00791 }
00792 
00793 int fei::FEI_Impl::setRHSScalars(int numScalars,
00794           const int* IDs,
00795           const double* scalars)
00796 {
00797   rhsScalars_.resize(rhsIDs_.size());
00798 
00799   for(int i=0; i<numScalars; ++i) {
00800     std::vector<int>::const_iterator
00801       iter = std::lower_bound(rhsIDs_.begin(), rhsIDs_.end(), IDs[i]);
00802     if (iter == rhsIDs_.end() || *iter != IDs[i]) {
00803       continue;
00804     }
00805 
00806     unsigned index = iter - rhsIDs_.begin();
00807     rhsScalars_[index] = scalars[i];
00808   }
00809 
00810   rhsScalarsSet_ = true;
00811 
00812   return(0);
00813 }
00814 
00815 int fei::FEI_Impl::loadComplete(bool applyBCs,
00816                                 bool globalAssemble)
00817 {
00818   if (!newData_) {
00819     return(0);
00820   }
00821 
00822   if (linSys_.get() == NULL) {
00823     FEI_CERR << "fei::FEI_Impl::loadComplete: loadComplete can not be called"
00824    << " until after initComplete has been called."<<FEI_ENDL;
00825     return(-1);
00826   }
00827 
00828   if (solveType_ == FEI_AGGREGATE_SUM) {
00829     for(unsigned i=0; i<A_.size(); ++i) {
00830   CHK_ERR( A_[i]->gatherFromOverlap() );
00831     }
00832 
00833     for(unsigned j=0; j<b_.size(); ++j) {
00834       CHK_ERR( b_[j]->gatherFromOverlap() );
00835     }
00836   
00837     CHK_ERR( aggregateSystem() );
00838   }
00839 
00840   CHK_ERR( linSys_->loadComplete(applyBCs, globalAssemble) );
00841 
00842   if (b_.size() > 1) {
00843     int rhs_counter = 0;
00844     for(unsigned i=0; i<b_.size(); ++i) {
00845       FEI_OSTRINGSTREAM osstr;
00846       osstr << "rhs_" << rhs_counter++;
00847       CHK_ERR( linSys_->putAttribute(osstr.str().c_str(), b_[i].get()) );
00848     }
00849   }
00850 
00851   newData_ = false;
00852 
00853   return(0);
00854 }
00855 
00856 int fei::FEI_Impl::aggregateSystem()
00857 {
00858   if (wrapper_[0].get() == NULL) {
00859     ERReturn(-1);
00860   }
00861 
00862   if (wrapper_[0].get() != NULL) {
00863     CHK_ERR( aggregateSystem_LinSysCore() );
00864   }
00865 
00866   return(0);
00867 }
00868 
00869 int fei::FEI_Impl::aggregateSystem_LinSysCore()
00870 {
00871   fei::SharedPtr<LinearSystemCore> lsc = wrapper_[0]->getLinearSystemCore();
00872   if (lsc.get() == NULL) return(-1);
00873 
00874   if (soln_fei_matrix_ == NULL) {
00875     soln_fei_matrix_ = new Data();
00876 
00877     CHK_ERR( lsc->copyOutMatrix(1.0, *soln_fei_matrix_) );
00878   }
00879 
00880   if (soln_fei_vector_ == NULL) {
00881     soln_fei_vector_ = new Data();
00882 
00883     CHK_ERR( lsc->setRHSID(rhsIDs_[0]) );
00884     CHK_ERR( lsc->copyOutRHSVector(1.0, *soln_fei_vector_) );
00885   }
00886 
00887   Data tmp, tmpv;
00888   for(unsigned i=0; i<matrixIDs_.size(); ++i) {
00889     fei::SharedPtr<LinearSystemCore> lsci = wrapper_[i]->getLinearSystemCore();
00890     if (lsci.get() == NULL) return(-1);
00891 
00892     if (i==0) {
00893       CHK_ERR( lsci->copyInMatrix(matScalars_[i], *soln_fei_matrix_) );
00894     }
00895     else {
00896       CHK_ERR( lsci->getMatrixPtr(tmp) );
00897       CHK_ERR( lsc->sumInMatrix(matScalars_[i], tmp) );
00898     }
00899   }
00900 
00901   int last_mat = matrixIDs_.size() - 1;
00902   for(unsigned i=0; i<rhsIDs_.size(); ++i) {
00903     fei::SharedPtr<LinearSystemCore> lsci = (int)i<last_mat ?
00904       wrapper_[i]->getLinearSystemCore() : wrapper_[last_mat]->getLinearSystemCore();
00905 
00906     if (i==0) {
00907       CHK_ERR( lsci->copyInRHSVector(rhsScalars_[i], *soln_fei_vector_) );
00908     }
00909     else {
00910       CHK_ERR( lsci->setRHSID(rhsIDs_[i]) );
00911       CHK_ERR( lsci->getRHSVectorPtr(tmpv) );
00912       CHK_ERR( lsc->sumInRHSVector(rhsScalars_[i], tmpv) );
00913     }
00914   }
00915 
00916   return(0);
00917 }
00918 
00919 int fei::FEI_Impl::residualNorm(int whichNorm,
00920          int numFields,
00921          int* fieldIDs,
00922          double* norms)
00923 {
00924   CHK_ERR( loadComplete() );
00925 
00926   fei::SharedPtr<fei::VectorSpace> rowspace =
00927     matGraph_->getRowSpace();
00928   int numLocalEqns = rowspace->getNumIndices_Owned();
00929 
00930   std::vector<double> residValues(numLocalEqns);
00931   double* residValuesPtr = &residValues[0];
00932 
00933   std::vector<int> globalEqnOffsets;
00934   rowspace->getGlobalIndexOffsets(globalEqnOffsets);
00935   int firstLocalOffset = globalEqnOffsets[localProc_];
00936 
00937   if (wrapper_[0].get() == NULL) {
00938     fei::SharedPtr<fei::Vector> r = factory_[0]->createVector(rowspace);
00939 
00940     //form r = A*x
00941     CHK_ERR( A_[index_current_]->multiply(x_.get(), r.get()) );
00942 
00943     //form r = b - A*x
00944     CHK_ERR( r->update(1.0, b_[index_current_rhs_row_].get(), -1.0) );
00945 
00946     //now put the values from r into the residValues array.
00947     for(int ii=0; ii<numLocalEqns; ++ii) {
00948       int index = firstLocalOffset+ii;
00949       CHK_ERR( r->copyOut(1, &index, &(residValuesPtr[ii]) ) );
00950     }
00951   }
00952   else {
00953     fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
00954     if (linSysCore.get() != NULL) {
00955       CHK_ERR( linSysCore->formResidual(residValuesPtr, numLocalEqns) );
00956     }
00957     else {
00958       FEI_CERR << "FEI_Impl::residualNorm: warning: residualNorm not implemented"
00959      << " for FiniteElementData."<<FEI_ENDL;
00960       int offset = 0;
00961       for(int ii=0; ii<numFields; ++ii) {
00962   int fieldSize = rowspace->getFieldSize(fieldIDs[ii]);
00963   for(int jj=0; jj<fieldSize; ++jj) {
00964     norms[offset++] = -99.9;
00965   }
00966       }
00967       return(0);
00968     }
00969   }
00970 
00971   int numLocalNodes = rowspace->getNumOwnedIDs(nodeIDType_);
00972   std::vector<int> nodeIDs(numLocalNodes);
00973   int* nodeIDsPtr = &nodeIDs[0];
00974   int check;
00975   CHK_ERR( rowspace->getOwnedIDs(nodeIDType_, numLocalNodes,
00976           nodeIDsPtr, check) );
00977   std::vector<int> indices(numLocalEqns);
00978   int* indicesPtr = &indices[0];
00979 
00980   std::vector<double> tmpNorms(numFields, 0.0);
00981 
00982   for(int i=0; i<numFields; ++i) {
00983     int fieldSize = rowspace->getFieldSize(fieldIDs[i]);
00984 
00985     CHK_ERR( rowspace->getGlobalIndices(numLocalNodes, nodeIDsPtr,
00986           nodeIDType_, fieldIDs[i],
00987           indicesPtr) );
00988 
00989     double tmp = 0.0;
00990     for(int j=0; j<fieldSize*numLocalNodes; ++j) {
00991       if (indicesPtr[j] < 0) {
00992   continue;
00993       }
00994 
00995       double val = std::abs(residValuesPtr[indicesPtr[j]-firstLocalOffset]);
00996       switch(whichNorm) {
00997       case 0:
00998   if (val > tmp) tmp = val;
00999   break;
01000       case 1:
01001   tmp += val;
01002   break;
01003       case 2:
01004   tmp += val*val;
01005   break;
01006       default:
01007   FEI_CERR << "FEI_Impl::residualNorm: whichNorm="<<whichNorm<<" not recognized"
01008        << FEI_ENDL;
01009   return(-1);
01010       }
01011     }
01012     tmpNorms[i] = tmp;
01013   }
01014 
01015   std::vector<double> normsArray(numFields);
01016   for(int i=0; i<numFields; ++i) {
01017     normsArray[i] = norms[i];
01018   }
01019 
01020   switch(whichNorm) {
01021   case 0:
01022     CHK_ERR( fei::GlobalMax(comm_, tmpNorms, normsArray) );
01023     break;
01024   default:
01025     CHK_ERR( fei::GlobalSum(comm_, tmpNorms, normsArray) );
01026   }
01027 
01028   for(int i=0; i<numFields; ++i) {
01029     norms[i] = normsArray[i];
01030   }
01031 
01032   if (whichNorm == 2) {
01033     for(int ii=0; ii<numFields; ++ii) {
01034       norms[ii] = std::sqrt(norms[ii]);
01035     }
01036   }
01037 
01038   return(0);
01039 }
01040 
01041 int fei::FEI_Impl::solve(int& status)
01042 {
01043   CHK_ERR( loadComplete() );
01044 
01045   fei::SharedPtr<fei::Solver> solver = factory_[0]->createSolver();
01046 
01047   std::vector<std::string> stdstrings;
01048   fei::utils::char_ptrs_to_strings(numParams_, paramStrings_, stdstrings);
01049   fei::ParameterSet params;
01050   fei::utils::parse_strings(stdstrings, " ", params);
01051 
01052   int err = solver->solve(linSys_.get(),
01053         NULL,//preconditioningMatrix
01054         params, iterations_, status);
01055 
01056   CHK_ERR( x_->scatterToOverlap() );
01057 
01058   return(err);
01059 }
01060 
01061 int fei::FEI_Impl::iterations(int& itersTaken) const
01062 {
01063   itersTaken = iterations_;
01064 
01065   return(0);
01066 }
01067 
01068 int fei::FEI_Impl::version(const char*& versionString)
01069 {
01070   versionString = fei::utils::version();
01071 
01072   return(0);
01073 }
01074 
01075 int fei::FEI_Impl::cumulative_cpu_times(double& initTime,
01076            double& loadTime,
01077            double& solveTime,
01078            double& solnReturnTime)
01079 {
01080   initTime = initTime_;
01081   loadTime = loadTime_;
01082   solveTime = solveTime_;
01083   solnReturnTime = solnReturnTime_;
01084 
01085   return(0);
01086 }
01087 
01088 int fei::FEI_Impl::getBlockNodeSolution(GlobalID elemBlockID,  
01089            int numNodes, 
01090            const GlobalID *nodeIDs, 
01091            int *offsets,
01092            double *results)
01093 {
01094   (void)elemBlockID;
01095   return ( getNodalSolution(numNodes, nodeIDs, offsets, results) );
01096 }
01097 
01098 int fei::FEI_Impl::getNodalSolution(int numNodes,
01099              const GlobalID* nodeIDs,
01100              int* offsets,
01101              double* results)
01102 {
01103   int j;
01104   int offset = 0;
01105   std::vector<int> fieldIDs;
01106   for(int i=0; i<numNodes; ++i) {
01107     offsets[i] = offset;
01108 
01109     GlobalID nodeID = nodeIDs[i];
01110 
01111     int numFields = rowSpace_->getNumFields(nodeIDType_, nodeID);
01112     iwork_.resize( numFields*2 );
01113     int* fieldSizes = &iwork_[0];
01114     rowSpace_->getFields(nodeIDType_, nodeID, fieldIDs);
01115 
01116     int numDOF = 0;
01117     for(j=0; j<numFields; ++j) {
01118       fieldSizes[j] = rowSpace_->getFieldSize(fieldIDs[j]);
01119       numDOF += fieldSizes[j];
01120     }
01121 
01122     if (!rowSpace_->isLocal(nodeIDType_, nodeID)) {
01123       offset += numDOF;
01124       continue;
01125     }
01126 
01127     int results_offset = offset;
01128     for(j=0; j<numFields; ++j) {
01129       CHK_ERR( x_->copyOutFieldData(fieldIDs[j], nodeIDType_,
01130             1, &nodeID,
01131             &(results[results_offset])));
01132 
01133       results_offset += fieldSizes[j];
01134     }
01135 
01136     offset += numDOF;
01137   }
01138 
01139   offsets[numNodes] = offset;
01140 
01141   return(0);
01142 }
01143 
01144 int fei::FEI_Impl::getBlockFieldNodeSolution(GlobalID elemBlockID,
01145                                   int fieldID,
01146                                   int numNodes, 
01147                                   const GlobalID *nodeIDs, 
01148                                   double *results)
01149 {
01150   throw std::runtime_error("FEI_Impl::getBlockFieldNodeSolution not implemented.");
01151   return(0);
01152 }
01153 
01154 int fei::FEI_Impl::getBlockElemSolution(GlobalID elemBlockID,  
01155            int numElems, 
01156            const GlobalID *elemIDs,
01157            int& numElemDOFPerElement,
01158            double *results)
01159 {
01160   std::map<int,int>::const_iterator b_iter = block_dof_per_elem_.find(elemBlockID);
01161   if (b_iter == block_dof_per_elem_.end()) return(-1);
01162   numElemDOFPerElement = (*b_iter).second;
01163 
01164   fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(elemBlockID);
01165   if (block==NULL) {
01166     FEI_OSTRINGSTREAM osstr;
01167     osstr<< "fei::FEI_Impl::getBlockElemSolution ERROR, elemBlockID "
01168    << elemBlockID << " not valid.";
01169     throw std::runtime_error(osstr.str());
01170   }
01171 
01172   fei::Pattern* pattern = block->getRowPattern();
01173 
01174   int numIDs = pattern->getNumIDs();
01175   const int* idTypes = pattern->getIDTypes();
01176   const int* fIDs= pattern->getFieldIDs();
01177   const int* numFieldsPerID = pattern->getNumFieldsPerID();
01178 
01179   std::vector<int> fieldIDs;
01180   int foffset = 0;
01181   int i;
01182   for(i=0; i<numIDs; ++i) {
01183     if (idTypes[i] == elemIDType_) {
01184       for(int j=0; j<numFieldsPerID[i]; ++j) {
01185   fieldIDs.push_back(fIDs[foffset++]);
01186       }
01187       break;
01188     }
01189     foffset += numFieldsPerID[i];
01190   }
01191 
01192   if (fieldIDs.size() < 1) {
01193     ERReturn(-1);
01194   }
01195 
01196   int offset = 0;
01197   for(i=0; i<numElems; ++i) {
01198     foffset = offset;
01199     for(size_t j=0; j<fieldIDs.size(); ++j) {
01200       int fieldSize;
01201       getFieldSize(fieldIDs[j], fieldSize);
01202 
01203       CHK_ERR( x_->copyOutFieldData(fieldIDs[j], elemIDType_, 1, &(elemIDs[i]),
01204             &(results[foffset])) );
01205       foffset += fieldSize;
01206     }
01207 
01208     offset += numElemDOFPerElement;
01209   }
01210 
01211   return(0);
01212 }
01213 
01214 int fei::FEI_Impl::getNumCRMultipliers(int& numMultCRs)
01215 {
01216   numMultCRs = rowSpace_->getNumOwnedAndSharedIDs(constraintIDType_);
01217   return(0);
01218 }
01219 
01220 int fei::FEI_Impl::getCRMultIDList(int numMultCRs, int* multIDs)
01221 {
01222   int checkNum;
01223   CHK_ERR( rowSpace_->getOwnedAndSharedIDs(constraintIDType_, numMultCRs,
01224           multIDs, checkNum) );
01225   return(0);
01226 }
01227 
01228 int fei::FEI_Impl::getCRMultipliers(int numCRs,
01229              const int* CRIDs,
01230              double *multipliers)
01231 {
01232   iwork_.resize(numCRs);
01233 
01234   for(int i=0; i<numCRs; ++i) {
01235     CHK_ERR( rowSpace_->getGlobalIndex(constraintIDType_, CRIDs[i], iwork_[i]));
01236   }
01237 
01238   CHK_ERR( x_->copyOut(numCRs, &iwork_[0], multipliers) );
01239 
01240   return(0);
01241 }
01242 
01243 int fei::FEI_Impl::putBlockNodeSolution(GlobalID elemBlockID, 
01244                              int numNodes, 
01245                              const GlobalID *nodeIDs, 
01246                              const int *offsets,
01247                              const double *estimates)
01248 {
01249   throw std::runtime_error("FEI_Impl::putBlockNodeSolution not implemented.");
01250   return(0);
01251 }
01252 
01253 int fei::FEI_Impl::putBlockFieldNodeSolution(GlobalID elemBlockID, 
01254                                   int fieldID, 
01255                                   int numNodes, 
01256                                   const GlobalID *nodeIDs, 
01257                                   const double *estimates)
01258 {
01259   throw std::runtime_error("FEI_Impl::putBlockFieldNodeSolution not implemented.");
01260   return(0);
01261 }
01262 
01263 int fei::FEI_Impl::putBlockElemSolution(GlobalID elemBlockID,  
01264                              int numElems, 
01265                              const GlobalID *elemIDs, 
01266                              int dofPerElem,
01267                              const double *estimates)
01268 {
01269   throw std::runtime_error("FEI_Impl::putBlockElemSolution not implemented.");
01270   return(0);
01271 }
01272 
01273 int fei::FEI_Impl::putCRMultipliers(int numMultCRs, 
01274                          const int* CRIDs,
01275                          const double* multEstimates)
01276 {
01277   throw std::runtime_error("FEI_Impl::putCRMultipliers not implemented.");
01278   return(0);
01279 }
01280 
01281 int fei::FEI_Impl::getBlockNodeIDList(GlobalID elemBlockID,
01282                            int numNodes,
01283                            GlobalID *nodeIDs)
01284 {
01285   if (!nodeset_filled_ || elemBlockID != nodeset_blockid_) {
01286     CHK_ERR( fillNodeset(elemBlockID) );
01287   }
01288 
01289   fei::copySetToArray(nodeset_, numNodes, nodeIDs);
01290   return(0);
01291 }
01292 
01293 int fei::FEI_Impl::fillNodeset(int blockID) const
01294 {
01295   if (nodeset_filled_ && blockID == nodeset_blockid_) {
01296     return(0);
01297   }
01298 
01299   fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(blockID);
01300   if (block==NULL) {
01301     FEI_OSTRINGSTREAM osstr;
01302     osstr<< "fei::FEI_Impl::fillNodeset ERROR, blockID "
01303    << blockID << " not valid.";
01304     throw std::runtime_error(osstr.str());
01305   }
01306 
01307   std::vector<fei::Record<int>*>& nodes = block->getRowConnectivities();
01308 
01309   nodeset_.clear();
01310 
01311   fei::Record<int>** nodesPtr = &nodes[0];
01312   for(unsigned i=0; i<nodes.size(); ++i) {
01313     nodeset_.insert(nodesPtr[i]->getID());
01314   }
01315 
01316   nodeset_filled_ = true;
01317   nodeset_blockid_ = blockID;
01318 
01319   return(0);
01320 }
01321 
01322 int fei::FEI_Impl::getBlockElemIDList(GlobalID elemBlockID, 
01323                           int numElems, 
01324                           GlobalID* elemIDs)
01325 {
01326   fei::ConnectivityBlock* block = matGraph_->getConnectivityBlock(elemBlockID);
01327   if (block==NULL) {
01328     FEI_OSTRINGSTREAM osstr;
01329     osstr<< "fei::FEI_Impl::getBlockElemIDList ERROR, elemBlockID "
01330    << elemBlockID << " not valid.";
01331     throw std::runtime_error(osstr.str());
01332   }
01333 
01334   std::map<int,int>& elemIDSet = block->getConnectivityIDs();
01335 
01336   fei::copyKeysToArray(elemIDSet, numElems, elemIDs);
01337 
01338   return(0);
01339 }
01340 
01341 int fei::FEI_Impl::getNumSolnParams(GlobalID nodeID, int& numSolnParams) const
01342 {
01343   numSolnParams = rowSpace_->getNumDegreesOfFreedom( nodeIDType_, nodeID);
01344   return(0);
01345 }
01346 
01347 int fei::FEI_Impl::getNumElemBlocks(int& numElemBlocks) const
01348 {
01349   numElemBlocks = matGraph_->getNumConnectivityBlocks();
01350   return(0);
01351 }
01352 
01353 int fei::FEI_Impl::getNumBlockActNodes(GlobalID blockID, int& numNodes) const
01354 {
01355   if (!nodeset_filled_ || blockID != nodeset_blockid_) {
01356     CHK_ERR( fillNodeset(blockID) );
01357   }
01358 
01359   numNodes = nodeset_.size();
01360   return(0);
01361 }
01362 
01363 int fei::FEI_Impl::getNumBlockActEqns(GlobalID blockID, int& numEqns) const
01364 {
01365   throw std::runtime_error("fei::FEI_Impl::getNumBlockActEqns not implemented.");
01366   return(0);
01367 }
01368 
01369 int fei::FEI_Impl::getNumNodesPerElement(GlobalID blockID, int& nodesPerElem) const
01370 {
01371   nodesPerElem = matGraph_->getNumIDsPerConnectivityList(blockID);
01372   return(0);
01373 }
01374 
01375 int fei::FEI_Impl::getNumEqnsPerElement(GlobalID blockID, int& numEqns) const
01376 {
01377   numEqns = matGraph_->getConnectivityNumIndices(blockID);
01378   return(0);
01379 }
01380 
01381 int fei::FEI_Impl::getNumBlockElements(GlobalID blockID, int& numElems) const
01382 {
01383   const fei::ConnectivityBlock* cblock =
01384     matGraph_->getConnectivityBlock(blockID);
01385   if (cblock != NULL) {
01386     numElems = cblock->getConnectivityIDs().size();
01387   }
01388   else {
01389     numElems = 0;
01390   }
01391   return(0);
01392 }
01393 
01394 int fei::FEI_Impl::getNumBlockElemDOF(GlobalID blockID, int& DOFPerElem) const
01395 {
01396   std::map<int,int>::const_iterator b_iter = block_dof_per_elem_.find(blockID);
01397   if (b_iter == block_dof_per_elem_.end()) DOFPerElem = 0;
01398   else DOFPerElem = (*b_iter).second;
01399 
01400   return(0);
01401 }
01402 
01403 int fei::FEI_Impl::getParameters(int& numParams, char**& paramStrings)
01404 {
01405   numParams = numParams_;
01406   paramStrings = paramStrings_;
01407   return(0);
01408 }
01409 
01410 int fei::FEI_Impl::getFieldSize(int fieldID, int& numScalars)
01411 {
01412   numScalars = rowSpace_->getFieldSize(fieldID);
01413   return(0);
01414 }
01415 
01416 int fei::FEI_Impl::getEqnNumbers(GlobalID ID,
01417           int idType, 
01418           int fieldID,
01419           int& numEqns,
01420           int* eqnNumbers)
01421 {
01422   numEqns = rowSpace_->getFieldSize(fieldID);
01423   CHK_ERR( rowSpace_->getGlobalIndices(1, &ID, idType, fieldID, eqnNumbers) );
01424   return(0);
01425 }
01426 
01427 int fei::FEI_Impl::getNodalFieldSolution(int fieldID,
01428             int numNodes,
01429             const GlobalID* nodeIDs,
01430             double* results)
01431 {
01432   CHK_ERR( x_->copyOutFieldData(fieldID, nodeIDType_, numNodes,
01433         nodeIDs, results) );
01434   return(0);
01435 }
01436 
01437 int fei::FEI_Impl::getNumLocalNodes(int& numNodes)
01438 {
01439   numNodes = rowSpace_->getNumOwnedAndSharedIDs(nodeIDType_);
01440   return(0);
01441 }
01442 
01443 int fei::FEI_Impl::getLocalNodeIDList(int& numNodes,
01444                GlobalID* nodeIDs,
01445                int lenNodeIDs)
01446 {
01447   CHK_ERR( rowSpace_->getOwnedAndSharedIDs(nodeIDType_, lenNodeIDs, nodeIDs, numNodes) );
01448   return(0);
01449 }
01450 
01451 int fei::FEI_Impl::putNodalFieldData(int fieldID,
01452               int numNodes,
01453               const GlobalID* nodeIDs,
01454               const double* nodeData)
01455 {
01456   int err = 0;
01457   if (fieldID < 0) {
01458     bool data_passed = false;
01459     if (wrapper_[0].get() != NULL) {
01460       std::vector<int> numbers(numNodes);
01461       for(int i=0; i<numNodes; ++i) {
01462   err = rowSpace_->getGlobalBlkIndex(nodeIDType_, nodeIDs[i], numbers[i]);
01463   if (err != 0) {
01464     FEI_CERR << "fei::FEI_Impl::putNodalFieldData ERROR, nodeID "
01465          << nodeIDs[i] << " not found."<<FEI_ENDL;
01466     ERReturn(-1);
01467   }
01468       }
01469 
01470       int fieldSize = 0;
01471       try {
01472   fieldSize = rowSpace_->getFieldSize(fieldID);
01473       }
01474       catch (std::runtime_error& exc) {
01475   FEI_CERR << "fei::FEI_Impl::putNodalFieldData ERROR: " <<exc.what()<<FEI_ENDL;
01476   ERReturn(-1);
01477       }
01478 
01479       fei::SharedPtr<LinearSystemCore> linSysCore = wrapper_[0]->getLinearSystemCore();
01480       if (linSysCore.get() != NULL) {
01481   linSysCore->putNodalFieldData(fieldID, fieldSize, 
01482               &numbers[0], numNodes, nodeData);
01483   data_passed = true;
01484       }
01485       else {
01486   //If we enter this block, we're probably dealing with a FiniteElementData
01487   //instance.
01488   fei::SharedPtr<FiniteElementData> fedata = wrapper_[0]->getFiniteElementData();
01489   if (fedata.get() != NULL) {
01490     fedata->putNodalFieldData(fieldID, fieldSize, numNodes,
01491             &numbers[0], nodeData);
01492     data_passed = true;
01493   }
01494       }
01495     }
01496 
01497     if (!data_passed) {
01498       //If we get to here and data_passed is false, wrapper_[0] is probably NULL.
01499       if (wrapper_[0].get() == NULL) {
01500   fei::SharedPtr<fei::Vector> dataVector =factory_[0]->createVector(matGraph_);
01501 
01502   CHK_ERR( dataVector->copyInFieldData(fieldID, nodeIDType_,
01503                numNodes, nodeIDs, nodeData) );
01504   if (fieldID == -3) {
01505     CHK_ERR( linSys_->putAttribute("coordinates", dataVector.get()) );
01506   }
01507   else {
01508     FEI_OSTRINGSTREAM osstr;
01509     osstr << "fieldID:" << fieldID;
01510     CHK_ERR( linSys_->putAttribute(osstr.str().c_str(), dataVector.get()) );
01511   }
01512       }
01513       else {
01514   FEI_CERR << "fei::FEI_Impl::putNodalFieldData ERROR, non-null LibraryWrapper"
01515        << " contains neither LinearSystemCore or FiniteElementData. " <<FEI_ENDL;
01516   ERReturn(-1);
01517       }
01518     }
01519     return(0);
01520   }
01521 
01522   CHK_ERR( x_->copyInFieldData(fieldID, nodeIDType_,
01523              numNodes, nodeIDs, nodeData));
01524   return(0);
01525 }

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