00001
00002
00003
00004
00005
00006
00007 #include <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 #include "Sacado_ParameterRegistration.hpp"
00011 #include "Albany_Utils.hpp"
00012 #include <typeinfo>
00013
00014 namespace LCM {
00015
00016 template<typename EvalT, typename Traits>
00017 YieldStrength<EvalT, Traits>::
00018 YieldStrength(Teuchos::ParameterList& p) :
00019 yieldStrength(p.get<std::string>("QP Variable Name"),
00020 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout"))
00021 {
00022 Teuchos::ParameterList* elmd_list =
00023 p.get<Teuchos::ParameterList*>("Parameter List");
00024
00025 Teuchos::RCP<PHX::DataLayout> vector_dl =
00026 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00027 std::vector<PHX::DataLayout::size_type> dims;
00028 vector_dl->dimensions(dims);
00029 numQPs = dims[1];
00030 numDims = dims[2];
00031
00032 Teuchos::RCP<ParamLib> paramLib =
00033 p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00034
00035 std::string type = elmd_list->get("Yield Strength Type", "Constant");
00036 if (type == "Constant") {
00037 is_constant = true;
00038 constant_value = elmd_list->get("Value", 1.0);
00039
00040
00041 new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00042 "Yield Strength", this, paramLib);
00043 }
00044 else if (type == "Truncated KL Expansion") {
00045 is_constant = false;
00046 PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>
00047 fx(p.get<std::string>("QP Coordinate Vector Name"), vector_dl);
00048 coordVec = fx;
00049 this->addDependentField(coordVec);
00050
00051 exp_rf_kl =
00052 Teuchos::rcp(new Stokhos::KL::ExponentialRandomField<MeshScalarT>(*elmd_list));
00053 int num_KL = exp_rf_kl->stochasticDimension();
00054
00055
00056 rv.resize(num_KL);
00057 for (int i=0; i<num_KL; i++) {
00058 std::string ss = Albany::strint("Yield Strength KL Random Variable",i);
00059 new Sacado::ParameterRegistration<EvalT, SPL_Traits>(ss, this, paramLib);
00060 rv[i] = elmd_list->get(ss, 0.0);
00061 }
00062 }
00063 else {
00064 TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00065 "Invalid yield strength type " << type);
00066 }
00067
00068
00069
00070
00071 if ( p.isType<std::string>("QP Temperature Name") ) {
00072 Teuchos::RCP<PHX::DataLayout> scalar_dl =
00073 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00074 PHX::MDField<ScalarT,Cell,QuadPoint>
00075 tmp(p.get<std::string>("QP Temperature Name"), scalar_dl);
00076 Temperature = tmp;
00077 this->addDependentField(Temperature);
00078 isThermoElastic = true;
00079 dYdT_value = elmd_list->get("dYdT Value", 0.0);
00080 refTemp = p.get<RealType>("Reference Temperature", 0.0);
00081 new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00082 "dYdT Value", this, paramLib);
00083 }
00084 else {
00085 isThermoElastic=false;
00086 dYdT_value=0.0;
00087 }
00088
00089 if ( p.isType<std::string>("Lattice Concentration Name") ) {
00090 Teuchos::RCP<PHX::DataLayout> scalar_dl =
00091 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout");
00092 PHX::MDField<ScalarT,Cell,QuadPoint>
00093 tmp(p.get<std::string>("Lattice Concentration Name"), scalar_dl);
00094 CL = tmp;
00095 this->addDependentField(CL);
00096 CLname = p.get<std::string>("Lattice Concentration Name")+"_old";
00097
00098
00099
00100
00101
00102
00103
00104 isDiffuseDeformation = true;
00105 zeta = elmd_list->get("zeta Value", 1.0);
00106 new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00107 "zeta Value", this, paramLib);
00108 }
00109 else {
00110 isDiffuseDeformation=false;
00111 zeta=1.0;
00112 }
00113
00114
00115 this->addEvaluatedField(yieldStrength);
00116 this->setName("Yield Strength"+PHX::TypeString<EvalT>::value);
00117 }
00118
00119
00120 template<typename EvalT, typename Traits>
00121 void YieldStrength<EvalT, Traits>::
00122 postRegistrationSetup(typename Traits::SetupData d,
00123 PHX::FieldManager<Traits>& fm)
00124 {
00125 this->utils.setFieldData(yieldStrength,fm);
00126 if (!is_constant) this->utils.setFieldData(coordVec,fm);
00127 if (isThermoElastic) this->utils.setFieldData(Temperature,fm);
00128 if (isDiffuseDeformation) this->utils.setFieldData(CL,fm);
00129
00130 }
00131
00132
00133 template<typename EvalT, typename Traits>
00134 void YieldStrength<EvalT, Traits>::
00135 evaluateFields(typename Traits::EvalData workset)
00136 {
00137 bool print = false;
00138
00139
00140 if (print)
00141 std::cout << " *** YieldStrength *** " << std::endl;
00142
00143 std::size_t numCells = workset.numCells;
00144
00145 if (is_constant) {
00146 for (std::size_t cell=0; cell < numCells; ++cell) {
00147 for (std::size_t qp=0; qp < numQPs; ++qp) {
00148 yieldStrength(cell,qp) = constant_value;
00149 }
00150 }
00151 }
00152 else {
00153 for (std::size_t cell=0; cell < numCells; ++cell) {
00154 for (std::size_t qp=0; qp < numQPs; ++qp) {
00155 Teuchos::Array<MeshScalarT> point(numDims);
00156 for (std::size_t i=0; i<numDims; i++)
00157 point[i] = Sacado::ScalarValue<MeshScalarT>::eval(coordVec(cell,qp,i));
00158 yieldStrength(cell,qp) = exp_rf_kl->evaluate(point, rv);
00159 }
00160 }
00161 }
00162 if (isThermoElastic) {
00163 for (std::size_t cell=0; cell < numCells; ++cell) {
00164 for (std::size_t qp=0; qp < numQPs; ++qp) {
00165 yieldStrength(cell,qp) -= dYdT_value * (Temperature(cell,qp) - refTemp);
00166
00167 if (print)
00168 {
00169 std::cout << " Y : " << yieldStrength(cell,qp) << std::endl;
00170 std::cout << " temp: " << Temperature(cell,qp) << std::endl;
00171 std::cout << " dYdT: " << dYdT_value << std::endl;
00172 std::cout << " refT: " << refTemp << std::endl;
00173 }
00174 }
00175 }
00176 }
00177 if (isDiffuseDeformation) {
00178
00179 Albany::MDArray CLold = (*workset.stateArrayPtr)[CLname];
00180
00181 for (std::size_t cell=0; cell < numCells; ++cell) {
00182 for (std::size_t qp=0; qp < numQPs; ++qp) {
00183
00184 yieldStrength(cell,qp) -= constant_value*(zeta-1.0)*(CL(cell,qp) -CLold(cell,qp) );
00185
00186 if (print)
00187 {
00188 std::cout << " Y : " << yieldStrength(cell,qp) << std::endl;
00189 std::cout << " CT : " << CT(cell,qp) << std::endl;
00190 std::cout << " zeta : " << zeta << std::endl;
00191 }
00192 }
00193 }
00194 }
00195 }
00196
00197
00198 template<typename EvalT,typename Traits>
00199 typename YieldStrength<EvalT,Traits>::ScalarT&
00200 YieldStrength<EvalT,Traits>::getValue(const std::string &n)
00201 {
00202 if (n == "Yield Strength")
00203 return constant_value;
00204 else if (n == "dYdT Value")
00205 return dYdT_value;
00206 else if (n == "zeta Value")
00207 return zeta;
00208 for (int i=0; i<rv.size(); i++) {
00209 if (n == Albany::strint("Yield Strength KL Random Variable",i))
00210 return rv[i];
00211 }
00212 TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00213 std::endl <<
00214 "Error! Logic error in getting paramter " << n
00215 << " in YieldStrength::getValue()" << std::endl);
00216 return constant_value;
00217 }
00218
00219
00220
00221 }
00222