Albany: a Trilinos-based PDE code: CapExplicit

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 "Teuchos_TestForException.hpp"
00007 #include "Phalanx_DataLayout.hpp"
00008 
00009 #include "Intrepid_FunctionSpaceTools.hpp"
00010 
00011 namespace LCM
00012 {
00013 
00014 //**********************************************************************
00015   template<typename EvalT, typename Traits>
00016   CapExplicit<EvalT, Traits>::
00017   CapExplicit(const Teuchos::ParameterList& p,
00018               const Teuchos::RCP<Albany::Layouts>& dl) :
00019     elasticModulus(p.get<std::string>("Elastic Modulus Name"),dl->qp_scalar),
00020     poissonsRatio(p.get<std::string>("Poissons Ratio Name"),dl->qp_scalar),
00021     strain(p.get<std::string>("Strain Name"),dl->qp_tensor),
00022     stress(p.get<std::string>("Stress Name"),dl->qp_tensor),
00023     backStress(p.get<std::string>("Back Stress Name"),dl->qp_tensor),
00024     capParameter(p.get<std::string>("Cap Parameter Name"),dl->qp_scalar),
00025     //friction(p.get<std::string>("Friction Name"),dl->qp_scalar),
00026     //dilatancy(p.get<std::string>("Dilatancy Name"),dl->qp_scalar),
00027     eqps(p.get<std::string>("Eqps Name"),dl->qp_scalar),
00028     //hardeningModulus(p.get<std::string>("Hardening Modulus Name"),dl->qp_scalar),
00029     volPlasticStrain(p.get<std::string>("Vol Plastic Strain Name"),dl->qp_scalar),
00030     A(p.get<RealType>("A Name")),
00031     B(p.get<RealType>("B Name")),
00032     C(p.get<RealType>("C Name")),
00033     theta(p.get<RealType>("Theta Name")),
00034     R(p.get<RealType>("R Name")),
00035     kappa0(p.get<RealType>("Kappa0 Name")),
00036     W(p.get<RealType>("W Name")),
00037     D1(p.get<RealType>("D1 Name")),
00038     D2(p.get<RealType>("D2 Name")),
00039     calpha(p.get<RealType>("Calpha Name")),
00040     psi(p.get<RealType>("Psi Name")),
00041     N(p.get<RealType>("N Name")),
00042     L(p.get<RealType>("L Name")),
00043     phi(p.get<RealType>("Phi Name")),
00044     Q(p.get<RealType>("Q Name"))
00045   {
00046     // Pull out numQPs and numDims from a Layout
00047     Teuchos::RCP<PHX::DataLayout> tensor_dl = p.get<
00048         Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout");
00049     std::vector<PHX::DataLayout::size_type> dims;
00050     tensor_dl->dimensions(dims);
00051     numQPs = dims[1];
00052     numDims = dims[2];
00053 
00054     this->addDependentField(elasticModulus);
00055     // PoissonRatio not used in 1D stress calc
00056     if (numDims > 1)
00057       this->addDependentField(poissonsRatio);
00058     this->addDependentField(strain);
00059 
00060     // state variable
00061     strainName = p.get<std::string>("Strain Name") + "_old";
00062     stressName = p.get<std::string>("Stress Name") + "_old";
00063     backStressName = p.get<std::string>("Back Stress Name") + "_old";
00064     capParameterName = p.get<std::string>("Cap Parameter Name") + "_old";
00065     eqpsName = p.get<std::string>("Eqps Name") + "_old";
00066     volPlasticStrainName = p.get<std::string>("Vol Plastic Strain Name")
00067         + "_old";
00068 
00069     // evaluated fields
00070     this->addEvaluatedField(stress);
00071     this->addEvaluatedField(backStress);
00072     this->addEvaluatedField(capParameter);
00073     //this->addEvaluatedField(friction);
00074     //this->addEvaluatedField(dilatancy);
00075     this->addEvaluatedField(eqps);
00076     //this->addEvaluatedField(hardeningModulus);
00077     this->addEvaluatedField(volPlasticStrain);
00078 
00079     this->setName("Stress" + PHX::TypeString<EvalT>::value);
00080 
00081     // initialize tensor
00082     I = Intrepid::eye<ScalarT>(numDims);
00083     id1 = Intrepid::identity_1<ScalarT>(numDims);
00084     id2 = Intrepid::identity_2<ScalarT>(numDims);
00085     id3 = Intrepid::identity_3<ScalarT>(numDims);
00086     Celastic = Intrepid::Tensor4<ScalarT>(numDims);
00087     compliance = Intrepid::Tensor4<ScalarT>(numDims);
00088     depsilon = Intrepid::Tensor<ScalarT>(numDims);
00089     sigmaN = Intrepid::Tensor<ScalarT>(numDims);
00090     strainN = Intrepid::Tensor<ScalarT>(numDims);
00091     sigmaVal = Intrepid::Tensor<ScalarT>(numDims);
00092     alphaVal = Intrepid::Tensor<ScalarT>(numDims);
00093     deps_plastic = Intrepid::Tensor<ScalarT>(numDims);
00094     sigmaTr = Intrepid::Tensor<ScalarT>(numDims);
00095     alphaTr = Intrepid::Tensor<ScalarT>(numDims);
00096     dfdsigma = Intrepid::Tensor<ScalarT>(numDims);
00097     dgdsigma = Intrepid::Tensor<ScalarT>(numDims);
00098     dfdalpha = Intrepid::Tensor<ScalarT>(numDims);
00099     halpha = Intrepid::Tensor<ScalarT>(numDims);
00100     dfdotCe = Intrepid::Tensor<ScalarT>(numDims);
00101     sigmaK = Intrepid::Tensor<ScalarT>(numDims);
00102     alphaK = Intrepid::Tensor<ScalarT>(numDims);
00103     dsigma = Intrepid::Tensor<ScalarT>(numDims);
00104     dev_plastic = Intrepid::Tensor<ScalarT>(numDims);
00105     xi = Intrepid::Tensor<ScalarT>(numDims);
00106     sN = Intrepid::Tensor<ScalarT>(numDims);
00107     s = Intrepid::Tensor<ScalarT>(numDims);
00108     strainCurrent = Intrepid::Tensor<ScalarT>(numDims);
00109     dJ3dsigma = Intrepid::Tensor<ScalarT>(numDims);
00110     eps_dev = Intrepid::Tensor<ScalarT>(numDims);
00111 
00112   }
00113 
00114 //**********************************************************************
00115   template<typename EvalT, typename Traits>
00116   void CapExplicit<EvalT, Traits>::postRegistrationSetup(
00117       typename Traits::SetupData d, PHX::FieldManager<Traits>& fm)
00118   {
00119     this->utils.setFieldData(elasticModulus, fm);
00120     if (numDims > 1)
00121       this->utils.setFieldData(poissonsRatio, fm);
00122     this->utils.setFieldData(stress, fm);
00123     this->utils.setFieldData(strain, fm);
00124     this->utils.setFieldData(backStress, fm);
00125     this->utils.setFieldData(capParameter, fm);
00126     //this->utils.setFieldData(friction, fm);
00127     //this->utils.setFieldData(dilatancy, fm);
00128     this->utils.setFieldData(eqps, fm);
00129     //this->utils.setFieldData(hardeningModulus, fm);
00130     this->utils.setFieldData(volPlasticStrain, fm);
00131 
00132   }
00133 
00134 //**********************************************************************
00135 
00136   template<typename EvalT, typename Traits>
00137   void CapExplicit<EvalT, Traits>::evaluateFields(
00138       typename Traits::EvalData workset)
00139   {
00140 
00141     // previous state
00142     Albany::MDArray strainold = (*workset.stateArrayPtr)[strainName];
00143     Albany::MDArray stressold = (*workset.stateArrayPtr)[stressName];
00144     Albany::MDArray backStressold = (*workset.stateArrayPtr)[backStressName];
00145     Albany::MDArray capParameterold = (*workset.stateArrayPtr)[capParameterName];
00146     Albany::MDArray eqpsold = (*workset.stateArrayPtr)[eqpsName];
00147     Albany::MDArray volPlasticStrainold =
00148         (*workset.stateArrayPtr)[volPlasticStrainName];
00149 
00150     ScalarT lame, mu, bulkModulus;
00151     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00152       for (std::size_t qp = 0; qp < numQPs; ++qp) {
00153         // local parameters
00154         lame = elasticModulus(cell, qp) * poissonsRatio(cell, qp)
00155             / (1.0 + poissonsRatio(cell, qp))
00156             / (1.0 - 2.0 * poissonsRatio(cell, qp));
00157         mu = elasticModulus(cell, qp) / 2.0
00158             / (1.0 + poissonsRatio(cell, qp));
00159         bulkModulus = lame + (2. / 3.) * mu;
00160 
00161         // elastic matrix
00162         Celastic = lame * id3 + mu * (id1 + id2);
00163 
00164         // elastic compliance tangent matrix
00165         compliance = (1. / bulkModulus / 9.) * id3
00166             + (1. / mu / 2.) * (0.5 * (id1 + id2)- (1. / 3.) * id3);
00167 
00168         // trial state
00169         Intrepid::Tensor<ScalarT> depsilon(3);
00170         for (std::size_t i = 0; i < numDims; ++i) {
00171           for (std::size_t j = 0; j < numDims; ++j) {
00172             depsilon(i, j) = strain(cell, qp, i, j) - strainold(cell, qp, i, j);
00173             strainN(i, j) = strainold(cell, qp, i, j);
00174             sigmaN(i,j) = stressold(cell,qp,i,j);
00175             alphaVal(i,j) = backStressold(cell,qp,i,j);
00176           }
00177         }
00178 
00179         sigmaVal = sigmaN + Intrepid::dotdot(Celastic, depsilon);
00180         ScalarT kappaVal = capParameterold(cell, qp);
00181 
00182         // initialize friction and dilatancy
00183         // (which will be updated only if plasticity occurs)
00184         //friction(cell, qp) = 0.0;
00185         //dilatancy(cell, qp) = 0.0;
00186         //hardeningModulus(cell, qp) = 0.0;
00187 
00188         // define generalized plastic hardening modulus H
00189         //ScalarT H(0.0), Htan(0.0);
00190 
00191         // define plastic strain increment, its two invariants: dev, and vol
00192         ScalarT deqps(0.0), devolps(0.0);
00193 
00194         // define a temporary tensor to store trial stress tensors
00195         sigmaTr = sigmaVal;
00196         // define a temporary tensor to store previous back stress
00197         alphaTr = alphaVal;
00198 
00199         // check yielding
00200         ScalarT f = compute_f(sigmaVal, alphaVal, kappaVal);
00201 
00202         // plastic correction
00203         ScalarT dgamma = 0.0;
00204         if (f > 0.0) {
00205           dfdsigma = compute_dfdsigma(sigmaN, alphaVal, kappaVal);
00206 
00207           dgdsigma = compute_dgdsigma(sigmaN, alphaVal, kappaVal);
00208 
00209           dfdalpha = -dfdsigma;
00210 
00211           ScalarT dfdkappa = compute_dfdkappa(sigmaN, alphaVal, kappaVal);
00212 
00213           ScalarT J2_alpha = 0.5 * Intrepid::dotdot(alphaVal, alphaVal);
00214 
00215           halpha = compute_halpha(dgdsigma, J2_alpha);
00216 
00217           ScalarT I1_dgdsigma = Intrepid::trace(dgdsigma);
00218 
00219           ScalarT dedkappa = compute_dedkappa(kappaVal);
00220 
00221           ScalarT hkappa;
00222           if (dedkappa != 0.0)
00223             hkappa = I1_dgdsigma / dedkappa;
00224           else
00225             hkappa = 0.0;
00226 
00227           ScalarT kai(0.0);
00228           kai = Intrepid::dotdot(dfdsigma, Intrepid::dotdot(Celastic, dgdsigma))
00229               - Intrepid::dotdot(dfdalpha, halpha) - dfdkappa * hkappa;
00230 
00231           //H = -Intrepid::dotdot(dfdalpha, halpha) - dfdkappa * hkappa;
00232 
00233           dfdotCe = Intrepid::dotdot(dfdsigma, Celastic);
00234 
00235           if (kai != 0.0)
00236             dgamma = Intrepid::dotdot(dfdotCe, depsilon) / kai;
00237           else
00238             dgamma = 0.0;
00239 
00240           // update
00241           sigmaVal -= dgamma * Intrepid::dotdot(Celastic, dgdsigma);
00242 
00243           alphaVal += dgamma * halpha;
00244 
00245           // restrictions on kappa, only allow monotonic decreasing (cap hardening)
00246           ScalarT dkappa = dgamma * hkappa;
00247           if (dkappa > 0) {
00248             dkappa = 0;
00249             //H = -Intrepid::dotdot(dfdalpha, halpha);
00250           }
00251 
00252           kappaVal += dkappa;
00253 
00254           // stress correction algorithm to avoid drifting from yield surface
00255           bool condition = false;
00256           int iteration = 0;
00257           int max_iteration = 20;
00258           RealType tolerance = 1.0e-10;
00259           while (condition == false) {
00260             f = compute_f(sigmaVal, alphaVal, kappaVal);
00261 
00262             dfdsigma = compute_dfdsigma(sigmaVal, alphaVal, kappaVal);
00263 
00264             dgdsigma = compute_dgdsigma(sigmaVal, alphaVal, kappaVal);
00265 
00266             dfdalpha = -dfdsigma;
00267 
00268             ScalarT dfdkappa = compute_dfdkappa(sigmaVal, alphaVal, kappaVal);
00269 
00270             J2_alpha = 0.5 * Intrepid::dotdot(alphaVal, alphaVal);
00271 
00272             halpha = compute_halpha(dgdsigma, J2_alpha);
00273 
00274             I1_dgdsigma = Intrepid::trace(dgdsigma);
00275 
00276             dedkappa = compute_dedkappa(kappaVal);
00277 
00278             if (dedkappa != 0)
00279               hkappa = I1_dgdsigma / dedkappa;
00280             else
00281               hkappa = 0;
00282 
00283             //generalized plastic hardening modulus
00284             //H = -Intrepid::dotdot(dfdalpha, halpha) - dfdkappa * hkappa;
00285 
00286             kai = Intrepid::dotdot(dfdsigma,
00287                 Intrepid::dotdot(Celastic, dgdsigma));
00288             kai = kai - Intrepid::dotdot(dfdalpha, halpha) - dfdkappa * hkappa;
00289 
00290             if (std::abs(f) < tolerance)
00291               break;
00292             if (iteration > max_iteration) {
00293               // output for debug
00294               //std::cout << "no stress correction after iteration = "
00295               //<< iteration << " yield function abs(f) = " << abs(f)
00296               //<< std::endl;
00297               break;
00298             }
00299 
00300             ScalarT delta_gamma;
00301             if (kai != 0)
00302               delta_gamma = f / kai;
00303             else
00304               delta_gamma = 0;
00305 
00306             // restrictions on kappa, only allow monotonic decreasing
00307             dkappa = delta_gamma * hkappa;
00308             if (dkappa > 0.0) {
00309               dkappa = 0.0;
00310               //H = -Intrepid::dotdot(dfdalpha, halpha);
00311             }
00312 
00313             // update
00314             sigmaK = sigmaVal
00315                 - delta_gamma * Intrepid::dotdot(Celastic, dgdsigma);
00316             alphaK = alphaVal + delta_gamma * halpha;
00317             ScalarT kappaK = kappaVal + dkappa;
00318 
00319             ScalarT fK = compute_f(sigmaK, alphaK, kappaK);
00320 
00321             if (std::abs(fK) > std::abs(f)) {
00322               // if the corrected stress is further away from yield surface,
00323               // then use normal correction
00324               ScalarT dfdotdf = Intrepid::dotdot(dfdsigma, dfdsigma);
00325               if (dfdotdf != 0)
00326                 delta_gamma = f / dfdotdf;
00327               else
00328                 delta_gamma = 0.0;
00329 
00330               sigmaK = sigmaVal - delta_gamma * dfdsigma;
00331               alphaK = alphaVal;
00332               kappaK = kappaVal;
00333 
00334               //H = 0.0;
00335             }
00336 
00337             sigmaVal = sigmaK;
00338             alphaVal = alphaK;
00339             kappaVal = kappaK;
00340 
00341             iteration++;
00342 
00343           } // end of stress correction
00344 
00345           // compute plastic strain increment
00346           // deps_plastic = compliance ( sigma_tr - sigma_(n+1));
00347           dsigma = sigmaTr - sigmaVal;
00348           deps_plastic = Intrepid::dotdot(compliance, dsigma);
00349 
00350           // compute its two invariants
00351           // devolps (volumetric) and deqps (deviatoric)
00352           devolps = Intrepid::trace(deps_plastic);
00353           dev_plastic = deps_plastic - (1. / 3.) * devolps * I;
00354           // use altenative definition, differ by constants
00355           deqps = std::sqrt(2) * Intrepid::norm(dev_plastic);
00356 
00357           // dilatancy
00358           //if (deqps != 0)
00359           //  dilatancy(cell, qp) = devolps / deqps;
00360           //else
00361           //  dilatancy(cell, qp) = 0.0;
00362 
00363           // previous p and tau
00364           //ScalarT pN(0.0), tauN(0.0);
00365           //xi = sigmaN - alphaTr;
00366           //pN = Intrepid::trace(xi);
00367           //pN = pN / 3.;
00368           //sN = xi - pN * I;
00369           //tauN = sqrt(1. / 2.) * Intrepid::norm(sN);
00370 
00371           // current p, and tau
00372           //ScalarT p(0.0), tau(0.0);
00373           //xi = sigmaVal - alphaVal;
00374           //p = Intrepid::trace(xi);
00375           //p = p / 3.;
00376           //s = xi - p * I;
00377           //tau = sqrt(1. / 2.) * Intrepid::norm(s);
00378 
00379           // difference
00380           //ScalarT dtau = tau - tauN;
00381           //ScalarT dp = p - pN;
00382 
00383           // friction coefficient by finite difference
00384           //if (dp != 0)
00385           //  friction(cell, qp) = dtau / dp;
00386           //else
00387           //  friction(cell, qp) = 0.0;
00388 
00389           // previous gamma(gamma)
00390           //ScalarT evol3 = Intrepid::trace(strainN);
00391           //evol3 = evol3 / 3.;
00392           //eps_dev = strainN - evol3 * I;
00393           //ScalarT gammaN = sqrt(2.) * Intrepid::norm(eps_dev);
00394 
00395           // current gamma(gamma)
00396           //strainCurrent = strainN + depsilon;
00397           //evol3 = Intrepid::trace(strainCurrent);
00398           //evol3 = evol3 / 3.;
00399           //eps_dev = strainCurrent - evol3 * I;
00400           //ScalarT gamma = sqrt(2.) * Intrepid::norm(eps_dev);
00401 
00402           // difference
00403           //ScalarT dGamma = gamma - gammaN;
00404           // tagent hardening modulus
00405           //if (dGamma != 0)
00406           //  Htan = dtau / dGamma;
00407 
00408           //if (std::abs(1. - Htan / mu) > 1.0e-10)
00409           //  hardeningModulus(cell, qp) = Htan / (1. - Htan / mu);
00410           //else
00411           //  hardeningModulus(cell, qp) = 0.0;
00412 
00413         } // end of plastic correction
00414 
00415         // update
00416         for (std::size_t i = 0; i < numDims; ++i) {
00417           for (std::size_t j = 0; j < numDims; ++j) {
00418             stress(cell, qp, i, j) = sigmaVal(i, j);
00419             backStress(cell, qp, i, j) = alphaVal(i, j);
00420           }
00421         }
00422 
00423         capParameter(cell, qp) = kappaVal;
00424         eqps(cell, qp) = eqpsold(cell, qp) + deqps;
00425         volPlasticStrain(cell, qp) = volPlasticStrainold(cell, qp) + devolps;
00426 
00427       } //loop over qps
00428 
00429     } //loop over cell
00430 
00431   } // end of evaluateFields
00432 
00433 //**********************************************************************
00434 // all local functions
00435   template<typename EvalT, typename Traits>
00436   typename CapExplicit<EvalT, Traits>::ScalarT
00437   CapExplicit<EvalT, Traits>::compute_f(
00438       Intrepid::Tensor<ScalarT> & sigma, Intrepid::Tensor<ScalarT> & alpha,
00439       ScalarT & kappa)
00440   {
00441 
00442     xi = sigma - alpha;
00443 
00444     ScalarT I1 = Intrepid::trace(xi), p = I1 / 3;
00445 
00446     s = xi - p * Intrepid::identity<ScalarT>(3);
00447 
00448     ScalarT J2 = 0.5 * Intrepid::dotdot(s, s);
00449 
00450     ScalarT J3 = Intrepid::det(s);
00451 
00452     ScalarT Gamma = 1.0;
00453     if (psi != 0 && J2 != 0)
00454       Gamma = 0.5
00455           * (1 - 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)
00456               + (1 + 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)) / psi);
00457 
00458     ScalarT Ff_I1 = A - C * std::exp(B * I1) - theta * I1;
00459 
00460     ScalarT Ff_kappa = A - C * std::exp(B * kappa) - theta * kappa;
00461 
00462     ScalarT X = kappa - R * Ff_kappa;
00463 
00464     ScalarT Fc = 1.0;
00465 
00466     if ((kappa - I1) > 0 && ((X - kappa) != 0))
00467       Fc = 1.0 - (I1 - kappa) * (I1 - kappa) / (X - kappa) / (X - kappa);
00468 
00469     return Gamma * Gamma * J2 - Fc * (Ff_I1 - N) * (Ff_I1 - N);
00470   }
00471 
00472   template<typename EvalT, typename Traits>
00473   Intrepid::Tensor<typename CapExplicit<EvalT, Traits>::ScalarT>
00474   CapExplicit<EvalT, Traits>::compute_dfdsigma(
00475     Intrepid::Tensor<ScalarT> & sigma,
00476     Intrepid::Tensor<ScalarT> & alpha, ScalarT & kappa)
00477   {
00478 
00479     xi = sigma - alpha;
00480 
00481     ScalarT I1 = Intrepid::trace(xi), p = I1 / 3;
00482 
00483     s = xi - p * Intrepid::identity<ScalarT>(3);
00484 
00485     ScalarT J2 = 0.5 * Intrepid::dotdot(s, s);
00486 
00487     ScalarT J3 = Intrepid::det(s);
00488 
00489     //dI1dsigma = I;
00490     //dJ2dsigma = s;
00491     for (int i = 0; i < 3; i++)
00492       for (int j = 0; j < 3; j++)
00493         dJ3dsigma(i, j) = s(i, j) * s(i, j) - 2 * J2 * I(i, j) / 3;
00494 
00495     ScalarT Ff_I1 = A - C * std::exp(B * I1) - theta * I1;
00496 
00497     ScalarT Ff_kappa = A - C * std::exp(B * kappa) - theta * kappa;
00498 
00499     ScalarT X = kappa - R * Ff_kappa;
00500 
00501     ScalarT Fc = 1.0;
00502 
00503     if ((kappa - I1) > 0)
00504       Fc = 1.0 - (I1 - kappa) * (I1 - kappa) / (X - kappa) / (X - kappa);
00505 
00506     ScalarT Gamma = 1.0;
00507     if (psi != 0 && J2 != 0)
00508       Gamma = 0.5
00509           * (1 - 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)
00510               + (1 + 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)) / psi);
00511 
00512     // derivatives
00513     ScalarT dFfdI1 = -(B * C * std::exp(B * I1) + theta);
00514 
00515     ScalarT dFcdI1 = 0.0;
00516     if ((kappa - I1) > 0 && ((X - kappa) != 0))
00517       dFcdI1 = -2.0 * (I1 - kappa) / (X - kappa) / (X - kappa);
00518 
00519     ScalarT dfdI1 = -(Ff_I1 - N) * (2 * Fc * dFfdI1 + (Ff_I1 - N) * dFcdI1);
00520 
00521     ScalarT dGammadJ2 = 0.0;
00522     if (J2 != 0)
00523       dGammadJ2 = 9.0 * std::sqrt(3.0) * J3 / std::pow(J2, 2.5) / 8.0
00524           * (1.0 - 1.0 / psi);
00525 
00526     ScalarT dfdJ2 = 2.0 * Gamma * dGammadJ2 * J2 + Gamma * Gamma;
00527 
00528     ScalarT dGammadJ3 = 0.0;
00529     if (J2 != 0)
00530       dGammadJ3 = -3.0 * std::sqrt(3.0) / std::pow(J2, 1.5) / 4.0
00531           * (1.0 - 1.0 / psi);
00532 
00533     ScalarT dfdJ3 = 2.0 * Gamma * dGammadJ3 * J2;
00534 
00535     dfdsigma = dfdI1 * I + dfdJ2 * s + dfdJ3 * dJ3dsigma;
00536 
00537     return dfdsigma;
00538   }
00539 
00540   template<typename EvalT, typename Traits>
00541   Intrepid::Tensor<typename CapExplicit<EvalT, Traits>::ScalarT>
00542   CapExplicit<EvalT, Traits>::compute_dgdsigma(
00543     Intrepid::Tensor<ScalarT> & sigma,
00544     Intrepid::Tensor<ScalarT> & alpha, ScalarT & kappa)
00545   {
00546 
00547     xi = sigma - alpha;
00548 
00549     ScalarT I1 = Intrepid::trace(xi), p = I1 / 3;
00550 
00551     s = xi - p * Intrepid::identity<ScalarT>(3);
00552 
00553     ScalarT J2 = 0.5 * Intrepid::dotdot(s, s);
00554 
00555     ScalarT J3 = Intrepid::det(s);
00556 
00557     //dJ2dsigma = s;
00558     //dJ3dsigma(3);
00559     for (int i = 0; i < 3; i++)
00560       for (int j = 0; j < 3; j++)
00561         dJ3dsigma(i, j) = s(i, j) * s(i, j) - 2 * J2 * I(i, j) / 3;
00562 
00563     ScalarT Ff_I1 = A - C * std::exp(L * I1) - phi * I1;
00564 
00565     ScalarT Ff_kappa = A - C * std::exp(L * kappa) - phi * kappa;
00566 
00567     ScalarT X = kappa - Q * Ff_kappa;
00568 
00569     ScalarT Fc = 1.0;
00570 
00571     if ((kappa - I1) > 0)
00572       Fc = 1.0 - (I1 - kappa) * (I1 - kappa) / (X - kappa) / (X - kappa);
00573 
00574     ScalarT Gamma = 1.0;
00575     if (psi != 0 && J2 != 0)
00576       Gamma = 0.5
00577           * (1 - 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)
00578               + (1 + 3.0 * std::sqrt(3.0) * J3 / 2 / std::pow(J2, 1.5)) / psi);
00579 
00580     // derivatives
00581     ScalarT dFfdI1 = -(L * C * std::exp(L * I1) + phi);
00582 
00583     ScalarT dFcdI1 = 0.0;
00584     if ((kappa - I1) > 0 && ((X - kappa) != 0))
00585       dFcdI1 = -2.0 * (I1 - kappa) / (X - kappa) / (X - kappa);
00586 
00587     ScalarT dfdI1 = -(Ff_I1 - N) * (2 * Fc * dFfdI1 + (Ff_I1 - N) * dFcdI1);
00588 
00589     ScalarT dGammadJ2 = 0.0;
00590     if (J2 != 0)
00591       dGammadJ2 = 9.0 * std::sqrt(3.0) * J3 / std::pow(J2, 2.5) / 8.0
00592           * (1.0 - 1.0 / psi);
00593 
00594     ScalarT dfdJ2 = 2.0 * Gamma * dGammadJ2 * J2 + Gamma * Gamma;
00595 
00596     ScalarT dGammadJ3 = 0.0;
00597     if (J2 != 0)
00598       dGammadJ3 = -3.0 * std::sqrt(3.0) / std::pow(J2, 1.5) / 4.0
00599           * (1.0 - 1.0 / psi);
00600 
00601     ScalarT dfdJ3 = 2.0 * Gamma * dGammadJ3 * J2;
00602 
00603     dgdsigma = dfdI1 * I + dfdJ2 * s + dfdJ3 * dJ3dsigma;
00604 
00605     return dgdsigma;
00606   }
00607 
00608   template<typename EvalT, typename Traits>
00609   typename CapExplicit<EvalT, Traits>::ScalarT
00610   CapExplicit<EvalT, Traits>::compute_dfdkappa(
00611     Intrepid::Tensor<ScalarT> & sigma, Intrepid::Tensor<ScalarT> & alpha,
00612       ScalarT & kappa)
00613   {
00614     ScalarT dfdkappa;
00615 
00616     xi = sigma - alpha;
00617 
00618     ScalarT I1 = Intrepid::trace(xi), p = I1 / 3.0;
00619 
00620     s = xi - p * Intrepid::identity<ScalarT>(3);
00621 
00622     ScalarT J2 = 0.5 * Intrepid::dotdot(s, s);
00623 
00624     ScalarT J3 = Intrepid::det(s);
00625 
00626     ScalarT Ff_I1 = A - C * std::exp(B * I1) - theta * I1;
00627 
00628     ScalarT Ff_kappa = A - C * std::exp(B * kappa) - theta * kappa;
00629 
00630     ScalarT X = kappa - R * Ff_kappa;
00631 
00632     ScalarT dFcdkappa = 0.0;
00633 
00634     if ((kappa - I1) > 0 && ((X - kappa) != 0)) {
00635       dFcdkappa = 2 * (I1 - kappa)
00636           * ((X - kappa)
00637               + R * (I1 - kappa) * (theta + B * C * std::exp(B * kappa)))
00638           / (X - kappa) / (X - kappa) / (X - kappa);
00639     }
00640 
00641     dfdkappa = -dFcdkappa * (Ff_I1 - N) * (Ff_I1 - N);
00642 
00643     return dfdkappa;
00644   }
00645   template<typename EvalT, typename Traits>
00646   typename CapExplicit<EvalT, Traits>::ScalarT CapExplicit<EvalT, Traits>::compute_Galpha(
00647       ScalarT & J2_alpha)
00648   {
00649     if (N != 0)
00650       return 1.0 - std::pow(J2_alpha, 0.5) / N;
00651     else
00652       return 0.0;
00653   }
00654 
00655   template<typename EvalT, typename Traits>
00656   Intrepid::Tensor<typename CapExplicit<EvalT, Traits>::ScalarT>
00657   CapExplicit<EvalT, Traits>::compute_halpha(Intrepid::Tensor<ScalarT> & dgdsigma,
00658       ScalarT & J2_alpha)
00659   {
00660 
00661     ScalarT Galpha = compute_Galpha(J2_alpha);
00662 
00663     ScalarT I1 = Intrepid::trace(dgdsigma), p = I1 / 3;
00664 
00665     s = dgdsigma - p * I;
00666 
00667     for (int i = 0; i < 3; i++) {
00668       for (int j = 0; j < 3; j++) {
00669         halpha(i, j) = calpha * Galpha * s(i, j);
00670       }
00671     }
00672 
00673     return halpha;
00674   }
00675 
00676   template<typename EvalT, typename Traits>
00677   typename CapExplicit<EvalT, Traits>::ScalarT
00678   CapExplicit<EvalT, Traits>::compute_dedkappa(
00679       ScalarT & kappa)
00680   {
00681     ScalarT Ff_kappa0 = A - C * std::exp(L * kappa0) - phi * kappa0;
00682 
00683     ScalarT X0 = kappa0 - Q * Ff_kappa0;
00684 
00685     ScalarT Ff_kappa = A - C * std::exp(L * kappa) - phi * kappa;
00686 
00687     ScalarT X = kappa - Q * Ff_kappa;
00688 
00689     ScalarT dedX = (D1 - 2 * D2 * (X - X0))
00690         * std::exp((D1 - D2 * (X - X0)) * (X - X0)) * W;
00691 
00692     ScalarT dXdkappa = 1 + Q * C * L * std::exp(L * kappa) + Q * phi;
00693 
00694     return dedX * dXdkappa;
00695   }
00696 
00697 //**********************************************************************
00698 }// end LCM
CapExplicit_Def.hpp
