Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009
00010 #include "Intrepid_FunctionSpaceTools.hpp"
00011
00012 namespace PHAL {
00013
00014
00015 template<typename EvalT, typename Traits>
00016 NSBodyForce<EvalT, Traits>::
00017 NSBodyForce(const Teuchos::ParameterList& p) :
00018 force(p.get<std::string>("Body Force Name"),
00019 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00020 haveHeat(p.get<bool>("Have Heat"))
00021 {
00022
00023 Teuchos::ParameterList* bf_list =
00024 p.get<Teuchos::ParameterList*>("Parameter List");
00025
00026 std::string type = bf_list->get("Type", "None");
00027 if (type == "None") {
00028 bf_type = NONE;
00029 }
00030 else if (type == "Constant") {
00031 bf_type = CONSTANT;
00032 rho = PHX::MDField<ScalarT,Cell,QuadPoint>(
00033 p.get<std::string>("Density QP Variable Name"),
00034 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00035 this->addDependentField(rho);
00036 }
00037 else if (type == "Boussinesq") {
00038 TEUCHOS_TEST_FOR_EXCEPTION(haveHeat == false, std::logic_error,
00039 std::endl <<
00040 "Error! Must enable heat equation for Boussinesq " <<
00041 "body force term!");
00042 bf_type = BOUSSINESQ;
00043 T = PHX::MDField<ScalarT,Cell,QuadPoint>(
00044 p.get<std::string>("Temperature QP Variable Name"),
00045 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00046 rho = PHX::MDField<ScalarT,Cell,QuadPoint>(
00047 p.get<std::string>("Density QP Variable Name"),
00048 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00049 beta = PHX::MDField<ScalarT,Cell,QuadPoint>(
00050 p.get<std::string>(
00051 "Volumetric Expansion Coefficient QP Variable Name"),
00052 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00053 this->addDependentField(rho);
00054 this->addDependentField(beta);
00055 this->addDependentField(T);
00056 }
00057
00058 this->addEvaluatedField(force);
00059
00060 Teuchos::RCP<PHX::DataLayout> vector_dl =
00061 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00062 std::vector<PHX::DataLayout::size_type> dims;
00063 vector_dl->dimensions(dims);
00064 numQPs = dims[1];
00065 numDims = dims[2];
00066
00067 if (bf_type == CONSTANT || bf_type == BOUSSINESQ) {
00068 if (bf_list->isType<Teuchos::Array<double> >("Gravity Vector"))
00069 gravity = bf_list->get<Teuchos::Array<double> >("Gravity Vector");
00070 else {
00071 gravity.resize(numDims);
00072 gravity[numDims-1] = -1.0;
00073 }
00074 }
00075 this->setName("NSBodyForce"+PHX::TypeString<EvalT>::value);
00076 }
00077
00078
00079 template<typename EvalT, typename Traits>
00080 void NSBodyForce<EvalT, Traits>::
00081 postRegistrationSetup(typename Traits::SetupData d,
00082 PHX::FieldManager<Traits>& fm)
00083 {
00084 if (bf_type == CONSTANT) {
00085 this->utils.setFieldData(rho,fm);
00086 }
00087 else if (bf_type == BOUSSINESQ) {
00088 this->utils.setFieldData(T,fm);
00089 this->utils.setFieldData(rho,fm);
00090 this->utils.setFieldData(beta,fm);
00091 }
00092
00093 this->utils.setFieldData(force,fm);
00094 }
00095
00096
00097 template<typename EvalT, typename Traits>
00098 void NSBodyForce<EvalT, Traits>::
00099 evaluateFields(typename Traits::EvalData workset)
00100 {
00101 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00102 for (std::size_t qp=0; qp < numQPs; ++qp) {
00103 for (std::size_t i=0; i < numDims; ++i) {
00104 if (bf_type == NONE)
00105 force(cell,qp,i) = 0.0;
00106 else if (bf_type == CONSTANT)
00107 force(cell,qp,i) = rho(cell,qp)*gravity[i];
00108 else if (bf_type == BOUSSINESQ)
00109 force(cell,qp,i) = rho(cell,qp)*T(cell,qp)*beta(cell,qp)*gravity[i];
00110 }
00111 }
00112 }
00113 }
00114
00115 }
00116
