C++ SymmTensor类代码示例

void
AxisymmetricRZ::computeStrain( const unsigned qp,
                               const SymmTensor & total_strain_old,
                               SymmTensor & total_strain_new,
                               SymmTensor & strain_increment )
{
  strain_increment.xx() = _grad_disp_r[qp](0);
  strain_increment.yy() = _grad_disp_z[qp](1);
  strain_increment.zz() = (_solid_model.q_point(qp)(0) != 0.0 ? _disp_r[qp]/_solid_model.q_point(qp)(0) : 0.0);
  strain_increment.xy() = 0.5*(_grad_disp_r[qp](1) + _grad_disp_z[qp](0));
  if (_large_strain)
  {
    strain_increment.xx() += 0.5*(_grad_disp_r[qp](0)*_grad_disp_r[qp](0) +
                                  _grad_disp_z[qp](0)*_grad_disp_z[qp](0));
    strain_increment.yy() += 0.5*(_grad_disp_r[qp](1)*_grad_disp_r[qp](1) +
                                  _grad_disp_z[qp](1)*_grad_disp_z[qp](1));
    strain_increment.zz() += 0.5*(strain_increment.zz()*strain_increment.zz());
    strain_increment.xy() += 0.5*(_grad_disp_r[qp](0)*_grad_disp_r[qp](1) +
                                  _grad_disp_z[qp](0)*_grad_disp_z[qp](1));
  }

  total_strain_new = strain_increment;

  strain_increment -= total_strain_old;
}

bool
ConstitutiveModel::applyThermalStrain(unsigned qp,
                                      SymmTensor & strain_increment,
                                      SymmTensor & d_strain_dT)
{
  if (_has_temp && _t_step != 0)
  {
    Real inc_thermal_strain;
    Real d_thermal_strain_d_temp;

    Real old_temp;
    if (_t_step == 1 && _has_stress_free_temp)
      old_temp = _stress_free_temp;
    else
      old_temp = _temperature_old[qp];

    Real current_temp = _temperature[qp];

    Real delta_t = current_temp - old_temp;

    Real alpha = _alpha;

    if (_alpha_function)
    {
      Point p;
      Real alpha_current_temp = _alpha_function->value(current_temp,p);
      Real alpha_old_temp = _alpha_function->value(old_temp,p);

      if (_mean_alpha_function)
      {
        Real small(1e-6);

        Real numerator = alpha_current_temp * (current_temp - _ref_temp) - alpha_old_temp * (old_temp - _ref_temp);
        Real denominator = 1.0 + alpha_old_temp * (old_temp - _ref_temp);
        if (denominator < small)
          mooseError("Denominator too small in thermal strain calculation");
        inc_thermal_strain = numerator / denominator;
        d_thermal_strain_d_temp = alpha_current_temp * (current_temp - _ref_temp);
      }
      else
      {
        inc_thermal_strain = delta_t * 0.5 * (alpha_current_temp + alpha_old_temp);
        d_thermal_strain_d_temp = alpha_current_temp;
      }
    }
    else
    {
      inc_thermal_strain = delta_t * alpha;
      d_thermal_strain_d_temp = alpha;
    }

    strain_increment.addDiag(-inc_thermal_strain);
    d_strain_dT.addDiag(-d_thermal_strain_d_temp);

  }

  bool modified = true;
  return modified;
}

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;
}

Real
StressDivergence::computeQpJacobian()
{
  Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
  RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();

  Real jacobian = 0.0;
  // B^T_i * C * B_j
  jacobian += _Jacobian_mult[_qp].stiffness(
      _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j

  if (_volumetric_locking_correction)
  {
    // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
    // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j +  Bvol^T_i * C * B_j + B^T_i * C * Bvol_j

    // Bvol^T_i * C * Bvol_j
    jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
                (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;

    // B^T_i * C * Bvol_j
    jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
                (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;

    // Bvol^T_i * C * B_j
    SymmTensor phi;
    if (_component == 0)
    {
      phi.xx() = _grad_phi[_j][_qp](0);
      phi.xy() = _grad_phi[_j][_qp](1);
      phi.xz() = _grad_phi[_j][_qp](2);
    }
    else if (_component == 1)
    {
      phi.yy() = _grad_phi[_j][_qp](1);
      phi.xy() = _grad_phi[_j][_qp](0);
      phi.yz() = _grad_phi[_j][_qp](2);
    }
    else if (_component == 2)
    {
      phi.zz() = _grad_phi[_j][_qp](2);
      phi.xz() = _grad_phi[_j][_qp](0);
      phi.yz() = _grad_phi[_j][_qp](1);
    }

    SymmTensor tmp(_Jacobian_mult[_qp] * phi);

    jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
                (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
  }

  if (_dt > 0)
    return jacobian * (1 + _alpha + _zeta / _dt);
  else
    return jacobian;
}

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);
  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));
  }

  total_strain_new = strain_increment;

  strain_increment -= total_strain_old;
}

Real
vonMisesStress(const SymmTensor & symm_stress)
{
  Real value = std::pow(symm_stress.xx() - symm_stress.yy(), 2) +
               std::pow(symm_stress.yy() - symm_stress.zz(), 2) +
               std::pow(symm_stress.zz() - symm_stress.xx(), 2) +
               6 * (std::pow(symm_stress.xy(), 2) + std::pow(symm_stress.yz(), 2) +
                    std::pow(symm_stress.zx(), 2));
  return std::sqrt(value / 2.0);
}

Real
StressDivergenceRSpherical::computeQpOffDiagJacobian(unsigned int jvar)
{

  if ( _temp_coupled && jvar == _temp_var )
  {
    SymmTensor test;
    test.xx() = _grad_test[_i][_qp](0);
    test.yy() = _test[_i][_qp] / _q_point[_qp](0);
    test.zz() = test.yy();
    return _d_stress_dT[_qp].doubleContraction(test) * _phi[_j][_qp];
  }

  return 0;
}

bool PoroElasticity::formPermeabilityTensor (SymmTensor& K,
                                             const Vectors&,
                                             const FiniteElement&,
                                             const Vec3& X) const
{
  const PoroMaterial* pmat = dynamic_cast<const PoroMaterial*>(material);
  if (!pmat) return false;

  Vec3 permeability = pmat->getPermeability(X);

  K.zero();
  for (size_t i = 1; i <= K.dim(); i++)
    K(i,i) = permeability(i);

  return true;
}

Real
MaterialTensorAux::principalValue( const SymmTensor & tensor, unsigned int index )
{
  ColumnMajorMatrix eval(3,1);
  ColumnMajorMatrix evec(3,3);
  tensor.columnMajorMatrix().eigen(eval, evec);
  // Eigen computes low to high.  We want high first.
  int i = -index + 2;
  return eval(i);
}

Real
calcPrincipleValues(const SymmTensor & symm_tensor, unsigned int index, RealVectorValue & direction)
{
    ColumnMajorMatrix eval(3,1);
    ColumnMajorMatrix evec(3,3);
    symm_tensor.columnMajorMatrix().eigen(eval, evec);
    // Eigen computes low to high.  We want high first.
    direction(0) = evec(0,index);
    direction(1) = evec(1,index);
    direction(2) = evec(2,index);
    return eval(index);
}

Real
StressDivergenceRZ::computeQpOffDiagJacobian(unsigned int jvar)
{

  if ( _rdisp_coupled && jvar == _rdisp_var )
  {
    return calculateJacobian( _component, 0 );
  }
  else if ( _zdisp_coupled && jvar == _zdisp_var )
  {
    return calculateJacobian( _component, 1 );
  }
  else if ( _temp_coupled && jvar == _temp_var )
  {
    SymmTensor test;
    if (_component == 0)
    {
      test.xx() = _grad_test[_i][_qp](0);
      test.xy() = 0.5*_grad_test[_i][_qp](1);
      test.zz() = _test[_i][_qp] / _q_point[_qp](0);
    }
    else
    {
      test.xy() = 0.5*_grad_test[_i][_qp](0);
      test.yy() = _grad_test[_i][_qp](1);
    }
    return _d_stress_dT[_qp].doubleContraction(test) * _phi[_j][_qp];
  }

  return 0;
}

