fei_FEDataFilter.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 <limits>
00012 #include <cmath>
00013 
00014 #include <fei_defs.h>
00015 
00016 #include <fei_CommUtils.hpp>
00017 #include <fei_TemplateUtils.hpp>
00018 #include <snl_fei_Constraint.hpp>
00019 typedef snl_fei::Constraint<GlobalID> ConstraintType;
00020 
00021 #include <fei_LibraryWrapper.hpp>
00022 #include <SNL_FEI_Structure.hpp>
00023 #include <fei_FiniteElementData.hpp>
00024 #include <fei_Lookup.hpp>
00025 #include <FEI_Implementation.hpp>
00026 #include <fei_EqnCommMgr.hpp>
00027 #include <fei_EqnBuffer.hpp>
00028 #include <fei_NodeDatabase.hpp>
00029 #include <fei_NodeCommMgr.hpp>
00030 #include <fei_ProcEqns.hpp>
00031 #include <fei_BlockDescriptor.hpp>
00032 #include <fei_ConnectivityTable.hpp>
00033 #include <snl_fei_Utils.hpp>
00034 
00035 #include <fei_FEDataFilter.hpp>
00036 
00037 #undef fei_file
00038 #define fei_file "FEDataFilter.cpp"
00039 #include <fei_ErrMacros.hpp>
00040 
00041 #define ASSEMBLE_PUT 0
00042 #define ASSEMBLE_SUM 1
00043 
00044 //------------------------------------------------------------------------------
00045 void convert_eqns_to_nodenumbers_and_offsets(const NodeDatabase& nodeDB,
00046                                              int numEqns,
00047                                              const int* eqns,
00048                                              std::vector<int>& nodeNumbers,
00049                                              std::vector<int>& dofOffsets)
00050 {
00051   nodeNumbers.resize(numEqns);
00052   dofOffsets.resize(numEqns);
00053 
00054   for(int i=0; i<numEqns; ++i) {
00055     const NodeDescriptor* nodePtr = NULL;
00056     int err = nodeDB.getNodeWithEqn(eqns[i], nodePtr);
00057     if (err < 0) {
00058       nodeNumbers[i] = -1;
00059       dofOffsets[i] = -1;
00060       continue;
00061     }
00062 
00063     nodeNumbers[i] = nodePtr->getNodeNumber();
00064 
00065     int firstNodeEqn = nodePtr->getFieldEqnNumbers()[0];
00066     dofOffsets[i] = eqns[i] - firstNodeEqn;
00067   }
00068 }
00069 
00070 //------------------------------------------------------------------------------
00071 void convert_field_and_nodes_to_eqns(const NodeDatabase& nodeDB,
00072                                      int fieldID, int fieldSize,
00073                                      int numNodes, const GlobalID* nodeIDs,
00074                                      std::vector<int>& eqns)
00075 {
00076   eqns.assign(numNodes*fieldSize, -1);
00077 
00078   size_t offset = 0;
00079   for(int i=0; i<numNodes; ++i) {
00080     const NodeDescriptor* node = NULL;
00081     int err = nodeDB.getNodeWithID(nodeIDs[i], node);
00082     if (err < 0) {
00083       offset += fieldSize;
00084       continue;
00085     }
00086 
00087     int eqn = 0;
00088     node->getFieldEqnNumber(fieldID, eqn);
00089     for(int j=0; j<fieldSize; ++j) {
00090       eqns[offset++] = eqn+j;
00091     }
00092   }
00093 }
00094 
00095 //------------------------------------------------------------------------------
00096 FEDataFilter::FEDataFilter(FEI_Implementation* owner,
00097                            MPI_Comm comm,
00098                            SNL_FEI_Structure* probStruct,
00099                            LibraryWrapper* wrapper,
00100                            int masterRank)
00101  : Filter(probStruct),
00102    wrapper_(wrapper),
00103    feData_(NULL),
00104    useLookup_(true),
00105    internalFei_(0),
00106    newData_(false),
00107    localStartRow_(0),
00108    localEndRow_(0),
00109    numGlobalEqns_(0),
00110    reducedStartRow_(0),
00111    reducedEndRow_(0),
00112    numReducedRows_(0),
00113    iterations_(0),
00114    numRHSs_(0),
00115    currentRHS_(0),
00116    rhsIDs_(),
00117    outputLevel_(0),
00118    comm_(comm),
00119    masterRank_(masterRank),
00120    problemStructure_(probStruct),
00121    penCRIDs_(),
00122    rowIndices_(),
00123    rowColOffsets_(0),
00124    colIndices_(0),
00125    eqnCommMgr_(NULL),
00126    eqnCommMgr_put_(NULL),
00127    maxElemRows_(0),
00128    eStiff_(NULL),
00129    eStiff1D_(NULL),
00130    eLoad_(NULL),
00131    numRegularElems_(0),
00132    constraintBlocks_(0, 16),
00133    constraintNodeOffsets_(),
00134    packedFieldSizes_()
00135 {
00136   localRank_ = fei::localProc(comm_);
00137   numProcs_ = fei::numProcs(comm_);
00138 
00139   internalFei_ = 0;
00140 
00141   numRHSs_ = 1;
00142   rhsIDs_.resize(numRHSs_);
00143   rhsIDs_[0] = 0;
00144 
00145   eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00146   createEqnCommMgr_put();
00147 
00148   if (wrapper_->haveFiniteElementData()) {
00149     feData_ = wrapper_->getFiniteElementData();
00150   }
00151   else {
00152     FEI_CERR << "FEDataFilter::FEDataFilter ERROR, must be constructed with a "
00153          << "FiniteElementData interface. Aborting." << FEI_ENDL;
00154 #ifndef FEI_SER
00155     MPI_Abort(comm_, -1);
00156 #else
00157     abort();
00158 #endif
00159   }
00160 
00161   //We need to get the parameters from the owning FEI_Implementation, if we've
00162   //been given a non-NULL FEI_Implementation...
00163   if (owner != NULL) {
00164     int numParams = -1;
00165     char** paramStrings = NULL;
00166     int err = owner->getParameters(numParams, paramStrings);
00167 
00168     //Now let's pass them into our own parameter-handling mechanism.
00169     err = parameters(numParams, paramStrings);
00170     if (err != 0) {
00171       FEI_CERR << "FEDataFilter::FEDataFilter ERROR, parameters failed." << FEI_ENDL;
00172       MPI_Abort(comm_, -1);
00173     }
00174   }
00175 
00176   return;
00177 }
00178 
00179 //------------------------------------------------------------------------------
00180 FEDataFilter::FEDataFilter(const FEDataFilter& src)
00181  : Filter(NULL),
00182    wrapper_(NULL),
00183    feData_(NULL),
00184    useLookup_(true),
00185    internalFei_(0),
00186    newData_(false),
00187    localStartRow_(0),
00188    localEndRow_(0),
00189    numGlobalEqns_(0),
00190    reducedStartRow_(0),
00191    reducedEndRow_(0),
00192    numReducedRows_(0),
00193    iterations_(0),
00194    numRHSs_(0),
00195    currentRHS_(0),
00196    rhsIDs_(),
00197    outputLevel_(0),
00198    comm_(0),
00199    masterRank_(0),
00200    problemStructure_(NULL),
00201    penCRIDs_(),
00202    rowIndices_(),
00203    rowColOffsets_(0),
00204    colIndices_(0),
00205    eqnCommMgr_(NULL),
00206    eqnCommMgr_put_(NULL),
00207    maxElemRows_(0),
00208    eStiff_(NULL),
00209    eStiff1D_(NULL),
00210    eLoad_(NULL),
00211    numRegularElems_(0),
00212    constraintBlocks_(0, 16),
00213    constraintNodeOffsets_(),
00214    packedFieldSizes_()
00215 {
00216 }
00217 
00218 //------------------------------------------------------------------------------
00219 FEDataFilter::~FEDataFilter() {
00220 //
00221 //  Destructor function. Free allocated memory, etc.
00222 //
00223   numRHSs_ = 0;
00224 
00225   delete eqnCommMgr_;
00226   delete eqnCommMgr_put_;
00227 
00228   delete [] eStiff_;
00229   delete [] eStiff1D_;
00230   delete [] eLoad_;
00231 }
00232 
00233 //------------------------------------------------------------------------------
00234 int FEDataFilter::initialize()
00235 {
00236 // Determine final sparsity pattern for setting the structure of the
00237 // underlying sparse matrix.
00238 //
00239   debugOutput("#  initialize");
00240 
00241   // now, obtain the global equation info, such as how many equations there
00242   // are globally, and what the local starting and ending row-numbers are.
00243 
00244   // let's also get the number of global nodes, and a first-local-node-number.
00245   // node-number is a globally 0-based number we are assigning to nodes.
00246   // node-numbers are contiguous on a processor -- i.e., a processor owns a
00247   // contiguous block of node-numbers. This provides an easier-to-work-with
00248   // node numbering than the application-supplied nodeIDs which may not be
00249   // assumed to be contiguous or 0-based, or anything else.
00250 
00251   std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
00252   localStartRow_ = eqnOffsets[localRank_];
00253   localEndRow_ = eqnOffsets[localRank_+1]-1;
00254   numGlobalEqns_ = eqnOffsets[numProcs_];
00255 
00256   //--------------------------------------------------------------------------
00257   //  ----- end active equation calculations -----
00258 
00259   if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
00260   eqnCommMgr_ = NULL;
00261   if (eqnCommMgr_put_ != NULL) delete eqnCommMgr_put_;
00262   eqnCommMgr_put_ = NULL;
00263 
00264   eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00265   if (eqnCommMgr_ == NULL) ERReturn(-1);
00266 
00267   int err = createEqnCommMgr_put();
00268   if (err != 0) ERReturn(err);
00269 
00270   //(we need to set the number of RHSs in the eqn comm manager)
00271   eqnCommMgr_->setNumRHSs(numRHSs_);
00272 
00273   //let's let the underlying linear system know what the global offsets are.
00274   //While we're dealing with global offsets, we'll also calculate the starting
00275   //and ending 'reduced' rows, etc.
00276   CHK_ERR( initLinSysCore() );
00277 
00278   return(FEI_SUCCESS);
00279 }
00280 
00281 //------------------------------------------------------------------------------
00282 int FEDataFilter::createEqnCommMgr_put()
00283 {
00284   if (eqnCommMgr_put_ != NULL) return(0);
00285 
00286   eqnCommMgr_put_  = eqnCommMgr_->deepCopy();
00287   if (eqnCommMgr_put_ == NULL) ERReturn(-1);
00288 
00289   eqnCommMgr_put_->resetCoefs();
00290   eqnCommMgr_put_->accumulate_ = false;
00291   return(0);
00292 }
00293 
00294 //==============================================================================
00295 int FEDataFilter::initLinSysCore()
00296 {
00297   try {
00298 
00299   int err = wrapper_->getFiniteElementData()->setLookup(*problemStructure_);
00300 
00301   if (err != 0) {
00302     useLookup_ = false;
00303   }
00304 
00305   reducedStartRow_ = localStartRow_;
00306   reducedEndRow_ = localEndRow_;
00307 
00308   int numElemBlocks = problemStructure_->getNumElemBlocks();
00309   NodeDatabase& nodeDB     = problemStructure_->getNodeDatabase();
00310   NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
00311 
00312   int numNodes = nodeDB.getNumNodeDescriptors();
00313   int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
00314                          nodeCommMgr.getLocalNodeIDs().size();
00315   numNodes -= numRemoteNodes;
00316 
00317   int numSharedNodes = nodeCommMgr.getNumSharedNodes();
00318 
00319   std::vector<int> numElemsPerBlock(numElemBlocks);
00320   std::vector<int> numNodesPerElem(numElemBlocks);
00321   std::vector<int> elemMatrixSizePerBlock(numElemBlocks);
00322 
00323   for(int blk=0; blk<numElemBlocks; blk++) {
00324     BlockDescriptor* block = NULL;
00325     CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
00326 
00327     numElemsPerBlock[blk] = block->getNumElements();
00328     numNodesPerElem[blk]  = block->getNumNodesPerElement();
00329 
00330     int* fieldsPerNode = block->fieldsPerNodePtr();
00331     int** fieldIDsTable = block->fieldIDsTablePtr();
00332 
00333     elemMatrixSizePerBlock[blk] = 0;
00334 
00335     for(int nn=0; nn<numNodesPerElem[blk]; nn++) {
00336       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00337       
00338       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00339         elemMatrixSizePerBlock[blk] +=
00340           problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00341       }
00342     }
00343   }
00344 
00345   //Now we need to run the penalty constraint records and figure out how many
00346   //extra "element-blocks" to describe. (A penalty constraint will be treated 
00347   //exactly like an element.) So first, we need to figure out how many different
00348   //sizes of constraint connectivities there are, because the constraints with
00349   //the same numbers of constrained nodes will be grouped together in blocks.
00350 
00351   if (problemStructure_==NULL) {
00352     FEI_COUT << "problemStructrue_ NULL"<<FEI_ENDL;
00353     ERReturn(-1);
00354   }
00355 
00356   std::map<GlobalID,ConstraintType*>::const_iterator
00357     cr_iter = problemStructure_->getPenConstRecords().begin(),
00358     cr_end  = problemStructure_->getPenConstRecords().end();
00359 
00360   //constraintBlocks will be a sorted list with each "block-id" being the
00361   //num-nodes-per-constraint for constraints in that block.
00362 
00363   //numConstraintsPerBlock is the same length as constraintBlocks
00364   std::vector<int> numConstraintsPerBlock;
00365   std::vector<int> numDofPerConstraint;
00366   penCRIDs_.resize(problemStructure_->getNumPenConstRecords());
00367 
00368   int counter = 0;
00369   while(cr_iter != cr_end) {
00370     penCRIDs_[counter++] = (*cr_iter).first;
00371     ConstraintType& cr = *((*cr_iter).second);
00372     int nNodes = cr.getMasters().size();
00373 
00374     int insertPoint = -1;
00375     int offset = snl_fei::binarySearch(nNodes, constraintBlocks_, insertPoint);
00376 
00377     int nodeOffset = 0;
00378     if (offset < 0) {
00379       constraintBlocks_.insert(constraintBlocks_.begin()+insertPoint, nNodes);
00380       numConstraintsPerBlock.insert(numConstraintsPerBlock.begin()+insertPoint, 1);
00381       numDofPerConstraint.insert(numDofPerConstraint.begin()+insertPoint, 0);
00382 
00383       if (insertPoint > 0) {
00384         nodeOffset = constraintNodeOffsets_[insertPoint-1] +
00385            constraintBlocks_[insertPoint-1];
00386       }
00387       constraintNodeOffsets_.insert(constraintNodeOffsets_.begin()+insertPoint, nodeOffset);
00388       offset = insertPoint;
00389     }
00390     else {
00391       numConstraintsPerBlock[offset]++;
00392       ++cr_iter;
00393       continue;
00394     }
00395 
00396     std::vector<int>& fieldIDsvec = cr.getMasterFieldIDs();
00397     int* fieldIDs = &fieldIDsvec[0];
00398     for(int k=0; k<nNodes; ++k) {
00399       int fieldSize = problemStructure_->getFieldSize(fieldIDs[k]);
00400       packedFieldSizes_.insert(packedFieldSizes_.begin()+nodeOffset+k, fieldSize);
00401       numDofPerConstraint[offset] += fieldSize;
00402     }
00403     ++cr_iter;
00404   }
00405 
00406   //now combine the elem-block info with the penalty-constraint info.
00407   int numBlocksTotal = numElemBlocks + constraintBlocks_.size();
00408   for(size_t i=0; i<constraintBlocks_.size(); ++i) {
00409     numElemsPerBlock.push_back(numConstraintsPerBlock[i]);
00410     numNodesPerElem.push_back(constraintBlocks_[i]);
00411     elemMatrixSizePerBlock.push_back(numDofPerConstraint[i]);
00412   }
00413 
00414   int numMultCRs = problemStructure_->getNumMultConstRecords();
00415 
00416   CHK_ERR( feData_->describeStructure(numBlocksTotal,
00417                                       &numElemsPerBlock[0],
00418                                       &numNodesPerElem[0],
00419                                       &elemMatrixSizePerBlock[0],
00420                                       numNodes,
00421                                       numSharedNodes,
00422                                       numMultCRs) );
00423 
00424   numRegularElems_ = 0;
00425   std::vector<int> numDofPerNode;
00426 
00427   for(int i=0; i<numElemBlocks; i++) {
00428     BlockDescriptor* block = NULL;
00429     CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
00430 
00431     if (block->getNumElements() == 0) continue;
00432 
00433     ConnectivityTable& ctbl =
00434       problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
00435 
00436     std::vector<int> cNodeList(block->getNumNodesPerElement());
00437 
00438     int* fieldsPerNode = block->fieldsPerNodePtr();
00439     int** fieldIDsTable = block->fieldIDsTablePtr();
00440 
00441     numDofPerNode.resize(0);
00442     for(int nn=0; nn<numNodesPerElem[i]; nn++) {
00443       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00444       numDofPerNode.push_back(0);
00445       int indx = numDofPerNode.size()-1;
00446 
00447       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00448         numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00449       }
00450     }
00451 
00452     int nodesPerElement = block->getNumNodesPerElement();
00453     NodeDescriptor** elemConn = &((*ctbl.elem_conn_ptrs)[0]);
00454     int offset = 0;
00455     int numElems = block->getNumElements();
00456     numRegularElems_ += numElems;
00457     for(int j=0; j<numElems; j++) {
00458 
00459       for(int k=0; k<nodesPerElement; k++) {
00460         NodeDescriptor* node = elemConn[offset++];
00461         cNodeList[k] = node->getNodeNumber();
00462       }
00463 
00464       CHK_ERR( feData_->setConnectivity(i, ctbl.elemNumbers[j],
00465                                         block->getNumNodesPerElement(),
00466                                         &cNodeList[0],
00467                                         &numDofPerNode[0]) );
00468     }
00469   }
00470 
00471   std::vector<int> nodeNumbers;
00472   cr_iter = problemStructure_->getPenConstRecords().begin();
00473   int i = 0;
00474   while(cr_iter != cr_end) {
00475     ConstraintType& cr = *((*cr_iter).second);
00476     std::vector<GlobalID>& nodeIDsvec = cr.getMasters();
00477     GlobalID* nodeIDs = &nodeIDsvec[0];
00478     int nNodes = cr.getMasters().size();
00479     int index = snl_fei::binarySearch(nNodes, constraintBlocks_);
00480     if (index < 0) {
00481       ERReturn(-1);
00482     }
00483 
00484     int blockNum = numElemBlocks + index;
00485 
00486     nodeNumbers.resize(nNodes);
00487 
00488     for(int k=0; k<nNodes; ++k) {
00489       const NodeDescriptor& node = Filter::findNodeDescriptor(nodeIDs[k]);
00490       nodeNumbers[k] = node.getNodeNumber();
00491     }
00492 
00493     int offset = constraintNodeOffsets_[index];
00494     CHK_ERR( feData_->setConnectivity(blockNum, numRegularElems_+i++,
00495                                       nNodes, &nodeNumbers[0],
00496                                       &packedFieldSizes_[offset]) );
00497     ++cr_iter;
00498   }
00499 
00500   }
00501   catch(std::runtime_error& exc) {
00502     FEI_CERR << exc.what() << FEI_ENDL;
00503     ERReturn(-1);
00504   }
00505 
00506   return(FEI_SUCCESS);
00507 }
00508 
00509 //------------------------------------------------------------------------------
00510 int FEDataFilter::resetSystem(double s)
00511 {
00512   //
00513   //  This puts the value s throughout both the matrix and the vector.
00514   //
00515   if (Filter::logStream() != NULL) {
00516     (*logStream()) << "FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
00517   }
00518 
00519   CHK_ERR( feData_->reset() );
00520  
00521   debugOutput("#FEDataFilter leaving resetSystem");
00522 
00523   return(FEI_SUCCESS);
00524 }
00525 
00526 //------------------------------------------------------------------------------
00527 int FEDataFilter::deleteMultCRs()
00528 {
00529 
00530   debugOutput("#FEDataFilter::deleteMultCRs");
00531 
00532   int err = feData_->deleteConstraints();
00533 
00534   debugOutput("#FEDataFilter leaving deleteMultCRs");
00535 
00536   return(err);
00537 }
00538 
00539 //------------------------------------------------------------------------------
00540 int FEDataFilter::resetTheMatrix(double s)
00541 {
00542   //FiniteElementData implementations can't currently reset the matrix without
00543   //resetting the rhs vector too. 
00544   return(FEI_SUCCESS);
00545 }
00546 
00547 //------------------------------------------------------------------------------
00548 int FEDataFilter::resetTheRHSVector(double s)
00549 {
00550   //FiniteElementData implementations can't currently reset the rhs vector
00551   //without resetting the matrix too.
00552   return(FEI_SUCCESS);
00553 }
00554 
00555 //------------------------------------------------------------------------------
00556 int FEDataFilter::resetMatrix(double s)
00557 {
00558   //
00559   //  This puts the value s throughout both the matrix and the vector.
00560   //
00561 
00562   debugOutput("FEI: resetMatrix");
00563 
00564   CHK_ERR( resetTheMatrix(s) );
00565 
00566   eqnCommMgr_->resetCoefs();
00567 
00568   debugOutput("#FEDataFilter leaving resetMatrix");
00569 
00570   return(FEI_SUCCESS);
00571 }
00572 
00573 //------------------------------------------------------------------------------
00574 int FEDataFilter::resetRHSVector(double s)
00575 {
00576   //
00577   //  This puts the value s throughout the rhs vector.
00578   //
00579 
00580   debugOutput("FEI: resetRHSVector");
00581 
00582   CHK_ERR( resetTheRHSVector(s) );
00583 
00584   eqnCommMgr_->resetCoefs();
00585 
00586   debugOutput("#FEDataFilter leaving resetRHSVector");
00587 
00588   return(FEI_SUCCESS);
00589 }
00590 
00591 //------------------------------------------------------------------------------
00592 int FEDataFilter::resetInitialGuess(double s)
00593 {
00594   //
00595   //  This puts the value s throughout the initial guess (solution) vector.
00596   //
00597   if (Filter::logStream() != NULL) {
00598     FEI_OSTREAM& os = *logStream();
00599     os << "FEI: resetInitialGuess" << FEI_ENDL;
00600     os << "#value to which initial guess is to be set" << FEI_ENDL;
00601     os << s << FEI_ENDL;
00602   }
00603 
00604   //Actually, the FiniteElementData doesn't currently allow us to alter
00605   //values in any initial guess or solution vector.
00606 
00607   debugOutput("#FEDataFilter leaving resetInitialGuess");
00608 
00609   return(FEI_SUCCESS);
00610 }
00611 
00612 //------------------------------------------------------------------------------
00613 int FEDataFilter::loadNodeBCs(int numNodes,
00614                               const GlobalID *nodeIDs,
00615                               int fieldID,
00616                               const int* offsetsIntoField,
00617                               const double* prescribedValues)
00618 {
00619   //
00620   //  load boundary condition information for a given set of nodes
00621   //
00622   int size = problemStructure_->getFieldSize(fieldID);
00623   if (size < 1) {
00624     FEI_CERR << "FEI Warning: loadNodeBCs called for fieldID "<<fieldID
00625          <<", which was defined with size "<<size<<" (should be positive)."<<FEI_ENDL;
00626     return(0);
00627   }
00628 
00629    if (Filter::logStream() != NULL) {
00630     (*logStream())<<"FEI: loadNodeBCs"<<FEI_ENDL
00631                      <<"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
00632                      <<"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
00633                      <<"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
00634     (*logStream())<<"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
00635 
00636     for(int j=0; j<numNodes; j++) {
00637       int nodeID = nodeIDs[j];
00638       (*logStream())<<nodeID<<"  ";
00639       (*logStream())<< offsetsIntoField[j]<<" ";
00640       (*logStream())<< prescribedValues[j]<<FEI_ENDL;
00641     }
00642    }
00643 
00644    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00645 
00646    std::vector<int> essEqns(numNodes);
00647    std::vector<double> alpha(numNodes);
00648    std::vector<double> gamma(numNodes);
00649 
00650    for(int i=0; i<numNodes; ++i) {
00651      NodeDescriptor* node = 0;
00652      nodeDB.getNodeWithID(nodeIDs[i], node);
00653      int eqn = -1;
00654      if (!node->getFieldEqnNumber(fieldID, eqn)) {
00655        ERReturn(-1);
00656      }
00657 
00658      essEqns[i] = eqn + offsetsIntoField[i];
00659      gamma[i] = prescribedValues[i];
00660      alpha[i] = 1.0;
00661    }
00662 
00663    if (essEqns.size() > 0) {
00664       CHK_ERR( enforceEssentialBCs(&essEqns[0], &alpha[0],
00665                                    &gamma[0], essEqns.size()) );
00666    }
00667 
00668    return(FEI_SUCCESS);
00669 }
00670 
00671 //------------------------------------------------------------------------------
00672 int FEDataFilter::loadElemBCs(int numElems,
00673                               const GlobalID *elemIDs,
00674                               int fieldID,
00675                               const double *const *alpha,
00676                               const double *const *beta,
00677                               const double *const *gamma)
00678 {
00679    return(-1);
00680 }
00681 
00682 //------------------------------------------------------------------------------
00683 void FEDataFilter::allocElemStuff()
00684 {
00685    int nb = problemStructure_->getNumElemBlocks();
00686 
00687    for(int i=0; i<nb; i++) {
00688      BlockDescriptor* block = NULL;
00689      int err = problemStructure_->getBlockDescriptor_index(i, block);
00690      if (err) voidERReturn;
00691 
00692       int numEqns = block->getNumEqnsPerElement();
00693       if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
00694    }
00695 
00696    eStiff_ = new double*[maxElemRows_];
00697    eStiff1D_ = new double[maxElemRows_*maxElemRows_];
00698 
00699    if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
00700 
00701    for(int r=0; r<maxElemRows_; r++) {
00702       eStiff_[r] = eStiff1D_ + r*maxElemRows_;
00703    }
00704 
00705    eLoad_ = new double[maxElemRows_];
00706 
00707    if (eLoad_ == NULL) voidERReturn
00708 }
00709 
00710 //------------------------------------------------------------------------------
00711 int FEDataFilter::sumInElem(GlobalID elemBlockID,
00712                         GlobalID elemID,
00713                         const GlobalID* elemConn,
00714                         const double* const* elemStiffness,
00715                         const double* elemLoad,
00716                         int elemFormat)
00717 {
00718   if (Filter::logStream() != NULL) {
00719     (*logStream()) << "FEI: sumInElem" << FEI_ENDL <<"# elemBlockID " << FEI_ENDL
00720                       << static_cast<int>(elemBlockID) << FEI_ENDL
00721                       << "# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
00722     BlockDescriptor* block = NULL;
00723     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
00724     int numNodes = block->getNumNodesPerElement();
00725     (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
00726     (*logStream()) << "#connected nodes" << FEI_ENDL;
00727     for(int i=0; i<numNodes; ++i) {
00728       GlobalID nodeID = elemConn[i];
00729       (*logStream())<<static_cast<int>(nodeID)<<" ";
00730     }
00731     (*logStream())<<FEI_ENDL;
00732   }
00733 
00734   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
00735                           elemLoad, elemFormat));
00736 }
00737 
00738 //------------------------------------------------------------------------------
00739 int FEDataFilter::sumInElemMatrix(GlobalID elemBlockID,
00740                               GlobalID elemID,
00741                               const GlobalID* elemConn,
00742                               const double* const* elemStiffness,
00743                               int elemFormat)
00744 {
00745   if (Filter::logStream() != NULL) {
00746     (*logStream()) << "FEI: sumInElemMatrix"<<FEI_ENDL
00747                       << "#elemBlockID" << FEI_ENDL << static_cast<int>(elemBlockID)
00748                       << "# elemID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
00749     BlockDescriptor* block = NULL;
00750     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
00751     int numNodes = block->getNumNodesPerElement();
00752     (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
00753     (*logStream()) << "#connected nodes" << FEI_ENDL;
00754     for(int i=0; i<numNodes; ++i) {
00755       GlobalID nodeID = elemConn[i];
00756       (*logStream())<<static_cast<int>(nodeID)<<" ";
00757     }
00758     (*logStream())<<FEI_ENDL;
00759   }
00760 
00761   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
00762                           NULL, elemFormat));
00763 }
00764 
00765 //------------------------------------------------------------------------------
00766 int FEDataFilter::sumInElemRHS(GlobalID elemBlockID,
00767                            GlobalID elemID,
00768                            const GlobalID* elemConn,
00769                            const double* elemLoad)
00770 {
00771   if (Filter::logStream() != NULL) {
00772     (*logStream()) << "FEI: sumInElemRHS"<<FEI_ENDL<<"# elemBlockID " << FEI_ENDL
00773                       <<static_cast<int>(elemBlockID)
00774                       << "# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
00775     BlockDescriptor* block = NULL;
00776     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
00777     int numNodes = block->getNumNodesPerElement();
00778     (*logStream()) << "#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
00779     (*logStream()) << "#connected nodes" << FEI_ENDL;
00780     for(int i=0; i<numNodes; ++i) {
00781       GlobalID nodeID = elemConn[i];
00782       (*logStream())<<static_cast<int>(nodeID)<<" ";
00783     }
00784     (*logStream())<<FEI_ENDL;
00785   }
00786 
00787   return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
00788                           elemLoad, -1));
00789 }
00790 
00791 //------------------------------------------------------------------------------
00792 int FEDataFilter::generalElemInput(GlobalID elemBlockID,
00793                         GlobalID elemID,
00794                         const GlobalID* elemConn,
00795                         const double* const* elemStiffness,
00796                         const double* elemLoad,
00797                         int elemFormat)
00798 {
00799   (void)elemConn;
00800   return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
00801                           elemFormat) );
00802 }
00803 
00804 //------------------------------------------------------------------------------
00805 int FEDataFilter::generalElemInput(GlobalID elemBlockID,
00806                         GlobalID elemID,
00807                         const double* const* elemStiffness,
00808                         const double* elemLoad,
00809                         int elemFormat)
00810 {
00811   //first get the block-descriptor for this elemBlockID...
00812 
00813   BlockDescriptor* block = NULL;
00814   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
00815 
00816   //now allocate our local stiffness/load copy if we haven't already.
00817 
00818   if (maxElemRows_ <= 0) allocElemStuff();
00819 
00820   int numElemRows = block->getNumEqnsPerElement();
00821 
00822   //an std::vector.resize operation is free if the size is either shrinking or
00823   //staying the same.
00824 
00825   const double* const* stiff = NULL;
00826   if (elemStiffness != NULL) stiff = elemStiffness;
00827 
00828   const double* load = NULL;
00829   if (elemLoad != NULL) load = elemLoad;
00830 
00831   //we'll make a local dense copy of the element stiffness array
00832   //if the stiffness array was passed in using one of the "weird"
00833   //element formats, AND if the stiffness array is non-null.
00834   if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
00835     Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
00836     stiff = eStiff_;
00837   }
00838 
00839   if (stiff != NULL || load != NULL) newData_ = true;
00840 
00841   if (Filter::logStream() != NULL) {
00842     if (stiff != NULL) {
00843       (*logStream())
00844         << "#numElemRows"<< FEI_ENDL << numElemRows << FEI_ENDL
00845         << "#elem-stiff (after being copied into dense-row format)"
00846         << FEI_ENDL;
00847       for(int i=0; i<numElemRows; i++) {
00848         const double* stiff_i = stiff[i];
00849         for(int j=0; j<numElemRows; j++) {
00850           (*logStream()) << stiff_i[j] << " ";
00851         }
00852         (*logStream()) << FEI_ENDL;
00853       }
00854     }
00855 
00856     if (load != NULL) {
00857       (*logStream()) << "#elem-load" << FEI_ENDL;
00858       for(int i=0; i<numElemRows; i++) {
00859         (*logStream()) << load[i] << " ";
00860       }
00861       (*logStream())<<FEI_ENDL;
00862     }
00863   }
00864 
00865   //Now we'll proceed to gather the stuff we need to pass the stiffness
00866   //data through to the FiniteElementData interface...
00867 
00868   int blockNumber = problemStructure_->getIndexOfBlock(elemBlockID);
00869 
00870   ConnectivityTable& connTable = problemStructure_->
00871     getBlockConnectivity(elemBlockID);
00872 
00873   std::map<GlobalID,int>::iterator
00874     iter = connTable.elemIDs.find(elemID);
00875   if (iter == connTable.elemIDs.end()) {
00876     ERReturn(-1);
00877   }
00878 
00879   int elemIndex = iter->second;
00880 
00881   int elemNumber = connTable.elemNumbers[elemIndex];
00882 
00883   int numNodes = block->getNumNodesPerElement();
00884   int* fieldsPerNode = block->fieldsPerNodePtr();
00885   int** fieldIDsTable = block->fieldIDsTablePtr();
00886 
00887   int numDistinctFields = block->getNumDistinctFields();
00888   int fieldSize = 0;
00889   if (numDistinctFields == 1) {
00890     fieldSize = problemStructure_->getFieldSize(fieldIDsTable[0][0]);
00891   }
00892 
00893   static std::vector<int> iwork;
00894   iwork.resize(2*numNodes);
00895 
00896   int* dofsPerNode = &iwork[0];
00897   for(int i=0; i<numNodes; ++i) {
00898     dofsPerNode[i] = 0;
00899   }
00900 
00901   int* nodeNumbers = dofsPerNode+numNodes;
00902 
00903   NodeDescriptor** elemNodes =
00904     &((*connTable.elem_conn_ptrs)[elemIndex*numNodes]);
00905 
00906   for(int nn=0; nn<numNodes; nn++) {
00907     NodeDescriptor* node = elemNodes[nn];
00908     nodeNumbers[nn] = node->getNodeNumber();
00909 
00910     if (numDistinctFields == 1) {
00911       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00912         dofsPerNode[nn] += fieldSize;
00913       }
00914     }
00915     else {
00916       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00917         dofsPerNode[nn] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00918       }
00919     }
00920   }
00921 
00922   if (stiff != NULL) {
00923     CHK_ERR( feData_->setElemMatrix(blockNumber, elemNumber, numNodes,
00924                                     nodeNumbers, dofsPerNode, stiff) );
00925   }
00926 
00927   if (load != NULL) {
00928     CHK_ERR( feData_->setElemVector(blockNumber, elemNumber, numNodes,
00929                                     nodeNumbers, dofsPerNode, load) );
00930   }
00931 
00932   return(FEI_SUCCESS);
00933 }
00934 
00935 //------------------------------------------------------------------------------
00936 int FEDataFilter::putIntoRHS(int IDType,
00937                           int fieldID,
00938                           int numIDs,
00939                           const GlobalID* IDs,
00940                           const double* rhsEntries)
00941 {
00942   int fieldSize = problemStructure_->getFieldSize(fieldID);
00943 
00944   rowIndices_.resize(fieldSize*numIDs);
00945   int checkNumEqns;
00946 
00947   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
00948                                             checkNumEqns,
00949                                             &rowIndices_[0]));
00950   if (checkNumEqns != numIDs*fieldSize) {
00951     ERReturn(-1);
00952   }
00953 
00954   CHK_ERR( exchangeRemoteEquations() );
00955 
00956   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
00957 
00958   return(0);
00959 }
00960 
00961 //------------------------------------------------------------------------------
00962 int FEDataFilter::sumIntoRHS(int IDType,
00963                           int fieldID,
00964                           int numIDs,
00965                           const GlobalID* IDs,
00966                           const double* rhsEntries)
00967 {
00968   int fieldSize = problemStructure_->getFieldSize(fieldID);
00969 
00970   rowIndices_.resize(fieldSize*numIDs);
00971   int checkNumEqns;
00972 
00973   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
00974                                             checkNumEqns,
00975                                             &rowIndices_[0]));
00976   if (checkNumEqns != numIDs*fieldSize) {
00977     ERReturn(-1);
00978   }
00979 
00980   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
00981 
00982   return(0);
00983 }
00984 
00985 //------------------------------------------------------------------------------
00986 int FEDataFilter::sumIntoMatrixDiagonal(int  IDType,
00987                              int  fieldID,
00988                              int  numIDs,
00989                              const GlobalID*  IDs,
00990                              const double*  coefficients)
00991 {
00992   int fieldSize = problemStructure_->getFieldSize(fieldID);
00993   const NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00994 
00995   std::vector<int> eqns;
00996   convert_field_and_nodes_to_eqns(nodeDB, fieldID, fieldSize, numIDs, IDs, eqns);
00997 
00998   std::vector<int> nodeNumbers, dofOffsets;
00999   convert_eqns_to_nodenumbers_and_offsets(nodeDB, eqns.size(), &eqns[0],
01000                                       nodeNumbers, dofOffsets);
01001 
01002   std::vector<int> ones(nodeNumbers.size(), 1);
01003 
01004   CHK_ERR( feData_->sumIntoMatrix(nodeNumbers.size(), &nodeNumbers[0], &dofOffsets[0],
01005                                   &ones[0], &nodeNumbers[0], &dofOffsets[0], coefficients));
01006   return( 0 );
01007 }
01008 
01009 //------------------------------------------------------------------------------
01010 int FEDataFilter::enforceEssentialBCs(const int* eqns, 
01011                                       const double* alpha,
01012                                       const double* gamma, 
01013                                       int numEqns)
01014 {
01015   std::vector<double> values;
01016   std::vector<int> nodeNumbers;
01017   std::vector<int> dofOffsets;
01018 
01019   for(int i=0; i<numEqns; i++) {
01020     int reducedEqn = -1;
01021     bool isSlave = problemStructure_->
01022       translateToReducedEqn(eqns[i], reducedEqn);
01023     if (isSlave) continue;
01024 
01025     int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqns[i]);
01026 
01027     nodeNumbers.push_back(nodeNumber);
01028 
01029     const NodeDescriptor* node = NULL;
01030     CHK_ERR( problemStructure_->getNodeDatabase().
01031              getNodeWithNumber(nodeNumber, node));
01032 
01033     int firstEqn = node->getFieldEqnNumbers()[0];
01034     dofOffsets.push_back(eqns[i] - firstEqn);
01035 
01036     values.push_back(gamma[i]/alpha[i]);
01037   }
01038 
01039   CHK_ERR( feData_->setDirichletBCs(nodeNumbers.size(),
01040                                     &nodeNumbers[0], &dofOffsets[0], &values[0]) );
01041 
01042   newData_ = true;
01043 
01044   return(FEI_SUCCESS);
01045 }
01046 
01047 //------------------------------------------------------------------------------
01048 int FEDataFilter::loadFEDataMultCR(int CRID,
01049                                int numCRNodes,
01050                                const GlobalID* CRNodes, 
01051                                const int* CRFields,
01052                                const double* CRWeights,
01053                                double CRValue)
01054 {
01055   if (Filter::logStream() != NULL) {
01056     FEI_OSTREAM& os = *logStream();
01057     os<<"FEI: loadCRMult"<<FEI_ENDL;
01058     os<<"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
01059     os<<"#CRNodes:"<<FEI_ENDL;
01060     int i;
01061     for(i=0; i<numCRNodes; ++i) {
01062       GlobalID nodeID = CRNodes[i];
01063       os << static_cast<int>(nodeID) << " ";
01064     }
01065     os << FEI_ENDL << "#fields:"<<FEI_ENDL;
01066     for(i=0; i<numCRNodes; ++i) os << CRFields[i] << " ";
01067     os << FEI_ENDL << "#field-sizes:"<<FEI_ENDL;
01068     for(i=0; i<numCRNodes; ++i) {
01069       int size = problemStructure_->getFieldSize(CRFields[i]);
01070       os << size << " ";
01071     }
01072     os << FEI_ENDL<<"#weights:"<<FEI_ENDL;
01073     int offset = 0;
01074     for(i=0; i<numCRNodes; ++i) {
01075       int size = problemStructure_->getFieldSize(CRFields[i]);
01076       for(int j=0; j<size; ++j) {
01077         os << CRWeights[offset++] << " ";
01078       }
01079     }
01080     os << FEI_ENDL<<"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
01081   }
01082 
01083   if (numCRNodes <= 0) return(0);
01084 
01085   std::vector<int> nodeNumbers;
01086   std::vector<int> dofOffsets;
01087   std::vector<double> weights;
01088 
01089   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01090 
01091   double fei_min = std::numeric_limits<double>::min();
01092 
01093   int offset = 0;
01094   for(int i=0; i<numCRNodes; i++) {
01095     NodeDescriptor* node = NULL;
01096     CHK_ERR( nodeDB.getNodeWithID(CRNodes[i], node) );
01097 
01098     int firstEqn = node->getFieldEqnNumbers()[0];
01099     int fieldEqn = -1;
01100     bool hasField = node->getFieldEqnNumber(CRFields[i], fieldEqn);
01101     if (!hasField) ERReturn(-1);
01102 
01103     int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
01104 
01105     for(int f=0; f<fieldSize; f++) {
01106       double weight = CRWeights[offset++];
01107       if (std::abs(weight) > fei_min) {
01108         nodeNumbers.push_back(node->getNodeNumber());
01109         dofOffsets.push_back((fieldEqn+f)-firstEqn);
01110         weights.push_back(weight);
01111       }
01112     }
01113   }
01114 
01115   CHK_ERR( feData_->setMultiplierCR(CRID, nodeNumbers.size(),
01116                                     &nodeNumbers[0], &dofOffsets[0],
01117                                     &weights[0], CRValue) );
01118   newData_ = true;
01119 
01120   return(0);
01121 }
01122 
01123 //------------------------------------------------------------------------------
01124 int FEDataFilter::loadFEDataPenCR(int CRID,
01125                               int numCRNodes,
01126                               const GlobalID* CRNodes, 
01127                               const int* CRFields,
01128                               const double* CRWeights,
01129                               double CRValue, 
01130                               double penValue)
01131 {
01132   if (numCRNodes <= 0) return(0);
01133 
01134   std::vector<int> nodeNumbers;
01135   std::vector<int> dofsPerNode;
01136   std::vector<double> weights;
01137 
01138   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01139 
01140   int offset = 0;
01141   for(int i=0; i<numCRNodes; i++) {
01142     NodeDescriptor* node = NULL;
01143     CHK_ERR( nodeDB.getNodeWithID(CRNodes[i], node) );
01144 
01145     int fieldEqn = -1;
01146     bool hasField = node->getFieldEqnNumber(CRFields[i], fieldEqn);
01147     if (!hasField) ERReturn(-1);
01148 
01149     int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
01150 
01151     nodeNumbers.push_back(node->getNodeNumber());
01152     dofsPerNode.push_back(fieldSize);
01153 
01154     for(int f=0; f<fieldSize; f++) {
01155       double weight = CRWeights[offset++];
01156       weights.push_back(weight);
01157     }
01158   }
01159 
01160   std::vector<double*> matrixCoefs(weights.size());
01161   std::vector<double> rhsCoefs(weights.size());
01162   offset = 0;
01163   for(size_t i=0; i<weights.size(); ++i) {
01164     double* coefPtr = new double[weights.size()];
01165     for(size_t j=0; j<weights.size(); ++j) {
01166       coefPtr[j] = weights[i]*weights[j]*penValue;
01167     }
01168     matrixCoefs[i] = coefPtr;
01169     rhsCoefs[i] = weights[i]*penValue*CRValue;
01170   }
01171 
01172   int crIndex = snl_fei::binarySearch(CRID, penCRIDs_);
01173 
01174   int index = snl_fei::binarySearch(numCRNodes, constraintBlocks_);
01175 
01176   int blockNum = problemStructure_->getNumElemBlocks() + index;
01177   int elemNum = numRegularElems_ + crIndex;
01178 
01179   CHK_ERR( feData_->setElemMatrix(blockNum, elemNum,
01180                                   nodeNumbers.size(),
01181                                   &nodeNumbers[0],
01182                                   &dofsPerNode[0],
01183                                   &matrixCoefs[0]) );
01184 
01185   CHK_ERR( feData_->setElemVector(blockNum, elemNum, nodeNumbers.size(),
01186                                   &nodeNumbers[0], &dofsPerNode[0], &rhsCoefs[0]) );
01187 
01188   newData_ = true;
01189 
01190   for(size_t i=0; i<weights.size(); ++i) {
01191     delete [] matrixCoefs[i];
01192   }
01193 
01194   return(0);
01195 }
01196 
01197 //------------------------------------------------------------------------------
01198 int FEDataFilter::loadCRMult(int CRID, 
01199                          int numCRNodes,
01200                          const GlobalID* CRNodes, 
01201                          const int* CRFields,
01202                          const double* CRWeights,
01203                          double CRValue)
01204 {
01205 //
01206 // Load Lagrange multiplier constraint relation data
01207 //
01208 
01209   //Give the constraint data to the underlying solver using this special function...
01210   CHK_ERR( loadFEDataMultCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
01211                             CRValue) );
01212 
01213   return(FEI_SUCCESS);
01214 }
01215 
01216 //------------------------------------------------------------------------------
01217 int FEDataFilter::loadCRPen(int CRID, 
01218                         int numCRNodes, 
01219                         const GlobalID* CRNodes,
01220                         const int* CRFields,
01221                         const double* CRWeights,
01222                         double CRValue,
01223                         double penValue)
01224 {
01225 //
01226 // Load penalty constraint relation data
01227 //
01228 
01229    debugOutput("FEI: loadCRPen");
01230 
01231    //Give the constraint data to the underlying solver using this special function...
01232    CHK_ERR( loadFEDataPenCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
01233                             CRValue, penValue) );
01234 
01235    return(FEI_SUCCESS);
01236 }
01237 
01238 //------------------------------------------------------------------------------
01239 int FEDataFilter::parameters(int numParams, const char *const* paramStrings)
01240 {
01241 //
01242 // this function takes parameters for setting internal things like solver
01243 // and preconditioner choice, etc.
01244 //
01245    if (numParams == 0 || paramStrings == NULL) {
01246       debugOutput("#FEDataFilter::parameters --- no parameters.");
01247    }
01248    else {
01249 
01250       snl_fei::getIntParamValue("outputLevel",numParams, paramStrings,
01251                                 outputLevel_);
01252 
01253       snl_fei::getIntParamValue("internalFei",numParams,paramStrings,
01254                                 internalFei_);
01255 
01256       if (Filter::logStream() != NULL) {
01257         (*logStream())<<"#FEDataFilter::parameters"<<FEI_ENDL
01258                          <<"# --- numParams: "<< numParams<<FEI_ENDL;
01259          for(int i=0; i<numParams; i++){
01260            (*logStream())<<"#------ paramStrings["<<i<<"]: "
01261                             <<paramStrings[i]<<FEI_ENDL;
01262          }
01263       }
01264    }
01265 
01266    CHK_ERR( Filter::parameters(numParams, paramStrings) );
01267 
01268    debugOutput("#FEDataFilter leaving parameters function");
01269 
01270    return(FEI_SUCCESS);
01271 }
01272 
01273 //------------------------------------------------------------------------------
01274 int FEDataFilter::loadComplete()
01275 {
01276   debugOutput("#FEDataFilter calling FEData matrixLoadComplete");
01277 
01278   CHK_ERR( feData_->loadComplete() );
01279 
01280   newData_ = false;
01281 
01282   return(0);
01283 }
01284 
01285 //------------------------------------------------------------------------------
01286 int FEDataFilter::residualNorm(int whichNorm, int numFields,
01287                            int* fieldIDs, double* norms, double& residTime)
01288 {
01289 //
01290 //This function can do 3 kinds of norms: infinity-norm (denoted
01291 //by whichNorm==0), 1-norm and 2-norm.
01292 //
01293    debugOutput("FEI: residualNorm");
01294 
01295    CHK_ERR( loadComplete() );
01296 
01297    //for now, FiniteElementData doesn't do residual calculations.
01298 
01299    int fdbNumFields = problemStructure_->getNumFields();
01300    const int* fdbFieldIDs = problemStructure_->getFieldIDsPtr();
01301 
01302    int i;
01303 
01304    //Since we don't calculate actual residual norms, we'll fill the user's
01305    //array with norm data that is obviously not real norm data.
01306    int offset = 0;
01307    i = 0;
01308    while(offset < numFields && i < fdbNumFields) {
01309      if (fdbFieldIDs[i] >= 0) {
01310        fieldIDs[offset++] = fdbFieldIDs[i];
01311      }
01312      ++i;
01313    }
01314    for(i=0; i<numFields; ++i) {
01315       norms[i] = -99.9;
01316    }
01317 
01318    //fill out the end of the array with garbage in case the user-provided
01319    //array is longer than the list of fields we have in fieldDB.
01320    for(i=offset; i<numFields; ++i) {
01321       fieldIDs[i] = -99;
01322    }
01323 
01324    return(FEI_SUCCESS);
01325 }
01326 
01327 //------------------------------------------------------------------------------
01328 int FEDataFilter::formResidual(double* residValues, int numLocalEqns)
01329 {
01330   //FiniteElementData implementations can't currently do residuals.
01331   return(FEI_SUCCESS);
01332 }
01333 
01334 //------------------------------------------------------------------------------
01335 int FEDataFilter::solve(int& status, double& sTime) {
01336 
01337    debugOutput("FEI: solve");
01338 
01339    CHK_ERR( loadComplete() );
01340 
01341    debugOutput("#FEDataFilter in solve, calling launchSolver...");
01342  
01343    double start = MPI_Wtime();
01344 
01345    CHK_ERR( feData_->launchSolver(status, iterations_) );
01346 
01347    sTime = MPI_Wtime() - start;
01348 
01349    debugOutput("#FEDataFilter... back from solver");
01350  
01351    //now unpack the locally-owned shared entries of the solution vector into
01352    //the eqn-comm-mgr data structures.
01353    CHK_ERR( unpackSolution() );
01354 
01355    debugOutput("#FEDataFilter leaving solve");
01356 
01357    if (status != 0) return(1);
01358    else return(FEI_SUCCESS);
01359 }
01360 
01361 //------------------------------------------------------------------------------
01362 int FEDataFilter::setNumRHSVectors(int numRHSs, int* rhsIDs){
01363 
01364    if (numRHSs < 0) {
01365       FEI_CERR << "FEDataFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
01366       ERReturn(-1);
01367    }
01368 
01369    numRHSs_ = numRHSs;
01370 
01371    rhsIDs_.resize(numRHSs_);
01372    for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
01373 
01374   //(we need to set the number of RHSs in the eqn comm manager)
01375   eqnCommMgr_->setNumRHSs(numRHSs_);
01376 
01377    return(FEI_SUCCESS);
01378 }
01379 
01380 //------------------------------------------------------------------------------
01381 int FEDataFilter::setCurrentRHS(int rhsID)
01382 {
01383    std::vector<int>::iterator iter =
01384        std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
01385 
01386    if (iter == rhsIDs_.end()) ERReturn(-1)
01387    
01388    currentRHS_ = iter - rhsIDs_.begin();
01389 
01390    return(FEI_SUCCESS);
01391 }
01392 
01393 //------------------------------------------------------------------------------
01394 int FEDataFilter::giveToMatrix(int numPtRows, const int* ptRows,
01395                            int numPtCols, const int* ptCols,
01396                            const double* const* values, int mode)
01397 {
01398   //This isn't going to be fast... I need to optimize the whole structure
01399   //of code that's associated with passing data to FiniteElementData.
01400 
01401   std::vector<int> rowNodeNumbers, rowDofOffsets, colNodeNumbers, colDofOffsets;
01402   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01403   int i;
01404 
01405   //First, we have to get nodeNumbers and dofOffsets for each of the
01406   //row-numbers and col-numbers.
01407 
01408   for(i=0; i<numPtRows; i++) {
01409     int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
01410     if (nodeNumber < 0) ERReturn(-1);
01411     int fieldID = problemStructure_->getAssociatedFieldID(ptRows[i]);
01412     if (fieldID < 0) ERReturn(-1);
01413     const NodeDescriptor* node = NULL;
01414     CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
01415     int firstEqn = node->getFieldEqnNumbers()[0];
01416 
01417     rowNodeNumbers.push_back(nodeNumber);
01418     rowDofOffsets.push_back(ptRows[i] - firstEqn);
01419   }
01420 
01421   for(i=0; i<numPtCols; i++) {
01422     int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
01423     if (nodeNumber < 0) ERReturn(-1);
01424     int fieldID = problemStructure_->getAssociatedFieldID(ptCols[i]);
01425     if (fieldID < 0) ERReturn(-1);
01426     const NodeDescriptor* node = NULL;
01427     CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
01428     int firstEqn = node->getFieldEqnNumbers()[0];
01429 
01430     colNodeNumbers.push_back(nodeNumber);
01431     colDofOffsets.push_back(ptCols[i] - firstEqn);
01432   }
01433 
01434   //now we have to flatten the colNodeNumbers and colDofOffsets out into
01435   //an array of length numPtRows*numPtCols, where the nodeNumbers and
01436   //dofOffsets are repeated 'numPtRows' times.
01437 
01438   int len = numPtRows*numPtCols;
01439   std::vector<int> allColNodeNumbers(len), allColDofOffsets(len);
01440   std::vector<double> allValues(len);
01441 
01442   int offset = 0;
01443   for(i=0; i<numPtRows; i++) {
01444     for(int j=0; j<numPtCols; j++) {
01445       allColNodeNumbers[offset] = colNodeNumbers[j];
01446       allColDofOffsets[offset] = colDofOffsets[j];
01447       allValues[offset++] = values[i][j];
01448     }
01449   }
01450 
01451   //while we're at it, let's make an array with numPtCols replicated in it
01452   //'numPtRows' times.
01453   std::vector<int> numColsPerRow(numPtRows, numPtCols);
01454 
01455   //now we're ready to hand this stuff off to the FiniteElementData
01456   //instantiation.
01457 
01458   CHK_ERR( feData_->sumIntoMatrix(numPtRows,
01459                                   &rowNodeNumbers[0],
01460                                   &rowDofOffsets[0],
01461                                   &numColsPerRow[0],
01462                                   &allColNodeNumbers[0],
01463                                   &allColDofOffsets[0],
01464                                   &allValues[0]) );
01465 
01466   return(FEI_SUCCESS);
01467 }
01468 
01469 //------------------------------------------------------------------------------
01470 int FEDataFilter::giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
01471                                        int numPtCols, const int* ptCols,
01472                                        const double* const* values, int mode)
01473 {
01474   //This isn't going to be fast... I need to optimize the whole structure
01475   //of code that's associated with passing data to FiniteElementData.
01476 
01477   std::vector<int> rowNodeNumbers, rowDofOffsets, colNodeNumbers, colDofOffsets;
01478   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01479   int i;
01480 
01481   //First, we have to get nodeNumbers and dofOffsets for each of the
01482   //row-numbers and col-numbers.
01483 
01484   for(i=0; i<numPtRows; i++) {
01485     int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
01486     if (nodeNumber < 0) ERReturn(-1);
01487     int fieldID = problemStructure_->getAssociatedFieldID(ptRows[i]);
01488     if (fieldID < 0) ERReturn(-1);
01489     const NodeDescriptor* node = NULL;
01490     CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
01491     int firstEqn = node->getFieldEqnNumbers()[0];
01492 
01493     rowNodeNumbers.push_back(nodeNumber);
01494     rowDofOffsets.push_back(ptRows[i] - firstEqn);
01495   }
01496 
01497   for(i=0; i<numPtCols; i++) {
01498     int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
01499     if (nodeNumber < 0) ERReturn(-1);
01500     int fieldID = problemStructure_->getAssociatedFieldID(ptCols[i]);
01501     if (fieldID < 0) ERReturn(-1);
01502     const NodeDescriptor* node = NULL;
01503     CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
01504     int firstEqn = node->getFieldEqnNumbers()[0];
01505 
01506     colNodeNumbers.push_back(nodeNumber);
01507     colDofOffsets.push_back(ptCols[i] - firstEqn);
01508   }
01509 
01510   //now we have to flatten the colNodeNumbers and colDofOffsets out into
01511   //an array of length numPtRows*numPtCols, where the nodeNumbers and
01512   //dofOffsets are repeated 'numPtRows' times.
01513 
01514   int len = numPtRows*numPtCols;
01515   std::vector<int> allColNodeNumbers(len), allColDofOffsets(len);
01516   std::vector<double> allValues(len);
01517 
01518   int offset = 0;
01519   for(i=0; i<numPtRows; i++) {
01520     for(int j=0; j<numPtCols; j++) {
01521       allColNodeNumbers[offset] = colNodeNumbers[j];
01522       allColDofOffsets[offset] = colDofOffsets[j];
01523       allValues[offset++] = values[i][j];
01524     }
01525   }
01526 
01527   //while we're at it, let's make an array with numPtCols replicated in it
01528   //'numPtRows' times.
01529   std::vector<int> numColsPerRow(numPtRows, numPtCols);
01530 
01531   //now we're ready to hand this stuff off to the FiniteElementData
01532   //instantiation.
01533 
01534   CHK_ERR( feData_->sumIntoMatrix(numPtRows,
01535                                   &rowNodeNumbers[0],
01536                                   &rowDofOffsets[0],
01537                                   &numColsPerRow[0],
01538                                   &allColNodeNumbers[0],
01539                                   &allColDofOffsets[0],
01540                                   &allValues[0]) );
01541 
01542   return(FEI_SUCCESS);
01543 }
01544 
01545 //------------------------------------------------------------------------------
01546 int FEDataFilter::getFromMatrix(int numPtRows, const int* ptRows,
01547                             const int* rowColOffsets, const int* ptCols,
01548                             int numColsPerRow, double** values)
01549 {
01550   return(-1);
01551 
01552 }
01553 
01554 //------------------------------------------------------------------------------
01555 int FEDataFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
01556 {
01557   ERReturn(-1);
01558 }
01559 
01560 //------------------------------------------------------------------------------
01561 int FEDataFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
01562 {
01563   ERReturn(-1);
01564 }
01565 
01566 //------------------------------------------------------------------------------
01567 int FEDataFilter::giveToRHS(int num, const double* values,
01568                         const int* indices, int mode)
01569 {
01570   std::vector<int> workspace(num*2);
01571   int* rowNodeNumbers = &workspace[0];
01572   int* rowDofOffsets  = rowNodeNumbers+num;
01573   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01574 
01575   for(int i=0; i<num; ++i) {
01576     const NodeDescriptor* nodeptr = 0;
01577     int err = nodeDB.getNodeWithEqn(indices[i], nodeptr);
01578     if (err < 0) { 
01579         rowNodeNumbers[i] = -1;
01580         rowDofOffsets[i] = -1;
01581         continue;
01582     }
01583 
01584     rowNodeNumbers[i] = nodeptr->getNodeNumber();
01585 
01586     int firstEqn = nodeptr->getFieldEqnNumbers()[0];
01587 
01588     rowDofOffsets[i] = indices[i] - firstEqn;
01589   }
01590 
01591   if (mode == ASSEMBLE_SUM) {
01592     CHK_ERR( feData_->sumIntoRHSVector(num,
01593                                        rowNodeNumbers,
01594                                        rowDofOffsets,
01595                                        values) );
01596   }
01597   else {
01598     CHK_ERR( feData_->putIntoRHSVector(num,
01599                                        rowNodeNumbers,
01600                                        rowDofOffsets,
01601                                        values) );
01602   }
01603 
01604   return(FEI_SUCCESS);
01605 }
01606 
01607 //------------------------------------------------------------------------------
01608 int FEDataFilter::giveToLocalReducedRHS(int num, const double* values,
01609                                     const int* indices, int mode)
01610 {
01611   std::vector<int> rowNodeNumbers, rowDofOffsets;
01612   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01613   int i;
01614 
01615   for(i=0; i<num; i++) {
01616     int nodeNumber = problemStructure_->getAssociatedNodeNumber(indices[i]);
01617     if (nodeNumber < 0) ERReturn(-1);
01618 
01619     int fieldID = problemStructure_->getAssociatedFieldID(indices[i]);
01620     if (fieldID < 0) ERReturn(-1);
01621 
01622     const NodeDescriptor* node = NULL;
01623     CHK_ERR( nodeDB.getNodeWithNumber(nodeNumber, node) );
01624 
01625     int firstEqn = node->getFieldEqnNumbers()[0];
01626 
01627     rowNodeNumbers.push_back(nodeNumber);
01628     rowDofOffsets.push_back(indices[i] - firstEqn);
01629   }
01630 
01631   if (mode == ASSEMBLE_SUM) {
01632     CHK_ERR( feData_->sumIntoRHSVector(rowNodeNumbers.size(),
01633                                        &rowNodeNumbers[0],
01634                                        &rowDofOffsets[0],
01635                                        values) );
01636   }
01637   else {
01638     CHK_ERR( feData_->putIntoRHSVector(rowNodeNumbers.size(),
01639                                        &rowNodeNumbers[0],
01640                                        &rowDofOffsets[0],
01641                                        values) );
01642   }
01643 
01644   return(FEI_SUCCESS);
01645 }
01646 
01647 //------------------------------------------------------------------------------
01648 int FEDataFilter::getFromRHS(int num, double* values, const int* indices)
01649 {
01650   return(-1);
01651 }
01652 
01653 //------------------------------------------------------------------------------
01654 int FEDataFilter::getEqnSolnEntry(int eqnNumber, double& solnValue)
01655 {
01656   //This function's task is to retrieve the solution-value for a global
01657   //equation-number. eqnNumber may or may not be a slave-equation, and may or
01658   //may not be locally owned. If it is not locally owned, it should at least
01659   //be shared.
01660   //return 0 if the solution is successfully retrieved, otherwise return 1.
01661   //
01662 
01663   if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
01664     //Dig into the eqn-comm-mgr for the shared-remote solution value.
01665     CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
01666   }
01667   else {
01668     //It's local, simply get the solution from the assembled linear system.
01669     CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
01670   }
01671   return(0);
01672 }
01673 
01674 //------------------------------------------------------------------------------
01675 int FEDataFilter::getSharedRemoteSolnEntry(int eqnNumber, double& solnValue)
01676 {
01677   std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
01678   double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
01679 
01680   int index = snl_fei::binarySearch(eqnNumber, remoteEqnNumbers);
01681   if (index < 0) {
01682     FEI_CERR << "FEDataFilter::getSharedRemoteSolnEntry: ERROR, eqn "
01683          << eqnNumber << " not found." << FEI_ENDL;
01684     ERReturn(-1);
01685   }
01686   solnValue = remoteSoln[index];
01687   return(0);
01688 }
01689 
01690 //------------------------------------------------------------------------------
01691 int FEDataFilter::getReducedSolnEntry(int eqnNumber, double& solnValue)
01692 {
01693   //We may safely assume that this function is called with 'eqnNumber' that is
01694   //local in the underlying assembled linear system. i.e., it isn't a slave-
01695   //equation, it isn't remotely owned, etc.
01696   //
01697 
01698   int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqnNumber);
01699 
01700   //if nodeNumber < 0, it probably means we're trying to look up the
01701   //node for a lagrange-multiplier (which doesn't exist). In that
01702   //case, we're just going to ignore the request and return for now...
01703   if (nodeNumber < 0) {solnValue = -999.99; return(FEI_SUCCESS);}
01704 
01705   const NodeDescriptor* node = NULL;
01706   CHK_ERR( problemStructure_->getNodeDatabase().
01707            getNodeWithNumber(nodeNumber, node));
01708 
01709   int eqn = problemStructure_->translateFromReducedEqn(eqnNumber);
01710   int firstEqn = node->getFieldEqnNumbers()[0];
01711   int dofOffset = eqn - firstEqn;
01712 
01713   bool fetiHasNode = true;
01714   GlobalID nodeID = node->getGlobalNodeID();
01715   NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
01716   std::vector<GlobalID>& shNodeIDs = nodeCommMgr.getSharedNodeIDs();
01717   int shIndex = snl_fei::binarySearch(nodeID, shNodeIDs);
01718   if (shIndex >= 0) {
01719     if (!(problemStructure_->isInLocalElement(nodeNumber)) ) fetiHasNode = false;
01720   }
01721 
01722   if (fetiHasNode) {
01723     int err = feData_->getSolnEntry(nodeNumber, dofOffset, solnValue);
01724     if (err != 0) {
01725       FEI_CERR << "FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
01726            << " (nodeID " << node->getGlobalNodeID() << "), dofOffset "<<dofOffset
01727            << " couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
01728       ERReturn(-1);
01729     }
01730   }
01731 
01732   return(FEI_SUCCESS);
01733 }
01734 
01735 //------------------------------------------------------------------------------
01736 int FEDataFilter::unpackSolution()
01737 {
01738   //
01739   //This function should be called after the solver has returned,
01740   //and we know that there is a solution in the underlying vector.
01741   //This function ensures that any locally-owned shared solution values are
01742   //available on the sharing processors.
01743   //
01744   if (Filter::logStream() != NULL) {
01745     (*logStream())<< "#  entering unpackSolution, outputLevel: "
01746                      <<outputLevel_<<FEI_ENDL;
01747   }
01748 
01749   //what we need to do is as follows.
01750   //The eqn comm mgr has a list of what it calls 'recv eqns'. These are
01751   //equations that we own, for which we received contributions from other
01752   //processors. The solution values corresponding to these equations need
01753   //to be made available to those remote contributing processors.
01754 
01755    int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
01756    std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
01757 
01758    for(int i=0; i<numRecvEqns; i++) {
01759      int eqn = recvEqnNumbers[i];
01760 
01761      if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
01762        FEI_CERR << "FEDataFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
01763              << ") out of local range." << FEI_ENDL;
01764        MPI_Abort(comm_, -1);
01765      }
01766 
01767      double solnValue = 0.0;
01768 
01769      CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
01770 
01771      eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
01772    }
01773 
01774    eqnCommMgr_->exchangeSoln();
01775 
01776   debugOutput("#FEDataFilter leaving unpackSolution");
01777   return(FEI_SUCCESS);
01778 }
01779              
01780 //------------------------------------------------------------------------------
01781 void FEDataFilter::  setEqnCommMgr(EqnCommMgr* eqnCommMgr)
01782 {
01783   delete eqnCommMgr_;
01784   eqnCommMgr_ = eqnCommMgr;
01785 }
01786 
01787 //------------------------------------------------------------------------------
01788 int FEDataFilter::getBlockNodeSolution(GlobalID elemBlockID,  
01789                                        int numNodes, 
01790                                        const GlobalID *nodeIDs, 
01791                                        int *offsets,
01792                                        double *results)
01793 {        
01794    debugOutput("FEI: getBlockNodeSolution");
01795 
01796    int numActiveNodes = problemStructure_->getNumActiveNodes();
01797    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01798 
01799    if (numActiveNodes <= 0) return(0);
01800 
01801    int numSolnParams = 0;
01802 
01803    BlockDescriptor* block = NULL;
01804    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01805 
01806    //Traverse the node list, checking if nodes are associated with this block.
01807    //If so, put its 'answers' in the results list.
01808 
01809    int offset = 0;
01810    for(int i=0; i<numActiveNodes; i++) {
01811      NodeDescriptor* node_i = NULL;
01812      CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
01813 
01814       if (offset == numNodes) break;
01815 
01816       GlobalID nodeID = nodeIDs[offset];
01817 
01818       //first let's set the offset at which this node's solution coefs start.
01819       offsets[offset++] = numSolnParams;
01820 
01821       NodeDescriptor* node = NULL;
01822       int err = 0;
01823       //Obtain the NodeDescriptor of nodeID in the activeNodes list...
01824       //Don't call the getActiveNodeDesc_ID function unless we have to.
01825 
01826       if (nodeID == node_i->getGlobalNodeID()) {
01827         node = node_i;
01828       }
01829       else {
01830          err = nodeDB.getNodeWithID(nodeID, node);
01831       }
01832 
01833       //ok. If err is not 0, meaning nodeID is NOT in the
01834       //activeNodes list, then skip to the next loop iteration.
01835 
01836       if (err != 0) {
01837         continue;
01838       }
01839 
01840       int numFields = node->getNumFields();
01841       const int* fieldIDs = node->getFieldIDList();
01842 
01843       for(int j=0; j<numFields; j++) {
01844         if (block->containsField(fieldIDs[j])) {
01845           int size = problemStructure_->getFieldSize(fieldIDs[j]);
01846           if (size < 1) {
01847             continue;
01848           }
01849 
01850           int thisEqn = -1;
01851           node->getFieldEqnNumber(fieldIDs[j], thisEqn);
01852 
01853           double answer;
01854           for(int k=0; k<size; k++) {
01855             CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
01856               results[numSolnParams++] = answer;
01857           }
01858         }
01859       }//for(j<numFields)loop
01860    }
01861 
01862    offsets[numNodes] = numSolnParams;
01863 
01864    return(FEI_SUCCESS);
01865 }
01866 
01867 //------------------------------------------------------------------------------
01868 int FEDataFilter::getNodalSolution(int numNodes, 
01869                                    const GlobalID *nodeIDs, 
01870                                    int *offsets,
01871                                    double *results)
01872 {        
01873   debugOutput("FEI: getNodalSolution");
01874 
01875   int numActiveNodes = problemStructure_->getNumActiveNodes();
01876   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01877 
01878   if (numActiveNodes <= 0) return(0);
01879 
01880   int numSolnParams = 0;
01881 
01882   //Traverse the node list, checking if nodes are local.
01883   //If so, put 'answers' in the results list.
01884 
01885   int offset = 0;
01886   for(int i=0; i<numActiveNodes; i++) {
01887     NodeDescriptor* node_i = NULL;
01888     CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
01889 
01890     if (offset == numNodes) break;
01891 
01892     GlobalID nodeID = nodeIDs[offset];
01893 
01894     //first let's set the offset at which this node's solution coefs start.
01895     offsets[offset++] = numSolnParams;
01896 
01897     NodeDescriptor* node = NULL;
01898     int err = 0;
01899     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
01900     //Don't call the getNodeWithID function unless we have to.
01901 
01902     if (nodeID == node_i->getGlobalNodeID()) {
01903       node = node_i;
01904     }
01905     else {
01906       err = nodeDB.getNodeWithID(nodeID, node);
01907     }
01908 
01909     //ok. If err is not 0, meaning nodeID is NOT in the
01910     //activeNodes list, then skip to the next loop iteration.
01911 
01912     if (err != 0) {
01913       continue;
01914     }
01915 
01916     int numFields = node->getNumFields();
01917     const int* fieldIDs = node->getFieldIDList();
01918 
01919     for(int j=0; j<numFields; j++) {
01920       int size = problemStructure_->getFieldSize(fieldIDs[j]);
01921       if (size < 1) {
01922         continue;
01923       }
01924 
01925       int thisEqn = -1;
01926       node->getFieldEqnNumber(fieldIDs[j], thisEqn);
01927 
01928       double answer;
01929       for(int k=0; k<size; k++) {
01930         CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
01931           results[numSolnParams++] = answer;
01932       }
01933     }//for(j<numFields)loop
01934   }
01935 
01936   offsets[numNodes] = numSolnParams;
01937 
01938   return(FEI_SUCCESS);
01939 }
01940 
01941 //------------------------------------------------------------------------------
01942 int FEDataFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
01943                                         int fieldID,
01944                                         int numNodes, 
01945                                         const GlobalID *nodeIDs, 
01946                                         double *results)
01947 {
01948   debugOutput("FEI: getBlockFieldNodeSolution");
01949 
01950   int numActiveNodes = problemStructure_->getNumActiveNodes();
01951   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01952 
01953   if (numActiveNodes <= 0) return(0);
01954 
01955   BlockDescriptor* block = NULL;
01956   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01957 
01958   int fieldSize = problemStructure_->getFieldSize(fieldID);
01959   if (fieldSize <= 0) ERReturn(-1);
01960 
01961   if (!block->containsField(fieldID)) {
01962     FEI_CERR << "FEDataFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
01963          << " not contained in element-block " << static_cast<int>(elemBlockID) << FEI_ENDL;
01964     return(1);
01965   }
01966 
01967    //Traverse the node list, checking if nodes are associated with this block.
01968    //If so, put the answers in the results list.
01969 
01970    for(int i=0; i<numNodes; i++) {
01971      NodeDescriptor* node_i = NULL;
01972      CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
01973 
01974      GlobalID nodeID = nodeIDs[i];
01975 
01976      NodeDescriptor* node = NULL;
01977      int err = 0;
01978      //Obtain the NodeDescriptor of nodeID in the activeNodes list...
01979      //Don't call the getActiveNodeDesc_ID function unless we have to.
01980 
01981      if (nodeID == node_i->getGlobalNodeID()) {
01982        node = node_i;
01983      }
01984      else {
01985        err = nodeDB.getNodeWithID(nodeID, node);
01986      }
01987 
01988      //ok. If err is not 0, meaning nodeID is NOT in the
01989      //activeNodes list, then skip to the next loop iteration.
01990 
01991      if (err != 0) {
01992        continue;
01993      }
01994 
01995      int eqnNumber = -1;
01996      bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
01997      if (!hasField) continue;
01998 
01999      int offset = fieldSize*i;
02000      for(int j=0; j<fieldSize; j++) {
02001        double answer = 0.0;
02002        CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
02003        results[offset+j] = answer;
02004      }
02005    }
02006 
02007    return(FEI_SUCCESS);
02008 }
02009 
02010 //------------------------------------------------------------------------------
02011 int FEDataFilter::getNodalFieldSolution(int fieldID,
02012                                         int numNodes, 
02013                                         const GlobalID *nodeIDs, 
02014                                         double *results)
02015 {
02016   debugOutput("FEI: getNodalFieldSolution");
02017 
02018   int numActiveNodes = problemStructure_->getNumActiveNodes();
02019   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02020 
02021   if (numActiveNodes <= 0) return(0);
02022 
02023   if (problemStructure_->numSlaveEquations() != 0) {
02024     FEI_CERR << "FEDataFilter::getEqnSolnEntry ERROR FETI-support is not currently"
02025          << " compatible with the FEI's constraint reduction." << FEI_ENDL;
02026     ERReturn(-1);
02027   }
02028 
02029   int fieldSize = problemStructure_->getFieldSize(fieldID);
02030   if (fieldSize <= 0) {
02031     ERReturn(-1);
02032   }
02033 
02034   NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
02035 
02036   //Traverse the node list, checking if nodes have the specified field.
02037   //If so, put the answers in the results list.
02038 
02039   for(int i=0; i<numNodes; i++) {
02040     NodeDescriptor* node_i = NULL;
02041     CHK_ERR( nodeDB.getNodeAtIndex(i, node_i) );
02042 
02043     GlobalID nodeID = nodeIDs[i];
02044 
02045     NodeDescriptor* node = NULL;
02046     int err = 0;
02047     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
02048     //Don't call the getNodeWithID function unless we have to.
02049 
02050     if (nodeID == node_i->getGlobalNodeID()) {
02051       node = node_i;
02052     }
02053     else {
02054       err = nodeDB.getNodeWithID(nodeID, node);
02055     }
02056 
02057     //ok. If err is not 0, meaning nodeID is NOT in the
02058     //activeNodes list, then skip to the next loop iteration.
02059 
02060     if (err != 0) {
02061       continue;
02062     }
02063 
02064     int nodeNumber = node->getNodeNumber();
02065 
02066     int eqnNumber = -1;
02067     bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
02068 
02069     //If this node doesn't have the specified field, then skip to the
02070     //next loop iteration.
02071     if (!hasField) continue;
02072 
02073     std::vector<GlobalID>& shNodeIDs = nodeCommMgr.getSharedNodeIDs();
02074     int shIndex = snl_fei::binarySearch(nodeID, &shNodeIDs[0], shNodeIDs.size());
02075     if (shIndex > -1) {
02076       if (!(problemStructure_->isInLocalElement(nodeNumber))) continue;
02077     }
02078 
02079     int firstEqn = node->getFieldEqnNumbers()[0];
02080     int dofOffset = eqnNumber-firstEqn;
02081 
02082     int offset = fieldSize*i;
02083 
02084     for(int j=0; j<fieldSize; j++) {
02085       if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
02086         CHK_ERR( getSharedRemoteSolnEntry(eqnNumber+j, results[offset+j]) );
02087         continue;
02088       }
02089 
02090       err = feData_->getSolnEntry(nodeNumber, dofOffset+j, results[offset+j]);
02091       if (err != 0) {
02092         FEI_CERR << "FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
02093              << " (nodeID " << nodeID << "), dofOffset "<<dofOffset
02094              << " couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
02095         ERReturn(-1);
02096       }
02097     }
02098   }
02099 
02100   return(FEI_SUCCESS);
02101 }
02102 
02103 //------------------------------------------------------------------------------
02104 int FEDataFilter::putBlockNodeSolution(GlobalID elemBlockID,
02105                                    int numNodes, 
02106                                    const GlobalID *nodeIDs, 
02107                                    const int *offsets,
02108                                    const double *estimates) {
02109         
02110    debugOutput("FEI: putBlockNodeSolution");
02111 
02112    int numActiveNodes = problemStructure_->getNumActiveNodes();
02113 
02114    if (numActiveNodes <= 0) return(0);
02115 
02116    BlockDescriptor* block = NULL;
02117    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
02118 
02119    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02120 
02121    //traverse the node list, checking for nodes associated with this block
02122    //when an associated node is found, put its 'answers' into the linear system.
02123 
02124    for(int i=0; i<numNodes; i++) {
02125      NodeDescriptor* node = NULL;
02126      int err = nodeDB.getNodeWithID(nodeIDs[i], node);
02127 
02128       if (err != 0) continue;
02129    
02130       if (!node->containedInBlock(elemBlockID)) continue;
02131 
02132       if (node->getOwnerProc() != localRank_) continue;
02133 
02134       int numFields = node->getNumFields();
02135       const int* fieldIDs = node->getFieldIDList();
02136       const int* fieldEqnNumbers = node->getFieldEqnNumbers();
02137 
02138       if (fieldEqnNumbers[0] < localStartRow_ ||
02139           fieldEqnNumbers[0] > localEndRow_) continue;
02140 
02141       int offs = offsets[i];
02142 
02143       for(int j=0; j<numFields; j++) {
02144          int size = problemStructure_->getFieldSize(fieldIDs[j]);
02145 
02146          if (block->containsField(fieldIDs[j])) {
02147             for(int k=0; k<size; k++) {
02148                int reducedEqn;
02149                problemStructure_->
02150                  translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
02151 
02152 //                if (useLinSysCore_) {
02153 //                   CHK_ERR( lsc_->putInitialGuess(&reducedEqn,
02154 //                                             &estimates[offs+k], 1) );
02155 //                }
02156             }
02157          }
02158          offs += size;
02159       }//for(j<numFields)loop
02160    }
02161 
02162    return(FEI_SUCCESS);
02163 }
02164 
02165 //------------------------------------------------------------------------------
02166 int FEDataFilter::putBlockFieldNodeSolution(GlobalID elemBlockID, 
02167                                         int fieldID, 
02168                                         int numNodes, 
02169                                         const GlobalID *nodeIDs, 
02170                                         const double *estimates)
02171 {
02172    debugOutput("FEI: putBlockFieldNodeSolution");
02173 
02174    BlockDescriptor* block = NULL;
02175    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
02176    if (!block->containsField(fieldID)) return(1);
02177 
02178    int fieldSize = problemStructure_->getFieldSize(fieldID);
02179    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02180 
02181    //if we have a negative fieldID, we'll need a list of length numNodes,
02182    //in which to put nodeNumbers for passing to the solver... 
02183 
02184    std::vector<int> numbers(numNodes);
02185 
02186    //if we have a fieldID >= 0, then our numbers list will hold equation numbers
02187    //and we'll need fieldSize*numNodes of them.
02188 
02189    std::vector<double> data;
02190 
02191    if (fieldID >= 0) {
02192      if (fieldSize < 1) {
02193        FEI_CERR << "FEI Warning, putBlockFieldNodeSolution called for field "
02194             << fieldID<<", which has size "<<fieldSize<<FEI_ENDL;
02195        return(0);
02196      }
02197      try {
02198      numbers.resize(numNodes*fieldSize);
02199      data.resize(numNodes*fieldSize);
02200      }
02201      catch(std::runtime_error& exc) {
02202        FEI_CERR << exc.what()<<FEI_ENDL;
02203        ERReturn(-1);
02204      }
02205    }
02206 
02207    int count = 0;
02208 
02209    for(int i=0; i<numNodes; i++) {
02210      NodeDescriptor* node = NULL;
02211      CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
02212 
02213       if (fieldID < 0) numbers[count++] = node->getNodeNumber();
02214       else {
02215          int eqn = -1;
02216          if (node->getFieldEqnNumber(fieldID, eqn)) {
02217            if (eqn >= localStartRow_ && eqn <= localEndRow_) {
02218              for(int j=0; j<fieldSize; j++) { 
02219                data[count] = estimates[i*fieldSize + j];
02220                problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
02221              }
02222            }
02223          }
02224       }
02225    }
02226  
02227    if (fieldID < 0) {
02228      CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize, 
02229                                          numNodes, &numbers[0],
02230                                          estimates));
02231    }
02232 
02233    return(FEI_SUCCESS);
02234 }
02235 
02236 //------------------------------------------------------------------------------
02237 int FEDataFilter::getBlockElemSolution(GlobalID elemBlockID,
02238                                    int numElems, 
02239                                    const GlobalID *elemIDs,
02240                                    int& numElemDOFPerElement,
02241                                    double *results)
02242 {
02243 //
02244 //  return the elemental solution parameters associated with a
02245 //  particular block of elements
02246 //
02247    debugOutput("FEI: getBlockElemSolution");
02248 
02249    BlockDescriptor* block = NULL;
02250    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
02251 
02252    std::map<GlobalID,int>& elemIDList = problemStructure_->
02253                           getBlockConnectivity(elemBlockID).elemIDs;
02254 
02255    int len = block->getNumElements();
02256 
02257    //if the user is only asking for a subset of element-solutions, shrink len.
02258    if (len > numElems) len = numElems;
02259 
02260    numElemDOFPerElement = block->getNumElemDOFPerElement();
02261    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
02262    double answer;
02263 
02264 
02265    if (numElemDOFPerElement <= 0) return(0);
02266 
02267    std::map<GlobalID,int>::const_iterator
02268      elemid_end = elemIDList.end(),
02269      elemid_itr = elemIDList.begin();
02270 
02271    for(int i=0; i<len; i++) {
02272       int index = i;
02273 
02274       //if the user-supplied elemIDs are out of order, we need the index of
02275       //the location of this element.
02276       if (elemid_itr->first != elemIDs[i]) {
02277          index = elemid_itr->second;
02278       }
02279 
02280       if (index < 0) continue;
02281 
02282       int offset = i*numElemDOFPerElement;
02283 
02284       for(int j=0; j<numElemDOFPerElement; j++) {
02285          int eqn = elemDOFEqnNumbers[index] + j;
02286 
02287          CHK_ERR( getEqnSolnEntry(eqn, answer) )
02288 
02289          results[offset+j] = answer;
02290       }
02291 
02292       ++elemid_itr;
02293    }
02294 
02295    return(FEI_SUCCESS);
02296 }
02297 
02298 //------------------------------------------------------------------------------
02299 int FEDataFilter::putBlockElemSolution(GlobalID elemBlockID,
02300                                    int numElems,
02301                                    const GlobalID *elemIDs,
02302                                    int dofPerElem,
02303                                    const double *estimates)
02304 {
02305    debugOutput("FEI: putBlockElemSolution");
02306 
02307    BlockDescriptor* block = NULL;
02308    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
02309 
02310    std::map<GlobalID,int>& elemIDList = problemStructure_->
02311                           getBlockConnectivity(elemBlockID).elemIDs;
02312 
02313    int len = block->getNumElements();
02314    if (len > numElems) len = numElems;
02315 
02316    int DOFPerElement = block->getNumElemDOFPerElement();
02317    if (DOFPerElement != dofPerElem) {
02318      FEI_CERR << "FEI ERROR, putBlockElemSolution called with bad 'dofPerElem' ("
02319           <<dofPerElem<<"), block "<<elemBlockID<<" should have dofPerElem=="
02320           <<DOFPerElement<<FEI_ENDL;
02321      ERReturn(-1);
02322    }
02323 
02324    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
02325 
02326    if (DOFPerElement <= 0) return(0);
02327 
02328    std::map<GlobalID,int>::const_iterator
02329      elemid_end = elemIDList.end(),
02330      elemid_itr = elemIDList.begin();
02331 
02332    for(int i=0; i<len; i++) {
02333       int index = i;
02334       if (elemid_itr->first != elemIDs[i]) {
02335          index = elemid_itr->second;
02336       }
02337 
02338       if (index < 0) continue;
02339 
02340       for(int j=0; j<DOFPerElement; j++) {
02341          int reducedEqn;
02342          problemStructure_->
02343            translateToReducedEqn(elemDOFEqnNumbers[i] + j, reducedEqn);
02344 //         double soln = estimates[i*DOFPerElement + j];
02345 
02346 //       if (useLinSysCore_) {
02347 //            CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
02348 //          }
02349       }
02350 
02351       ++elemid_itr;
02352    }
02353 
02354    return(FEI_SUCCESS);
02355 }
02356 
02357 //------------------------------------------------------------------------------
02358 int FEDataFilter::getCRMultipliers(int numCRs,
02359                                    const int* CRIDs,
02360                                    double* multipliers)
02361 {
02362   for(int i=0; i<numCRs; i++) {
02363     //temporarily, FETI's getMultiplierSoln method isn't implemented.
02364     //CHK_ERR( feData_->getMultiplierSoln(CRIDs[i], multipliers[i]) );
02365     multipliers[i] = -999.99;
02366   }
02367 
02368   return(-1);
02369 }
02370 
02371 //------------------------------------------------------------------------------
02372 int FEDataFilter::putCRMultipliers(int numMultCRs,
02373                                const int* CRIDs,
02374                                const double *multEstimates)
02375 {
02376   debugOutput("FEI: putCRMultipliers");
02377 
02378   return(FEI_SUCCESS);
02379 }
02380 
02381 //------------------------------------------------------------------------------
02382 int FEDataFilter::putNodalFieldData(int fieldID,
02383                                 int numNodes,
02384                                 const GlobalID* nodeIDs,
02385                                 const double* nodeData)
02386 {
02387   debugOutput("FEI: putNodalFieldData");
02388 
02389   int fieldSize = problemStructure_->getFieldSize(fieldID);
02390   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02391 
02392   std::vector<int> nodeNumbers(numNodes);
02393 
02394   for(int i=0; i<numNodes; i++) {
02395     NodeDescriptor* node = NULL;
02396     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
02397 
02398     int nodeNumber = node->getNodeNumber();
02399     if (nodeNumber < 0) {
02400       GlobalID nodeID = nodeIDs[i];
02401       FEI_CERR << "FEDataFilter::putNodalFieldData ERROR, node with ID " 
02402            << static_cast<int>(nodeID) << " doesn't have an associated nodeNumber "
02403            << "assigned. putNodalFieldData shouldn't be called until after the "
02404            << "initComplete method has been called." << FEI_ENDL;
02405       ERReturn(-1);
02406     }
02407 
02408     nodeNumbers[i] = nodeNumber;
02409   }
02410 
02411   CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
02412                                       numNodes, &nodeNumbers[0],
02413                                       nodeData) );
02414 
02415   return(0);
02416 }
02417 
02418 //------------------------------------------------------------------------------
02419 int FEDataFilter::putNodalFieldSolution(int fieldID,
02420                                 int numNodes,
02421                                 const GlobalID* nodeIDs,
02422                                 const double* nodeData)
02423 {
02424   debugOutput("FEI: putNodalFieldSolution");
02425 
02426   if (fieldID < 0) {
02427     return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
02428   }
02429 
02430   int fieldSize = problemStructure_->getFieldSize(fieldID);
02431   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02432 
02433   std::vector<int> eqnNumbers(fieldSize);
02434 
02435   for(int i=0; i<numNodes; i++) {
02436     NodeDescriptor* node = NULL;
02437     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
02438 
02439     int eqn = -1;
02440     bool hasField = node->getFieldEqnNumber(fieldID, eqn);
02441     if (!hasField) continue;
02442 
02443   }
02444 
02445   return(0);
02446 }
02447 
02448 //------------------------------------------------------------------------------
02449 int FEDataFilter::assembleEqns(int numRows, 
02450                                int numCols,
02451                                const int* rowNumbers,
02452                                const int* colIndices,
02453                                const double* const* coefs,
02454                                bool structurallySymmetric,
02455                                int mode)
02456 {
02457   if (numRows == 0) return(FEI_SUCCESS);
02458 
02459   const int* indPtr = colIndices;
02460   for(int i=0; i<numRows; i++) {
02461     int row = rowNumbers[i];
02462 
02463     const double* coefPtr = coefs[i];
02464 
02465     CHK_ERR(giveToMatrix(1, &row, numCols, indPtr, &coefPtr, mode));
02466 
02467     if (!structurallySymmetric) indPtr += numCols;
02468   }
02469 
02470   return(FEI_SUCCESS);
02471 }
02472 
02473 //------------------------------------------------------------------------------
02474 int FEDataFilter::assembleRHS(int numValues,
02475                               const int* indices,
02476                               const double* coefs,
02477                               int mode) {
02478 //
02479 //This function hands the data off to the routine that finally
02480 //sticks it into the RHS vector.
02481 //
02482   if (problemStructure_->numSlaveEquations() == 0) {
02483     CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
02484     return(FEI_SUCCESS);
02485   }
02486 
02487   for(int i = 0; i < numValues; i++) {
02488     int eqn = indices[i];
02489 
02490     CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
02491   }
02492 
02493   return(FEI_SUCCESS);
02494 }
02495 
02496 //==============================================================================
02497 void FEDataFilter::debugOutput(const char* mesg)
02498 {
02499   if (Filter::logStream() != NULL) {
02500     (*logStream()) << mesg << FEI_ENDL;
02501    }
02502 }

Generated on Wed May 12 21:30:41 2010 for FEI by  doxygen 1.4.7