chenold.cpp

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

#include <math.h>
#include "chen.h"
#include "global.h"
#include "genfile.h"
#include "intpoints.h"
#include "vecttens.h"
#include "alias.h"
#include "stochdriver.h"
#include "matrix.h"

chen::chen (void)
{
  fyc=0.0;
  fyt=0.0;
  fybc=0.0;
  fc=0.0;
  ft=0.0;
  fbc=0.0;
  epsu=0.0;
  
  ay=0.0;  au=0.0;
  ky=0.0;  ku=0.0;
  alpha=0.0;  beta=0.0;
  
  //  hardening is switched off
  //hard=0;
  //  hardening is switched on
  hard=1;
}

chen::~chen (void)
{

}

/**
   function reads material parameters of the model
   
   21.2.2005
*/
void chen::read (FILE *in)
{
  fscanf (in,"%lf %lf %lf %lf %lf %lf %lf",&fyc,&fyt,&fybc,&fc,&ft,&fbc,&epsu);
  sra.read (in);
}


/**
   function returns elastic stiffness matrix
   
   @param d - elastic stiffness matrix
   
   4.8.2001
*/
void chen::matstiff (matrix &d,long ipp,long ido)
{
  if (Mp->stmat==initial_stiff){
    //  initial elastic matrix
    Mm->elmatstiff (d,ipp);
  }
  if (Mp->stmat==tangent_stiff){
    //  tangent stiffness matrix
    //fprintf (stderr,"\n\n tangent stiffness matrix is not implemented yet\n");
    
    matrix ad(d.m,d.n);
    Mm->elmatstiff (ad,ipp);
    tangentstiff (ad,d,ipp,ido);
  }
}


void chen::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
{
  long ncomp=Mm->ip[ipp].ncompstr;
  double denom,gamma;
  vector str,av(d.m),q(1);
  matrix sig(3,3),am(d.m,d.n);
  
  gamma=Mm->ip[ipp].eqother[ido+ncomp];
  if (gamma<1.0e-10){
    copym (d,td);
  }
  else{
    
    allocv (ncomp,str);
    
    Mm->givestress (0,ipp,str);
    vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
    deryieldfdsigma (sig,q,sig);
    tensor_vector (str,sig,Mm->ip[ipp].ssst,stress);
    
    if (Mm->ip[ipp].ssst==planestress){
      vector auxstr(3);
      auxstr[0]=str[0];auxstr[1]=str[1];auxstr[2]=str[2];
      destrv (str);
      allocv (d.m,str);
      str[0]=auxstr[0];str[1]=auxstr[1];str[2]=auxstr[2];
    }
    
    mxv (d,str,av);
    scprd (av,str,denom);
    
    
    q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
    denom+= plasmodscalar(q);
    
    if (fabs(denom)<1.0e-10){
      copym (d,td);
    }
    else{
      vxv (str,str,am);
      mxm (d,am,td);
      mxm (td,d,am);
      
      cmulm (1.0/denom,am);
      
      subm (d,am,td);
    }
  }
  
}

/**
   function evaluates yield function for given stresses
   
   @param sig - stresses
   @param q - internal parameters (hardening)
   
   4.8.2001
*/
double chen::yieldfunction (matrix &sig,vector &q)
{
  double f,invar,j2,zone,kappa;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  j2=second_invar(dev);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  //  compression-compression
  if(invar<0.0 && zone<0.0)
    {
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      alpha = (au-ay)/(ku-ky);
      beta = (ay*ku-au*ky)/(ku-ky);
      
      kappa = sqrt(ky + q[0]*(ku-ky));
      if (kappa>sqrt(ku))
	kappa=sqrt(ku);
      
      if (hard==0){
	//  yield function without hardening
	f = j2 + ay/3.0*invar-ky;
      }
      if (hard==1){
	//  yield function with hardening
	f = j2 + beta/3.0*invar-kappa*kappa*(1.0-alpha/3.0*invar);
      }
    }
  //  otherwise
  else
    {
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      alpha = (au-ay)/(ku-ky);
      beta = (ay*ku-au*ky)/(ku-ky);
      
      kappa = sqrt(ky + q[0]*(ku-ky));
      if (kappa>sqrt(ku))
	kappa=sqrt(ku);
      
      if (hard==0){
	//  yield function without hardening
	f = j2 - (invar*invar)/6.0+ay/3.0*invar-ky;
      }
      if (hard==1){
	//  yield function with hardening
	f = j2 - (invar*invar)/6.0 + beta/3.0*invar-kappa*kappa*(1.0-alpha/3.0*invar);
      }
      
    }
  //fprintf (Out,"\n a0  %le",a0);
  //fprintf (Out,"\n tau02  %le",tau02);
  
  //  f = tensornorm (dev) - (sqrt(2.0/3.0)*fs+q[0]);
  return f;
}


