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