boermat.cpp

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

#include <math.h>
#include "boermat.h"
#include "global.h"
#include "genfile.h"
#include "intpoints.h"
#include "matrix.h"




#define nijac 20
#define limit 1.0e-8

/**
  This constructor inializes attributes to zero values; exponent n
  is initialized with Groen constant.
*/
boermat::boermat (void)
{
  phi=0.0;  c=0.0;  psi=0.0;
  n=1.0/0.229; // Groen constant

}

/**
  This destructor is only for the formal purposes.
*/
boermat::~boermat (void)
{

}


/**
  This function reads material parameters from the opened text file given
  by the parameter in. Then it computes material constants alpha, alpha1, alpha2,
  beta and beta1

  @param in - pointer to the opned text file

*/
void boermat::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf",&phi,&c,&psi);
  sra.read (in);
  
  a=pow((3.0+sin(phi))/(3.0-sin(phi)),n);
  delta=(a-1.0)/(a+1.0);

  alpha=6.0*sin(phi)/(3.0-sin(phi));
  alpha1=alpha*pow(1.0-delta,1.0/n);
  alpha2=6.0*sin(psi)/(3.0-sin(psi));

  beta=6.0*cos(phi)/(3.0-sin(phi));
  beta1=beta*pow(1.0-delta,1.0/n);

}

/**
   This function computes the value of yield functions.

   @param sig - stress tensor

   @retval The function returns value of yield function for the given stress tensor

   25.3.2002
*/
double boermat::yieldfunction (matrix &sig)
{
  double i1s,i3s,j2s,a1,a2,f,sin3m;
  matrix dev(3,3);

  i1s = first_invar (sig)/3.0;

  deviator (sig,dev);

//  i3s = third_invar (sig);
  i3s = third_invar (dev);
  j2s = second_invar (dev);

  if (fabs(j2s)>Mp->zero)
    sin3m = - 3.0*sqrt(3.0)/2.0 * i3s/sqrt(j2s*j2s*j2s);
  else
    sin3m = 0.0;

  a1 = sqrt(3.0*j2s) * pow((1.0-delta*sin3m), 1.0/n);
  a2 = alpha1 * i1s;

  f = a1 + a2 - beta1 * c;

  return f;
}

/**
   This function computes derivatives of as-th yield function
   with respect of vector sigma.

   @param sig - stress tensor
   @param dfds - %matrix where the resulting derivatives are stored


   4.1.2002
*/
void boermat::deryieldfsigma (matrix &sig,matrix &dfds)
{
  double j2s,j3s,sin3m, a1, a2;
  matrix dev(3,3);
  matrix da1ds(3,3), da2ds(3,3), da3ds(3,3), dj3sds(3,3);

  deviator (sig,dev);

  j2s = second_invar (dev);
//  i3s = third_invar (sig);
  j3s = third_invar (dev);
  if (fabs(j2s)>Mp->zero)
    sin3m = - 3.0*sqrt(3.0)/2.0 * j3s/sqrt(j2s*j2s*j2s);
  else
    sin3m = 0.0;
  a1 = sqrt(3.0*j2s);
  a2 = pow(1.0-delta*sin3m, 1.0/n);
  // derivative d_a1/d_sigma
  copym(dev,da1ds);
  cmulm(0.5*sqrt(3.0/j2s), da1ds);
  // derivative d_a2/d_sigma
  //  d I3/d sigma_ij
  dj3sds[0][0]=(sig[1][1]*sig[2][2]-sig[1][2]*sig[2][1])*(2.0/3.0); // d/d sigma_11
  dj3sds[0][1]=sig[1][2]*sig[2][0]-sig[1][0]*sig[2][2]; // d/d sigma_12
  dj3sds[0][2]=sig[1][0]*sig[2][1]-sig[1][1]*sig[2][0]; // d/d sigma_13
  dj3sds[1][0]=sig[0][2]*sig[2][1]-sig[0][1]*sig[2][2]; // d/d sigma_21
  dj3sds[1][1]=(sig[0][0]*sig[2][2]-sig[0][2]*sig[2][0])*(2.0/3.0); // d/d sigma_22
  dj3sds[1][2]=sig[0][1]*sig[2][0]-sig[0][0]*sig[2][1]; // d/d sigma_23
  dj3sds[2][0]=sig[0][1]*sig[1][2]-sig[0][2]*sig[1][1]; // d/d sigma_31
  dj3sds[2][1]=sig[0][2]*sig[1][0]-sig[1][2]*sig[0][0]; // d/d sigma_32
  dj3sds[2][2]=(sig[0][0]*sig[1][1]-sig[0][1]*sig[1][0])*(2.0/3.0); // d/d sigma_33
  cmulm(sqrt(j2s*j2s*j2s), dj3sds);
  copym(dev, da2ds);
  cmulm(3.0/2.0*sqrt(j2s)*j3s, da2ds);
  subm(dj3sds, da2ds, da2ds);
  cmulm(3.0*sqrt(3.0)/2.0/(j2s*j2s*j2s), da2ds);
  cmulm(delta/n*pow(1-delta*sin3m, 1-1.0/n), da2ds);
  // derivative d_a3/d_sigma
  fillm(0.0, da3ds);
  da3ds[0][0] = da3ds[1][1] = da3ds[2][2] = 1.0/3.0*alpha1;

  addm(da1ds, da2ds, dfds);
  addm(dfds, da3ds, dfds);

/*
  c=-1.0*sqrt(3.0*j2s)*(1.0/n)*pow(1.0-delta*sin3m,1.0/n-1)*delta;
  cmulm (c,dfds);

  c=sqrt(3.0)/2.0/sqrt(j2s)*pow(1.0-delta*sin3m,n);
  cmulm (c,dev);

  addm (dev,dfds,dfds);

  dfds[0][0]+=alpha1/3.0;
  dfds[1][1]+=alpha1/3.0;
  dfds[2][2]+=alpha1/3.0;*/

}

