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