elastisomat.cpp

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

#include "elastisomat.h"
#include "matrix.h"
#include "vector.h"
#include "stochdriver.h"
#include "global.h"
#include "intpoints.h"

elastisomat::elastisomat (void)
{
  e = 0.0;  nu = 0.0;
}
elastisomat::~elastisomat (void)
{

}

/**
   function reads material parameters
   
   @param in - input stream
*/
void elastisomat::read (XFILE *in)
{
  xfscanf (in,"%lf %lf",&e,&nu);
}

void elastisomat::matstiff (matrix &d,strastrestate ssst)
{
  switch (ssst){
  case bar:{
    matstiff_bar (d);
    break;
  }
  case plbeam:{
    matstiff_plbeam (d);
    break;
  }
  case spacebeam:{
    matstiff_spacebeam (d);
    break;
  }
  case planestress:{
    matstiff_plstress (d);
    break;
  }
  case planestrain:{
    matstiff_plstrain (d);
    break;
  }
  case plate:{
    matstiff_plate (d);
    break;
  }
  case axisymm:{
    matstiff_axi (d);
    break;
  }
  case spacestress:{
    matstiff_spacestr (d);
    break;
  }
  default:{
    fprintf (stderr,"\n unknown number of components of stress tensor is required");
    fprintf (stderr,"\n in function elastisomat::matstiff (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

void elastisomat::elmatstiff (matrix &d,strastrestate ssst)
{
  switch (ssst){
  case bar:{
    matstiff_bar (d);
    break;
  }
  case plbeam:{
    matstiff_plbeam (d);
    break;
  }
  case spacebeam:{
    matstiff_spacebeam (d);
    break;
  }
  case planestress:{
    matstiff_plstress (d);
    break;
  }
  case planestrain:{
    matstiff_plstrain (d);
    break;
  }
  case plate:{
    matstiff_plate (d);
    break;
  }
  case axisymm:{
    matstiff_axi (d);
    break;
  }
  case spacestress:{
    matstiff_spacestr (d);
    break;
  }
  default:{
    fprintf (stderr,"\n unknown number of components of stress tensor is required");
    fprintf (stderr,"\n in function elastisomat::elmatstiff (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function creates stiffness matrix of the elastic
   isotropic material for bar elements
   
   d - stiffness matrix of the material
   
   11.9.2001
*/
void elastisomat::matstiff_bar (matrix &d)
{
  d[0][0] = e;
}

/**
   function creates stiffness matrix of the elastic
   isotropic material for plane beam elements
   
   @param d - stiffness matrix of the material
   
   11.9.2001
*/
void elastisomat::matstiff_plbeam (matrix &d)
{
  d[0][0] = e;
  d[1][1] = e/2.0/(1.0+nu);
  d[2][2] = e;
}

/**
   function creates stiffness matrix of the elastic
   isotropic material for plane beam elements
   
   @param d - stiffness matrix of the material
   
   14.3.2003
*/
void elastisomat::matstiff_spacebeam (matrix &d)
{
  d[0][0] = e;
  d[1][1] = e/2.0/(1.0+nu);
  d[2][2] = e/2.0/(1.0+nu);
  d[3][3] = e/2.0/(1.0+nu);
  d[4][4] = e;
  d[5][5] = e;
}
/**
   function creates stiffness matrix of the elastic
   isotropic material for 2D problems (plane stress)

   @param d - stiffness matrix of the material

   19.7.2001
*/
void elastisomat::matstiff_plstress (matrix &d)
{
  double c;
  
  fillm(0.0,d);

  c = e/(1.0-nu*nu);
  
  d[0][0] = c;     d[0][1] = c*nu;  d[0][2] = 0.0;
  d[1][0] = c*nu;  d[1][1] = c;     d[1][2] = 0.0;
  d[2][0] = 0.0;   d[2][1] = 0.0;   d[2][2] = e/2.0/(1.0+nu);
}

/**
   function creates stiffness matrix of the elastic
   isotropic material for 2D problems (plane strain)

   @param d - stiffness matrix of the material

   19.7.2001
*/
void elastisomat::matstiff_plstrain (matrix &d)
{
  double c;
  
  fillm(0.0,d);

  c = e/(1.0+nu)/(1.0-2.0*nu);
  
  d[0][0] = c*(1.0-nu);   d[0][1] = c*nu;         d[0][2] = 0.0;
  d[1][0] = c*nu;         d[1][1] = c*(1.0-nu);   d[1][2] = 0.0;
  d[2][0] = 0.0;          d[2][1] = 0.0;          d[2][2] = e/2.0/(1.0+nu);

  if (d.m > 3)
    {
      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];
    }
}

void elastisomat::matstiff_axi (matrix &d)
{
  double g,s;
  
  fillm(0.0,d);
  
  g = e/2.0/(1.0+nu);
  s = e/(1.0+nu)/(1.0-2.0*nu);
  
  d[0][0]=s*(1-nu);  d[0][1]=s*nu;     d[0][2]=d[0][1];
  d[1][0]=d[0][1];   d[1][1]=d[0][0];  d[1][2]=d[0][1];
  d[2][0]=d[0][1];   d[2][1]=d[0][1];  d[2][2]=d[0][0];

  d[3][3]=g;

}

/**
   function creates stiffness matrix of the elastic
   isotropic material for plate elements
   
   @param d - stiffness matrix
   
   19.7.2001
*/
void elastisomat::matstiff_plate (matrix &d)
{
  double c,g;
  
  fillm(0.0,d);

  c = e/12.0/(1.0-nu*nu);
  g = e/2.0/(1.0+nu);
  
  d[0][0]=c;        d[0][1]=c*nu;  d[0][2]=0.0;
  d[1][0]=d[0][1];  d[1][1]=c;     d[1][2]=0.0;
  d[2][0]=0.0;      d[2][1]=0.0;   d[2][2]=g/12.0;
  
  d[3][3]=g;  d[4][4]=g;
}

/**
   function creates stiffness matrix of the elastic
   isotropic material for 3D problems
   
   @param d - stiffness matrix of the material

   19.7.2001
*/
void elastisomat::matstiff_spacestr (matrix &d)
{
  double g,s;
  
  fillm(0.0,d);
  
  g = e/2.0/(1.0+nu);
  s = e/(1.0+nu)/(1.0-2.0*nu);
  
  d[0][0]=s*(1-nu);  d[0][1]=s*nu;     d[0][2]=s*nu;
  d[1][0]=d[0][1];   d[1][1]=d[0][0];  d[1][2]=d[0][1];
  d[2][0]=d[0][1];   d[2][1]=d[0][1];  d[2][2]=d[0][0];

  d[3][3]=g;         d[4][4]=g;        d[5][5]=g;
}





void elastisomat::matcompl (matrix &c,strastrestate ssst)
{
  switch (ssst){
  case bar:{
    matcompl_bar (c);
    break;
  }
  case plbeam:{
    matcompl_plbeam (c);
    break;
  }
  case planestress:{
    matcompl_plstress (c);
    break;
  }
  case planestrain:{
    matcompl_plstrain (c);
    break;
  }
  case axisymm:{
    matcompl_axi (c);
    break;
  }
  case spacestress:{
    matcompl_spacestr (c);
    break;
  }
  default:{
    fprintf (stderr,"\n unknown number of components of stress tensor is required");
    fprintf (stderr,"\n in function elastisomat::matcompl (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function creates compliance matrix of the elastic
   isotropic material for bar elements
   
   c - compliance matrix of the material
   
   5.11.2002
*/
void elastisomat::matcompl_bar (matrix &c)
{
  c[0][0] = 1.0/e;
}

/**
   function creates compliance matrix of the elastic
   isotropic material for plane beam elements
   
   @param c - complinace matrix of the material
   
   5.11.2002
*/
void elastisomat::matcompl_plbeam (matrix &c)
{
  c[0][0] = 1.0/e;
  c[1][1] = 2.0*(1.0+nu)/e;
  c[2][2] = 1.0/e;
}

/**
   function creates compliance matrix of the elastic
   isotropic material for 2D problems (plane stress)

   @param c - compliance matrix of the material

   5.11.2002
*/
void elastisomat::matcompl_plstress (matrix &c)
{
  fillm(0.0,c);
  
  c[0][0] =  1.0/e;     c[0][1] = -1.0*nu/e;  c[0][2] = 0.0;
  c[1][0] = -1.0*nu/e;  c[1][1] =  1.0/e;     c[1][2] = 0.0;
  c[2][0] = 0.0;        c[2][1] = 0.0;        c[2][2] = 2.0*(1.0+nu)/e;
}

/**
   function creates compliance matrix of the elastic
   isotropic material for 2D problems (plane strain)

   @param c - compliance matrix of the material

   5.11.2002
*/
void elastisomat::matcompl_plstrain (matrix &c)
{
  double g;
  
  fillm(0.0,c);
  
  g = (1.0+nu)/e;
  
  c[0][0] = g*(1.0-nu);   c[0][1] = -1.0*g*nu;    c[0][2] = 0.0;
  c[1][0] = -1.0*g*nu;    c[1][1] = g*(1.0-nu);   c[1][2] = 0.0;
  c[2][0] = 0.0;          c[2][1] = 0.0;          c[2][2] = 2.0*g;
}


/**
   function creates compliance matrix of the elastic
   isotropic material for axisymmetric problems
   
   @param c - compliance matrix of the material

   8.4.2005, JK
*/
void elastisomat::matcompl_axi (matrix &c)
{
  double g;
  fillm(0.0,c);
  
  g = 2.0*(1.0+nu)/e;
  
  c[0][0]=1.0/e;     c[0][1]=-1.0*nu/e;  c[0][2]=c[0][1];    c[0][3]=0.0;
  c[1][0]=c[0][1];   c[1][1]=c[0][0];    c[1][2]=c[0][1];    c[1][3]=0.0;
  c[2][0]=c[0][1];   c[2][1]=c[0][1];    c[2][2]=c[0][0];    c[2][3]=0.0;

  c[3][0]=c[0][3];   c[3][1]=c[1][3];    c[3][2]=c[2][3];    c[3][3]=g;
}

/**
   function creates compliance matrix of the elastic
   isotropic material for 3D problems
   
   @param c - compliance matrix of the material

   5.11.2002
*/
void elastisomat::matcompl_spacestr (matrix &c)
{
  double g;
  fillm(0.0,c);
  
  g = 2.0*(1.0+nu)/e;
  
  c[0][0]=1.0/e;     c[0][1]=-1.0*nu/e;  c[0][2]=c[0][1];
  c[1][0]=c[0][1];   c[1][1]=c[0][0];    c[1][2]=c[0][1];
  c[2][0]=c[0][1];   c[2][1]=c[0][1];    c[2][2]=c[0][0];

  c[3][3]=g;         c[4][4]=g;          c[5][5]=g;
}

void elastisomat::nlstresses (long ipp)
{
  long i, n = Mm->ip[ipp].ncompstr;
  vector eps(n),sig(n);
  strastrestate ssst = Mm->ip[ipp].ssst;
  matrix d(n,n);
  
  //  initial values
  for (i=0;i<n;i++){
    eps[i]=Mm->ip[ipp].strain[i];
  }
  
  matstiff (d,ssst);
  mxv (d,eps,sig);
  if (Mm->ip[ipp].ssst == planestress)
    Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (eps[0]+eps[1]);

  for (i=0;i<n;i++){
    Mm->ip[ipp].stress[i]=sig[i];
  }
  
}

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