camclay.cpp

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

  strastrestate ssst = Mm->ip[ipp].ssst;
  long ncompstr = Mm->ip[ipp].ncompstr;
  long i;

  vector_tensor (eps, epst, ssst, strain);
  epsv = epst[0][0] + epst[1][1] + epst[2][2];
  // initial value of specific volume under small pressure p_1
  v_kappa1 = Mm->ic[ipp][0];
  // initial reference presure
  p1 = Mm->ic[ipp][1];
  // initial preconsolidation pressure
  pc_0  = Mm->ic[ipp][2];
  if (q[0] == 0.0)
  // original specific volume before any new loading
  {
    Mm->ip[ipp].eqother[ido+ncompstr+2] = v_pc0     = v_kappa1 - kappa*log(pc_0/p1);
    Mm->ip[ipp].eqother[ido+ncompstr+3] = v_lambda1 = v_pc0 + lambda*log(pc_0/p1);
    p_ini = 0.0;
    for (i=0; i< ncompstr; i++)
      p_ini += Mm->ic[ipp][i+3];
    p_ini /= 3.0;
    Mm->ip[ipp].eqother[ido+ncompstr+4] = v_ini    = v_pc0 + kappa*log(pc_0/p_ini);
    q[0]=pc_0;
    return;
  }
  else
  {
    v_pc0     = Mm->ip[ipp].eqother[ido+ncompstr+2];
    v_lambda1 = Mm->ip[ipp].eqother[ido+ncompstr+3];
    v_ini     = Mm->ip[ipp].eqother[ido+ncompstr+4];
  }
  i1s = first_invar (sig)/3.0;
  // vk = v_ini*(1+epsv-epsv0)+kappa*log(i1s/p1);
  vk = v_ini*(1+epsv)+kappa*log(i1s/p1);
  pc_new = p1 * exp((vk-v_lambda1)/(kappa-lambda));
  q[0] = pc_new;
  return;
}


/**
  This function computes material stiffnes matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number
  @param ido - index of internal variables for given material in the ipp other array

*/
void camclay::matstiff (matrix &d, long ipp,long ido)
{
  double zero=1.0e-20;

  if (Mp->stmat==1)
  {
    //  initial elastic matrix
    Mm->elmatstiff (d,ipp);
  }
  if (Mp->stmat==2)
  {
    //  consitent tangent stiffness matrix
    matrix e(d.m, d.n), de(d.m, d.n), ddfdst(6,6), tmp(d.m, d.n), ddfds(d.m,d.n);
    long n = Mm->ip[ipp].ncompstr, i;    
    double gamma = Mm->ip[ipp].other[ido+n];
    for (i = 0; i < d.m; i++)
      e[i][i] = 1.0;
    Mm->elmatstiff (de,ipp);
    dderyieldfsigma (ddfdst);
    //  matrix representation of the fourth order tensor
    tensor4_matrix (ddfds,ddfdst,Mm->ip[ipp].ssst);
    cmulm(gamma, de, d);
    mxm(d, ddfds, tmp);
    addm(e, tmp, d);
    invm(d, tmp,zero);
    mxm(tmp, de, d);
  }
}


