mohrc.cpp

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

{  
  return 0.0;
}


/**
   Function updates hardening/softening parameter in case that hardening/softening rule is adopted

   @param ipp   - integration point pointer
   @param epsp  - plastic strain components
   @param q     - hardening parameters

   4.11.2003
*/
void mohrcoulomb::updateq(long ipp, vector &epsp, vector &q)
{
  return;
}

/**
   function checks ordering of principal stresses and vertex singularity

   @param psig  - principal stress vector sorted in this way psig[0] < psig[1] < psig[2]

   @retval 12 - in case that the psig[0] < psig[1] condition is not statisfied
   @retval 23 - in case that the psig[1] < psig[2] condition is not statisfied
   @retval 1  - in case that the vertex singularity occurs i.e. tensile strength is exceeded
   @retval 0  - in other cases

   4.11.2003
*/
long mohrcoulomb::checkpsig(vector &psig)
{
  double s1, s2, sigma, tau, maxsigt;

  // vertex return
  s1 = psig[0];
  s2 = psig[2];
  sigma   =  (s1 + s2) / 2.0;
  tau     = -(s1 - s2) / (2.0 * cos(phi));
  maxsigt =  c / tan(phi);
  if ((phi != 0.0) && (tau < ((sigma - maxsigt)/tan(phi))))
    return 1;

  // double vector return
  if (psig[0] > psig[1])
    return 12;
  if (psig[1] > psig[2])
    return 23;

  // normal return
  return 0;
}



/**
   function checks ordering of principal stresses and vertex singularity

   @param psig  - principal stress vector sorted in this way psig[0] < psig[1] < psig[2]

   @retval It returns index of zero psig

   4.11.2003
*/
long mohrcoulomb::checkzeropsig(vector &psig)
{
  long i;

  for (i=0; i<3; i++)
  {
    if (psig[i] == 0.0)
      return i;
  }
  return -1;  
}



/**
   Function returns stresses on sufrace of plasticity
   multisurface cutting plane method at principal stresses is used.

   gamma,epsp and q will be replaced by new values

   @param ipp - integration point pointer
   @param gamma - consistency parameter
   @param epsn - total strain components
   @param epsp - plastic strain components
   @param q - hardening parameters
   @param mu - active yield surfaces indicator
   @param ni - maximum number of iterations
   @param err - required error
   
   4.8.2001
*/
void mohrcoulomb::mc_msurf_cp(long ipp, double &gamma, vector &epsn, vector &epsp, vector &q, long mu,long ni,double err)
{
  long i,j,ncomp=epsn.n,idem,zps;
  double e,nu;
  long stat[2], nas;
  vector epsa(ncomp), sig(ncomp), epsptr(ncomp), depsp(ncomp);
  vector psig(3),peps(3),pdepsp(3),pepsa(3),dpepsp(3), f(2);
  matrix d(ncomp,ncomp),sigt(3,3),epsat(3,3),epspt(3,3);
  matrix pd(3,3),pvect(3,3);
  matrix dfds, dgds, am, cpm, hcpm;
  vector dgamma, ff;   

  // Initialization
  // checking for the correct type of the elastic material
  idem = Mm->ip[ipp].gemid();
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\nError - The Mohr-Coulomb material is combined with\n");
    fprintf(stderr, "        the nonsupported elastic material.\n");
    fprintf(stderr, "        Only elisomat type is supported\n");
    exit(0);
  }

  e  = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
  // elastic stiffness matrix computation
  Mm->elmatstiff(d, ipp);
  pelmatstiff (ipp, pd, e, nu);
  // computation complementary values for plain stress/strain
  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];
  }
