drprag.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 "drprag.h"
#include "global.h"
#include "vecttens.h"
#include "intpoints.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"



#define nijac 20
#define limit 1.0e-8

/**
  This constructor inializes attributes to zero values.
*/
drprag::drprag (void)
{
  phi=0.0;  c=0.0;  psi=0.0;
  theta = 0.0; clim = 0.0;
}

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

}


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

  @param in - pointer to the opned text file

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

  alpha=6.0*sin(phi)/(3.0-sin(phi));
  alpha1=6.0*sin(psi)/(3.0-sin(psi));

  beta=6.0*cos(phi)/(3.0-sin(phi));
}

/**
   function computes the value of yield functions

   @param sig - stress components
   @param q   - %vector of hardening parameters

   25.3.2002
*/
double drprag::yieldfunction (matrix &sig, vector &q)
{
  double i1s,j2s,f;
  matrix dev(3,3);

  i1s = first_invar (sig)/3.0;

  deviator (sig,dev);

  j2s = second_invar (dev);

  f = sqrt(3.0*j2s) + alpha*i1s - beta*cohesion(q);

  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 drprag::deryieldfsigma (matrix &sig,matrix &dfds)
{
  double k;
  matrix dev(3,3);

  deviator (sig,dev);
  normedtensor (dev,dfds);
  k=sqrt(6.0)/2.0;
  cmulm (k,dfds);

  dfds[0][0]+=alpha/3.0;
  dfds[1][1]+=alpha/3.0;
  dfds[2][2]+=alpha/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 drprag::derpotsigma (matrix &sig,matrix &dgds)
{
  double k;
  matrix dev(3,3);

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

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

/**
  This function computes value of cohesion which depends on the
  hardening parameter.

  @param qtr - %vector of hardening parameters

  @retval This function returns value of the cohesion coresponding with given hardening parameter.
*/
double drprag::cohesion(vector &qtr)
{
  double tgt = tan(theta);
  double cq = c + tgt*qtr[0];
  if (((cq < clim) && (c < clim)) ||
      ((cq > clim) && (c > clim)))
    return cq;
  return clim;
}

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

   @param qtr - %vector of the hardening parameters
   @param dfds - %matrix where the resulting derivatives are stored


   4.1.2002
*/
void drprag::deryieldfq(vector &qtr, vector &dfq)
{
  double tgt = tan(theta);
  double cq = c + tgt*qtr[0];
  if (((cq < clim) && (c < clim)) ||
      ((cq > clim) && (c > clim)))
    dfq[0] = -beta*tgt;
  else
    dfq[0] = 0.0;
  return;
}

/**
   This function computes derivatives of hardening paramters
   with respect of consistency parameter gamma.

   @param dqdg - %matrix where the resulting derivatives are stored


   4.1.2002
*/
void drprag::der_q_gamma(vector &dqdg)
{
  dqdg[0] = sqrt(1 + 2/9*alpha1*alpha1);
}

/**
  This function computes plastic modulus.

  @param qtr -%vector of hardening parameters

  @retval The function returns value of the plastic modulus.
*/
double drprag::plasmodscalar(vector &qtr)
{
  double ret;
  vector dfq(qtr.n);
  vector dqg(qtr.n);

  deryieldfq(qtr, dfq);
  der_q_gamma(dqg);
  scprd(dfq, dqg, ret);
//  ret = 0.0;

  return -ret;
}

/**
  This function computes new value of teh hardening parameter q.

  @param ipp  - integration point pointer
  @param epsp - %vector of the reached plastic strains
  @param q    - %vector of the hardening parameters

*/
void drprag::updateq(long ipp, vector &epsp, vector &q)
{
  matrix epst(3,3);
  strastrestate ssst=Mm->ip[ipp].ssst;

  vector_tensor (epsp, epst, ssst, strain);
  cumulstrain(epst, q[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 drprag::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 drprag::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].eqother[ido+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];
}

/**
  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 drprag::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];
}

/**
  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 drprag::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 drprag::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 drprag::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 drprag::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;
    }
    case 3:{
      theta=val[i];
      break;
    }
    case 4:{
      clim=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }

  alpha=6.0*sin(phi)/(3.0-sin(phi));
  alpha1=6.0*sin(psi)/(3.0-sin(psi));

  beta=6.0*cos(phi)/(3.0-sin(phi));
}