FEI Version of the Day
fei_Filter.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_sstream.hpp"
00045 
00046 #include "fei_CommUtils.hpp"
00047 #include "fei_TemplateUtils.hpp"
00048 
00049 #include "fei_defs.h"
00050 #include "fei_NodeDescriptor.hpp"
00051 #include "fei_NodeDatabase.hpp"
00052 #include "fei_BlockDescriptor.hpp"
00053 #include "SNL_FEI_Structure.hpp"
00054 #include "snl_fei_Utils.hpp"
00055 #include "fei_Filter.hpp"
00056 
00057 #include <cmath>
00058 #include <algorithm>
00059 
00060 #undef fei_file
00061 #define fei_file "fei_Filter.cpp"
00062 
00063 #include "fei_ErrMacros.hpp"
00064 
00065 //------------------------------------------------------------------------------
00066 Filter::Filter(SNL_FEI_Structure* probStruct)
00067   : problemStructure_(probStruct),
00068     logInput_(false),
00069     logInputStream_(NULL),
00070     outputLevel_(0)
00071 {
00072 }
00073 
00074 //------------------------------------------------------------------------------
00075 Filter::~Filter()
00076 {}
00077 
00078 //------------------------------------------------------------------------------
00079 void Filter::setLogStream(FEI_OSTREAM* logstrm)
00080 {
00081   logInputStream_ = logstrm;
00082 }
00083 
00084 //------------------------------------------------------------------------------
00085 FEI_OSTREAM* Filter::logStream()
00086 {
00087   return( logInputStream_ );
00088 }
00089 
00090 //------------------------------------------------------------------------------
00091 void Filter::copyStiffness(const double* const* elemStiff,
00092          int numRows, int elemFormat,
00093          double** copy)
00094 {
00095   //
00096   //Unpack the element stiffness array in elemStiff into a full dense
00097   //'copy'.
00098   //
00099   int i, j;
00100   const double* elStiff_i = NULL;
00101 
00102   switch (elemFormat) {
00103 
00104   case FEI_DENSE_ROW:
00105     for (i = 0; i < numRows; i++) {
00106       elStiff_i = elemStiff[i];
00107       for (j = 0; j < numRows; j++) {
00108   copy[i][j] = elStiff_i[j];
00109       }
00110     }
00111     break;
00112 
00113   case FEI_UPPER_SYMM_ROW:
00114     for (i = 0; i < numRows; i++) {
00115       elStiff_i = elemStiff[i];
00116       int jcol=0;
00117       for (j = i; j < numRows; j++) {
00118   copy[i][j] = elStiff_i[jcol++];
00119   copy[j][i] = copy[i][j];
00120       }
00121     }
00122     break;
00123 
00124   case FEI_LOWER_SYMM_ROW:
00125     for (i = 0; i < numRows; i++) {
00126       elStiff_i = elemStiff[i];
00127       for (j = 0; j <=i; j++) {
00128   copy[i][j] = elStiff_i[j];
00129   copy[j][i] = copy[i][j];
00130       }
00131     }
00132     break;
00133 
00134   case FEI_DENSE_COL:
00135     for (i = 0; i < numRows; i++) {
00136       elStiff_i = elemStiff[i];
00137       for (j = 0; j < numRows; j++) {
00138   copy[j][i] = elStiff_i[j];
00139       }
00140     }
00141     break;
00142 
00143   case FEI_UPPER_SYMM_COL:
00144     for (i = 0; i < numRows; i++) {
00145       elStiff_i = elemStiff[i];
00146       for (j = 0; j <= i; j++) {
00147   copy[i][j] = elStiff_i[j];
00148   copy[j][i] = copy[i][j];
00149       }
00150     }
00151     break;
00152 
00153   case FEI_LOWER_SYMM_COL:
00154     for (i = 0; i < numRows; i++) {
00155       elStiff_i = elemStiff[i];
00156       int jcol=0;
00157       for (j = i; j < numRows; j++) {
00158   copy[i][j] = elStiff_i[jcol++];
00159   copy[j][i] = copy[i][j];
00160       }
00161     }
00162     break;
00163 
00164   default:
00165     throw std::runtime_error("copyStiffness ERROR, unrecognized elem-format");
00166   }
00167 }
00168 
00169 //------------------------------------------------------------------------------
00170 const NodeDescriptor* Filter::findNode(GlobalID nodeID) const {
00171 //
00172 //This function returns a NodeDescriptor ptr, may return NULL.
00173 //
00174   const NodeDescriptor* node = NULL;
00175   problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
00176   return node;
00177 }
00178 
00179 //------------------------------------------------------------------------------
00180 const NodeDescriptor& Filter::findNodeDescriptor(GlobalID nodeID) const {
00181 //
00182 //This function returns a NodeDescriptor reference if nodeID is an active node.
00183 //
00184   const NodeDescriptor* node = NULL;
00185   int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
00186 
00187   if (err != 0) {
00188     fei::console_out() << "ERROR, Filter::findNodeDescriptor unable to find node "
00189    << static_cast<int>(nodeID) << FEI_ENDL;
00190     std::abort();
00191   }
00192 
00193   return( *node );
00194 }
00195 //------------------------------------------------------------------------------
00196 int Filter::calculateResidualNorms(int whichNorm, int numFields,
00197            int* fieldIDs, double* norms,
00198            std::vector<double>& residValues)
00199 {
00200   std::vector<double> normsArray(numFields, 0.0);
00201 
00202   std::fill(fieldIDs, fieldIDs+numFields, -999);
00203 
00204   std::vector<double> tmpNorms(numFields);
00205   double* tmpNormsPtr = &tmpNorms[0];
00206 
00207   double* residPtr = &(residValues[0]);
00208 
00209   const std::vector<int>& pfieldIDs = problemStructure_->getFieldIDs();
00210   int numDBFields = pfieldIDs.size();
00211   std::vector<int>::const_iterator
00212     f_iter = pfieldIDs.begin(),
00213     f_end = pfieldIDs.end();
00214 
00215   int DBFieldSize = 0;
00216 
00217   int offset = 0;
00218   for(; f_iter != f_end; ++f_iter) {
00219 
00220     if (offset == 0) DBFieldSize = problemStructure_->getFieldSize(*f_iter);
00221 
00222     if (*f_iter > -1) {
00223       if (offset < numFields) {
00224         fieldIDs[offset] = *f_iter;
00225         tmpNormsPtr[offset++] = 0.0;
00226       }
00227     }
00228   }
00229 
00230   int reducedStartRow = problemStructure_->getFirstReducedEqn();
00231   int reducedEndRow   = problemStructure_->getLastReducedEqn();
00232 
00233   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
00234   int numNodes = nodeDB.getNumNodeDescriptors();
00235 
00236   bool haveSlaves = problemStructure_->numSlaveEquations() > 0;
00237 
00238   for(int i=0; i<numNodes; i++) {
00239     const NodeDescriptor* node = NULL;
00240     nodeDB.getNodeAtIndex(i, node);
00241 
00242     if (node==NULL || node->getOwnerProc() != localRank_) continue;
00243 
00244     const int* fieldIDList = node->getFieldIDList();
00245     const int* fieldEqnNums = node->getFieldEqnNumbers();
00246     int numNodeFields = node->getNumFields();
00247 
00248     for(int j=0; j<numNodeFields; j++) {
00249       int fIndex = 0;
00250       int fSize = DBFieldSize;
00251 
00252       if (numDBFields > 1) {
00253         fIndex = fei::binarySearch(fieldIDList[j], fieldIDs, numFields);
00254         if (fIndex < 0) return(-1);
00255         fSize = problemStructure_->getFieldSize(fieldIDList[j]);
00256         if (fSize < 0) return(-1);
00257       }
00258 
00259       for(int k=0; k<fSize; k++) {
00260         int eqn = fieldEqnNums[j]+k;
00261 
00262         if (haveSlaves) {
00263           if (problemStructure_->isSlaveEqn(eqn)) continue;
00264           int reducedEqn;
00265           problemStructure_->translateToReducedEqn(eqn, reducedEqn);
00266 
00267           if (reducedEqn < reducedStartRow || reducedEqn > reducedEndRow) {
00268             continue;
00269           }
00270           eqn = reducedEqn;
00271         }
00272 
00273         double rval = residPtr[eqn - reducedStartRow];
00274 
00275         switch(whichNorm) {
00276           case 0:
00277             if (tmpNormsPtr[fIndex] < std::abs(rval)) tmpNormsPtr[fIndex] = rval;
00278             break;
00279           case 1:
00280             tmpNormsPtr[fIndex] += std::abs(rval);
00281             break;
00282           case 2:
00283             tmpNormsPtr[fIndex] += rval*rval;
00284             break;
00285           default:
00286             FEI_COUT << "Filter::residualNorm: ERROR, whichNorm="<<whichNorm
00287               << " not recognized." << FEI_ENDL;
00288             return(-1);
00289         }
00290       }
00291     }
00292   }
00293 
00294   //so at this point we have the local norms. We now need to perform the
00295   //global max or sum, depending on which norm it is.
00296 
00297   MPI_Comm comm = problemStructure_->getCommunicator();
00298 
00299   if (whichNorm != 0) {
00300     CHK_ERR( fei::GlobalSum(comm, tmpNorms, normsArray) );
00301   }
00302   else {
00303     CHK_ERR( fei::GlobalMax(comm, tmpNorms, normsArray) );
00304   }
00305 
00306   for(int i=0; i<numFields; ++i) {
00307     norms[i] = normsArray[i];
00308   }
00309 
00310   if (whichNorm == 2) {
00311     for(int i=0; i<numFields; ++i) norms[i] = std::sqrt(norms[i]);
00312   }
00313 
00314   return(0);
00315 }
00316 
00317 //------------------------------------------------------------------------------
00318 int Filter::parameters(int numParams, const char *const* paramStrings)
00319 {
00320   if (numParams == 0 || paramStrings == NULL) return(0);
00321 
00322   const char* param = snl_fei::getParam("outputLevel",numParams,paramStrings);
00323 
00324   if ( param != NULL){
00325     std::string str(&(param[11]));
00326     FEI_ISTRINGSTREAM isstr(str);
00327     isstr >> outputLevel_;
00328   }
00329 
00330   return(0);
00331 }
00332 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends