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   std::map<int,int>& fieldDB = problemStructure_->getFieldDatabase();
00175   int numDBFields = fieldDB.size();
00176   std::map<int,int>::const_iterator
00177     db_iter = fieldDB.begin(),
00178     db_end = fieldDB.end();
00179 
00180   int DBFieldSize = 0;
00181 
00182   int offset = 0;
00183   for(; db_iter != db_end; ++db_iter) {
00184     const std::pair<const int,int>& dbpair = *db_iter;
00185 
00186     if (offset == 0) DBFieldSize = dbpair.second;
00187 
00188     if (dbpair.first > -1) {
00189       if (offset < numFields) {
00190   fieldIDs[offset] = dbpair.first;
00191   tmpNormsPtr[offset++] = 0.0;
00192       }
00193     }
00194   }
00195 
00196   int reducedStartRow = problemStructure_->getFirstReducedEqn();
00197   int reducedEndRow   = problemStructure_->getLastReducedEqn();
00198 
00199   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00200   int numNodes = nodeDB.getNumNodeDescriptors();
00201 
00202   bool haveSlaves = problemStructure_->numSlaveEquations() > 0;
00203 
00204   for(int i=0; i<numNodes; i++) {
00205     const NodeDescriptor* node = NULL;
00206     nodeDB.getNodeAtIndex(i, node);
00207 
00208     if (node==NULL || node->getOwnerProc() != localRank_) continue;
00209 
00210     const int* fieldIDList = node->getFieldIDList();
00211     const int* fieldEqnNums = node->getFieldEqnNumbers();
00212     int numNodeFields = node->getNumFields();
00213 
00214     for(int j=0; j<numNodeFields; j++) {
00215       int fIndex = 0;
00216       int fSize = DBFieldSize;
00217 
00218       if (numDBFields > 1) {
00219         fIndex = fei::binarySearch(fieldIDList[j], fieldIDs, numFields);
00220         if (fIndex < 0) return(-1);
00221         fSize = problemStructure_->getFieldSize(fieldIDList[j]);
00222         if (fSize < 0) return(-1);
00223       }
00224 
00225       for(int k=0; k<fSize; k++) {
00226         int eqn = fieldEqnNums[j]+k;
00227 
00228         if (haveSlaves) {
00229           if (problemStructure_->isSlaveEqn(eqn)) continue;
00230           int reducedEqn;
00231           problemStructure_->translateToReducedEqn(eqn, reducedEqn);
00232 
00233           if (reducedEqn < reducedStartRow || reducedEqn > reducedEndRow) {
00234             continue;
00235           }
00236           eqn = reducedEqn;
00237         }
00238 
00239         double rval = residPtr[eqn - reducedStartRow];
00240 
00241         switch(whichNorm) {
00242           case 0:
00243             if (tmpNormsPtr[fIndex] < std::abs(rval)) tmpNormsPtr[fIndex] = rval;
00244             break;
00245           case 1:
00246             tmpNormsPtr[fIndex] += std::abs(rval);
00247             break;
00248           case 2:
00249             tmpNormsPtr[fIndex] += rval*rval;
00250             break;
00251           default:
00252             FEI_COUT << "Filter::residualNorm: ERROR, whichNorm="<<whichNorm
00253               << " not recognized." << FEI_ENDL;
00254             return(-1);
00255         }
00256       }
00257     }
00258   }
00259 
00260   //so at this point we have the local norms. We now need to perform the
00261   //global max or sum, depending on which norm it is.
00262 
00263   MPI_Comm comm = problemStructure_->getCommunicator();
00264 
00265   if (whichNorm != 0) {
00266     CHK_ERR( fei::GlobalSum(comm, tmpNorms, normsArray) );
00267   }
00268   else {
00269     CHK_ERR( fei::GlobalMax(comm, tmpNorms, normsArray) );
00270   }
00271 
00272   for(int i=0; i<numFields; ++i) {
00273     norms[i] = normsArray[i];
00274   }
00275 
00276   if (whichNorm == 2) {
00277     for(int i=0; i<numFields; ++i) norms[i] = std::sqrt(norms[i]);
00278   }
00279 
00280   return(0);
00281 }
00282 
00283 //------------------------------------------------------------------------------
00284 int Filter::parameters(int numParams, const char *const* paramStrings)
00285 {
00286   if (numParams == 0 || paramStrings == NULL) return(0);
00287 
00288   const char* param = snl_fei::getParam("outputLevel",numParams,paramStrings);
00289 
00290   if ( param != NULL){
00291     std::string str(&(param[11]));
00292     FEI_ISTRINGSTREAM isstr(str);
00293     isstr >> outputLevel_;
00294   }
00295 
00296   return(0);
00297 }
00298 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends