C++ SymmTensor::trace方法代码示例

Real
volumetricStrain(const SymmTensor & symm_strain)
{
  Real value = symm_strain.trace();
  value += symm_strain.xx() * symm_strain.yy() + symm_strain.yy() * symm_strain.zz() +
           symm_strain.zz() * symm_strain.xx() +
           symm_strain.xx() * symm_strain.yy() * symm_strain.zz();
  return value;
}

void
CardiacHolzapfel2009Material::computeQpStressProperties(const SymmTensor &C, const SymmTensor & /*E*/)
{
  const SymmTensor CC(square(C));
  // invariants (We will not need all of them. However, defining them avoids to forget any.)
  const Real I1(C.trace());
  /*const Real I2(0.5*(I1*I1-CC.trace()));*/
  /*const Real I3(_J[_qp]);               */
  const Real I4f(_Ef[_qp]*(C*_Ef[_qp]));
  const Real I4s(_Es[_qp]*(C*_Es[_qp]));
  /*const Real I4n(_En[_qp]*(C*_En[_qp])); */
  /*const Real I5f(_Ef[_qp]*(CC*_Ef[_qp]));*/
  /*const Real I5s(_Es[_qp]*(CC*_Es[_qp]));*/
  /*const Real I5n(_En[_qp]*(CC*_En[_qp]));*/
  const Real I8fs(_Ef[_qp]*(C*_Es[_qp]));
  /*const Real I8fn(_Ef[_qp]*(C*_En[_qp]));*/
  /*const Real I8sn(_Es[_qp]*(C*_En[_qp]));*/

  // the following will be needed in the stress as well as in the energy and stress_derivative
  const Real  i_term(   _p[A1 ]*std::exp(_p[B1 ]*(I1 -3)        ) );
  const Real  f_term( 2*_p[Af ]*std::exp(_p[Bf ]*(I4f-1)*(I4f-1)) );
  const Real  s_term( 2*_p[As ]*std::exp(_p[Bs ]*(I4s-1)*(I4s-1)) );
  const Real fs_term(   _p[Afs]*std::exp(_p[Bfs]* I8fs  * I8fs  ) );

  // elastic energy contribution
  _W[_qp] =  i_term             /(2*_p[B1 ])
         + ( f_term - 2*_p[Af ])/(2*_p[Bf ])
         + ( s_term - 2*_p[As ])/(2*_p[Bs ])
         + (fs_term - 2*_p[Afs])/(2*_p[Bfs]);

  // tensors for constructing stress and stress_derivative
  const SymmTensor EfEf(kron(_Ef[_qp]));
  const SymmTensor EsEs(kron(_Es[_qp]));
  const SymmTensor EfEs(kronSym(_Ef[_qp],_Es[_qp]));

  // stress
  _stress[_qp] = scaledID(i_term) + EfEf*(I4f-1)*f_term + EsEs*(I4s-1)*s_term + EfEs*I8fs*fs_term;

  // stress derivative                          /* fancy lambda function syntax makes things much easier here */
  _stress_derivative[_qp].fill_from_minor_iter( [&](const unsigned int M,
                                                    const unsigned int N,
                                                    const unsigned int P,
                                                    const unsigned int Q) -> Real { return                          _p[B1 ]  *  i_term *  _id(M,N)  * _id(P,Q)
                                                                                           + (1 + (I4f-1)*(I4f-1)*2*_p[Bf ]) *  f_term * EfEf(M,N) * EfEf(P,Q)
                                                                                           + (1 + (I4f-1)*(I4f-1)*2*_p[Bf ]) *  s_term * EsEs(M,N) * EsEs(P,Q)
                                                                                           + (1 +  I8fs          *2*_p[Bfs]) * fs_term * EfEs(M,N) * EfEs(P,Q); } );
}

void
PLC_LSH::computeLSH( const SymmTensor & strain_increment,
                     SymmTensor & plastic_strain_increment,
                     SymmTensor & stress_new )
{

  // compute deviatoric trial stress
  SymmTensor dev_trial_stress(stress_new);
  dev_trial_stress.addDiag( -stress_new.trace()/3 );

  // effective trial stress
  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);

  // determine if yield condition is satisfied
  Real yield_condition = effective_trial_stress - _hardening_variable_old[_qp] - _yield_stress;

  _hardening_variable[_qp] = _hardening_variable_old[_qp];
  _plastic_strain[_qp] = _plastic_strain_old[_qp];

  if (yield_condition > 0)  // then use newton iteration to determine effective plastic strain increment
  {
    unsigned int it = 0;
    Real plastic_residual = 0;
    Real norm_plas_residual = 10;
    Real first_norm_plas_residual = 10;
    Real scalar_plastic_strain_increment = 0;

    while (it < _max_its &&
          norm_plas_residual > _absolute_tolerance &&
          (norm_plas_residual/first_norm_plas_residual) > _relative_tolerance)
    {
      plastic_residual = effective_trial_stress - (3. * _shear_modulus * scalar_plastic_strain_increment) -
        _hardening_variable[_qp] - _yield_stress;
      norm_plas_residual = std::abs(plastic_residual);
      if (it == 0)
      {
        first_norm_plas_residual = norm_plas_residual;
      }

      scalar_plastic_strain_increment += plastic_residual / (3. * _shear_modulus + _hardening_constant);

      _hardening_variable[_qp] = _hardening_variable_old[_qp] + (_hardening_constant * scalar_plastic_strain_increment);

      if (_output_iteration_info == true)
      {
        _console
          << "pls_it="    << it
          << " trl_strs=" << effective_trial_stress
          << " del_p="    << scalar_plastic_strain_increment
          << " harden="   << _hardening_variable[_qp]
          << " rel_res="  << norm_plas_residual/first_norm_plas_residual
          << " rel_tol="  << _relative_tolerance
          << " abs_res="  << norm_plas_residual
          << " abs_tol="  << _absolute_tolerance
          << std::endl;
      }

      ++it;

    }

    if (it == _max_its &&
       norm_plas_residual > _absolute_tolerance &&
       (norm_plas_residual/first_norm_plas_residual) > _relative_tolerance)
    {
      mooseError("Max sub-newton iteration hit during plasticity increment solve!");
    }

    if (effective_trial_stress < 0.01)
    {
      effective_trial_stress = 0.01;
    }
    plastic_strain_increment = dev_trial_stress;
    plastic_strain_increment *= (1.5*scalar_plastic_strain_increment/effective_trial_stress);

    SymmTensor elastic_strain_increment(strain_increment);
    elastic_strain_increment -= plastic_strain_increment;

    // compute stress increment
    stress_new = *elasticityTensor() * elastic_strain_increment;

    // update stress and plastic strain
    stress_new += _stress_old;
    _plastic_strain[_qp] += plastic_strain_increment;

  } // end of if statement

}

void
ReturnMappingModel::computeStress(const Elem & /*current_elem*/,
                                  const SymmElasticityTensor & elasticityTensor,
                                  const SymmTensor & stress_old,
                                  SymmTensor & strain_increment,
                                  SymmTensor & stress_new,
                                  SymmTensor & inelastic_strain_increment)
{
  // compute deviatoric trial stress
  SymmTensor dev_trial_stress(stress_new);
  dev_trial_stress.addDiag(-dev_trial_stress.trace() / 3.0);

  // compute effective trial stress
  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);

  // compute effective strain increment
  SymmTensor dev_strain_increment(strain_increment);
  dev_strain_increment.addDiag(-strain_increment.trace() / 3.0);
  _effective_strain_increment = dev_strain_increment.doubleContraction(dev_strain_increment);
  _effective_strain_increment = std::sqrt(2.0 / 3.0 * _effective_strain_increment);

  const SymmIsotropicElasticityTensor * iso_e_t =
      dynamic_cast<const SymmIsotropicElasticityTensor *>(&elasticityTensor);
  if (!iso_e_t)
    mooseError("Models derived from ReturnMappingModel require a SymmIsotropicElasticityTensor");
  _three_shear_modulus = 3.0 * iso_e_t->shearModulus();

  computeStressInitialize(effective_trial_stress, elasticityTensor);

  Real scalar;
  returnMappingSolve(effective_trial_stress, scalar, _console);

  // compute inelastic and elastic strain increments
  if (_legacy_return_mapping)
  {
    if (effective_trial_stress < 0.01)
      effective_trial_stress = 0.01;

    inelastic_strain_increment = dev_trial_stress;
    inelastic_strain_increment *= (1.5 * scalar / effective_trial_stress);
  }
  else
  {
    if (scalar != 0.0)
      inelastic_strain_increment = dev_trial_stress * (1.5 * scalar / effective_trial_stress);
    else
      inelastic_strain_increment = 0.0;
  }

  strain_increment -= inelastic_strain_increment;
  _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + scalar;

  // compute stress increment
  stress_new = elasticityTensor * strain_increment;

  // update stress
  stress_new += stress_old;

  computeStressFinalize(inelastic_strain_increment);
}

bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc,
                                          const SymmTensor& epsil, double* Phi,
                                          SymmTensor* sigma, Matrix* dSdE,
                                          bool postProc, bool printElm) const
{
  PROFILE3("FractureEl::evalStress");

  unsigned short int a = 0, b = 0;

  // Define a Lambda-function to set up the isotropic constitutive matrix
  auto&& setIsotropic = [this,a,b](Matrix& C, double lambda, double mu) mutable
  {
    for (a = 1; a <= C.rows(); a++)
      if (a > nsd)
        C(a,a) = mu;
      else
      {
        C(a,a) = 2.0*mu;
        for (b = 1; b <= nsd; b++)
          C(a,b) += lambda;
      }
  };

  // Define some material constants
  double trEps = epsil.trace();
  double C0 = trEps >= -epsZ ? Gc*lambda : lambda;
  double Cp = Gc*mu;

  if (trEps >= -epsZ && trEps <= epsZ)
  {
    // No strains, stress free configuration
    Phi[0] = 0.0;
    if (postProc)
      Phi[1] = Phi[2] = Phi[3] = 0.0;
    if (sigma)
      *sigma = 0.0;
    if (dSdE)
      setIsotropic(*dSdE,C0,Cp);
    return true;
  }

  // Calculate principal strains and the associated directions
  Vec3 eps;
  std::vector<SymmTensor> M(nsd,SymmTensor(nsd));
  {
    PROFILE4("Tensor::principal");
    if (!epsil.principal(eps,M.data()))
      return false;
  }

  // Split the strain tensor into positive and negative parts
  SymmTensor ePos(nsd), eNeg(nsd);
  for (a = 0; a < nsd; a++)
    if (eps[a] > 0.0)
      ePos += eps[a]*M[a];
    else if (eps[a] < 0.0)
      eNeg += eps[a]*M[a];

  if (sigma)
  {
    // Evaluate the stress tensor
    *sigma = C0*trEps;
    *sigma += 2.0*mu*(Gc*ePos + eNeg);
  }

  // Evaluate the tensile energy
  Phi[0] = mu*(ePos*ePos).trace();
  if (trEps > 0.0) Phi[0] += 0.5*lambda*trEps*trEps;
  if (postProc)
  {
    // Evaluate the compressive energy
    Phi[1] = mu*(eNeg*eNeg).trace();
    if (trEps < 0.0) Phi[1] += 0.5*lambda*trEps*trEps;
    // Evaluate the total strain energy
    Phi[2] = Gc*Phi[0] + Phi[1];
    // Evaluate the bulk energy
    Phi[3] = Gc*(Phi[0] + Phi[1]);
  }
  else if (sigmaC > 0.0) // Evaluate the crack driving function
    Phi[0] = this->MieheCrit56(eps,lambda,mu);

#if INT_DEBUG > 4
  std::cout <<"eps_p = "<< eps <<"\n";
  for (a = 0; a < nsd; a++)
    std::cout <<"M("<< 1+a <<") =\n"<< M[a];
  std::cout <<"ePos =\n"<< ePos <<"eNeg =\n"<< eNeg;
  if (sigma) std::cout <<"sigma =\n"<< *sigma;
  std::cout <<"Phi = "<< Phi[0];
  if (postProc) std::cout <<" "<< Phi[1] <<" "<< Phi[2] <<" "<< Phi[3];
  std::cout << std::endl;
#else
  if (printElm)
  {
    std::cout <<"g(c) = "<< Gc
              <<"\nepsilon =\n"<< epsil <<"eps_p = "<< eps
              <<"\nePos =\n"<< ePos <<"eNeg =\n"<< eNeg;
    if (sigma) std::cout <<"sigma =\n"<< *sigma;
    std::cout <<"Phi = "<< Phi[0];
    if (postProc) std::cout <<" "<< Phi[1] <<" "<< Phi[2] <<" "<< Phi[3];
    std::cout << std::endl;
//.........这里部分代码省略.........

void
Linear::computeStrain( const unsigned qp,
                       const SymmTensor & total_strain_old,
                       SymmTensor & total_strain_new,
                       SymmTensor & strain_increment )
{
  strain_increment.xx( _grad_disp_x[qp](0) );
  strain_increment.yy( _grad_disp_y[qp](1) );
  strain_increment.zz( _grad_disp_z[qp](2) );
  strain_increment.xy( 0.5*(_grad_disp_x[qp](1)+_grad_disp_y[qp](0)) );
  strain_increment.yz( 0.5*(_grad_disp_y[qp](2)+_grad_disp_z[qp](1)) );
  strain_increment.zx( 0.5*(_grad_disp_z[qp](0)+_grad_disp_x[qp](2)) );
  if (_large_strain)
  {
    strain_increment.xx() += 0.5*(_grad_disp_x[qp](0)*_grad_disp_x[qp](0) +
                                  _grad_disp_y[qp](0)*_grad_disp_y[qp](0) +
                                  _grad_disp_z[qp](0)*_grad_disp_z[qp](0));
    strain_increment.yy() += 0.5*(_grad_disp_x[qp](1)*_grad_disp_x[qp](1) +
                                  _grad_disp_y[qp](1)*_grad_disp_y[qp](1) +
                                  _grad_disp_z[qp](1)*_grad_disp_z[qp](1));
    strain_increment.zz() += 0.5*(_grad_disp_x[qp](2)*_grad_disp_x[qp](2) +
                                  _grad_disp_y[qp](2)*_grad_disp_y[qp](2) +
                                  _grad_disp_z[qp](2)*_grad_disp_z[qp](2));
    strain_increment.xy() += 0.5*(_grad_disp_x[qp](0)*_grad_disp_x[qp](1) +
                                  _grad_disp_y[qp](0)*_grad_disp_y[qp](1) +
                                  _grad_disp_z[qp](0)*_grad_disp_z[qp](1));
    strain_increment.yz() += 0.5*(_grad_disp_x[qp](1)*_grad_disp_x[qp](2) +
                                  _grad_disp_y[qp](1)*_grad_disp_y[qp](2) +
                                  _grad_disp_z[qp](1)*_grad_disp_z[qp](2));
    strain_increment.zx() += 0.5*(_grad_disp_x[qp](2)*_grad_disp_x[qp](0) +
                                  _grad_disp_y[qp](2)*_grad_disp_y[qp](0) +
                                  _grad_disp_z[qp](2)*_grad_disp_z[qp](0));
  }

  if (_volumetric_locking_correction)
  {
    // volumetric locking correction - averaging the volumertic strain over the element
    Real volumetric_strain = 0.0;
    Real volume = 0.0;
    for (unsigned int qp_loop = 0; qp_loop < _solid_model.qrule()->n_points(); ++qp_loop)
    {
      volumetric_strain += (_grad_disp_x[qp_loop](0) + _grad_disp_y[qp_loop](1) + _grad_disp_z[qp_loop](2)) / 3.0 * _solid_model.JxW(qp_loop);

      volume += _solid_model.JxW(qp_loop);

      if (_large_strain)
      {
        volumetric_strain += 0.5 * (_grad_disp_x[qp](0) * _grad_disp_x[qp](0) +
                                    _grad_disp_y[qp](0) * _grad_disp_y[qp](0) +
                                    _grad_disp_z[qp](0) * _grad_disp_z[qp](0)) / 3.0 * _solid_model.JxW(qp_loop);
        volumetric_strain += 0.5 * (_grad_disp_x[qp](1) * _grad_disp_x[qp](1) +
                                    _grad_disp_y[qp](1) * _grad_disp_y[qp](1) +
                                    _grad_disp_z[qp](1) * _grad_disp_z[qp](1)) / 3.0 * _solid_model.JxW(qp_loop);
        volumetric_strain += 0.5 * (_grad_disp_x[qp](2) * _grad_disp_x[qp](2) +
                                    _grad_disp_y[qp](2) * _grad_disp_y[qp](2) +
                                    _grad_disp_z[qp](2) * _grad_disp_z[qp](2)) / 3.0 * _solid_model.JxW(qp_loop);
      }
    }

    volumetric_strain /= volume; // average volumetric strain

    // strain increment at _qp
    Real trace = strain_increment.trace();
    strain_increment.xx() += volumetric_strain - trace / 3.0;
    strain_increment.yy() += volumetric_strain - trace / 3.0;
    strain_increment.zz() += volumetric_strain - trace / 3.0;
  }

  total_strain_new = strain_increment;
  strain_increment -= total_strain_old;
}

Real
firstInvariant(const SymmTensor & symm_tensor)
{
    return symm_tensor.trace();
}

void
PlaneStrain::computeStrain( const unsigned qp,
                               const SymmTensor & total_strain_old,
                               SymmTensor & total_strain_new,
                               SymmTensor & strain_increment )
{
  strain_increment.xx() = _grad_disp_x[qp](0);
  strain_increment.yy() = _grad_disp_y[qp](1);

  if (_have_strain_zz)
    strain_increment.zz() = _strain_zz[qp];
  else if (_have_scalar_strain_zz && _scalar_strain_zz.size()>0)
    strain_increment.zz() = _scalar_strain_zz[0];
  else
    strain_increment.zz() = 0;

  strain_increment.xy() = 0.5*(_grad_disp_x[qp](1) + _grad_disp_y[qp](0));
  strain_increment.yz() = 0;
  strain_increment.zx() = 0;
  if (_large_strain)
  {
    strain_increment.xx() += 0.5*(_grad_disp_x[qp](0)*_grad_disp_x[qp](0) +
                                  _grad_disp_y[qp](0)*_grad_disp_y[qp](0));
    strain_increment.yy() += 0.5*(_grad_disp_x[qp](1)*_grad_disp_x[qp](1) +
                                  _grad_disp_y[qp](1)*_grad_disp_y[qp](1));
    strain_increment.xy() += 0.5*(_grad_disp_x[qp](0)*_grad_disp_x[qp](1) +
                                  _grad_disp_y[qp](0)*_grad_disp_y[qp](1));
  }

  if (_volumetric_locking_correction)
  {
    // volumetric locking correction
    Real volumetric_strain = 0.0;
    Real volume = 0.0;
    Real dim = 3.0;
    for (unsigned int qp_loop = 0; qp_loop < _solid_model.qrule()->n_points(); ++qp_loop)
    {
      if (_have_strain_zz)
        volumetric_strain += (_grad_disp_x[qp_loop](0) + _grad_disp_y[qp_loop](1) + _strain_zz[qp_loop]) / dim * _solid_model.JxW(qp_loop);
      else if (_have_scalar_strain_zz && _scalar_strain_zz.size() > 0)
        volumetric_strain += (_grad_disp_x[qp_loop](0) + _grad_disp_y[qp_loop](1) + _scalar_strain_zz[0]) / dim * _solid_model.JxW(qp_loop);
      else
        volumetric_strain += (_grad_disp_x[qp_loop](0) + _grad_disp_y[qp_loop](1)) / dim * _solid_model.JxW(qp_loop);

      volume += _solid_model.JxW(qp_loop);

      if (_large_strain)
      {
        volumetric_strain += 0.5 * (_grad_disp_x[qp_loop](0) * _grad_disp_x[qp_loop](0) +
                                    _grad_disp_y[qp_loop](0) * _grad_disp_y[qp_loop](0)) / dim * _solid_model.JxW(qp_loop);
        volumetric_strain += 0.5 * (_grad_disp_x[qp_loop](1) * _grad_disp_x[qp_loop](1) +
                                    _grad_disp_y[qp_loop](1) * _grad_disp_y[qp_loop](1)) / dim * _solid_model.JxW(qp_loop);
      }
    }

    volumetric_strain /= volume; // average volumetric strain

    // strain increment at _qp
    Real trace = strain_increment.trace();
    strain_increment.xx() += volumetric_strain - trace / dim;
    strain_increment.yy() += volumetric_strain - trace / dim;
    strain_increment.zz() += volumetric_strain - trace / dim;
  }

  total_strain_new = strain_increment;

  strain_increment -= total_strain_old;
}

void
ReturnMappingModel::computeStress(const Elem & /*current_elem*/, unsigned qp,
                                  const SymmElasticityTensor & elasticityTensor,
                                  const SymmTensor & stress_old,
                                  SymmTensor & strain_increment,
                                  SymmTensor & stress_new,
                                  SymmTensor & inelastic_strain_increment)
{
  // compute deviatoric trial stress
  SymmTensor dev_trial_stress(stress_new);
  dev_trial_stress.addDiag(-dev_trial_stress.trace()/3.0);

  // compute effective trial stress
  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);

  // compute effective strain increment
  SymmTensor dev_strain_increment(strain_increment);
  dev_strain_increment.addDiag(-strain_increment.trace()/3.0);
  _effective_strain_increment = dev_strain_increment.doubleContraction(dev_strain_increment);
  _effective_strain_increment = std::sqrt(2.0/3.0 * _effective_strain_increment);

  computeStressInitialize(qp, effective_trial_stress, elasticityTensor);

  // Use Newton sub-iteration to determine inelastic strain increment

  Real scalar = 0;
  unsigned int it = 0;
  Real residual = 10;
  Real norm_residual = 10;
  Real first_norm_residual = 10;

  std::string iter_output;

  while (it < _max_its &&
        norm_residual > _absolute_tolerance &&
        (norm_residual/first_norm_residual) > _relative_tolerance)
  {
    iterationInitialize(qp, scalar);

    residual = computeResidual(qp, effective_trial_stress, scalar);
    norm_residual = std::abs(residual);
    if (it == 0)
    {
      first_norm_residual = norm_residual;
      if (first_norm_residual == 0)
      {
        first_norm_residual = 1;
      }
    }

    scalar -= residual / computeDerivative(qp, effective_trial_stress, scalar);

    if (_output_iteration_info == true ||
        _output_iteration_info_on_error == true)
    {
        iter_output = "In the element " + Moose::stringify(_current_elem->id()) +
                         + " and the qp point " + Moose::stringify(qp) + ": \n" +
                         + " iteration = " + Moose::stringify(it ) + "\n" +
                         + " effective trial stress = " + Moose::stringify(effective_trial_stress) + "\n" +
                         + " scalar effective inelastic strain = " + Moose::stringify(scalar) +"\n" +
                         + " relative residual = " + Moose::stringify(norm_residual/first_norm_residual) + "\n" +
                         + " relative tolerance = " + Moose::stringify(_relative_tolerance) + "\n" +
                         + " absolute residual = " + Moose::stringify(norm_residual) + "\n" +
                         + " absolute tolerance = " + Moose::stringify(_absolute_tolerance) + "\n";
      }
    iterationFinalize(qp, scalar);
    ++it;
  }

  if (_output_iteration_info)
    _console << iter_output;

  if (it == _max_its &&
     norm_residual > _absolute_tolerance &&
     (norm_residual/first_norm_residual) > _relative_tolerance)
  {
    if (_output_iteration_info_on_error)
    {
      Moose::err << iter_output;
    }
    mooseError("Exceeded maximum iterations in ReturnMappingModel solve for material: " << _name << ".  Rerun with  'output_iteration_info_on_error = true' for more information.");
  }

  // compute inelastic and elastic strain increments (avoid potential divide by zero - how should this be done)?
  if (effective_trial_stress < 0.01)
  {
    effective_trial_stress = 0.01;
  }

  inelastic_strain_increment = dev_trial_stress;
  inelastic_strain_increment *= (1.5*scalar/effective_trial_stress);

  strain_increment -= inelastic_strain_increment;

  // compute stress increment
  stress_new = elasticityTensor * strain_increment;

  // update stress
  stress_new += stress_old;
//.........这里部分代码省略.........

void
ReturnMappingModel::computeStress( const Elem & /*current_elem*/, unsigned qp,
                                   const SymmElasticityTensor & elasticityTensor,
                                   const SymmTensor & stress_old,
                                   SymmTensor & strain_increment,
                                   SymmTensor & stress_new,
                                   SymmTensor & inelastic_strain_increment )
{

  // compute deviatoric trial stress
  SymmTensor dev_trial_stress(stress_new);
  dev_trial_stress.addDiag( -dev_trial_stress.trace()/3.0 );

  // compute effective trial stress
  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);

  // compute effective strain increment
  SymmTensor dev_strain_increment(strain_increment);
  dev_strain_increment.addDiag( -strain_increment.trace()/3.0);
  _effective_strain_increment = dev_strain_increment.doubleContraction(dev_strain_increment);
  _effective_strain_increment = std::sqrt(2.0/3.0 * _effective_strain_increment);

  computeStressInitialize(qp, effective_trial_stress, elasticityTensor);

  // Use Newton sub-iteration to determine inelastic strain increment

  Real scalar = 0;
  unsigned int it = 0;
  Real residual = 10;
  Real norm_residual = 10;
  Real first_norm_residual = 10;

  std::stringstream iter_output;

  while (it < _max_its &&
        norm_residual > _absolute_tolerance &&
        (norm_residual/first_norm_residual) > _relative_tolerance)
  {
    iterationInitialize( qp, scalar );

    residual = computeResidual(qp, effective_trial_stress, scalar);
    norm_residual = std::abs(residual);
    if (it == 0)
    {
      first_norm_residual = norm_residual;
      if (first_norm_residual == 0)
      {
        first_norm_residual = 1;
      }
    }

    scalar -= residual / computeDerivative(qp, effective_trial_stress, scalar);

    if (_output_iteration_info == true ||
        _output_iteration_info_on_error == true)
    {
      iter_output
        << " it="       << it
        << " trl_strs=" << effective_trial_stress
        << " scalar="   << scalar
        << " rel_res="  << norm_residual/first_norm_residual
        << " rel_tol="  << _relative_tolerance
        << " abs_res="  << norm_residual
        << " abs_tol="  << _absolute_tolerance
        << std::endl;
    }

    iterationFinalize( qp, scalar );

    ++it;
  }

  if (_output_iteration_info)
    _console << iter_output.str();


  if (it == _max_its &&
     norm_residual > _absolute_tolerance &&
     (norm_residual/first_norm_residual) > _relative_tolerance)
  {
    if (_output_iteration_info_on_error)
    {
      Moose::err << iter_output.str();
    }
    mooseError("Max sub-newton iteration hit during nonlinear constitutive model solve! (" << _name << ")");
  }

  // compute inelastic and elastic strain increments (avoid potential divide by zero - how should this be done)?
  if (effective_trial_stress < 0.01)
  {
    effective_trial_stress = 0.01;
  }

  inelastic_strain_increment = dev_trial_stress;
  inelastic_strain_increment *= (1.5*scalar/effective_trial_stress);

  strain_increment -= inelastic_strain_increment;

  // compute stress increment
//.........这里部分代码省略.........

Real
MaterialTensorAux::getTensorQuantity(const SymmTensor & tensor,
                                     const MTA_ENUM quantity,
                                     const MooseEnum & quantity_moose_enum,
                                     const int index,
                                     const Point * curr_point,
                                     const Point * p1,
                                     const Point * p2)
{
  Real value(0);
  if (quantity == MTA_COMPONENT)
  {
    value = tensor.component(index);
  }
  else if ( quantity == MTA_VONMISES )
  {
    value = std::sqrt(0.5*(
                           std::pow(tensor.xx() - tensor.yy(), 2) +
                           std::pow(tensor.yy() - tensor.zz(), 2) +
                           std::pow(tensor.zz() - tensor.xx(), 2) + 6 * (
                           std::pow(tensor.xy(), 2) +
                           std::pow(tensor.yz(), 2) +
                           std::pow(tensor.zx(), 2))));
  }
  else if ( quantity == MTA_PLASTICSTRAINMAG )
  {
    value = std::sqrt(2.0/3.0 * tensor.doubleContraction(tensor));
  }
  else if ( quantity == MTA_HYDROSTATIC )
  {
    value = tensor.trace()/3.0;
  }
  else if ( quantity == MTA_HOOP )
  {
    // This is the location of the stress.  A vector from the cylindrical axis to this point will define the x' axis.
    Point p0( *curr_point );

    // The vector p1 + t*(p2-p1) defines the cylindrical axis.  The point along this
    // axis closest to p0 is found by the following for t:
    const Point p2p1( *p2 - *p1 );
    const Point p2p0( *p2 - p0 );
    const Point p1p0( *p1 - p0 );
    const Real t( -(p1p0*p2p1)/p2p1.size_sq() );
    // The nearest point on the cylindrical axis to p0 is p.
    const Point p( *p1 + t * p2p1 );
    Point xp( p0 - p );
    xp /= xp.size();
    Point yp( p2p1/p2p1.size() );
    Point zp( xp.cross( yp ));
    //
    // The following works but does more than we need
    //
//    // Rotation matrix R
//    ColumnMajorMatrix R(3,3);
//    // Fill with direction cosines
//    R(0,0) = xp(0);
//    R(1,0) = xp(1);
//    R(2,0) = xp(2);
//    R(0,1) = yp(0);
//    R(1,1) = yp(1);
//    R(2,1) = yp(2);
//    R(0,2) = zp(0);
//    R(1,2) = zp(1);
//    R(2,2) = zp(2);
//    // Rotate
//    ColumnMajorMatrix tensor( _tensor[_qp].columnMajorMatrix() );
//    ColumnMajorMatrix tensorp( R.transpose() * ( tensor * R ));
//    // Hoop stress is the zz stress
//    value = tensorp(2,2);
    //
    // Instead, tensorp(2,2) = R^T(2,:)*tensor*R(:,2)
    //
    const Real zp0( zp(0) );
    const Real zp1( zp(1) );
    const Real zp2( zp(2) );
    value = (zp0*tensor(0,0)+zp1*tensor(1,0)+zp2*tensor(2,0))*zp0 +
            (zp0*tensor(0,1)+zp1*tensor(1,1)+zp2*tensor(2,1))*zp1 +
            (zp0*tensor(0,2)+zp1*tensor(1,2)+zp2*tensor(2,2))*zp2;
  }
  else if ( quantity == MTA_RADIAL )
  {
    // This is the location of the stress.  A vector from the cylindrical axis to this point will define the x' axis
    // which is the radial direction in which we want the stress.
    Point p0( *curr_point );

    // The vector p1 + t*(p2-p1) defines the cylindrical axis.  The point along this
    // axis closest to p0 is found by the following for t:
    const Point p2p1( *p2 - *p1 );
    const Point p2p0( *p2 - p0 );
    const Point p1p0( *p1 - p0 );
    const Real t( -(p1p0*p2p1)/p2p1.size_sq() );
    // The nearest point on the cylindrical axis to p0 is p.
    const Point p( *p1 + t * p2p1 );
    Point xp( p0 - p );
    xp /= xp.size();
    const Real xp0( xp(0) );
    const Real xp1( xp(1) );
    const Real xp2( xp(2) );
    value = (xp0*tensor(0,0)+xp1*tensor(1,0)+xp2*tensor(2,0))*xp0 +
            (xp0*tensor(0,1)+xp1*tensor(1,1)+xp2*tensor(2,1))*xp1 +
//.........这里部分代码省略.........