void
Nonlinear::computeStrainIncrement( const ColumnMajorMatrix & Fhat,
                                   SymmTensor & strain_increment )
{

  //
  // A calculation of the strain at the mid-interval is probably more
  // accurate (second vs. first order).  This would require the
  // incremental deformation gradient at the mid-step, which we
  // currently don't have.  We would then have to calculate the
  // rotation for the whole step.
  //
  //
  // We are looking for:
  //     log( Uhat )
  //  =  log( sqrt( Fhat^T*Fhat ) )
  //  =  log( sqrt( Chat ) )
  // A Taylor series expansion gives:
  //     ( Chat - 0.25 * Chat^T*Chat - 0.75 * I )
  //  =  ( - 0.25 * Chat^T*Chat + Chat - 0.75 * I )
  //  =  ( (0.25*Chat - 0.75*I) * (Chat - I) )
  //  =  ( B * A )
  //    B
  //  = 0.25*Chat - 0.75*I
  //  = 0.25*(Chat - I) - 0.5*I
  //  = 0.25*A - 0.5*I
  //

  const Real Uxx = Fhat(0,0);
  const Real Uxy = Fhat(0,1);
  const Real Uxz = Fhat(0,2);
  const Real Uyx = Fhat(1,0);
  const Real Uyy = Fhat(1,1);
  const Real Uyz = Fhat(1,2);
  const Real Uzx = Fhat(2,0);
  const Real Uzy = Fhat(2,1);
  const Real Uzz = Fhat(2,2);

  const Real Axx = Uxx*Uxx + Uyx*Uyx + Uzx*Uzx - 1.0;
  const Real Axy = Uxx*Uxy + Uyx*Uyy + Uzx*Uzy;
  const Real Axz = Uxx*Uxz + Uyx*Uyz + Uzx*Uzz;
  const Real Ayy = Uxy*Uxy + Uyy*Uyy + Uzy*Uzy - 1.0;
  const Real Ayz = Uxy*Uxz + Uyy*Uyz + Uzy*Uzz;
  const Real Azz = Uxz*Uxz + Uyz*Uyz + Uzz*Uzz - 1.0;

  const Real Bxx = 0.25 * Axx - 0.5;
  const Real Bxy = 0.25 * Axy;
  const Real Bxz = 0.25 * Axz;
  const Real Byy = 0.25 * Ayy - 0.5;
  const Real Byz = 0.25 * Ayz;
  const Real Bzz = 0.25 * Azz - 0.5;

  strain_increment.xx( -(Bxx*Axx + Bxy*Axy + Bxz*Axz) );
  strain_increment.xy( -(Bxx*Axy + Bxy*Ayy + Bxz*Ayz) );
  strain_increment.zx( -(Bxx*Axz + Bxy*Ayz + Bxz*Azz) );
  strain_increment.yy( -(Bxy*Axy + Byy*Ayy + Byz*Ayz) );
  strain_increment.yz( -(Bxy*Axz + Byy*Ayz + Byz*Azz) );
  strain_increment.zz( -(Bxz*Axz + Byz*Ayz + Bzz*Azz) );

}

ColumnMajorMatrix
CrackFrontDefinition::rotateToCrackFrontCoords(const SymmTensor tensor, const unsigned int node_index) const
{
  ColumnMajorMatrix tensor_CMM;
  tensor_CMM(0,0) = tensor.xx();
  tensor_CMM(0,1) = tensor.xy();
  tensor_CMM(0,2) = tensor.xz();
  tensor_CMM(1,0) = tensor.xy();
  tensor_CMM(1,1) = tensor.yy();
  tensor_CMM(1,2) = tensor.yz();
  tensor_CMM(2,0) = tensor.xz();
  tensor_CMM(2,1) = tensor.yz();
  tensor_CMM(2,2) = tensor.zz();

  ColumnMajorMatrix tmp = _rot_matrix[node_index] * tensor_CMM;
  ColumnMajorMatrix rotT = _rot_matrix[node_index].transpose();
  ColumnMajorMatrix rotated_tensor = tmp * rotT;

  return rotated_tensor;
}

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); } );
}

bool Elasticity::kinematics (const Vector& eV,
			     const Vector& N, const Matrix& dNdX, double r,
			     Matrix& B, Tensor&, SymmTensor& eps) const
{
  // Evaluate the strain-displacement matrix, B
  if (axiSymmetry)
  {
    if (!this->formBmatrix(B,N,dNdX,r))
      return false;
  }
  else
  {
    if (!this->formBmatrix(B,dNdX))
      return false;
  }

  if (eV.empty() || eps.dim() == 0)
    return true;

  // Evaluate the strains
  return B.multiply(eV,eps); // eps = B*eV
}