chen.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

/**
   function evaluates derivatives of yield function
   with respect to internal variable q

   @param sig - stress tensor (stored in 3x3 %matrix)
   @param q - internal variables
   @param dfdq - derivatives of yield function with respect to internal variable q (stored in 1x1 %matrix)
   
   JK, 20.4.2005
*/
void chen::deryieldfdq (matrix &sig,vector &q,vector &dfdq)
{
  double invar,j2,zone;
  matrix dev(3,3);
  
  //  computation of stress deviator
  deviator (sig,dev);
  //  first invariant of the stress tensor
  invar=first_invar (sig);
  //  second invariant of the stress deviator
  j2=second_invar(dev);

  fillv (0.0,dfdq);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  if (state==1){
    //  compression-compression
    if(invar<0.0 && zone<0.0){
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      
      if (q[0]<1.0)
	dfdq[0] = (au-ay)*invar/3.0 - (ku-ky);
      else
	dfdq[0]=0.0;
    }
    //  otherwise
    else{
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      
      if (q[0]<1.0)
	dfdq[0] = (au-ay)*invar/3.0 - (ku-ky);
      else
	dfdq[0]=0.0;
    }
  }
  if (state==2){
    ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
    au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
    ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
    ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
    
    if (q[0]<1.0)
      dfdq[0] = (au-ay)*invar/3.0 - (ku-ky);
    else
      dfdq[0]=0.0;
  }
  if (state==3){
    ay = (fyc-fyt)/2.0;
    au = (fc-ft)/2.0;
    ky = fyc*fyt/6.0;
    ku = fc*ft/6.0;
    
    if (q[0]<1.0)
      dfdq[0] = (au-ay)*invar/3.0 - (ku-ky);
    else
      dfdq[0]=0.0;
  }
}

/**
   function evaluates second derivatives of yield function
   with respect to stress components and internal variables
   
   @param sig - stress tensor (stored in 3x3 %matrix)
   @param q - hardening parameters
   @param dfdsdq - derivatives of yield function with respect to stress components (stored in 6x1 %matrix)
   
   JK, 20.4.2005
*/
void chen::deryieldfdsigmadq (matrix &sig,vector &q,matrix &dfdsdq)
{
  double invar,j2,zone;
  matrix dev(3,3);
  
  //  computation of stress deviator
  deviator (sig,dev);
  //  first invariant of the stress tensor
  invar=first_invar (sig);
  //  second invariant of the stress deviator
  j2=second_invar(dev);
  
  fillm (0.0,dfdsdq);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  if (state==1){
    //compression-compression
    if(invar<0.0 && zone<0.0){
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      
      if (q[0]<1.0){
	dfdsdq[0][0]=(au-ay)/3.0;
	dfdsdq[1][0]=(au-ay)/3.0;
	dfdsdq[2][0]=(au-ay)/3.0;
      }

      dfdsdq[3][0]=0.0;
      dfdsdq[4][0]=0.0;
      dfdsdq[5][0]=0.0;
    }
    //otherwise
    else{ 
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      
      if (q[0]<1.0){
	dfdsdq[0][0]=(au-ay)/3.0;
	dfdsdq[1][0]=(au-ay)/3.0;
	dfdsdq[2][0]=(au-ay)/3.0;
      }
      
      dfdsdq[3][0]=0.0;
      dfdsdq[4][0]=0.0;
      dfdsdq[5][0]=0.0;
    }
  }
  if (state==2){
    ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
    au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
    ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
    ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
    
    if (q[0]<1.0){
      dfdsdq[0][0]=(au-ay)/3.0;
      dfdsdq[1][0]=(au-ay)/3.0;
      dfdsdq[2][0]=(au-ay)/3.0;
    }
    
    dfdsdq[3][0]=0.0;
    dfdsdq[4][0]=0.0;
    dfdsdq[5][0]=0.0;
  }
  if (state==3){
    ay = (fyc-fyt)/2.0;
    au = (fc-ft)/2.0;
    ky = fyc*fyt/6.0;
    ku = fc*ft/6.0;
    
    if (q[0]<1.0){
      dfdsdq[0][0]=(au-ay)/3.0;
      dfdsdq[1][0]=(au-ay)/3.0;
      dfdsdq[2][0]=(au-ay)/3.0;
    }
    
    dfdsdq[3][0]=0.0;
    dfdsdq[4][0]=0.0;
    dfdsdq[5][0]=0.0;
  }
  
}



/**
   function computes derivatives of hardening function with respect to stress components
   
   @param sig - stress components (stored in 3x3 %matrix)
   @param q - internal parameters
   @param dhds - derivatives of hardening functions with respect to stresses (stored in 6 x ncomphard %matrix)
   
   JK, 17.2.2007
*/
void chen::dhdsigma (matrix &sigt,vector &q,matrix &dhds)
{
  matrix dfds(3,3),dfdsds(6,6);

  //  first derivatives of yield function with respect to stresses
  deryieldfdsigma (sigt,q,dfds);
  //  second derivatives of yield function with respect to stresses
  deryieldfdsigmadsigma (sigt,dfdsds);
  //  derivatives of hardening function with respect to stresses
  hs.dhdsigma (sigt,dfds,dfdsds,dhds);
  
  if (q[0]>=1.0)
    fillm (0.0,dhds);
}