/**
   function evaluates derivatives of yield function
   with respect of stress components
   4.8.2001
*/
void chen::deryieldfdsigma (matrix &sig,vector &q,matrix &dfds)
{
  double invar,j2,zone,kappa;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  j2=second_invar(dev);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  //  compression-compression
  if(invar<0.0 && zone<0.0)
    {
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      alpha = (au-ay)/(ku-ky);
      beta = (ay*ku-au*ky)/(ku-ky);
      
      kappa = sqrt(ky + q[0]*(ku-ky));
      if (kappa>sqrt(ku))
	kappa=sqrt(ku);
      
      if (hard==0){
	dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
	dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
	dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + beta/3.0+alpha/3.0*kappa*kappa;
      }
      if (hard==1){
	dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0;
	dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0;
	dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0;
      }

      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];
    }
  
  
  
  //  otherwise
  else
    { 
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      alpha = (au-ay)/(ku-ky);
      beta = (ay*ku-au*ky)/(ku-ky);

      kappa = sqrt(ky + q[0]*(ku-ky));
      if (kappa>sqrt(ku))
	kappa=sqrt(ku);
      
      if (hard==0){
	dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
	dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + beta/3.0+alpha/3.0*kappa*kappa;
	dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + beta/3.0+alpha/3.0*kappa*kappa;
      }
      if (hard==1){
	dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0;
	dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0;
	dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0;
      }
      
      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];

      dfds[0][0]-=1.0/3.0*invar;
      dfds[1][1]-=1.0/3.0*invar;
      dfds[2][2]-=1.0/3.0*invar;

      //fprintf (Out,"\n derivace: jsme v otherwise");
    }
  
}

/**
   function evaluates derivatives of yield function
   with respect of stress components
   4.8.2001
*/
void chen::deryieldfdsigma_old (matrix &sig,matrix &dfds)
{
  double invar,j2,zone;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  j2=second_invar(dev);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  //compression-compression
  if(invar<0.0 && zone<0.0)
    {
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);

      dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2])+ay/3.0;
      dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2])+ay/3.0;
      dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1])+ay/3.0;

      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];
    }
  
  
  
  //otherwise
  else
    { 
      ay = (fyc-fyt)/2.0;

      dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2])+ay/3.0;
      dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2])+ay/3.0;
      dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1])+ay/3.0;

      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];

      dfds[0][0]-=1.0/3.0*invar;
      dfds[1][1]-=1.0/3.0*invar;
      dfds[2][2]-=1.0/3.0*invar;

      //fprintf (Out,"\n derivace: jsme v otherwise");
    }
  
}


/**
   function computes second derivatives of yiled function with
   respect to stresses
*/
void chen::deryieldfdsigmadsigma (matrix &sig,matrix &dfdsds)
{
  double invar,j2,zone;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  j2=second_invar(dev);
  
  fillm (0.0,dfdsds);
  
  zone = sqrt(j2)+invar/sqrt(3.0);

  //  compression-compression
  if(invar<0.0 && zone<0.0){
    dfdsds[0][0] =  2.0/3.0;
    dfdsds[0][1] = -1.0/3.0;
    dfdsds[0][2] = -1.0/3.0;
    
    dfdsds[1][0] = -1.0/3.0;
    dfdsds[1][1] =  2.0/3.0;
    dfdsds[1][2] = -1.0/3.0;

    dfdsds[2][0] = -1.0/3.0;
    dfdsds[2][1] = -1.0/3.0;
    dfdsds[2][2] =  2.0/3.0;

    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
  }
  
  //  otherwise
  else{ 
    dfdsds[0][0] =  1.0/3.0;
    dfdsds[0][1] = -2.0/3.0;
    dfdsds[0][2] = -2.0/3.0;
    
    dfdsds[1][0] = -2.0/3.0;
    dfdsds[1][1] =  1.0/3.0;
    dfdsds[1][2] = -2.0/3.0;
    
    dfdsds[2][0] = -2.0/3.0;
    dfdsds[2][1] = -2.0/3.0;
    dfdsds[2][2] =  1.0/3.0;

    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
  }

}

/**
   function evaluates derivatives of yield function
   with respect of stress components
   4.8.2001
*/
void chen::deryieldfdsigma_old_old (matrix &sig,matrix &dfds)
{
  double invar,j2,a0;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  normedtensor (dev,dfds);    
  j2=second_invar(dev);
  
  //compression-compression
  if(invar<0.0 && (sqrt(j2)+invar/sqrt(3.0))<0.0)
    {
      a0 = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      dfds[0][0]+=a0/3.0;
      dfds[1][1]+=a0/3.0;
      dfds[2][2]+=a0/3.0;
    }
  
  
  
  //otherwise
  else
    { 
      a0 = (fyc-fyt)/2.0;
      dfds[0][0]+=a0/3.0-1.0/3.0*invar;
      dfds[1][1]+=a0/3.0-1.0/3.0*invar;
      dfds[2][2]+=a0/3.0-1.0/3.0*invar;
    }
  
}

/**
   function evaluates derivatives of yield function
   with respect of internal variable q
   27.10.2001
*/
void chen::deryieldfdq (matrix &sig,vector &q,vector &dfdq)
{
  double invar,j2,zone,kappa;
  matrix dev(3,3);
  
  if (hard==0){
    nullv (dfdq.a,dfdq.n);
  }
  if (hard==1){
    deviator (sig,dev);
    invar=first_invar (sig);
    j2=second_invar(dev);