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