mohrcparab.cpp

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

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "mohrcparab.h"
#include "global.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"
#include "intpoints.h"


#define nijac 20
#define limit 1.0e-8

/**
  This constructor inializes attributes to zero values.
*/
mohrcoulombpar::mohrcoulombpar (void)
{
  phi=0.0;  c=0.0;  psi=0.0;
  alphaphi = 0.0; betaphi = 0.0; sigcphi = 0.0;
  alphapsi = 0.0; betapsi = 0.0; sigcpsi = 0.0;
}

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

}


/**
  This function reads material parameters from the opened text file given
  by the parameter in. Internal parameters alpha, beta and sigc
  for phi and psi are evaluated.

  @param in - pointer to the opned text file

*/
void mohrcoulombpar::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf",&phi,&c,&psi);
  sra.read (in);

  alphaphi = sin(phi)/sqrt(9+3*sin(phi)*sin(phi));
  betaphi  = alphaphi / sqrt(3*(3*c*c - alphaphi*alphaphi));
  sigcphi  = sqrt(3*(c*c - alphaphi*alphaphi/3.0));

  alphapsi = sin(psi)/sqrt(9+3*sin(psi)*sin(psi));
  betapsi  = alphapsi / sqrt(3*(3*c*c - alphapsi*alphapsi));
  sigcpsi  = sqrt(3*(c*c - alphapsi*alphapsi/3.0));
}

/**
   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 mohrcoulombpar::yieldfunction (matrix &sig)
{
  double f, tmp;
  double j2, j1;
  matrix dev(3,3);

  deviator(sig, dev);
  j1 = first_invar(dev);
  j2 = second_invar(dev);
  tmp = 3 * j2 + sqrt(3.0) * betaphi * sigcphi * j1;
  f = sqrt(tmp) - sigcphi;

  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 mohrcoulombpar::deryieldfsigma (matrix &sig,matrix &dfds)
{
  double j2, j1, k1, k2;
  matrix dev(3,3);

  deviator(sig, dev);
  j1 = first_invar(dev);
  j2 = second_invar(dev);
  k1 = 1.0 / (sqrt(j2 + sqrt(3.0)*betaphi*sigcphi*j1));
  k2 = sqrt(3.0) * betaphi * sigcphi;
  dfds[0][0] = k1 * (3.0 * dev[0][0] + k2);
  dfds[0][1] = k1 * (3.0 * dev[0][1]);
  dfds[0][2] = k1 * (3.0 * dev[0][2]);
  dfds[1][0] = k1 * (3.0 * dev[1][0]);
  dfds[1][1] = k1 * (3.0 * dev[1][1] + k2);
  dfds[1][2] = k1 * (3.0 * dev[1][2]);
  dfds[2][0] = k1 * (3.0 * dev[2][0]);
  dfds[2][1] = k1 * (3.0 * dev[2][1]);
  dfds[2][2] = k1 * (3.0 * dev[2][2] + k2);
  return;
}

/**
   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 mohrcoulombpar::derplaspotsigma (matrix &sig,matrix &dfds)
{
  double j2, j1, k1, k2;
  matrix dev(3,3);

  deviator(sig, dev);
  j1 = first_invar(dev);
  j2 = second_invar(dev);
  k1 = 1.0 / (sqrt(j2 + sqrt(3.0)*betapsi*sigcpsi*j1));
  k2 = sqrt(3.0) * betapsi * sigcpsi;
  dfds[0][0] = k1 * (3.0 * dev[0][0] + k2);
  dfds[0][1] = k1 * (3.0 * dev[0][1]);
  dfds[0][2] = k1 * (3.0 * dev[0][2]);
  dfds[1][0] = k1 * (3.0 * dev[1][0]);
  dfds[1][1] = k1 * (3.0 * dev[1][1] + k2);
  dfds[1][2] = k1 * (3.0 * dev[1][2]);
  dfds[2][0] = k1 * (3.0 * dev[2][0]);
  dfds[2][1] = k1 * (3.0 * dev[2][1]);
  dfds[2][2] = k1 * (3.0 * dev[2][2] + k2);
  return;
}

/**
  This function computes material stiffnes matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number

*/
void mohrcoulombpar::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
    //  nedodelano
    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 mohrcoulombpar::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;
}


/**
  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 mohrcoulombpar::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 mohrcoulombpar::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 mohrcoulombpar::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 mohrcoulombpar::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 mohrcoulombpar::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      phi=val[0];
      break;
    }
    case 1:{
      c=val[1];
      break;
    }
    case 2:{
      psi=val[2];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  alphaphi = sin(phi)/sqrt(9+3*sin(phi)*sin(phi));
  betaphi  = alphaphi / sqrt(3*(3*c*c - alphaphi*alphaphi));
  sigcphi  = sqrt(3*(c*c - alphaphi*alphaphi/3.0));

  alphapsi = sin(psi)/sqrt(9+3*sin(psi)*sin(psi));
  betapsi  = alphapsi / sqrt(3*(3*c*c - alphapsi*alphapsi));
  sigcpsi  = sqrt(3*(c*c - alphapsi*alphapsi/3.0));
}