/**
   This function computes derivatives of plastic potential function
   with respect of vector sigma.

   @param sig - stress tensor
   @param dgds - %matrix where the resulting derivatives are stored
*/
void boermat::derpotsigma (matrix &sig,matrix &dgds)
{
  double c;
  matrix dev(3,3);

  deviator (sig,dev);
  normedtensor (dev,dgds);
  c=sqrt(6.0)/2.0;
  cmulm (c,dgds);

  dgds[0][0]+=alpha1/3.0;
  dgds[1][1]+=alpha1/3.0;
  dgds[2][2]+=alpha1/3.0;
}

/**
  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 boermat::matstiff (matrix &d, long ipp, long ido)
{
  if (Mp->stmat==0)
  {
    //  initial elastic matrix
    Mm->elmatstiff (d,ipp);
  }
  if (Mp->stmat==1)
  {
    //  tangent stiffness matrix
    //  not completed
    Mm->elmatstiff (d,ipp);
  }
}

/**
  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 im  - index of material type for given ip
  @param ido - index of internal variables for given material in the ipp other array
*/
void boermat::nlstresses (long ipp, long im, long ido)
  //
{
  long i,ni,nc=Mm->ip[ipp].ncompstr;
  double gamma,err;
  vector epsn(nc),epsp(nc),q(0);

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

  //  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<nc;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+nc]=gamma;
  /*  for (i=0; i<nc; i++){
      Mm->ip[ipp].other[ido+i]=epsp[i];
      }
      Mm->ip[ipp].other[ido+nc]=gamma;*/
}

/**
  This function computes stresses at given integration point ipp,
  depending on the reached averaged nonlocal 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 im  - index of material type for given ip
  @param ido - index of internal variables for given material in the ipp other array
*/
void boermat::nonloc_nlstresses (long ipp,long im, long ido)
  //
{
  long i,ni,nc=Mm->ip[ipp].ncompstr;
  double gamma,err;
  vector epsn(nc),epsp(nc),q(0);

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

  //  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<nc;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+nc]=gamma;
}

/**
  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 im  - index of material type for given ip
  @param ido - index of internal variables for given material in the ipp other array

*/
void boermat::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 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 boermat::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

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

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

  return gamma;
}

void boermat::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      phi=val[i];
      break;
    }
    case 1:{
      c=val[i];
      break;
    }
    case 2:{
      psi=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  a=pow((3.0+sin(phi))/(3.0-sin(phi)),n);
  delta=(a-1.0)/(a+1.0);

  alpha=6.0*sin(phi)/(3.0-sin(phi));
  alpha1=alpha*pow(1.0-delta,1.0/n);
  alpha2=6.0*sin(psi)/(3.0-sin(psi));

  beta=6.0*cos(phi)/(3.0-sin(phi));
  beta1=beta*pow(1.0-delta,1.0/n);
  
}