/**
   function computes derivatives of hardening function with respect to internal parameters

   @param sig - stress components (stored in 3x3 %matrix)
   @param q - internal parameters
   @param dhdq - derivatives of hardening functions with respect to internal parameters (stored in ncomphard x ncomphard %matrix)
   
   JK, 18.2.2007
*/
void chen::dhdqpar (matrix &sigt,vector &q,matrix &dhdq)
{
  matrix dfds(3,3),dfdsdq(6,1);

  //  first derivatives of yield function with respect to stresses
  deryieldfdsigma (sigt,q,dfds);
  //  second derivatives of yield function with respect to stresses
  deryieldfdsigmadq (sigt,q,dfdsdq);
  //  derivatives of hardening function with respect to hardening parameters
  hs.dhdqpar (sigt,dfds,dfdsdq,dhdq);

  if (q[0]>=1.0)
    fillm (0.0,dhdq);
}

/**
   function computes derivatives of hardening function with respect to consistency parameter

   @param dhdg - derivatives of hardening function with respect to consistency parameter stored in ncomphard x 1 %vector
   
   JK, 18.2.2007
*/
void chen::dhdgamma (vector &dhdg)
{
  //  derivatives of hardening function with respect to consistency parameter
  hs.dhdgamma (dhdg);
}

/**
   function computes values of hardening function
   
   @param sig - stress components (stored in 3x3 %matrix)
   @param q - internal parameters stored in %vector
   @param h - values of hardening function stored in %vector
   
   JK, 18.2.2007
*/
void chen::hvalues (matrix &sigt,vector &q,vector &h)
{
  matrix dfds(3,3),dfdsds(6,6);

  //  first derivatives of yield function with respect to stresses
  deryieldfdsigma (sigt,q,dfds);
  //  values of hardening function
  hs.hvalues (sigt,dfds,h);
}



/**
   function computes true stresses from attained strains
   
   @param ipp - number of integration point
   @param im - type of material
   @param ido - index in the array other
   
   21.2.2005
*/
void chen::nlstresses (long ipp, long im, long ido)
{
  long i,ni,n;//nhard;
  double gamma,err;
  
  //  number of strain/stress components
  n = Mm->ip[ipp].ncompstr;
  
  vector epsn(n),epsp(n),q(1);
  
  //  initial values
  for (i=0;i<n;i++){
    //  new total strains
    epsn[i]=Mm->ip[ipp].strain[i];
    //  attained equilibriated plastic starins
    epsp[i]=Mm->ip[ipp].eqother[ido+i];
  }
  //  consistency parameter
  gamma = Mm->ip[ipp].eqother[ido+n];
  
  //if (gamma>0.0){
  //printf ("\n ipp %ld  gamma %le",ipp,gamma);
  //}
  
  //  hardening parameter
  q[0]=Mm->ip[ipp].eqother[ido+n+1];

  //  stress return algorithm
  if (sra.give_tsra () == cp){
    ni=sra.give_ni ();
    err=sra.give_err ();
    
    //Mm->newton_stress_return (ipp,im,ido,gamma,epsn,epsp,q,ni,err);

    //Mm->newton_stress_return_2 (ipp,im,ido,gamma,epsn,epsp,q,ni,err,epslim);
    Mm->newton_stress_return_2 (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
    
    //Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
  }
  else{
    fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }
  
  //fprintf (Out,"\n q %le         gamma %le",q[0],gamma);
  
  //  new data storage
  for (i=0;i<n;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+n]=gamma;
  Mm->ip[ipp].other[ido+n+1]=q[0];
  
}

void chen::nonloc_nlstresses (long ipp, long im, long ido)
{
  long i,ni, n = Mm->ip[ipp].ncompstr;
  double gamma,err;
  vector epsn(n),epsp(n),q(1);

  //  initial values
  for (i=0;i<n;i++){
    epsn[i]=Mm->ip[ipp].strain[i];
    epsp[i]=Mm->ip[ipp].nonloc[i];
  }
  gamma = Mm->ip[ipp].eqother[ido+n];
  q[0] = Mm->ip[ipp].eqother[ido+n+1];

  //  stress return algorithm
  if (sra.give_tsra () == cp){
    ni=sra.give_ni ();
    err=sra.give_err ();
    
    Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
  }
  else{
    fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }
  
  //  new data storage
  for (i=0;i<n;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+n]=gamma;
  Mm->ip[ipp].other[ido+n+1]=q[0];
}

/**
   function updates values of the array eqother
   
   @param ipp - number of integration point
*/
void chen::updateval (long ipp, long im, long ido)
{
  long i,n = Mm->givencompeqother(ipp, im);

  for (i=0;i<n;i++){
    Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
  }
}

/**
   function changes material parameters in stochastic/fuzzy computations
   
   21.2.2005
*/
void chen::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.nba;i++){
    switch (atm.atrib[i]){
    case 0:{
      fyc=val[i];
      break;
    }
    case 1:{
      fyt=val[i];
      break;
    }
    case 2:{
      fybc=val[i];
      break;
    }
    case 3:{
      fc=val[i];
      break;
    }
    case 4:{
      ft=val[i];
      break;
    }
    case 5:{
      fbc=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  //  hardening/softening
  hs.changeparam (atm,val);
  
}

double chen::give_consparam (long ipp,long ido)
{
  long ncompstr;
  double gamma;
  
  ncompstr=Mm->ip[ipp].ncompstr;
  gamma = Mm->ip[ipp].eqother[ido+ncompstr];
  
  return gamma;
}

long chen::give_num_interparam ()
{
  return 1;
}

void chen::give_interparam (long ipp,long ido,vector &q)
{
  long ncompstr=Mm->ip[ipp].ncompstr;
  
  q[0]=Mm->ip[ipp].eqother[ido+ncompstr+1];
}