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