Albany: a Trilinos-based PDE code: IPtoNodalField

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 <fstream>
00008 #include <Teuchos_TestForException.hpp>
00009 #include "Albany_Utils.hpp"
00010 #include "Adapt_NodalDataBlock.hpp"
00011 
00012 namespace LCM
00013 {
00014 
00015 //------------------------------------------------------------------------------
00016 template<typename EvalT, typename Traits>
00017 IPtoNodalFieldBase<EvalT, Traits>::
00018 IPtoNodalFieldBase(Teuchos::ParameterList& p,
00019       const Teuchos::RCP<Albany::Layouts>& dl) :
00020   weights_("Weights", dl->qp_scalar),
00021   nodal_weights_name_("nodal_weights")
00022 {
00023 
00025   Teuchos::ParameterList* plist = 
00026     p.get<Teuchos::ParameterList*>("Parameter List");
00027   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00028     this->getValidIPtoNodalFieldParameters();
00029   plist->validateParameters(*reflist,0);
00030 
00031   output_to_exodus_ = plist->get<bool>("Output to File", true);
00032 
00034   Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00035   Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00036   Teuchos::RCP<PHX::DataLayout> cell_dl = dl->cell_scalar;
00037   Teuchos::RCP<PHX::DataLayout> node_dl = dl->node_qp_vector;
00038   Teuchos::RCP<PHX::DataLayout> vert_vector_dl = dl->vertices_vector;
00039   num_pts_ = vector_dl->dimension(1);
00040   num_dims_ = vector_dl->dimension(2);
00041   num_nodes_ = node_dl->dimension(1);
00042   num_vertices_ = vert_vector_dl->dimension(2);
00043 
00045   this->p_state_mgr_ = p.get< Albany::StateManager* >("State Manager Ptr");
00046 
00047   // register the nodal weights
00048   this->addDependentField(weights_);  
00049   this->p_state_mgr_->registerStateVariable(nodal_weights_name_,
00050                                             dl->node_node_scalar,
00051                                             dl->dummy, "all", 
00052                                             "scalar", 0.0, false,
00053                                             true);
00054 
00055   // loop over the number of fields and register
00056   number_of_fields_ = plist->get<int>("Number of Fields", 0);
00057 
00058   // resize field vectors
00059   ip_field_names_.resize(number_of_fields_);
00060   ip_field_layouts_.resize(number_of_fields_);
00061   nodal_field_names_.resize(number_of_fields_);
00062   ip_fields_.resize(number_of_fields_);
00063 
00064   for (int field(0); field < number_of_fields_; ++field) {
00065     ip_field_names_[field] = plist->get<std::string>(Albany::strint("IP Field Name", field));
00066     ip_field_layouts_[field] = plist->get<std::string>(Albany::strint("IP Field Layout", field));
00067     nodal_field_names_[field] = "nodal_" + ip_field_names_[field];
00068 
00069     if (ip_field_layouts_[field] == "Scalar") {
00070       PHX::MDField<ScalarT> s(ip_field_names_[field],dl->qp_scalar);
00071       ip_fields_[field] = s;
00072     } else if (ip_field_layouts_[field] == "Vector") {
00073       PHX::MDField<ScalarT> v(ip_field_names_[field],dl->qp_vector);
00074       ip_fields_[field] = v;
00075     } else if (ip_field_layouts_[field] == "Tensor") {
00076       PHX::MDField<ScalarT> t(ip_field_names_[field],dl->qp_tensor);
00077       ip_fields_[field] = t;
00078     } else {
00079       TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"Field Layout unknown");
00080     }
00081 
00082     this->addDependentField(ip_fields_[field]);
00083 
00084     if (ip_field_layouts_[field] == "Scalar" ) {
00085       this->p_state_mgr_->registerStateVariable(nodal_field_names_[field],
00086                                                 dl->node_node_scalar,
00087                                                 dl->dummy, "all", 
00088                                                 "scalar", 0.0, false,
00089                                                 output_to_exodus_);
00090     } else if (ip_field_layouts_[field] == "Vector" ) {
00091       this->p_state_mgr_->registerStateVariable(nodal_field_names_[field],
00092                                                 dl->node_node_vector,
00093                                                 dl->dummy, "all", 
00094                                                 "scalar", 0.0, false,
00095                                                 output_to_exodus_);
00096     } else if (ip_field_layouts_[field] == "Tensor" ) {
00097       this->p_state_mgr_->registerStateVariable(nodal_field_names_[field],
00098                                                 dl->node_node_tensor,
00099                                                 dl->dummy, "all", 
00100                                                 "scalar", 0.0, false,
00101                                                 output_to_exodus_);
00102     }
00103   }
00104 
00105   // Create field tag
00106   field_tag_ = 
00107     Teuchos::rcp(new PHX::Tag<ScalarT>("IP to Nodal Field", dl->dummy));
00108 
00109   this->addEvaluatedField(*field_tag_);
00110 }
00111 
00112 //------------------------------------------------------------------------------
00113 template<typename EvalT, typename Traits>
00114 void IPtoNodalFieldBase<EvalT, Traits>::
00115 postRegistrationSetup(typename Traits::SetupData d,
00116                       PHX::FieldManager<Traits>& fm)
00117 {
00118   this->utils.setFieldData(weights_,fm);
00119   for (int field(0); field < number_of_fields_; ++field) {
00120     this->utils.setFieldData(ip_fields_[field],fm);
00121   }
00122 }
00123 
00124 //------------------------------------------------------------------------------
00125 // Specialization: Residual
00126 //------------------------------------------------------------------------------
00127 //------------------------------------------------------------------------------
00128 template<typename Traits>
00129 IPtoNodalField<PHAL::AlbanyTraits::Residual, Traits>::
00130 IPtoNodalField(Teuchos::ParameterList& p, const Teuchos::RCP<Albany::Layouts>& dl) :
00131   IPtoNodalFieldBase<PHAL::AlbanyTraits::Residual, Traits>(p, dl)
00132 {
00133 }
00134 
00135 //------------------------------------------------------------------------------
00136 template<typename Traits>
00137 void IPtoNodalField<PHAL::AlbanyTraits::Residual, Traits>::
00138 preEvaluate(typename Traits::PreEvalData workset)
00139 {
00140   Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->p_state_mgr_->getStateInfoStruct()->getNodalDataBlock();
00141   node_data->initializeVectors(0.0);
00142 }
00143 
00144 //------------------------------------------------------------------------------
00145 template<typename Traits>
00146 void IPtoNodalField<PHAL::AlbanyTraits::Residual, Traits>::
00147 evaluateFields(typename Traits::EvalData workset)
00148 {
00149   // volume averaged field, store as nodal data that will be scattered
00150   // and summed
00151 
00152   // Get the node data block container
00153   Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->p_state_mgr_->getStateInfoStruct()->getNodalDataBlock();
00154   Teuchos::RCP<Epetra_Vector> data = node_data->getLocalNodeVec();
00155   Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >  wsElNodeID = workset.wsElNodeID;
00156   Teuchos::RCP<const Epetra_BlockMap> local_node_map = node_data->getLocalMap();
00157 
00158   int num_nodes = this->num_nodes_;
00159   int num_dims  = this->num_dims_;
00160   int num_pts   = this->num_pts_;
00161   int blocksize = node_data->getBlocksize();
00162 
00163   // deal with weights
00164   int  node_weight_offset;
00165   int  node_weight_ndofs;
00166   node_data->getNDofsAndOffset(this->nodal_weights_name_, node_weight_offset, node_weight_ndofs);
00167   for (int cell = 0; cell < workset.numCells; ++cell) {
00168     for (int node = 0; node < num_nodes; ++node) {
00169       int global_node = wsElNodeID[cell][node];
00170       int local_node = local_node_map->LID(global_node);
00171       if(local_node < 0) continue;
00172       for (int pt = 0; pt < num_pts; ++pt) {
00173         (*data)[local_node * blocksize + node_weight_offset] += this->weights_(cell,pt);
00174       }
00175     }
00176   }
00177   
00178   // deal with each of the fields
00179 
00180   for (int field(0); field < this->number_of_fields_; ++field) {
00181     int  node_var_offset;
00182     int  node_var_ndofs;
00183     node_data->getNDofsAndOffset(this->nodal_field_names_[field], node_var_offset, node_var_ndofs);
00184     for (int cell = 0; cell < workset.numCells; ++cell) {
00185       for (int node = 0; node < num_nodes; ++node) {
00186         int global_node = wsElNodeID[cell][node];
00187         int local_node = local_node_map->LID(global_node);
00188         if(local_node < 0) continue;
00189         for (int pt = 0; pt < num_pts; ++pt) {
00190           if (this->ip_field_layouts_[field] == "Scalar" ) {
00191             // save the scalar component
00192             (*data)[local_node * blocksize + node_var_offset] += 
00193               this->ip_fields_[field](cell,pt) * this->weights_(cell,pt);
00194           } else if (this->ip_field_layouts_[field] == "Vector" ) {
00195             for (int dim0 = 0; dim0 < num_dims; ++dim0) {
00196               // save the vector component
00197               (*data)[local_node * blocksize + node_var_offset + dim0] += 
00198                 this->ip_fields_[field](cell,pt,dim0) * this->weights_(cell,pt);
00199             }
00200           } else if (this->ip_field_layouts_[field] == "Tensor" ) {
00201             for (int dim0 = 0; dim0 < num_dims; ++dim0) {
00202               for (int dim1 = 0; dim1 < num_dims; ++dim1) {
00203                 // save the tensor component
00204                 (*data)[local_node * blocksize + node_var_offset + dim0*num_dims + dim1] += 
00205                   this->ip_fields_[field](cell,pt,dim0,dim1) * this->weights_(cell,pt);
00206               }
00207             }
00208           }
00209         }
00210       }    
00211     } // end cell loop
00212   } // end field loop
00213 } 
00214 //------------------------------------------------------------------------------
00215 template<typename Traits>
00216 void IPtoNodalField<PHAL::AlbanyTraits::Residual, Traits>::
00217 postEvaluate(typename Traits::PostEvalData workset)
00218 {
00219   // Note: we are in postEvaluate so all PEs call this
00220 
00221   // Get the node data block container
00222   Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->p_state_mgr_->getStateInfoStruct()->getNodalDataBlock();
00223   Teuchos::RCP<Epetra_Vector> data = node_data->getOverlapNodeVec();
00224   Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >  wsElNodeID = workset.wsElNodeID;
00225   Teuchos::RCP<const Epetra_BlockMap> overlap_node_map = node_data->getOverlapMap();
00226 
00227   // Build the exporter
00228   node_data->initializeExport();
00229 
00230   // do the export
00231   node_data->exportAddNodalDataBlock();
00232 
00233   int num_nodes = overlap_node_map->NumMyElements();
00234   int blocksize = node_data->getBlocksize();
00235 
00236   // get weight info
00237   int  node_weight_offset;
00238   int  node_weight_ndofs;
00239   node_data->getNDofsAndOffset(this->nodal_weights_name_, node_weight_offset, node_weight_ndofs);
00240 
00241   for (int field(0); field < this->number_of_fields_; ++field) {
00242     int  node_var_offset;
00243     int  node_var_ndofs;
00244     node_data->getNDofsAndOffset(this->nodal_field_names_[field], node_var_offset, node_var_ndofs);
00245 
00246     // all PEs divide the accumulated value(s) by the weights
00247     for (int overlap_node=0; overlap_node < num_nodes; ++overlap_node)
00248       for (int k=0; k < node_var_ndofs; ++k) 
00249         (*data)[overlap_node * blocksize + node_var_offset + k] /=
00250           (*data)[overlap_node * blocksize + node_weight_offset];
00251 
00252   }
00253   // Export the data from the local to overlapped decomposition
00254   // Divide the overlap field through by the weights
00255   // Store the overlapped vector data back in stk in the field "field_name"
00256   node_data->saveNodalDataState();
00257 }
00258 
00259 //------------------------------------------------------------------------------
00260 template<typename EvalT, typename Traits>
00261 Teuchos::RCP<const Teuchos::ParameterList>
00262 IPtoNodalFieldBase<EvalT,Traits>::getValidIPtoNodalFieldParameters() const
00263 {
00264   Teuchos::RCP<Teuchos::ParameterList> validPL =
00265       rcp(new Teuchos::ParameterList("Valid IPtoNodalField Params"));;
00266 
00267   validPL->set<std::string>("Name", "", "Name of field Evaluator");
00268   validPL->set<int>("Number of Fields", 0);
00269   validPL->set<std::string>("IP Field Name 0", "", "IP Field prefix");
00270   validPL->set<std::string>("IP Field Layout 0", "", "IP Field Layout: Scalar, Vector, or Tensor");
00271 
00272   validPL->set<std::string>("IP Field Name 1", "", "IP Field prefix");
00273   validPL->set<std::string>("IP Field Layout 1", "", "IP Field Layout: Scalar, Vector, or Tensor");
00274 
00275   validPL->set<std::string>("IP Field Name 2", "", "IP Field prefix");
00276   validPL->set<std::string>("IP Field Layout 2", "", "IP Field Layout: Scalar, Vector, or Tensor");
00277 
00278   validPL->set<std::string>("IP Field Name 3", "", "IP Field prefix");
00279   validPL->set<std::string>("IP Field Layout 3", "", "IP Field Layout: Scalar, Vector, or Tensor");
00280 
00281   validPL->set<std::string>("IP Field Name 4", "", "IP Field prefix");
00282   validPL->set<std::string>("IP Field Layout 4", "", "IP Field Layout: Scalar, Vector, or Tensor");
00283 
00284   validPL->set<std::string>("IP Field Name 5", "", "IP Field prefix");
00285   validPL->set<std::string>("IP Field Layout 5", "", "IP Field Layout: Scalar, Vector, or Tensor");
00286 
00287   validPL->set<std::string>("IP Field Name 6", "", "IP Field prefix");
00288   validPL->set<std::string>("IP Field Layout 6", "", "IP Field Layout: Scalar, Vector, or Tensor");
00289 
00290   validPL->set<std::string>("IP Field Name 7", "", "IP Field prefix");
00291   validPL->set<std::string>("IP Field Layout 7", "", "IP Field Layout: Scalar, Vector, or Tensor");
00292 
00293   validPL->set<std::string>("IP Field Name 8", "", "IP Field prefix");
00294   validPL->set<std::string>("IP Field Layout 8", "", "IP Field Layout: Scalar, Vector, or Tensor");
00295 
00296   validPL->set<std::string>("IP Field Name 9", "", "IP Field prefix");
00297   validPL->set<std::string>("IP Field Layout 9", "", "IP Field Layout: Scalar, Vector, or Tensor");
00298 
00299   validPL->set<bool>("Output to File", true, "Whether nodal field info should be output to a file");
00300   validPL->set<bool>("Generate Nodal Values", true, "Whether values at the nodes should be generated");
00301 
00302   return validPL;
00303 }
00304 
00305 //------------------------------------------------------------------------------
00306 }
00307
IPtoNodalField_Def.hpp
