j2flow.cpp

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

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

j2flow::j2flow (void)
{
  fs=0.0;  k=0.0;
}
j2flow::~j2flow (void)
{

}

void j2flow::read (XFILE *in)
{
  xfscanf (in,"%lf %lf",&fs,&k);
  sra.read (in);
}

/**
   function returns elastic stiffness matrix
   
   @param d - elastic stiffness matrix
   
   4.8.2001
*/
void j2flow::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
    
    matrix ad(d.m,d.n);
    Mm->elmatstiff (ad,ipp);
    tangentstiff (ad,d,ipp,ido);
    
    //if (Mm->ip[ipp].other[0]>0.0)
    //d[0][0]*=1.0e-4;
  }
}

void j2flow::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);
    deryieldfsigma (sig,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 j2flow::yieldfunction (matrix &sig,vector &q)
{
  double f;
  matrix dev(3,3);

  deviator (sig,dev);
  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 j2flow::deryieldfsigma (matrix &sig,matrix &dfds)
{
  deviator (sig,sig);
  normedtensor (sig,dfds);
}

/**
   function evaluates derivatives of yield function
   with respect of internal variable q
   27.10.2001
*/
void j2flow::deryieldfq (vector &dq)
{
  //dq[0]=k*(-1.0);
  dq[0]=-1.0;
}

void j2flow::plasmod (matrix &h)
{
  h[0][0]=k;
}

double j2flow::plasmodscalar(vector &qtr)
{
  double ret;
  vector dfq(qtr.n), hp(qtr.n);
  matrix h(qtr.n, qtr.n);

  deryieldfq(dfq);
  plasmod (h);
  mxv (h,dfq,hp);
  scprd (hp,dfq,ret);
  return ret;
}

void j2flow::updateq(double dgamma, vector &q)
{
  long j;
  vector dfq(q.n), hp(q.n);
  matrix h(q.n, q.n);

  deryieldfq(dfq);
  plasmod (h);
  mxv (h,dfq,hp);
  plasmod (h);
  for (j=0;j<q.n;j++)
    q[j]-=dgamma*hp[j];
}

void j2flow::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];

  //  new data storage
  /*
  for (i=0;i<n;i++){
    Mm->ip[ipp].other[i]=epsp[i];
  }

  Mm->ip[ipp].other[n]=gamma;
  Mm->ip[ipp].other[n+1]=q[0];
  */
}

void j2flow::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];
}

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



void j2flow::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      fs=val[i];
      break;
    }
    case 1:{
      k=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
}

double j2flow::give_consparam (long ipp,long ido)
{
  long ncompstr;
  double gamma;
  
  ncompstr=Mm->ip[ipp].ncompstr;
  gamma = Mm->ip[ipp].eqother[ido+ncompstr];
  
  return gamma;
}

long j2flow::give_num_interparam ()
{
  return 1;
}

void j2flow::give_interparam (long ipp,long ido,vector &q)
{
  long ncompstr=Mm->ip[ipp].ncompstr;
  
  q[0]=Mm->ip[ipp].eqother[ido+ncompstr+1];
}