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