/*  if (Mm->ip[ipp].ssst == planestress)
    epsn[3] = nu / (1.0 - nu) * (epsn[0]+epsn[1]);*/
  copyv(epsp, epsptr);

  // Computing principal directions - they do not change during stress return
  // elastic strain
  subv (epsn, epsp, epsa);
  // trial stress computation
  mxv (d, epsa, sig);
  // principal stresses and their directions
  vector_tensor (sig, sigt, Mm->ip[ipp].ssst, stress);
  princ_val (sigt, psig, pvect, nijac, limit, Mp->zero, 3);
  // transformation elastic strain into the computed principal directions
  vector_tensor (epsa, epsat, Mm->ip[ipp].ssst, strain);
  glmatrixtransf(epsat, pvect);
  peps[0] = epsat[0][0];
  peps[1] = epsat[1][1];
  peps[2] = epsat[2][2];
  if (Mm->ip[ipp].ssst == planestress)
  {
    zps = checkzeropsig(psig);
    fillrow(0.0, zps, pd);
    fillcol(0.0, zps, pd);
  }

 // Main iteration loop
 // number of yield surfaces is always 2, for every mu value
 for (i=0; i<ni; i++)
 {
    //  elastic strain
    subv (peps, pdepsp, pepsa);
    //  trial stress computation
    mxv (pd, pepsa, psig);
    //  yield function control
    yieldfunction(psig, mu, f);
    // detection of active yield surfaces
    nas=0;
    for (j=0; j<2; j++)
    {
      if (f[j] >= err)
      {
        stat[j] = 1;
        nas++;
      }
      else
        stat[j] = 0;
    }
    if (nas==0)  // no active yield surface
      break;
    
    // allocation of necessarry matrices and vectors        
    allocm (nas, 3, dfds); allocm (nas, 3, dgds);
    allocm (nas, 3, am);   allocm (nas, nas, cpm);
    allocm (nas, nas, hcpm);
    allocv (nas, dgamma);
    allocv (nas, ff);
    yieldfunction(psig, mu, stat, ff);
    // computing matrices of derivations
    dfdsigma(stat, mu, dfds);
    dgdsigma(stat, mu, dgds);
    plasmodscalar(ipp, q, mu, hcpm);
    // assembling resulting matrix for the dgamma solution
    mxm(dfds, pd, am);
    mxmt(am, dgds, cpm);
    addm(cpm, hcpm, cpm);
    // solution of the equation system
    gemp(cpm.a, dgamma.a, ff.a, nas, 1, Mp->zero, 1);
    // new plastic strain vector increments
    mtxv(dgds, dgamma, dpepsp);
    // compute new principal plastic strain increments and their transformed values into g.c.s.
    for (j=0; j<3; j++)
    {
      epspt[j][j]+=dpepsp[j];
      pdepsp[j]+=dpepsp[j];
    }
    for (j=0; j<nas; j++)
      gamma += dgamma[j];
    lgmatrixtransf(epspt, pvect);
    tensor_vector(depsp, epspt, Mm->ip[ipp].ssst, strain);
    addv(epsp, depsp, epsptr);
    // update hardening parameter
    updateq(ipp,epsptr, q);
    // deallocating matrices and vectors used in this loop
    destrm (dfds); destrm (dgds);
    destrm (am);   destrm (cpm);
    destrm (hcpm);
    destrv (dgamma);
    destrv(ff);
  }
  // updating plastic strains
  copyv(epsptr, epsp);
  //  elastic strain
  subv (epsn, epsp, epsa);
  //  trial stress computation
  mxv (d, epsa, sig);
  //  storing resulting stresses
  Mm->storestress (0, ipp, sig);
  return;
}

void mohrcoulomb::plasmodscalar(long ipp, vector &q, long mu, matrix &hcpm)
{
  fillm(0.0, hcpm);
}  

void mohrcoulomb::yieldfunction(vector &psig, long mu, long *stat, vector &f)
{
  long i = 0;
  double k1, k2, s1, s2;

  s1 = psig[0];
  s2 = psig[2];
  k1 = -1.0 + sin(phi);  k2 = 1.0 + sin(phi);

  if (stat[0])
  {
    s1 = psig[0];
    s2 = psig[2];
    f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
    i++;
  }
  if (stat[1])
  {
    switch (mu)
    {
      case 12 : // swap s1 and s2 at the condition s1 < s2 < s3
        s1 = psig[1];
        s2 = psig[2];
        f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
        break;
      case 23 : // swap s2 and s3 at the condition s1 < s2 < s3
        s1 = psig[0];
        s2 = psig[1];
        f[i] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
        break;
      case 1 : // tension strength cutoff
        f[i] = s1 + s2 - 2.0*c/tan(phi);
        break;
      default :
        break;
    }
  }
  return;
} 

void mohrcoulomb::yieldfunction(vector &psig, long mu, vector &f)
{
  double k1, k2, s1, s2;

  k1 = -1.0 + sin(phi);  k2 = 1.0 + sin(phi);
  s1 = psig[0];
  s2 = psig[2];

  f[0] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
  switch (mu)
  {
    case 12 : // swap s1 and s2 at the condition s1 < s2 < s3
      s1 = psig[1];
      s2 = psig[2];
      f[1] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
      break;
    case 23 : // swap s2 and s3 at the condition s1 < s2 < s3
      s1 = psig[0];
      s2 = psig[1];
      f[1] = k1*s1 + k2*s2 - 2.0*c*cos(phi);
      break;
    case 1 : // tension strength cutoff
      f[1] = s1 + s2 - 2.0*c/tan(phi);
      break;
    default :
      break;
  }
  return;
}
 
void mohrcoulomb::dfdsigma(long *stat, long mu, matrix &dfds)
{
  long i = 0;
  if (stat[0])
  {
    dfds[i][0] = -1.0 + sin(phi);
    dfds[i][1] =  0.0;
    dfds[i][2] =  1.0 + sin(phi);
    i++;
  }
  if (stat[1])
  {
    switch (mu)
    {
      case 12 : // swap s1 and s2 at the condition s1 < s2 < s3
        dfds[i][0] =  0.0;
        dfds[i][1] = -1.0 + sin(phi);
        dfds[i][2] =  1.0 + sin(phi);
        break;
      case 23 : // swap s2 and s3 at the condition s1 < s2 < s3
        dfds[i][0] = -1.0 + sin(phi);
        dfds[i][1] =  1.0 + sin(phi);
        dfds[i][2] =  0.0;
        break;
      case 1 : // tension strength cutoff
        dfds[i][0] = 1.0;
        dfds[i][1] = 0.0;
        dfds[i][2] = 1.0;
        break;
      default :
        break;
    }
  }
  return;
}

void mohrcoulomb::dgdsigma(long *stat, long mu, matrix &dgds)
{
  long i = 0;
  if (stat[0])
  {
    dgds[i][0] = -1.0 + sin(psi);
    dgds[i][1] =  0.0;
    dgds[i][2] =  1.0 + sin(psi);
    i++;
  }
  if (stat[1])
  {
    switch (mu)
    {
      case 12 : // swap s1 and s2 at the condition s1 < s2 < s3
        dgds[i][0] =  0.0;
        dgds[i][1] = -1.0 + sin(psi);
        dgds[i][2] =  1.0 + sin(psi);
        break;
      case 23 : // swap s2 and s3 at the condition s1 < s2 < s3
        dgds[i][0] = -1.0 + sin(psi);
        dgds[i][1] =  1.0 + sin(psi);
        dgds[i][2] =  0.0;
        break;
    case 1 : // tension strength cutoff
        dgds[i][0] = 1.0;
        dgds[i][1] = 0.0;
        dgds[i][2] = 1.0;
        break;
      default :
        break;
    }
  }
  return;
}


void mohrcoulomb::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 mohrcoulomb::giveirrstrains (long ipp, long ido, vector &epsp)
{
  long i;
  for (i=0;i<epsp.n;i++)
    epsp[i] = Mm->ip[ipp].eqother[ido+i];
}



double mohrcoulomb::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 mohrcoulomb::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__);
    }
    }
  }
}