/**
  This function computes stresses at given integration point ipp,
  depending on the reached strains.
  The cutting plane algorithm is used. The stress and the other attribute of
  given integration point is actualized.

  @param ipp - integration point number in the mechmat ip array.
  @param ido - index of internal variables for given material in the ipp other array
*/
void camclay::nlstresses (long ipp, long im, long ido)
  //
{
  long i,ni,ncomp=Mm->ip[ipp].ncompstr;
  long j;
  double gamma,err;
  vector epsn(ncomp),epsp(ncomp),q(1);
  vector epsa(ncomp),sig(ncomp),dfds(ncomp),dgds(ncomp);
  matrix d(ncomp,ncomp),sigt(3,3),dfdst(3,3),dgdst(3,3);
  matrix epst(3,3), epspt(3,3);
  matrix dev(3,3);
  double i1s, j2s, epsv, epsvp;

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


  //  stress return algorithm
  switch(sra.give_tsra ())
  {
    case cp :
      ni=sra.give_ni ();
      err=sra.give_err ();
      Mm->cutting_plane (ipp, im, ido, gamma, epsn, epsp, q, ni, err);
      break;
    case gsra :
      ni=sra.give_ni ();
      err=sra.give_err ();
      Mm->newton_stress_return (ipp, im, ido, gamma, epsn, epsp, q, ni, err);
      break;
    default :
      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<ncomp; i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];    
  }
  Mm->ip[ipp].other[ido+ncomp]=gamma;
  Mm->ip[ipp].other[ido+ncomp+1]=q[0];
  Mm->givestress(0, ipp, sig);
  for (j=0; j < ncomp; j++)
    sig[j] += Mm->ic[ipp][3+j];
  vector_tensor (sig,sigt,Mm->ip[ipp].ssst,strain);
  vector_tensor (epsn,epst,Mm->ip[ipp].ssst,strain);
  vector_tensor (epsp,epspt,Mm->ip[ipp].ssst,strain);
  deviator (sigt,dev);
  i1s = first_invar (sigt)/3.0;
  j2s = sqrt(second_invar (dev));
  epsv = first_invar (epst);
  epsvp = first_invar (epspt);
  Mm->ip[ipp].other[ido+ncomp+5] = i1s;  
  Mm->ip[ipp].other[ido+ncomp+6] = j2s;  
  Mm->ip[ipp].other[ido+ncomp+7] = epsv;  
  Mm->ip[ipp].other[ido+ncomp+8] = epsvp;  
}



/*
void camclay::nlstresses (long ipp, long ido)
  //
{
  long i,ni,ncomp=Mm->ip[ipp].ncompstr;
  long j,idem;
  double gamma,err;
  double f,denom,dgamma,nu;
  vector epsn(ncomp),epsp(ncomp),q(1);
  vector epsa(ncomp),sig(ncomp),dfds(ncomp),dgds(ncomp);
  matrix d(ncomp,ncomp),sigt(3,3),dfdst(3,3),dgdst(3,3);
  matrix epst(3,3), epspt(3,3);
  matrix dev(3,3);
  double i1s, j2s, epsv, epsvp;

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


  //  stress return algorithm
  if (sra.give_tsra () == cp){
    ni=sra.give_ni ();
    err=sra.give_err ();
 
    //  initialization
    dgamma=0.0;

    //  elastic stiffness matrix
    Mm->elmatstiff (d,ipp);
    if (Mm->ip[ipp].ssst == planestrain)
    {
       d[0][3] = d[0][1]; d[1][3] = d[1][0];
       d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
    }
    idem = Mm->ip[ipp].gemid();
    if (Mm->ip[ipp].ssst == planestress)
    {
      nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
      epsn[3] = nu / (1.0 - nu) * (epsn[0]+epsn[1]);
    }
    //  main iteration loop
    for (i=0;i<ni;i++){
      //  elastic strain
      subv (epsn,epsp,epsa);
      //  trial stress computation
      mxv (d,epsa,sig);
      for (j=0; j < ncomp; j++)
        sig[j] += Mm->ic[ipp][3+j];
      vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
      // updating hardening variables from previous step
      updateq(ipp, ido, epsn, epsp, sigt, q);
      // checking yield function
      f = yieldfunction (sigt,q);
    
      //if (f>0.0) printf ("\n plasticity at int. point %ld",ipp);
    
      if (f<err)  break;
      if (i==ni-1 && f>err)  
        fprintf (stderr,"\n yield surface was not reached");
      deryieldfsigma (sigt, q, dfdst);
      tensor_vector (dfds,dfdst,Mm->ip[ipp].ssst,stress);

      mxv (d,dfds,epsa);
      scprd (dfds,epsa,denom);
      denom += plasmodscalar(ipp, ido, sigt, q);
      //  new increment of consistency parameter
      dgamma = f/denom;
      //  new internal variables
      gamma+=dgamma;
      // updating plastic strains
      for (j=0;j<ncomp;j++){
        epsp[j]+=dgamma*dfds[j];
      }
    }
    
    //  elastic strain
    subv (epsn,epsp,epsa);
    //  trial stress computation
    mxv (d,epsa,sig);
    //  storing resulting stresses
    Mm->storestress (0,ipp,sig);
  }
  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<ncomp; i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];    
  }
  Mm->ip[ipp].other[ido+ncomp]=gamma;
  Mm->ip[ipp].other[ido+ncomp+1]=q[0];
  for (j=0; j < ncomp; j++)
    sig[j] += Mm->ic[ipp][3+j];
  vector_tensor (sig,sigt,Mm->ip[ipp].ssst,strain);
  vector_tensor (epsn,epst,Mm->ip[ipp].ssst,strain);
  vector_tensor (epsp,epspt,Mm->ip[ipp].ssst,strain);
  deviator (sigt,dev);
  i1s = first_invar (sigt)/3.0;
  j2s = sqrt(second_invar (dev));
  epsv = first_invar (epst);
  epsvp = first_invar (epspt);
  Mm->ip[ipp].other[ido+ncomp+5] = i1s;  
  Mm->ip[ipp].other[ido+ncomp+6] = j2s;  
  Mm->ip[ipp].other[ido+ncomp+7] = epsv;  
  Mm->ip[ipp].other[ido+ncomp+8] = epsvp;  
}*/

/**
  This function updates values in the other array reached in the previous equlibrium state to
  values reached in the new actual equilibrium state.

  @param ipp - integration point number in the mechmat ip array.
  @param ido - index of internal variables for given material in the ipp other array

*/
void camclay::updateval (long ipp,long ido)
{
  long i,n = Mm->ip[ipp].ncompstr;

  for (i=0;i<n;i++){
    Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
  }
  Mm->ip[ipp].eqother[ido+n]   = Mm->ip[ipp].other[ido+n];       // gamma
  Mm->ip[ipp].eqother[ido+n+1] = Mm->ip[ipp].other[ido+n+1];     // pc (hardening) 
//  Mm->ip[ipp].eqother[ido+n+2] = Mm->ip[ipp].other[ido+n+2];     // v_pc0
//  Mm->ip[ipp].eqother[ido+n+3] = Mm->ip[ipp].other[ido+n+3];     // v_lambda1
//  Mm->ip[ipp].eqother[ido+n+4] = Mm->ip[ipp].other[ido+n+4];     // v_ini
  Mm->ip[ipp].eqother[ido+n+5] = Mm->ip[ipp].other[ido+n+5];     // i1s;  
  Mm->ip[ipp].eqother[ido+n+6] = Mm->ip[ipp].other[ido+n+6];     // j2s;  
  Mm->ip[ipp].eqother[ido+n+7] = Mm->ip[ipp].other[ido+n+7];     // epsv;  
  Mm->ip[ipp].eqother[ido+n+8] = Mm->ip[ipp].other[ido+n+8];     // epsvp;  
}



/**
  Function returns irreversible plastic strains.

  @param ipp   - integration point number in the mechmat ip array.
  @param ido   - index of the first internal variable for given material in the ipp other array
  @param epsp  - %vector of irreversible strains
 
  Returns vector of irreversible strains via parameter epsp
*/
void camclay::giveirrstrains (long ipp, long ido, vector &epsp)
{
  long i;
  for (i=0;i<epsp.n;i++)
    epsp[i] = Mm->ip[ipp].eqother[ido+i];
}



/**
  This function extracts consistency parametr gamma for the reached equilibrium state
  from the integration point other array.

  @param ipp - integration point number in the mechmat ip array.
  @param ido - index of internal variables for given material in the ipp other array

  @retval The function returns value of consistency parameter.

*/
double camclay::give_consparam (long ipp,long ido)
{ 
  long ncompstr;
  double gamma;

  ncompstr=Mm->ip[ipp].ncompstr;
  gamma = Mm->ip[ipp].eqother[ido+ncompstr];

  return gamma;
}

void camclay::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      m=val[i];
      break;
    }
    case 1:{
      lambda=val[i];
      break;
    }
    case 2:{
      kappa=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
}