FEI Version of the Day
fei_Filter.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include "fei_sstream.hpp"
00010 
00011 #include "fei_CommUtils.hpp"
00012 #include "fei_TemplateUtils.hpp"
00013 
00014 #include "fei_defs.h"
00015 #include "fei_NodeDescriptor.hpp"
00016 #include "fei_NodeDatabase.hpp"
00017 #include "fei_BlockDescriptor.hpp"
00018 #include "SNL_FEI_Structure.hpp"
00019 #include "snl_fei_Utils.hpp"
00020 #include "fei_Filter.hpp"
00021 
00022 #include <cmath>
00023 #include <algorithm>
00024 
00025 #undef fei_file
00026 #define fei_file "fei_Filter.cpp"
00027 
00028 #include "fei_ErrMacros.hpp"
00029 
00030 //------------------------------------------------------------------------------
00031 Filter::Filter(SNL_FEI_Structure* probStruct)
00032   : problemStructure_(probStruct),
00033     logInput_(false),
00034     logInputStream_(NULL),
00035     outputLevel_(0)
00036 {
00037 }
00038 
00039 //------------------------------------------------------------------------------
00040 Filter::~Filter()
00041 {}
00042 
00043 //------------------------------------------------------------------------------
00044 void Filter::setLogStream(FEI_OSTREAM* logstrm)
00045 {
00046   logInputStream_ = logstrm;
00047 }
00048 
00049 //------------------------------------------------------------------------------
00050 FEI_OSTREAM* Filter::logStream()
00051 {
00052   return( logInputStream_ );
00053 }
00054 
00055 //------------------------------------------------------------------------------
00056 void Filter::copyStiffness(const double* const* elemStiff,
00057          int numRows, int elemFormat,
00058          double** copy)
00059 {
00060   //
00061   //Unpack the element stiffness array in elemStiff into a full dense
00062   //'copy'.
00063   //
00064   int i, j;
00065   const double* elStiff_i = NULL;
00066 
00067   switch (elemFormat) {
00068 
00069   case FEI_DENSE_ROW:
00070     for (i = 0; i < numRows; i++) {
00071       elStiff_i = elemStiff[i];
00072       for (j = 0; j < numRows; j++) {
00073   copy[i][j] = elStiff_i[j];
00074       }
00075     }
00076     break;
00077 
00078   case FEI_UPPER_SYMM_ROW:
00079     for (i = 0; i < numRows; i++) {
00080       elStiff_i = elemStiff[i];
00081       int jcol=0;
00082       for (j = i; j < numRows; j++) {
00083   copy[i][j] = elStiff_i[jcol++];
00084   copy[j][i] = copy[i][j];
00085       }
00086     }
00087     break;
00088 
00089   case FEI_LOWER_SYMM_ROW:
00090     for (i = 0; i < numRows; i++) {
00091       elStiff_i = elemStiff[i];
00092       for (j = 0; j <=i; j++) {
00093   copy[i][j] = elStiff_i[j];
00094   copy[j][i] = copy[i][j];
00095       }
00096     }
00097     break;
00098 
00099   case FEI_DENSE_COL:
00100     for (i = 0; i < numRows; i++) {
00101       elStiff_i = elemStiff[i];
00102       for (j = 0; j < numRows; j++) {
00103   copy[j][i] = elStiff_i[j];
00104       }
00105     }
00106     break;
00107 
00108   case FEI_UPPER_SYMM_COL:
00109     for (i = 0; i < numRows; i++) {
00110       elStiff_i = elemStiff[i];
00111       for (j = 0; j <= i; j++) {
00112   copy[i][j] = elStiff_i[j];
00113   copy[j][i] = copy[i][j];
00114       }
00115     }
00116     break;
00117 
00118   case FEI_LOWER_SYMM_COL:
00119     for (i = 0; i < numRows; i++) {
00120       elStiff_i = elemStiff[i];
00121       int jcol=0;
00122       for (j = i; j < numRows; j++) {
00123   copy[i][j] = elStiff_i[jcol++];
00124   copy[j][i] = copy[i][j];
00125       }
00126     }
00127     break;
00128 
00129   default:
00130     throw std::runtime_error("copyStiffness ERROR, unrecognized elem-format");
00131   }
00132 }
00133 
00134 //------------------------------------------------------------------------------
00135 const NodeDescriptor& Filter::findNodeDescriptor(GlobalID nodeID) const {
00136 //
00137 //This function returns a NodeDescriptor reference if nodeID is an active node.
00138 //
00139   const NodeDescriptor* node = NULL;
00140   int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
00141 
00142   if (err != 0) {
00143     FEI_CERR << "ERROR, Filter::findNodeDescriptor unable to find node "
00144    << static_cast<int>(nodeID) << FEI_ENDL;
00145     std::abort();
00146   }
00147 
00148   return( *node );
00149 }
00150 
00151 //------------------------------------------------------------------------------
00152 int Filter::calculateResidualNorms(int whichNorm, int numFields,
00153            int* fieldIDs, double* norms,
00154            std::vector<double>& residValues)
00155 {
00156   std::vector<double> normsArray(numFields, 0.0);
00157 
00158   std::fill(fieldIDs, fieldIDs+numFields, -999);
00159 
00160   std::vector<double> tmpNorms(numFields);
00161   double* tmpNormsPtr = &tmpNorms[0];
00162 
00163   double* residPtr = &(residValues[0]);
00164 
00165   std::map<int,int>& fieldDB = problemStructure_->getFieldDatabase();
00166   int numDBFields = fieldDB.size();
00167   std::map<int,int>::const_iterator
00168     db_iter = fieldDB.begin(),
00169     db_end = fieldDB.end();
00170 
00171   int DBFieldSize = 0;
00172 
00173   int offset = 0;
00174   for(; db_iter != db_end; ++db_iter) {
00175     const std::pair<const int,int>& dbpair = *db_iter;
00176 
00177     if (offset == 0) DBFieldSize = dbpair.second;
00178 
00179     if (dbpair.first > -1) {
00180       if (offset < numFields) {
00181   fieldIDs[offset] = dbpair.first;
00182   tmpNormsPtr[offset++] = 0.0;
00183       }
00184     }
00185   }
00186 
00187   int reducedStartRow = problemStructure_->getFirstReducedEqn();
00188   int reducedEndRow   = problemStructure_->getLastReducedEqn();
00189 
00190   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00191   int numNodes = nodeDB.getNumNodeDescriptors();
00192 
00193   bool haveSlaves = problemStructure_->numSlaveEquations() > 0;
00194 
00195   for(int i=0; i<numNodes; i++) {
00196     const NodeDescriptor* node = NULL;
00197     CHK_ERR( nodeDB.getNodeAtIndex(i, node) );
00198 
00199     if (node->getOwnerProc() != localRank_) continue;
00200 
00201     const int* fieldIDList = node->getFieldIDList();
00202     const int* fieldEqnNums = node->getFieldEqnNumbers();
00203     int numNodeFields = node->getNumFields();
00204 
00205     for(int j=0; j<numNodeFields; j++) {
00206       int fIndex = 0;
00207       int fSize = DBFieldSize;
00208 
00209       if (numDBFields > 1) {
00210   fIndex = fei::binarySearch(fieldIDList[j], fieldIDs, numFields);
00211   if (fIndex < 0) return(-1);
00212   fSize = problemStructure_->getFieldSize(fieldIDList[j]);
00213   if (fSize < 0) return(-1);
00214       }
00215 
00216       for(int k=0; k<fSize; k++) {
00217   int eqn = fieldEqnNums[j]+k;
00218 
00219   if (haveSlaves) {
00220     if (problemStructure_->isSlaveEqn(eqn)) continue;
00221     int reducedEqn;
00222     problemStructure_->translateToReducedEqn(eqn, reducedEqn);
00223 
00224     if (reducedEqn < reducedStartRow || reducedEqn > reducedEndRow) {
00225       continue;
00226     }
00227     eqn = reducedEqn;
00228   }
00229 
00230   double rval = residPtr[eqn - reducedStartRow];
00231 
00232   switch(whichNorm) {
00233   case 0:
00234     if (tmpNormsPtr[fIndex] < std::abs(rval)) tmpNormsPtr[fIndex] = rval;
00235     break;
00236   case 1:
00237     tmpNormsPtr[fIndex] += std::abs(rval);
00238     break;
00239   case 2:
00240     tmpNormsPtr[fIndex] += rval*rval;
00241     break;
00242   default:
00243     FEI_COUT << "Filter::residualNorm: ERROR, whichNorm="<<whichNorm
00244          << " not recognized." << FEI_ENDL;
00245     return(-1);
00246   }
00247       }
00248     }
00249   }
00250 
00251   //so at this point we have the local norms. We now need to perform the
00252   //global max or sum, depending on which norm it is.
00253 
00254   MPI_Comm comm = problemStructure_->getCommunicator();
00255 
00256   if (whichNorm != 0) {
00257     CHK_ERR( fei::GlobalSum(comm, tmpNorms, normsArray) );
00258   }
00259   else {
00260     CHK_ERR( fei::GlobalMax(comm, tmpNorms, normsArray) );
00261   }
00262 
00263   for(int i=0; i<numFields; ++i) {
00264     norms[i] = normsArray[i];
00265   }
00266 
00267   if (whichNorm == 2) {
00268     for(int i=0; i<numFields; ++i) norms[i] = std::sqrt(norms[i]);
00269   }
00270 
00271   return(0);
00272 }
00273 
00274 //------------------------------------------------------------------------------
00275 int Filter::parameters(int numParams, const char *const* paramStrings)
00276 {
00277   if (numParams == 0 || paramStrings == NULL) return(0);
00278 
00279   const char* param = snl_fei::getParam("outputLevel",numParams,paramStrings);
00280 
00281   if ( param != NULL){
00282     std::string str(&(param[11]));
00283     FEI_ISTRINGSTREAM isstr(str);
00284     isstr >> outputLevel_;
00285   }
00286 
00287   return(0);
00288 }
00289 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends