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 "Teuchos_TestForException.hpp" 00008 #include "Phalanx_DataLayout.hpp" 00009 00010 #include "Intrepid_FunctionSpaceTools.hpp" 00011 00012 namespace LCM { 00013 00014 //********************************************************************** 00015 template<typename EvalT, typename Traits> 00016 TotalStress<EvalT, Traits>:: 00017 TotalStress(const Teuchos::ParameterList& p) : 00018 effStress (p.get<std::string> ("Effective Stress Name"), 00019 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ), 00020 biotCoefficient (p.get<std::string> ("Biot Coefficient Name"), 00021 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 00022 porePressure (p.get<std::string> ("QP Variable Name"), 00023 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 00024 stress (p.get<std::string> ("Total Stress Name"), 00025 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ) 00026 { 00027 // Pull out numQPs and numDims from a Layout 00028 Teuchos::RCP<PHX::DataLayout> tensor_dl = 00029 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout"); 00030 std::vector<PHX::DataLayout::size_type> dims; 00031 tensor_dl->dimensions(dims); 00032 numQPs = dims[1]; 00033 numDims = dims[2]; 00034 00035 this->addDependentField(effStress); 00036 this->addDependentField(biotCoefficient); 00037 this->addDependentField(porePressure); 00038 00039 this->addEvaluatedField(stress); 00040 00041 this->setName("TotalStress"+PHX::TypeString<EvalT>::value); 00042 00043 } 00044 00045 //********************************************************************** 00046 template<typename EvalT, typename Traits> 00047 void TotalStress<EvalT, Traits>:: 00048 postRegistrationSetup(typename Traits::SetupData d, 00049 PHX::FieldManager<Traits>& fm) 00050 { 00051 this->utils.setFieldData(stress,fm); 00052 this->utils.setFieldData(effStress,fm); 00053 this->utils.setFieldData(biotCoefficient,fm); 00054 this->utils.setFieldData(porePressure,fm); 00055 } 00056 00057 //********************************************************************** 00058 template<typename EvalT, typename Traits> 00059 void TotalStress<EvalT, Traits>:: 00060 evaluateFields(typename Traits::EvalData workset) 00061 { 00062 00063 for (std::size_t cell=0; cell < workset.numCells; ++cell) { 00064 for (std::size_t qp=0; qp < numQPs; ++qp) { 00065 for (std::size_t dim=0; dim<numDims; ++ dim) { 00066 for (std::size_t j=0; j<numDims; ++ j) { 00067 stress(cell,qp,dim,j) = effStress(cell,qp, dim,j); 00068 } 00069 } 00070 } 00071 } 00072 00073 for (std::size_t cell=0; cell < workset.numCells; ++cell) { 00074 for (std::size_t qp=0; qp < numQPs; ++qp) { 00075 for (std::size_t dim=0; dim<numDims; ++ dim) { 00076 stress(cell,qp,dim,dim) -= biotCoefficient(cell,qp)* 00077 porePressure(cell,qp); 00078 } 00079 } 00080 } 00081 /* ScalarT lambda, mu; // B is the Biot coefficient 1 - K/K(s) where 00082 00083 switch (numDims) { 00084 case 1: 00085 Intrepid::FunctionSpaceTools::tensorMultiplyDataData<ScalarT>(stress, elasticModulus, strain); 00086 for (std::size_t cell=0; cell < workset.numCells; ++cell) { 00087 for (std::size_t qp=0; qp < numQPs; ++qp) { 00088 stress(cell, qp) = stress(cell, qp) - porePressure(cell,qp); 00089 } 00090 } 00091 break; 00092 case 2: 00093 // Compute Stress (with the plane strain assumption for now) 00094 for (std::size_t cell=0; cell < workset.numCells; ++cell) { 00095 for (std::size_t qp=0; qp < numQPs; ++qp) { 00096 lambda = ( elasticModulus(cell,qp) * poissonsRatio(cell,qp) ) / ( ( 1 + poissonsRatio(cell,qp) ) * ( 1 - 2 * poissonsRatio(cell,qp) ) ); 00097 mu = elasticModulus(cell,qp) / ( 2 * ( 1 + poissonsRatio(cell,qp) ) ); 00098 stress(cell,qp,0,0) = 2.0 * mu * ( strain(cell,qp,0,0) ) + lambda * ( strain(cell,qp,0,0) + strain(cell,qp,1,1) ) - biotCoefficient(cell,qp)*porePressure(cell,qp) ; 00099 stress(cell,qp,1,1) = 2.0 * mu * ( strain(cell,qp,1,1) ) + lambda * ( strain(cell,qp,0,0) + strain(cell,qp,1,1) ) - biotCoefficient(cell,qp)*porePressure(cell,qp); 00100 stress(cell,qp,0,1) = 2.0 * mu * ( strain(cell,qp,0,1) ); 00101 stress(cell,qp,1,0) = stress(cell,qp,0,1); 00102 } 00103 } 00104 break; 00105 case 3: 00106 // Compute Stress 00107 for (std::size_t cell=0; cell < workset.numCells; ++cell) { 00108 for (std::size_t qp=0; qp < numQPs; ++qp) { 00109 lambda = ( elasticModulus(cell,qp) * poissonsRatio(cell,qp) ) / ( ( 1 + poissonsRatio(cell,qp) ) * ( 1 - 2 * poissonsRatio(cell,qp) ) ); 00110 mu = elasticModulus(cell,qp) / ( 2 * ( 1 + poissonsRatio(cell,qp) ) ); 00111 stress(cell,qp,0,0) = 2.0 * mu * ( strain(cell,qp,0,0) ) + lambda * ( strain(cell,qp,0,0) + strain(cell,qp,1,1) + strain(cell,qp,2,2) ) 00112 - biotCoefficient(cell,qp)*porePressure(cell,qp); 00113 stress(cell,qp,1,1) = 2.0 * mu * ( strain(cell,qp,1,1) ) + lambda * ( strain(cell,qp,0,0) + strain(cell,qp,1,1) + strain(cell,qp,2,2) ) 00114 - biotCoefficient(cell,qp)*porePressure(cell,qp); 00115 stress(cell,qp,2,2) = 2.0 * mu * ( strain(cell,qp,2,2) ) + lambda * ( strain(cell,qp,0,0) + strain(cell,qp,1,1) + strain(cell,qp,2,2) ) 00116 - biotCoefficient(cell,qp)*porePressure(cell,qp); 00117 stress(cell,qp,0,1) = 2.0 * mu * ( strain(cell,qp,0,1) ); 00118 stress(cell,qp,1,2) = 2.0 * mu * ( strain(cell,qp,1,2) ); 00119 stress(cell,qp,2,0) = 2.0 * mu * ( strain(cell,qp,2,0) ); 00120 stress(cell,qp,1,0) = stress(cell,qp,0,1); 00121 stress(cell,qp,2,1) = stress(cell,qp,1,2); 00122 stress(cell,qp,0,2) = stress(cell,qp,2,0); 00123 } 00124 } 00125 break; 00126 } 00127 */ 00128 } 00129 00130 //********************************************************************** 00131 } 00132
1.7.1