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::findNode(GlobalID nodeID) const {
00136 //
00137 //This function returns a NodeDescriptor ptr, may return NULL.
00138 //
00139   const NodeDescriptor* node = NULL;
00140   problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
00141   return node;
00142 }
00143 
00144 //------------------------------------------------------------------------------
00145 const NodeDescriptor& Filter::findNodeDescriptor(GlobalID nodeID) const {
00146 //
00147 //This function returns a NodeDescriptor reference if nodeID is an active node.
00148 //
00149   const NodeDescriptor* node = NULL;
00150   int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
00151 
00152   if (err != 0) {
00153     fei::console_out() << "ERROR, Filter::findNodeDescriptor unable to find node "
00154    << static_cast<int>(nodeID) << FEI_ENDL;
00155     std::abort();
00156   }
00157 
00158   return( *node );
00159 }
00160 //------------------------------------------------------------------------------
00161 int Filter::calculateResidualNorms(int whichNorm, int numFields,
00162            int* fieldIDs, double* norms,
00163            std::vector<double>& residValues)
00164 {
00165   std::vector<double> normsArray(numFields, 0.0);
00166 
00167   std::fill(fieldIDs, fieldIDs+numFields, -999);
00168 
00169   std::vector<double> tmpNorms(numFields);
00170   double* tmpNormsPtr = &tmpNorms[0];
00171 
00172   double* residPtr = &(residValues[0]);
00173 
00174   const std::vector<int>& pfieldIDs = problemStructure_->getFieldIDs();
00175   int numDBFields = pfieldIDs.size();
00176   std::vector<int>::const_iterator
00177     f_iter = pfieldIDs.begin(),
00178     f_end = pfieldIDs.end();
00179 
00180   int DBFieldSize = 0;
00181 
00182   int offset = 0;
00183   for(; f_iter != f_end; ++f_iter) {
00184 
00185     if (offset == 0) DBFieldSize = problemStructure_->getFieldSize(*f_iter);
00186 
00187     if (*f_iter > -1) {
00188       if (offset < numFields) {
00189         fieldIDs[offset] = *f_iter;
00190         tmpNormsPtr[offset++] = 0.0;
00191       }
00192     }
00193   }
00194 
00195   int reducedStartRow = problemStructure_->getFirstReducedEqn();
00196   int reducedEndRow   = problemStructure_->getLastReducedEqn();
00197 
00198   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00199   int numNodes = nodeDB.getNumNodeDescriptors();
00200 
00201   bool haveSlaves = problemStructure_->numSlaveEquations() > 0;
00202 
00203   for(int i=0; i<numNodes; i++) {
00204     const NodeDescriptor* node = NULL;
00205     nodeDB.getNodeAtIndex(i, node);
00206 
00207     if (node==NULL || node->getOwnerProc() != localRank_) continue;
00208 
00209     const int* fieldIDList = node->getFieldIDList();
00210     const int* fieldEqnNums = node->getFieldEqnNumbers();
00211     int numNodeFields = node->getNumFields();
00212 
00213     for(int j=0; j<numNodeFields; j++) {
00214       int fIndex = 0;
00215       int fSize = DBFieldSize;
00216 
00217       if (numDBFields > 1) {
00218         fIndex = fei::binarySearch(fieldIDList[j], fieldIDs, numFields);
00219         if (fIndex < 0) return(-1);
00220         fSize = problemStructure_->getFieldSize(fieldIDList[j]);
00221         if (fSize < 0) return(-1);
00222       }
00223 
00224       for(int k=0; k<fSize; k++) {
00225         int eqn = fieldEqnNums[j]+k;
00226 
00227         if (haveSlaves) {
00228           if (problemStructure_->isSlaveEqn(eqn)) continue;
00229           int reducedEqn;
00230           problemStructure_->translateToReducedEqn(eqn, reducedEqn);
00231 
00232           if (reducedEqn < reducedStartRow || reducedEqn > reducedEndRow) {
00233             continue;
00234           }
00235           eqn = reducedEqn;
00236         }
00237 
00238         double rval = residPtr[eqn - reducedStartRow];
00239 
00240         switch(whichNorm) {
00241           case 0:
00242             if (tmpNormsPtr[fIndex] < std::abs(rval)) tmpNormsPtr[fIndex] = rval;
00243             break;
00244           case 1:
00245             tmpNormsPtr[fIndex] += std::abs(rval);
00246             break;
00247           case 2:
00248             tmpNormsPtr[fIndex] += rval*rval;
00249             break;
00250           default:
00251             FEI_COUT << "Filter::residualNorm: ERROR, whichNorm="<<whichNorm
00252               << " not recognized." << FEI_ENDL;
00253             return(-1);
00254         }
00255       }
00256     }
00257   }
00258 
00259   //so at this point we have the local norms. We now need to perform the
00260   //global max or sum, depending on which norm it is.
00261 
00262   MPI_Comm comm = problemStructure_->getCommunicator();
00263 
00264   if (whichNorm != 0) {
00265     CHK_ERR( fei::GlobalSum(comm, tmpNorms, normsArray) );
00266   }
00267   else {
00268     CHK_ERR( fei::GlobalMax(comm, tmpNorms, normsArray) );
00269   }
00270 
00271   for(int i=0; i<numFields; ++i) {
00272     norms[i] = normsArray[i];
00273   }
00274 
00275   if (whichNorm == 2) {
00276     for(int i=0; i<numFields; ++i) norms[i] = std::sqrt(norms[i]);
00277   }
00278 
00279   return(0);
00280 }
00281 
00282 //------------------------------------------------------------------------------
00283 int Filter::parameters(int numParams, const char *const* paramStrings)
00284 {
00285   if (numParams == 0 || paramStrings == NULL) return(0);
00286 
00287   const char* param = snl_fei::getParam("outputLevel",numParams,paramStrings);
00288 
00289   if ( param != NULL){
00290     std::string str(&(param[11]));
00291     FEI_ISTRINGSTREAM isstr(str);
00292     isstr >> outputLevel_;
00293   }
00294 
00295   return(0);
00296 }
00297 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends