Albany: a Trilinos-based PDE code: PHAL_Source

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <cmath>
00008 #include <sstream>
00009 #include <iostream>
00010 #include <fstream>
00011 #include <algorithm>
00012 
00013 #include "Sacado_ParameterAccessor.hpp"
00014 #include "Sacado_ParameterRegistration.hpp"
00015 #include "Teuchos_VerboseObject.hpp"
00016 #include "Albany_Utils.hpp"
00017 #include "Stokhos_KL_ExponentialRandomField.hpp"
00018 #include "Teuchos_Array.hpp"
00019 #include "Teuchos_TestForException.hpp"
00020 
00021 namespace PHAL {
00022 
00023 namespace Source_Functions {
00024 const double pi = 3.1415926535897932385;
00025 
00026 template <typename EvalT, typename Traits>
00027 class Source_Base {
00028 public :
00029   Source_Base(){};
00030   virtual ~Source_Base(){};
00031   virtual void EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p)                     = 0;
00032   virtual void DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p)                     = 0;
00033   virtual void FieldData      (PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm) = 0;
00034   virtual void evaluateFields (typename Traits::EvalData workset)                                         = 0;
00035 };
00036 
00038 
00039 template<typename EvalT, typename Traits>
00040 class Constant : 
00041     public Source_Base<EvalT,Traits>, 
00042     public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00043 public :
00044   typedef typename EvalT::ScalarT ScalarT;
00045   typedef typename EvalT::MeshScalarT MeshScalarT;
00046   static bool check_for_existance(Teuchos::ParameterList* source_list);
00047   Constant(Teuchos::ParameterList& p);
00048   virtual ~Constant(){}
00049   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00050              Teuchos::ParameterList& p);
00051   virtual void DependentFields(Source<EvalT,Traits> &source, 
00052              Teuchos::ParameterList& p);
00053   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00054        PHX::FieldManager<Traits>& fm);
00055   virtual void evaluateFields (typename Traits::EvalData workset);
00056   virtual ScalarT & getValue(const std::string &n) { return m_constant;};
00057 private :
00058   ScalarT     m_constant;
00059   std::size_t m_num_qp;
00060   Teuchos::ParameterList* m_source_list;
00061   PHX::MDField<ScalarT,Cell,Point> m_source;
00062 };
00063 
00064 template<typename EvalT,typename Traits>
00065 bool 
00066 Constant<EvalT,Traits>::
00067 check_for_existance(Teuchos::ParameterList* source_list)
00068 {
00069   const bool exists = source_list->getEntryPtr("Constant");
00070   return exists;
00071 }
00072 
00073 template<typename EvalT,typename Traits>
00074 Constant<EvalT,Traits>::
00075 Constant(Teuchos::ParameterList& p) {
00076   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00077   Teuchos::ParameterList& paramList = m_source_list->sublist("Constant");
00078   m_constant = paramList.get("Value", 0.0);
00079   // Add the factor as a Sacado-ized parameter
00080   Teuchos::RCP<ParamLib> paramLib = 
00081     p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00082   new Sacado::ParameterRegistration<EvalT, SPL_Traits> ("Constant Source Value",
00083               this, paramLib);
00084 }
00085 
00086 template<typename EvalT,typename Traits>
00087 void Constant<EvalT,Traits>::
00088 EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00089   Teuchos::RCP<PHX::DataLayout> dl = 
00090     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00091   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00092   m_source = f ;
00093   source.addEvaluatedField(m_source);
00094 }
00095 
00096 template<typename EvalT,typename Traits>
00097 void 
00098 Constant<EvalT,Traits>::
00099 DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00100 }
00101 
00102 template<typename EvalT,typename Traits>
00103 void 
00104 Constant<EvalT,Traits>::
00105 FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00106     PHX::FieldManager<Traits>& fm){
00107   utils.setFieldData(m_source, fm);
00108   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00109   m_source.dimensions(dims);
00110   m_num_qp = dims[1];
00111 }
00112 
00113 template<typename EvalT,typename Traits>
00114 void 
00115 Constant<EvalT,Traits>::
00116 evaluateFields(typename Traits::EvalData workset){
00117 
00118   // Loop over cells, quad points: compute Constant Source Term
00119   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00120     for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00121       m_source(cell, iqp) = m_constant;
00122   }
00123 }
00124 
00125 
00127 
00128 template<typename EvalT, typename Traits>
00129 class Table : 
00130     public Source_Base<EvalT,Traits>, 
00131     public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00132 public :
00133   typedef typename EvalT::ScalarT ScalarT;
00134   typedef typename EvalT::MeshScalarT MeshScalarT;
00135   static bool check_for_existance(Teuchos::ParameterList* source_list);
00136   Table(Teuchos::ParameterList& p);
00137   virtual ~Table(){}
00138   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00139              Teuchos::ParameterList& p);
00140   virtual void DependentFields(Source<EvalT,Traits> &source, 
00141              Teuchos::ParameterList& p);
00142   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00143        PHX::FieldManager<Traits>& fm);
00144   virtual void evaluateFields (typename Traits::EvalData workset);
00145   virtual ScalarT & getValue(const std::string &n) { return m_constant;};
00146 private :
00147   ScalarT     m_constant;
00148   std::size_t m_num_qp;
00149   Teuchos::ParameterList* m_source_list;
00150   PHX::MDField<ScalarT,Cell,Point> m_source;
00151   std::vector<double> time;
00152   std::vector<double> sourceval;
00153   int num_time_vals;
00154 };
00155 
00156 template<typename EvalT,typename Traits>
00157 bool 
00158 Table<EvalT,Traits>::
00159 check_for_existance(Teuchos::ParameterList* source_list)
00160 {
00161   const bool exists = source_list->getEntryPtr("Table");
00162   return exists;
00163 }
00164 
00165 template<typename EvalT,typename Traits>
00166 Table<EvalT,Traits>::
00167 Table(Teuchos::ParameterList& p) {
00168   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00169   Teuchos::ParameterList& paramList = m_source_list->sublist("Table");
00170   std::string filename = paramList.get("Filename", "missing");
00171 
00172   // open file
00173   std::ifstream inFile(&filename[0]); 
00174 
00175   TEUCHOS_TEST_FOR_EXCEPTION(!inFile, Teuchos::Exceptions::InvalidParameter, std::endl <<
00176          "Error! Cannot open tabular data file \"" << filename 
00177          << "\" in source table fill" << std::endl);
00178 
00179   // Count lines in file
00180   int array_size = std::count(std::istreambuf_iterator<char>(inFile), 
00181              std::istreambuf_iterator<char>(), '\n');
00182 
00183   // Allocate and fill arrays
00184   time.resize(array_size);
00185   sourceval.resize(array_size);
00186 
00187   // rewind file
00188   inFile.seekg(0);
00189 
00190   for(num_time_vals = 0; num_time_vals < array_size; num_time_vals++){
00191 
00192     if(inFile.eof()) break;
00193     inFile >> time[num_time_vals];
00194     if(inFile.eof()) break;
00195     inFile >> sourceval[num_time_vals];
00196 
00197 //std::cout << "time " << num_time_vals << " is " << time[num_time_vals] 
00198 // << " " << sourceval[num_time_vals] << std::endl;
00199   }
00200 
00201   inFile.close();
00202 
00203   m_constant = sourceval[0];
00204 
00205   // Add the factor as a Sacado-ized parameter
00206   Teuchos::RCP<ParamLib> paramLib = 
00207     p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00208   new Sacado::ParameterRegistration<EvalT, SPL_Traits> ("Table Source Value",
00209               this, paramLib);
00210 }
00211 
00212 template<typename EvalT,typename Traits>
00213 void Table<EvalT,Traits>::
00214 EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00215   Teuchos::RCP<PHX::DataLayout> dl = 
00216     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00217   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00218   m_source = f ;
00219   source.addEvaluatedField(m_source);
00220 }
00221 
00222 template<typename EvalT,typename Traits>
00223 void 
00224 Table<EvalT,Traits>::
00225 DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00226 }
00227 
00228 template<typename EvalT,typename Traits>
00229 void 
00230 Table<EvalT,Traits>::
00231 FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00232     PHX::FieldManager<Traits>& fm){
00233   utils.setFieldData(m_source, fm);
00234   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00235   m_source.dimensions(dims);
00236   m_num_qp = dims[1];
00237 }
00238 
00239 template<typename EvalT,typename Traits>
00240 void 
00241 Table<EvalT,Traits>::
00242 evaluateFields(typename Traits::EvalData workset){
00243 
00244   if(workset.current_time <= 0.0) // if time is uninitialized or zero, just take first value
00245 
00246     m_constant = sourceval[0];
00247 
00248   else { // Interpolate between time values
00249 
00250     bool found_it = false;
00251 
00252     for(int i = 0; i < num_time_vals - 1; i++) // Stride through time
00253 
00254       if(workset.current_time >= time[i] && workset.current_time <= time[i + 1] ){ // Have bracketed current time
00255 
00256         double s = (workset.current_time - time[i]) / (time[i + 1] - time[i]); // 0 \leq s \leq 1
00257 
00258         m_constant = sourceval[i] + s * (sourceval[i + 1] - sourceval[i]); // interp value corresponding to s
00259 
00260         found_it = true;
00261 
00262         break;
00263 
00264       }
00265 
00266     TEUCHOS_TEST_FOR_EXCEPTION(!found_it, Teuchos::Exceptions::InvalidParameter, std::endl <<
00267          "Error! Cannot locate the current time \"" << workset.current_time 
00268          << "\" in the time series data between the endpoints " << time[0]
00269           << " and " << time[num_time_vals - 1] << "." << std::endl);
00270   }
00271 
00272   // Loop over cells, quad points: compute Table Source Term
00273   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00274     for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00275       m_source(cell, iqp) = m_constant;
00276   }
00277 }
00278 
00280 
00281 
00282 template<typename EvalT, typename Traits>
00283 class Trigonometric : 
00284     public Source_Base<EvalT,Traits>, 
00285     public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00286 public :
00287   typedef typename EvalT::ScalarT ScalarT;
00288   typedef typename EvalT::MeshScalarT MeshScalarT;
00289   static bool check_for_existance(Teuchos::ParameterList* source_list);
00290   Trigonometric(Teuchos::ParameterList& p);
00291   virtual ~Trigonometric(){}
00292   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00293              Teuchos::ParameterList& p);
00294   virtual void DependentFields(Source<EvalT,Traits> &source, 
00295              Teuchos::ParameterList& p);
00296   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00297        PHX::FieldManager<Traits>& fm);
00298   virtual void evaluateFields (typename Traits::EvalData workset);
00299   virtual ScalarT & getValue(const std::string &n) { return m_constant;};
00300 private :
00301   ScalarT     m_constant; 
00302   std::size_t m_num_qp;
00303   std::size_t m_num_dim;
00304   Teuchos::ParameterList* m_source_list;
00305   PHX::MDField<ScalarT,Cell,Point> m_source;
00306   PHX::MDField<MeshScalarT,Cell,Point,Dim> coordVec;
00307 };
00308 
00309 template<typename EvalT,typename Traits>
00310 bool 
00311 Trigonometric<EvalT,Traits>::
00312 check_for_existance(Teuchos::ParameterList* source_list)
00313 {
00314   const bool exists = source_list->getEntryPtr("Trigonometric");
00315   return exists;
00316 }
00317 
00318 template<typename EvalT,typename Traits>
00319 Trigonometric<EvalT,Traits>::
00320 Trigonometric(Teuchos::ParameterList& p) {
00321   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00322   Teuchos::ParameterList& paramList = m_source_list->sublist("Trigonometric");
00323   m_constant = paramList.get("Value", 1.0); 
00324   // Add the factor as a Sacado-ized parameter
00325   Teuchos::RCP<ParamLib> paramLib = 
00326     p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00327   new Sacado::ParameterRegistration<EvalT, SPL_Traits> ("Trigonometric Source Value",
00328               this, paramLib);
00329 }
00330 
00331 template<typename EvalT,typename Traits>
00332 void Trigonometric<EvalT,Traits>::
00333 EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00334   Teuchos::RCP<PHX::DataLayout> dl = 
00335     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00336   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00337   m_source = f ;
00338   source.addEvaluatedField(m_source);
00339 }
00340 
00341 template<typename EvalT,typename Traits>
00342 void 
00343 Trigonometric<EvalT,Traits>::
00344 DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) 
00345 {
00346   Teuchos::RCP<PHX::DataLayout> scalar_qp = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00347   Teuchos::RCP<PHX::DataLayout> vector_qp = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");  
00348 
00349   PHX::MDField<MeshScalarT,Cell,Point,Dim> f0
00350     (p.get<std::string>("QP Coordinate Vector Name"),  vector_qp);
00351   coordVec = f0;
00352   source.addDependentField(coordVec);
00353 }
00354 
00355 template<typename EvalT,typename Traits>
00356 void 
00357 Trigonometric<EvalT,Traits>::
00358 FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00359     PHX::FieldManager<Traits>& fm){
00360   utils.setFieldData(m_source, fm);
00361   utils.setFieldData(coordVec,fm);
00362   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00363   m_source.dimensions(dims);
00364   coordVec.dimensions(dims); 
00365   m_num_qp = dims[1];
00366   m_num_dim = dims[2]; 
00367 }
00368 
00369 template<typename EvalT,typename Traits>
00370 void 
00371 Trigonometric<EvalT,Traits>::
00372 evaluateFields(typename Traits::EvalData workset){
00373 
00374   // Loop over cells, quad points: compute Trigonometric Source Term
00375   if (m_num_dim == 2) {
00376     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00377       for (std::size_t iqp=0; iqp<m_num_qp; iqp++) {
00378         MeshScalarT *X = &coordVec(cell,iqp,0); 
00379         m_source(cell, iqp) = 8.0*pi*pi*sin(2.0*pi*X[0])*cos(2.0*pi*X[1]);
00380       }
00381     }
00382   }
00383   else {
00384     std::cout << "Trigonometric source implemented only for 2D; setting f = 1 constant source." << std::endl; 
00385     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00386       for (std::size_t iqp=0; iqp<m_num_qp; iqp++) {
00387         m_source(cell, iqp) = m_constant;
00388       }
00389     }
00390   }
00391 }
00392 
00394 
00395 
00396 template<typename EvalT, typename Traits>
00397 class Quadratic : 
00398     public Source_Base<EvalT,Traits>, 
00399     public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00400 public :
00401   typedef typename EvalT::ScalarT ScalarT;
00402   typedef typename EvalT::MeshScalarT MeshScalarT;
00403   static bool check_for_existance(Teuchos::ParameterList* source_list);
00404   Quadratic(Teuchos::ParameterList& p);
00405   virtual ~Quadratic(){}
00406   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00407              Teuchos::ParameterList& p);
00408   virtual void DependentFields(Source<EvalT,Traits> &source, 
00409              Teuchos::ParameterList& p);
00410   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00411        PHX::FieldManager<Traits>& fm);
00412   virtual void evaluateFields (typename Traits::EvalData workset);
00413   virtual ScalarT & getValue(const std::string &n) { return m_factor;};
00414 private :
00415   ScalarT     m_factor;
00416   double      m_constant;
00417   std::size_t m_num_qp;
00418   Teuchos::ParameterList* m_source_list;
00419   PHX::MDField<ScalarT,Cell,Point>   m_source;
00420   PHX::MDField<ScalarT,Cell,Point>   m_baseField;
00421 };
00422 
00423 template<typename EvalT,typename Traits>
00424 bool 
00425 Quadratic<EvalT,Traits>::
00426 check_for_existance(Teuchos::ParameterList* source_list)
00427 {
00428   const bool exists = source_list->getEntryPtr("Quadratic");
00429   return exists;
00430 }
00431 
00432 template<typename EvalT,typename Traits>
00433 Quadratic<EvalT,Traits>::
00434 Quadratic(Teuchos::ParameterList& p) {
00435   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00436   Teuchos::ParameterList& paramList = m_source_list->sublist("Quadratic");
00437   m_factor   = paramList.get("Nonlinear Factor", 0.0);
00438   m_constant = paramList.get("Constant", false);
00439   // Add the factor as a Sacado-ized parameter
00440   Teuchos::RCP<ParamLib> paramLib = 
00441     p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00442   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (
00443     "Quadratic Nonlinear Factor", this, paramLib);
00444 }
00445 
00446 template<typename EvalT,typename Traits>
00447 void 
00448 Quadratic<EvalT,Traits>::
00449 EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00450   Teuchos::RCP<PHX::DataLayout> dl = 
00451     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00452   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00453   m_source = f ;
00454   source.addEvaluatedField(m_source);
00455 }
00456 template<typename EvalT,typename Traits>
00457 void 
00458 Quadratic<EvalT,Traits>::
00459 DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00460   Teuchos::RCP<PHX::DataLayout> dl = 
00461     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00462   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Variable Name"), dl);
00463   m_baseField = f;
00464   source.addDependentField(m_baseField);
00465 }
00466 template<typename EvalT,typename Traits>
00467 void 
00468 Quadratic<EvalT,Traits>::
00469 FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00470     PHX::FieldManager<Traits>& fm){
00471   utils.setFieldData(m_source,   fm);
00472   utils.setFieldData(m_baseField,fm);
00473   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00474   m_source.dimensions(dims);
00475   m_num_qp = dims[1];
00476 }
00477 
00478 template<typename EvalT,typename Traits>
00479 void 
00480 Quadratic<EvalT,Traits>::
00481 evaluateFields(typename Traits::EvalData workset){
00482 
00483   // Loop over cells, quad points: compute Quadratic Source Term
00484   if (m_constant) {
00485     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00486       for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00487         m_source(cell, iqp) = m_factor;
00488     }
00489   } 
00490   else {
00491     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00492       for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00493         m_source(cell, iqp) = 
00494     m_factor * m_baseField(cell,iqp) * m_baseField(cell,iqp);
00495     }
00496   }
00497 }
00498 
00500 
00501 template<typename EvalT, typename Traits>
00502 class TruncatedKL : 
00503     public Source_Base<EvalT,Traits>, 
00504     public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00505 public :
00506   typedef typename EvalT::ScalarT ScalarT;
00507   typedef typename EvalT::MeshScalarT MeshScalarT;
00508   static bool check_for_existance(Teuchos::ParameterList* source_list);
00509   TruncatedKL(Teuchos::ParameterList& p);
00510   virtual ~TruncatedKL(){}
00511   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00512              Teuchos::ParameterList& p);
00513   virtual void DependentFields(Source<EvalT,Traits> &source, 
00514              Teuchos::ParameterList& p);
00515   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00516        PHX::FieldManager<Traits>& fm);
00517   virtual void evaluateFields (typename Traits::EvalData workset);
00518   virtual ScalarT & getValue(const std::string &n);
00519 private :
00520   std::size_t m_num_qp;
00521   std::size_t m_num_dims;
00522   Teuchos::ParameterList* m_source_list;
00523   PHX::MDField<ScalarT,Cell,Point>   m_source;
00524   PHX::MDField<MeshScalarT,Cell,Point,Dim> m_coordVec;
00525   Teuchos::RCP< Stokhos::KL::ExponentialRandomField<MeshScalarT> > m_exp_rf_kl;
00526   Teuchos::Array<ScalarT> m_rv;
00527   Teuchos::Array<MeshScalarT> m_point;
00528   std::string param_name_base;
00529 };
00530 
00531 template<typename EvalT,typename Traits>
00532 bool 
00533 TruncatedKL<EvalT,Traits>::
00534 check_for_existance(Teuchos::ParameterList* source_list)
00535 {
00536   const bool exists = source_list->getEntryPtr("Truncated KL Expansion");
00537   return exists;
00538 }
00539 
00540 template<typename EvalT,typename Traits>
00541 TruncatedKL<EvalT,Traits>::
00542 TruncatedKL(Teuchos::ParameterList& p) {
00543   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00544   Teuchos::ParameterList& paramList = 
00545     m_source_list->sublist("Truncated KL Expansion");
00546   
00547   m_exp_rf_kl = 
00548       Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalarT>(paramList));
00549   int num_KL = m_exp_rf_kl->stochasticDimension();
00550 
00551   param_name_base = p.get<std::string>("Source Name") + " KL Random Variable";
00552 
00553   // Add KL random variables as Sacado-ized parameters
00554   m_rv.resize(num_KL);
00555   Teuchos::RCP<ParamLib> paramLib = 
00556     p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00557   for (int i=0; i<num_KL; i++) {
00558     std::string ss = Albany::strint(param_name_base,i);
00559     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(ss, this, paramLib);
00560     m_rv[i] = paramList.get(ss, 0.0);
00561   }
00562 }
00563 
00564 template<typename EvalT,typename Traits>
00565 void 
00566 TruncatedKL<EvalT,Traits>::
00567 EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00568   Teuchos::RCP<PHX::DataLayout> dl = 
00569     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00570   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00571   m_source = f ;
00572   source.addEvaluatedField(m_source);
00573 }
00574 
00575 template<typename EvalT,typename Traits>
00576 void 
00577 TruncatedKL<EvalT,Traits>::
00578 DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00579   Teuchos::RCP<PHX::DataLayout> vector_dl =
00580     p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00581   PHX::MDField<MeshScalarT,Cell,Point,Dim>
00582     fx(p.get<std::string>("QP Coordinate Vector Name"), vector_dl);
00583   m_coordVec = fx;
00584   source.addDependentField(m_coordVec);
00585 }
00586 
00587 template<typename EvalT,typename Traits>
00588 void 
00589 TruncatedKL<EvalT,Traits>::
00590 FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00591     PHX::FieldManager<Traits>& fm){
00592   utils.setFieldData(m_source,   fm);
00593   utils.setFieldData(m_coordVec,fm);
00594   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00595   m_coordVec.dimensions(dims);
00596   m_num_qp = dims[1];
00597   m_num_dims = dims[2];
00598   m_point.resize(m_num_dims);
00599 }
00600 
00601 template<typename EvalT,typename Traits>
00602 void 
00603 TruncatedKL<EvalT,Traits>::
00604 evaluateFields(typename Traits::EvalData workset){
00605 
00606   // Loop over cells, quad points: compute TruncatedKL Source Term
00607   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00608     for (std::size_t qp=0; qp<m_num_qp; qp++) {
00609       for (std::size_t i=0; i<m_num_dims; i++)
00610   m_point[i] = 
00611     Sacado::ScalarValue<MeshScalarT>::eval(m_coordVec(cell,qp,i));
00612         m_source(cell, qp) = m_exp_rf_kl->evaluate(m_point, m_rv);
00613     }
00614   }
00615 }
00616 
00617 template<typename EvalT,typename Traits>
00618 typename TruncatedKL<EvalT,Traits>::ScalarT& 
00619 TruncatedKL<EvalT,Traits>::
00620 getValue(const std::string &n) {
00621   for (int i=0; i<m_rv.size(); i++) {
00622     if (n == Albany::strint(param_name_base,i))
00623       return m_rv[i];
00624   }
00625   TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00626          std::endl <<
00627          "Error! Logic error in getting parameter " << n
00628          << " in TruncatedKL::getValue()" << std::endl);
00629 }
00630 
00632 
00633 template<typename EvalT, typename Traits>
00634 class NeutronFission : 
00635     public Source_Base<EvalT,Traits> {
00636 public :
00637   typedef typename EvalT::ScalarT ScalarT;
00638   typedef typename EvalT::MeshScalarT MeshScalarT;
00639 
00640   static bool check_for_existance(Teuchos::ParameterList* source_list) {
00641     return source_list->isSublist("Neutron Fission");
00642   }
00643 
00644   NeutronFission(Teuchos::ParameterList& p) {
00645     m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00646   }
00647 
00648   virtual ~NeutronFission() {}
00649 
00650   virtual void EvaluatedFields(Source<EvalT,Traits> &source, 
00651              Teuchos::ParameterList& p) {
00652     Teuchos::RCP<PHX::DataLayout> dl = 
00653       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00654     PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00655     m_source = f ;
00656     source.addEvaluatedField(m_source);
00657   }
00658 
00659   virtual void DependentFields(Source<EvalT,Traits> &source, 
00660              Teuchos::ParameterList& p) {
00661     Teuchos::RCP<PHX::DataLayout> dl = 
00662       p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00663     m_phi = PHX::MDField<ScalarT,Cell,Point>(
00664       p.get<std::string>("Neutron Flux Name"), dl);
00665     m_sigma_f = PHX::MDField<ScalarT,Cell,Point>(
00666       p.get<std::string>("Fission Cross Section Name"), dl);
00667     m_E_f = PHX::MDField<ScalarT,Cell,Point>(
00668       p.get<std::string>("Energy Released per Fission Name"), dl);
00669     source.addDependentField(m_phi);
00670     source.addDependentField(m_sigma_f);
00671     source.addDependentField(m_E_f);
00672   }
00673 
00674   virtual void FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, 
00675        PHX::FieldManager<Traits>& fm) {
00676     utils.setFieldData(m_source,   fm);
00677     utils.setFieldData(m_phi,fm);
00678     utils.setFieldData(m_sigma_f,fm);
00679     utils.setFieldData(m_E_f,fm);
00680     typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00681     m_source.dimensions(dims);
00682     m_num_qp = dims[1];
00683   }
00684 
00685   virtual void evaluateFields (typename Traits::EvalData workset) {
00686     for (std::size_t cell = 0; cell < workset.numCells; ++cell)
00687       for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00688         m_source(cell, iqp) = 
00689     m_phi(cell,iqp) * m_sigma_f(cell,iqp) * m_E_f(cell,iqp);
00690   }
00691 
00692 private :
00693   std::size_t m_num_qp;
00694   Teuchos::ParameterList* m_source_list;
00695   PHX::MDField<ScalarT,Cell,Point>   m_source;
00696   PHX::MDField<ScalarT,Cell,Point>   m_phi;
00697   PHX::MDField<ScalarT,Cell,Point>   m_sigma_f;
00698   PHX::MDField<ScalarT,Cell,Point>   m_E_f;
00699 };
00700 
00702 
00703 template<typename EvalT, typename Traits>
00704 class MVQuadratic : public Source_Base<EvalT,Traits>, public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00705 public :
00706   typedef typename EvalT::ScalarT ScalarT;
00707   typedef typename EvalT::MeshScalarT MeshScalarT;
00708   static bool check_for_existance(Teuchos::ParameterList* source_list);
00709   MVQuadratic(Teuchos::ParameterList& p);
00710   virtual ~MVQuadratic(){}
00711   virtual void EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
00712   virtual void DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
00713   virtual void FieldData      (PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm);
00714   virtual void evaluateFields (typename Traits::EvalData workset);
00715   virtual ScalarT & getValue(const std::string &n);
00716 private :
00717   std::vector<ScalarT>     m_factor;
00718   std::size_t m_num_qp;
00719   Teuchos::ParameterList* m_source_list;
00720   PHX::MDField<ScalarT,Cell,Point>   m_source;
00721   PHX::MDField<ScalarT,Cell,Point>   m_baseField;
00722 };
00723 
00724 template<typename EvalT,typename Traits>
00725 bool MVQuadratic<EvalT,Traits>::check_for_existance(Teuchos::ParameterList* source_list)
00726 {
00727   const bool exists = source_list->getEntryPtr("Multivariate Quadratic");
00728   return exists;
00729 }
00730 
00731 template<typename EvalT,typename Traits>
00732 MVQuadratic<EvalT,Traits>::MVQuadratic(Teuchos::ParameterList& p) {
00733   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00734   Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Quadratic");
00735   int num_vars = paramList.get("Dimension", 1);
00736   m_factor.resize(num_vars);
00737   Teuchos::RCP<ParamLib> paramLib = p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00738   for (int i=0; i<num_vars; i++) {
00739     m_factor[i] = paramList.get(Albany::strint("Nonlinear Factor",i), 0.0);
00740 
00741     // Add the factor as a Sacado-ized parameter
00742     new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00743       Albany::strint("Multivariate Quadratic Nonlinear Factor",i), this, paramLib);
00744   }
00745 }
00746 
00747 template<typename EvalT,typename Traits>
00748 void MVQuadratic<EvalT,Traits>::EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00749   //Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Quadratic");
00750   Teuchos::RCP<PHX::DataLayout> dl = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00751   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00752   m_source = f ;
00753   source.addEvaluatedField(m_source);
00754 }
00755 template<typename EvalT,typename Traits>
00756 void MVQuadratic<EvalT,Traits>::DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00757   //Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Quadratic");
00758   Teuchos::RCP<PHX::DataLayout> dl = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00759   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Variable Name"), dl);
00760   m_baseField = f;
00761   source.addDependentField(m_baseField);
00762 }
00763 template<typename EvalT,typename Traits>
00764 void MVQuadratic<EvalT,Traits>::FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm){
00765   utils.setFieldData(m_source,   fm);
00766   utils.setFieldData(m_baseField,fm);
00767   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00768   m_source.dimensions(dims);
00769   m_num_qp = dims[1];
00770 }
00771 
00772 template<typename EvalT,typename Traits>
00773 void MVQuadratic<EvalT,Traits>::evaluateFields(typename Traits::EvalData workset){
00774 
00775   ScalarT a = 0.0;
00776   for (unsigned int j=0; j<m_factor.size(); j++) {
00777     a += m_factor[j];
00778   }
00779   a /= static_cast<double>(m_factor.size());
00780 
00781   // Loop over cells, quad points: compute MVQuadratic Source Term
00782   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00783     for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00784       m_source(cell, iqp) = a * m_baseField(cell,iqp) * m_baseField(cell,iqp);
00785   }
00786 }
00787 
00788 template<typename EvalT,typename Traits>
00789 typename MVQuadratic<EvalT,Traits>::ScalarT& 
00790 MVQuadratic<EvalT,Traits>::getValue(const std::string &n)
00791 {
00792   for (unsigned int i=0; i<m_factor.size(); i++) {
00793     if (n == Albany::strint("Multivariate Quadratic Nonlinear Factor",i))
00794       return m_factor[i];
00795   }
00796   return m_factor[0];
00797 }
00798 
00799 template<typename EvalT, typename Traits>
00800 class MVExponential : public Source_Base<EvalT,Traits>, public Sacado::ParameterAccessor<EvalT, SPL_Traits> {
00801 public :
00802   typedef typename EvalT::ScalarT ScalarT;
00803   typedef typename EvalT::MeshScalarT MeshScalarT;
00804   static bool check_for_existance(Teuchos::ParameterList* source_list);
00805   MVExponential(Teuchos::ParameterList& p);
00806   virtual ~MVExponential(){}
00807   virtual void EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
00808   virtual void DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
00809   virtual void FieldData      (PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm);
00810   virtual void evaluateFields (typename Traits::EvalData workset);
00811   virtual ScalarT & getValue(const std::string &n);
00812 private :
00813   std::vector<ScalarT>     m_factor;
00814   std::size_t m_num_qp;
00815   Teuchos::ParameterList* m_source_list;
00816   PHX::MDField<ScalarT,Cell,Point>   m_source;
00817   PHX::MDField<ScalarT,Cell,Point>   m_baseField;
00818 };
00819 
00820 template<typename EvalT,typename Traits>
00821 bool MVExponential<EvalT,Traits>::check_for_existance(Teuchos::ParameterList* source_list)
00822 {
00823   const bool exists = source_list->getEntryPtr("Multivariate Exponential");
00824   return exists;
00825 }
00826 
00827 template<typename EvalT,typename Traits>
00828 MVExponential<EvalT,Traits>::MVExponential(Teuchos::ParameterList& p) {
00829   m_source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
00830   Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Exponential");
00831   int num_vars = paramList.get("Dimension", 1);
00832   m_factor.resize(num_vars);
00833   Teuchos::RCP<ParamLib> paramLib = p.get< Teuchos::RCP<ParamLib> > ("Parameter Library", Teuchos::null);
00834   for (int i=0; i<num_vars; i++) {
00835     m_factor[i] = paramList.get(Albany::strint("Nonlinear Factor",i), 0.0);
00836 
00837     // Add the factor as a Sacado-ized parameter
00838     new Sacado::ParameterRegistration<EvalT, SPL_Traits> (
00839       Albany::strint("Multivariate Exponential Nonlinear Factor",i), this, paramLib);
00840   }
00841 }
00842 
00843 template<typename EvalT,typename Traits>
00844 void MVExponential<EvalT,Traits>::EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00845   //Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Exponential");
00846   Teuchos::RCP<PHX::DataLayout> dl = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00847   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Source Name"), dl);
00848   m_source = f ;
00849   source.addEvaluatedField(m_source);
00850 }
00851 template<typename EvalT,typename Traits>
00852 void MVExponential<EvalT,Traits>::DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p) {
00853   //Teuchos::ParameterList& paramList = m_source_list->sublist("Multivariate Exponential");
00854   Teuchos::RCP<PHX::DataLayout> dl = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00855   PHX::MDField<ScalarT,Cell,Point> f(p.get<std::string>("Variable Name"), dl);
00856   m_baseField = f;
00857   source.addDependentField(m_baseField);
00858 }
00859 template<typename EvalT,typename Traits>
00860 void MVExponential<EvalT,Traits>::FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm){
00861   utils.setFieldData(m_source,   fm);
00862   utils.setFieldData(m_baseField,fm);
00863   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Node>::size_type > dims;
00864   m_source.dimensions(dims);
00865   m_num_qp = dims[1];
00866 }
00867 
00868 template<typename EvalT,typename Traits>
00869 void MVExponential<EvalT,Traits>::evaluateFields(typename Traits::EvalData workset){
00870 
00871   ScalarT a = 0.0;
00872   for (unsigned int j=0; j<m_factor.size(); j++) {
00873     a += m_factor[j];
00874   }
00875   a /= static_cast<double>(m_factor.size());
00876 
00877   // Loop over cells, quad points: compute MVExponential Source Term
00878   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00879     for (std::size_t iqp=0; iqp<m_num_qp; iqp++)
00880       m_source(cell, iqp) = a*std::exp(m_baseField(cell,iqp));
00881   }
00882 }
00883 
00884 template<typename EvalT,typename Traits>
00885 typename MVExponential<EvalT,Traits>::ScalarT& 
00886 MVExponential<EvalT,Traits>::getValue(const std::string &n)
00887 {
00888   for (unsigned int i=0; i<m_factor.size(); i++) {
00889     if (n == Albany::strint("Multivariate Exponential Nonlinear Factor",i))
00890       return m_factor[i];
00891   }
00892   return m_factor[0];
00893 }
00894 
00895 
00896 // This needs to be templated to get mesh derivatives
00897 template<typename EvalT>
00898 class Spatial_Base {
00899 public :
00900   virtual ~Spatial_Base(){}
00901   typedef typename EvalT::MeshScalarT MeshScalarT;
00902   virtual MeshScalarT evaluateFields  (const std::vector<MeshScalarT> &coords)=0;
00903 };
00904 
00905 template<typename EvalT>
00906 class Gaussian : public Spatial_Base<EvalT> {
00907 public :
00908   static bool check_for_existance(Teuchos::ParameterList &source_list);
00909   typedef typename EvalT::MeshScalarT MeshScalarT;
00910   Gaussian(Teuchos::ParameterList &source_list, std::size_t num);
00911   virtual ~Gaussian(){}
00912   virtual MeshScalarT evaluateFields(const std::vector<MeshScalarT> &coords);
00913 private :
00914   double      m_amplitude        ;
00915   double      m_radius           ;
00916   double      m_sigma_sq         ;
00917   double      m_sigma_pi         ;
00918   Teuchos::Array<double> m_centroid ;
00919 };
00920 
00921 template<typename EvalT>
00922 inline bool Gaussian<EvalT>::check_for_existance(Teuchos::ParameterList &source_list)
00923 {
00924   std::string g("Gaussian");
00925   bool exists = source_list.getEntryPtr("Type");
00926   if (exists) exists = g==source_list.get("Type",g);
00927   return exists;
00928 }
00929 
00930 template<typename EvalT>
00931 inline Gaussian<EvalT>::Gaussian(Teuchos::ParameterList &source_list, std::size_t num)
00932 {
00933   Teuchos::ParameterList& paramList = source_list.sublist("Spatial",true);
00934   m_amplitude = paramList.get("Amplitude",      1.0);
00935   m_radius    = paramList.get("Radius",         1.0);
00936   m_centroid  = source_list.get(Albany::strint("Center",num),m_centroid);
00937   m_sigma_sq = 1.0/(2.0*std::pow(m_radius, 2));
00938   const double pi = 3.1415926535897932385;
00939   m_sigma_pi = 1.0/(m_radius*std::sqrt(2*pi));
00940 }
00941 
00942 template<typename EvalT>
00943 typename EvalT::MeshScalarT Gaussian<EvalT>::
00944 evaluateFields(const std::vector<typename EvalT::MeshScalarT> &coords)
00945 {
00946   MeshScalarT exponent=0;
00947   const std::size_t nsd = coords.size();
00948   for (std::size_t i=0; i<nsd; ++i) {
00949     exponent += std::pow(m_centroid[i]-coords[i],2);
00950   }
00951   exponent *= m_sigma_sq;  
00952   MeshScalarT x(0.0);
00953   if (nsd==1)
00954     x = m_amplitude *          m_sigma_pi    * std::exp(-exponent);           
00955   else if (nsd==2)
00956     x = m_amplitude * std::pow(m_sigma_pi,2) * std::exp(-exponent);           
00957   else if (nsd==3)
00958     x = m_amplitude * std::pow(m_sigma_pi,3) * std::exp(-exponent);           
00959   return x;
00960 }
00961 
00962 
00963 class Wavelet_Base {
00964 public :
00965   virtual ~Wavelet_Base(){}
00966   virtual RealType evaluateFields  (const RealType time) = 0;
00967 };
00968 
00969 class Monotone : public Wavelet_Base {
00970 public :
00971   static bool check_for_existance(Teuchos::ParameterList &source_list);
00972   Monotone(Teuchos::ParameterList &source_list);
00973   virtual ~Monotone(){}
00974   virtual RealType evaluateFields  (const RealType time);
00975 private :
00976 };
00977 
00978 inline bool Monotone::check_for_existance(Teuchos::ParameterList &source_list)
00979 { 
00980   std::string g("Monotone");
00981   bool exists = source_list.getEntryPtr("Type");
00982   if (exists) exists = g==source_list.get("Type",g);
00983   return exists;
00984 }
00985 
00986 inline Monotone::Monotone(Teuchos::ParameterList &source_list)
00987 {
00988 }
00989 
00990 inline RealType Monotone::evaluateFields(const RealType time)
00991 { return 1.0; }
00992 
00993 
00994 template<typename EvalT, typename Traits>
00995 class PointSource : public Source_Base<EvalT,Traits> {
00996 public :
00997   typedef typename EvalT::ScalarT ScalarT;
00998   typedef typename EvalT::MeshScalarT MeshScalarT;
00999   static bool check_for_existance(Teuchos::ParameterList* source_list);
01000   PointSource(Teuchos::ParameterList* source_list);
01001   virtual ~PointSource();
01002   virtual void EvaluatedFields (Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
01003   virtual void DependentFields (Source<EvalT,Traits> &source, Teuchos::ParameterList& p);
01004   virtual void FieldData       (PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm);
01005   virtual void evaluateFields  (typename Traits::EvalData workset);
01006 private :
01007   std::size_t                 m_num_qp;
01008   std::size_t                 m_num_dim;
01009   Wavelet_Base               *m_wavelet;
01010   std::vector<Spatial_Base<EvalT> * >  m_spatials;
01011 
01012   Teuchos::ParameterList* m_source_list;
01013               PHX::MDField<ScalarT,Cell,Point>   m_pressure_source;
01014   PHX::MDField<MeshScalarT,Cell,Point,Dim> coordVec;
01015 };
01016 template<typename EvalT, typename Traits>
01017 bool PointSource<EvalT,Traits>::check_for_existance(Teuchos::ParameterList* source_list)
01018 {
01019   const bool exists = source_list->getEntryPtr("Point");
01020   return exists;
01021 }
01022 
01023 template<typename EvalT, typename Traits>
01024 PointSource<EvalT,Traits>::PointSource(Teuchos::ParameterList* source_list) :
01025   m_num_qp(0), m_wavelet(NULL), m_spatials(), m_source_list(source_list)
01026 {
01027   Teuchos::ParameterList& paramList = source_list->sublist("Point",true);
01028   const std::size_t num  = paramList.get("Number", 0);
01029   Teuchos::ParameterList& spatial_param = paramList.sublist("Spatial",true);
01030   if (Gaussian<EvalT>::check_for_existance(spatial_param)) {
01031     for (std::size_t i=0; i<num; ++i) {
01032       Gaussian<EvalT> *s = new Gaussian<EvalT>(paramList,i);
01033       m_spatials.push_back(s);
01034     }
01035   } else {
01036     TEUCHOS_TEST_FOR_EXCEPTION(m_spatials.empty(), std::logic_error,
01037                        "Point: Did not find a single spatial component.  Specify Gaussian.");
01038   }
01039   Teuchos::ParameterList& wavelet_param = paramList.sublist("Time Wavelet",true);
01040   if (Monotone::check_for_existance(wavelet_param)) {
01041     Monotone *monotone = new Monotone(wavelet_param);
01042     m_wavelet = monotone;
01043   } else {
01044     TEUCHOS_TEST_FOR_EXCEPTION(!m_wavelet, std::logic_error,
01045                        "Point: Did not find a single wavelet component. Specify Monotone.");
01046   }
01047 }
01048 template<typename EvalT, typename Traits>
01049 PointSource<EvalT,Traits>::~PointSource()
01050 {
01051   delete m_wavelet;
01052   for (std::size_t i=0; i<m_spatials.size(); ++i) {
01053     Spatial_Base<EvalT> *s = m_spatials[i];
01054     delete s;
01055   }
01056   m_spatials.clear();
01057 }
01058 
01059 template<typename EvalT, typename Traits>
01060 void PointSource<EvalT,Traits>::EvaluatedFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p)
01061 {
01062   Teuchos::RCP<PHX::DataLayout> scalar_qp = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
01063   PHX::MDField<ScalarT,Cell,Point> f0(p.get<std::string>("Pressure Source Name"), scalar_qp);
01064   m_pressure_source       = f0;
01065   source.addEvaluatedField(m_pressure_source);
01066 }
01067 
01068 template<typename EvalT, typename Traits>
01069 void PointSource<EvalT,Traits>::DependentFields(Source<EvalT,Traits> &source, Teuchos::ParameterList& p)
01070 {
01071   Teuchos::RCP<PHX::DataLayout> scalar_qp = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
01072   Teuchos::RCP<PHX::DataLayout> vector_qp = p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
01073 
01074   PHX::MDField<MeshScalarT,Cell,Point,Dim> f0
01075     (p.get<std::string>("QP Coordinate Vector Name"),  vector_qp);
01076   coordVec = f0;
01077   source.addDependentField(coordVec);
01078 }
01079 
01080 template<typename EvalT, typename Traits>
01081 void    PointSource<EvalT,Traits>::FieldData(PHX::EvaluatorUtilities<EvalT,Traits> &utils, PHX::FieldManager<Traits>& fm){
01082   utils.setFieldData(m_pressure_source, fm);
01083   utils.setFieldData(coordVec,     fm);
01084   typename std::vector< typename PHX::template MDField<ScalarT,Cell,Point,Dim>::size_type > dims;
01085   coordVec.dimensions(dims);
01086   m_num_qp  = dims[1];
01087   m_num_dim = dims[2];
01088 }
01089 
01090 template<typename EvalT, typename Traits>
01091 void PointSource<EvalT,Traits>::evaluateFields(typename Traits::EvalData workset)
01092 {
01093   const RealType time  = workset.current_time;
01094   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
01095     for (std::size_t iqp=0; iqp<m_num_qp; iqp++) {
01096       std::vector<MeshScalarT> coord;
01097       for (std::size_t i=0; i<m_num_dim; ++i) {
01098         const MeshScalarT  x  = coordVec(cell,iqp,i);
01099         coord.push_back(x);
01100       }
01101       ScalarT &p_s = m_pressure_source(cell,iqp);
01102       p_s = 0;
01103       const RealType wavelet = m_wavelet->evaluateFields(time);
01104       for (std::size_t i=0; i<m_spatials.size(); ++i) {
01105         const MeshScalarT spatial = m_spatials[i]->evaluateFields(coord);
01106         p_s += wavelet*spatial;
01107       }
01108     }
01109   }
01110 }
01111 
01112 
01113 }
01114 
01115 using namespace Source_Functions;
01116 
01117 template<typename EvalT, typename Traits>
01118 Source<EvalT, Traits>::Source(Teuchos::ParameterList& p)
01119 {
01120 
01121   Teuchos::ParameterList* source_list = p.get<Teuchos::ParameterList*>("Parameter List", NULL);
01122   if (Constant<EvalT,Traits>::check_for_existance(source_list)) {
01123     Constant<EvalT,Traits>    *q = new Constant<EvalT,Traits>(p);
01124     Source_Base<EvalT,Traits> *sb = q;
01125     m_sources.push_back(sb);
01126     this->setName("ConstantSource"+PHX::TypeString<EvalT>::value);
01127   }
01128   if (Table<EvalT,Traits>::check_for_existance(source_list)) {
01129     Table<EvalT,Traits>    *q = new Table<EvalT,Traits>(p);
01130     Source_Base<EvalT,Traits> *sb = q;
01131     m_sources.push_back(sb);
01132     this->setName("TableSource"+PHX::TypeString<EvalT>::value);
01133   }
01134   if (Trigonometric<EvalT,Traits>::check_for_existance(source_list)) {
01135     Trigonometric<EvalT,Traits>    *q = new Trigonometric<EvalT,Traits>(p);
01136     Source_Base<EvalT,Traits> *sb = q;
01137     m_sources.push_back(sb);
01138     this->setName("TrigonometricSource"+PHX::TypeString<EvalT>::value);
01139   }
01140   if (TruncatedKL<EvalT,Traits>::check_for_existance(source_list)) {
01141     TruncatedKL<EvalT,Traits>    *q = new TruncatedKL<EvalT,Traits>(p);
01142     Source_Base<EvalT,Traits> *sb = q;
01143     m_sources.push_back(sb);
01144     this->setName("TruncatedKLSource"+PHX::TypeString<EvalT>::value);
01145   }
01146   if (Quadratic<EvalT,Traits>::check_for_existance(source_list)) {
01147     Quadratic<EvalT,Traits>    *q = new Quadratic<EvalT,Traits>(p);
01148     Source_Base<EvalT,Traits> *sb = q;
01149     m_sources.push_back(sb);
01150     this->setName("QuadraticSource"+PHX::TypeString<EvalT>::value);
01151   }
01152   if (NeutronFission<EvalT,Traits>::check_for_existance(source_list)) {
01153     NeutronFission<EvalT,Traits>    *q = new NeutronFission<EvalT,Traits>(p);
01154     Source_Base<EvalT,Traits> *sb = q;
01155     m_sources.push_back(sb);
01156     this->setName("NeutronFissionSource"+PHX::TypeString<EvalT>::value);
01157   }
01158   if (MVQuadratic<EvalT,Traits>::check_for_existance(source_list)) {
01159     MVQuadratic<EvalT,Traits> *q = new MVQuadratic<EvalT,Traits>(p);
01160     Source_Base<EvalT,Traits> *sb = q;
01161     m_sources.push_back(sb);
01162     this->setName("MVQuadraticSource"+PHX::TypeString<EvalT>::value);
01163   }
01164   if (MVExponential<EvalT,Traits>::check_for_existance(source_list)) {
01165     MVExponential<EvalT,Traits> *q = new MVExponential<EvalT,Traits>(p);
01166     Source_Base<EvalT,Traits> *sb = q;
01167     m_sources.push_back(sb);
01168     this->setName("MVExponentialSource"+PHX::TypeString<EvalT>::value);
01169   }
01170   if (PointSource<EvalT,Traits>::check_for_existance(source_list)) {
01171     PointSource<EvalT,Traits>       *s = new PointSource<EvalT,Traits>(source_list);
01172     Source_Base<EvalT,Traits> *sb = s;
01173     m_sources.push_back(sb);
01174     this->setName("PointSource"+PHX::TypeString<EvalT>::value);
01175   }
01176   for (std::size_t i=0; i<m_sources.size(); ++i) {
01177     Source_Base<EvalT,Traits>* sb =  m_sources[i];
01178     sb->DependentFields(*this, p);
01179     sb->EvaluatedFields(*this, p);
01180   }
01181 
01182 }
01183 
01184 template<typename EvalT, typename Traits>
01185 Source<EvalT, Traits>::~Source() {
01186   for (std::size_t i=0; i<m_sources.size(); ++i) {
01187     Source_Base<EvalT,Traits>* sb =  m_sources[i];
01188     delete sb;
01189   }
01190   m_sources.clear();
01191 }
01192 
01193 //**********************************************************************
01194 template<typename EvalT, typename Traits>
01195 void Source<EvalT, Traits>::
01196 postRegistrationSetup(typename Traits::SetupData d,
01197                       PHX::FieldManager<Traits>& fm)
01198 {
01199   for (std::size_t i=0; i<m_sources.size(); ++i) {
01200     Source_Base<EvalT,Traits>* sb =  m_sources[i];
01201     sb->FieldData(this->utils, fm);
01202   }
01203 }
01204 
01205 //**********************************************************************
01206 template<typename EvalT, typename Traits>
01207 void Source<EvalT, Traits>::evaluateFields(typename Traits::EvalData workset)
01208 {
01209   for (std::size_t i=0; i<m_sources.size(); ++i) {
01210     Source_Base<EvalT,Traits>* sb =  m_sources[i];
01211     sb->evaluateFields(workset);
01212   }
01213   return;
01214 }
01215 
01216 //**********************************************************************
01217 }
01218
PHAL_Source_Def.hpp
