Albany: a Trilinos-based PDE code: QCAD_ResponseCenterOfMass

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 "Phalanx_DataLayout.hpp"
00010 #include "Teuchos_CommHelpers.hpp"
00011 #include "Phalanx.hpp"
00012 
00013 template<typename EvalT, typename Traits>
00014 QCAD::ResponseCenterOfMass<EvalT, Traits>::
00015 ResponseCenterOfMass(Teuchos::ParameterList& p,
00016          const Teuchos::RCP<Albany::Layouts>& dl) :
00017   coordVec("Coord Vec", dl->qp_vector),
00018   weights("Weights", dl->qp_scalar)
00019 {
00020   // get and validate Response parameter list
00021   Teuchos::ParameterList* plist = 
00022     p.get<Teuchos::ParameterList*>("Parameter List");
00023   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00024     this->getValidResponseParameters();
00025   plist->validateParameters(*reflist,0);
00026 
00027   // number of quad points per cell and dimension of space
00028   Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00029   Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00030   
00031   std::vector<PHX::DataLayout::size_type> dims;
00032   vector_dl->dimensions(dims);
00033   numQPs  = dims[1];
00034   numDims = dims[2];
00035 
00037   Teuchos::RCP<QCAD::MaterialDatabase> materialDB;
00038   Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem = 
00039     p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");
00040   if(paramsFromProblem != Teuchos::null)
00041     materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
00042   else materialDB = Teuchos::null;
00043   
00044   // User-specified parameters
00045   fieldName = plist->get<std::string>("Field Name");
00046   opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );
00047   
00048   // setup field
00049   PHX::MDField<ScalarT> f(fieldName, scalar_dl); field = f;
00050 
00051   // add dependent fields
00052   this->addDependentField(field);
00053   this->addDependentField(coordVec);
00054   this->addDependentField(weights);
00055   opRegion->addDependentFields(this);
00056 
00057   this->setName(fieldName+" Response Center of Mass"+PHX::TypeString<EvalT>::value);
00058 
00059   using PHX::MDALayout;
00060 
00061   // Setup scatter evaluator
00062   p.set("Stand-alone Evaluator", false);
00063   std::string local_response_name = 
00064     fieldName + " Local Response Center of Mass";
00065   std::string global_response_name = 
00066     fieldName + " Global Response Center of Mass";
00067   int worksetSize = scalar_dl->dimension(0);
00068   int responseSize = 4;
00069   Teuchos::RCP<PHX::DataLayout> local_response_layout =
00070     Teuchos::rcp(new MDALayout<Cell,Dim>(worksetSize, responseSize));
00071   Teuchos::RCP<PHX::DataLayout> global_response_layout =
00072     Teuchos::rcp(new MDALayout<Dim>(responseSize));
00073   PHX::Tag<ScalarT> local_response_tag(local_response_name, 
00074                local_response_layout);
00075   PHX::Tag<ScalarT> global_response_tag(global_response_name, 
00076           global_response_layout);
00077   p.set("Local Response Field Tag", local_response_tag);
00078   p.set("Global Response Field Tag", global_response_tag);
00079   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);
00080 }
00081 
00082 // **********************************************************************
00083 template<typename EvalT, typename Traits>
00084 void QCAD::ResponseCenterOfMass<EvalT, Traits>::
00085 postRegistrationSetup(typename Traits::SetupData d,
00086                       PHX::FieldManager<Traits>& fm)
00087 {
00088   this->utils.setFieldData(field,fm);
00089   this->utils.setFieldData(coordVec,fm);
00090   this->utils.setFieldData(weights,fm);
00091   opRegion->postRegistrationSetup(fm);
00092   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postRegistrationSetup(d,fm);
00093 }
00094 
00095 // **********************************************************************
00096 template<typename EvalT, typename Traits>
00097 void QCAD::ResponseCenterOfMass<EvalT, Traits>::
00098 preEvaluate(typename Traits::PreEvalData workset)
00099 {
00100   for (typename PHX::MDField<ScalarT>::size_type i=0; 
00101        i<this->global_response.size(); i++)
00102     this->global_response[i] = 0.0;
00103 
00104   // Do global initialization
00105   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::preEvaluate(workset);
00106 }
00107 
00108 // **********************************************************************
00109 template<typename EvalT, typename Traits>
00110 void QCAD::ResponseCenterOfMass<EvalT, Traits>::
00111 evaluateFields(typename Traits::EvalData workset)
00112 {
00113   // Zero out local response
00114   for (typename PHX::MDField<ScalarT>::size_type i=0; 
00115        i<this->local_response.size(); i++)
00116     this->local_response[i] = 0.0;
00117 
00118   ScalarT integral, moment;
00119 
00120   if(!opRegion->elementBlockIsInRegion(workset.EBName))
00121     return;
00122 
00123   for (std::size_t cell=0; cell < workset.numCells; ++cell) 
00124   {
00125     if(!opRegion->cellIsInRegion(cell)) continue;
00126 
00127     // Add to running total volume and mass moment
00128     for (std::size_t qp=0; qp < numQPs; ++qp) {
00129       integral = field(cell,qp) * weights(cell,qp);
00130       this->local_response(cell,3) += integral;
00131       this->global_response(3) += integral;
00132 
00133       for(std::size_t i=0; i<numDims && i<3; i++) {
00134   moment = field(cell,qp) * weights(cell,qp) * coordVec(cell,qp,i);
00135   this->local_response(cell,i) += moment;
00136   this->global_response(i) += moment;
00137       }
00138     }
00139 
00140   }
00141 
00142   // Do any local-scattering necessary
00143   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::evaluateFields(workset);
00144 }
00145 
00146 // **********************************************************************
00147 template<typename EvalT, typename Traits>
00148 void QCAD::ResponseCenterOfMass<EvalT, Traits>::
00149 postEvaluate(typename Traits::PostEvalData workset)
00150 {
00151   // Add contributions across processors
00152   Teuchos::RCP< Teuchos::ValueTypeSerializer<int,ScalarT> > serializer =
00153     workset.serializerManager.template getValue<EvalT>();
00154   Teuchos::reduceAll(
00155     *workset.comm, *serializer, Teuchos::REDUCE_SUM,
00156     this->global_response.size(), &this->global_response[0], 
00157     &this->global_response[0]);
00158 
00159   int iNormalizer = 3;
00160   if( fabs(this->global_response[iNormalizer]) > 1e-9 ) {
00161     for( int i=0; i < this->global_response.size(); i++) {
00162       if( i == iNormalizer ) continue;
00163       this->global_response[i] /= this->global_response[iNormalizer];
00164     }
00165     this->global_response[iNormalizer] = 1.0;
00166   }
00167 
00168   // Do global scattering
00169   PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postEvaluate(workset);
00170 }
00171 
00172 // **********************************************************************
00173 template<typename EvalT,typename Traits>
00174 Teuchos::RCP<const Teuchos::ParameterList>
00175 QCAD::ResponseCenterOfMass<EvalT,Traits>::getValidResponseParameters() const
00176 {
00177   Teuchos::RCP<Teuchos::ParameterList> validPL =
00178       rcp(new Teuchos::ParameterList("Valid ResponseCenterOfMass Params"));;
00179   Teuchos::RCP<const Teuchos::ParameterList> baseValidPL =
00180     PHAL::SeparableScatterScalarResponse<EvalT,Traits>::getValidResponseParameters();
00181   validPL->setParameters(*baseValidPL);
00182 
00183   Teuchos::RCP<const Teuchos::ParameterList> regionValidPL =
00184     QCAD::MeshRegion<EvalT,Traits>::getValidParameters();
00185   validPL->setParameters(*regionValidPL);
00186 
00187   validPL->set<std::string>("Name", "", "Name of response function");
00188   validPL->set<int>("Phalanx Graph Visualization Detail", 0, "Make dot file to visualize phalanx graph");
00189   validPL->set<std::string>("Field Name", "", "Scalar field from which to compute center of mass");
00190   validPL->set<std::string>("Description", "", "Description of this response used by post processors");
00191 
00192   return validPL;
00193 }
00194 
00195 // **********************************************************************
00196
QCAD_ResponseCenterOfMass_Def.hpp
