chen0407.cpp

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

    dfdsds[1][2] = -1.0/3.0;
    
    dfdsds[2][0] = -1.0/3.0;
    dfdsds[2][1] = -1.0/3.0;
    dfdsds[2][2] =  2.0/3.0;
    
    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
    
  }
  if (state==3){
    dfdsds[0][0] =  1.0/3.0;
    dfdsds[0][1] = -2.0/3.0;
    dfdsds[0][2] = -2.0/3.0;
    
    dfdsds[1][0] = -2.0/3.0;
    dfdsds[1][1] =  1.0/3.0;
    dfdsds[1][2] = -2.0/3.0;
    
    dfdsds[2][0] = -2.0/3.0;
    dfdsds[2][1] = -2.0/3.0;
    dfdsds[2][2] =  1.0/3.0;
    
    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
  }
}

/**
   function evaluates derivatives of yield function
   with respect to internal variable q
   
   JK, 20.4.2005
*/
void chen::deryieldfdq (matrix &sig,vector &q,vector &dfdq)
{
  double invar,j2,zone;
  matrix dev(3,3);
  
  deviator (sig,dev);
  invar=first_invar (sig);
  j2=second_invar(dev);
  
  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;
  }
}

/**

*/
void chen::deryieldfdsigmadq (matrix &sig,vector &q,matrix &dfdsdq)
{
  double invar,j2,zone;
  matrix dev(3,3);
  
  deviator (sig,dev);
  invar=first_invar (sig);
  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;
  }
  
}

void chen::deryieldfdqdq (matrix &dfdqdq)
{
  dfdqdq[0][0]=0.0;
}

void chen::updateq(long ipp,double dgamma,vector &epsp,matrix &sig,vector &q)
{

}


/**
   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];
  //  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->closest_point_proj (ipp,im,ido,gamma,epsn,epsp,q,ni,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);
    //Mm->newton_stress_return_3 (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.num;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;
    }
    case 8:{
      hp=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
}

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