• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

TvergaardHutchinson_Def.hpp

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 #include <Intrepid_MiniTensor.h>
00007 
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 
00011 namespace LCM{
00012   //**********************************************************************
00013   template<typename EvalT, typename Traits>
00014   TvergaardHutchinson<EvalT, Traits>::
00015   TvergaardHutchinson(const Teuchos::ParameterList& p,
00016                       const Teuchos::RCP<Albany::Layouts>& dl) :
00017     cubature        (p.get<Teuchos::RCP<Intrepid::Cubature<RealType> > >("Cubature")),
00018     intrepidBasis   (p.get<Teuchos::RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >("Intrepid Basis")),
00019     jump(p.get<std::string>("Vector Jump Name"),dl->qp_vector),
00020     currentBasis(p.get<std::string>("Current Basis Name"),dl->qp_tensor),
00021     cohesiveTraction(p.get<std::string>("Cohesive Traction Name"),dl->qp_vector),
00022     delta_1(p.get<RealType>("delta_1 Name")),
00023     delta_2(p.get<RealType>("delta_2 Name")),
00024     delta_c(p.get<RealType>("delta_c Name")),
00025     sigma_c(p.get<RealType>("sigma_c Name")),
00026     beta_0(p.get<RealType>("beta_0 Name")),
00027     beta_1(p.get<RealType>("beta_1 Name")),
00028     beta_2(p.get<RealType>("beta_2 Name"))
00029   {
00030     this->addDependentField(jump);
00031     this->addDependentField(currentBasis);
00032     this->addEvaluatedField(cohesiveTraction);
00033 
00034     this->setName("TvergaardHutchinson" + PHX::TypeString<EvalT>::value);
00035 
00036     // get dimension
00037     std::vector<PHX::DataLayout::size_type> dims;
00038     dl->qp_vector->dimensions(dims);
00039     worksetSize = dims[0];
00040     numDims = dims[2];
00041 
00042     numQPs = cubature->getNumPoints();
00043 
00044   }
00045 
00046   //**********************************************************************
00047   template<typename EvalT, typename Traits>
00048   void TvergaardHutchinson<EvalT, Traits>::
00049   postRegistrationSetup(typename Traits::SetupData d,
00050                         PHX::FieldManager<Traits>& fm)
00051   {
00052     this->utils.setFieldData(jump,fm);
00053     this->utils.setFieldData(currentBasis,fm);
00054     this->utils.setFieldData(cohesiveTraction,fm);
00055   }
00056   //**********************************************************************
00057   template<typename EvalT, typename Traits>
00058   void TvergaardHutchinson<EvalT, Traits>::
00059   evaluateFields(typename Traits::EvalData workset)
00060   {
00061 
00062     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00063       for (std::size_t pt = 0; pt < numQPs; ++pt) {
00064 
00065         //current basis vector
00066         Intrepid::Vector<ScalarT> g_0(3, &currentBasis(cell, pt, 0, 0));
00067         Intrepid::Vector<ScalarT> g_1(3, &currentBasis(cell, pt, 1, 0));
00068         Intrepid::Vector<ScalarT> n(3, &currentBasis(cell, pt, 2, 0));
00069 
00070         //construct orthogonal unit basis
00071         Intrepid::Vector<ScalarT> t_0(0,0,0), t_1(0,0,0);
00072         t_0 = g_0 / norm(g_0);
00073         t_1 = cross(n,t_0);
00074 
00075         //construct transformation matrix Q (2nd order tensor)
00076         Intrepid::Tensor<ScalarT> Q(3, Intrepid::ZEROS);
00077         // manually fill Q = [t_0; t_1; n];
00078         Q(0,0) = t_0(0); Q(1,0) = t_0(1);  Q(2,0) = t_0(2);
00079         Q(0,1) = t_1(0); Q(1,1) = t_1(1);  Q(2,1) = t_1(2);
00080         Q(0,2) = n(0);   Q(1,2) = n(1);    Q(2,2) = n(2);
00081 
00082         //global and local jump
00083         Intrepid::Vector<ScalarT> jump_global(3, &jump(cell, pt, 0));
00084         Intrepid::Vector<ScalarT> jump_local(3);
00085         jump_local = Intrepid::dot(Intrepid::transpose(Q), jump_global);
00086 
00087         // matrix beta that controls relative effect of shear and normal opening
00088         Intrepid::Tensor<ScalarT> beta(3, Intrepid::ZEROS);
00089         beta(0,0) = beta_0; beta(1,1) = beta_1; beta(2,2) = beta_2;
00090 
00091         // compute scalar effective jump
00092         ScalarT jump_eff, tmp2;
00093         Intrepid::Vector<ScalarT> tmp1;
00094         tmp1 = Intrepid::dot(beta,jump_local);
00095         tmp2 = Intrepid::dot(jump_local,tmp1);
00096 
00097         jump_eff =
00098             std::sqrt(Intrepid::dot(jump_local,Intrepid::dot(beta,jump_local)));
00099 
00100         // traction-separation law from Tvergaard-Hutchinson 1992
00101         ScalarT sigma_eff;
00102         //Sacado::ScalarValue<ScalarT>::eval
00103         if(jump_eff <= delta_1)
00104           sigma_eff = sigma_c * jump_eff / delta_1;
00105         else if(jump_eff > delta_1 && jump_eff <= delta_2)
00106           sigma_eff = sigma_c;
00107         else if(jump_eff > delta_2 && jump_eff <= delta_c)
00108           sigma_eff = sigma_c * (delta_c - jump_eff) / (delta_c - delta_2);
00109         else
00110           sigma_eff = 0.0;
00111 
00112         // construct traction vector
00113         Intrepid::Vector<ScalarT> traction_local(3);
00114         traction_local.clear();
00115         if(jump_eff != 0)
00116           traction_local = Intrepid::dot(beta,jump_local) * sigma_eff / jump_eff;
00117 
00118         // global traction vector
00119         Intrepid::Vector<ScalarT> traction_global(3);
00120         traction_global = Intrepid::dot(Q, traction_local);
00121 
00122         cohesiveTraction(cell,pt,0) = traction_global(0);
00123         cohesiveTraction(cell,pt,1) = traction_global(1);
00124         cohesiveTraction(cell,pt,2) = traction_global(2);
00125 
00126         // output for debug, I'll keep it for now until the implementation fully verified
00127         // Q.Chen
00128         //std::cout << "jump_eff " << Sacado::ScalarValue<ScalarT>::eval(jump_eff) << std::endl;
00129         //std::cout << "sigma_eff " << Sacado::ScalarValue<ScalarT>::eval(sigma_eff) << std::endl;
00130 
00131       }
00132     }
00133   }
00134   //**********************************************************************
00135 }

Generated on Wed Mar 26 2014 18:36:46 for Albany: a Trilinos-based PDE code by  doxygen 1.7.1