vecttens.cpp

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

#include "vecttens.h"
#include "matrix.h"
#include "vector.h"



/**
   function creates tensor components from %vector components
   
   @param v - %vector components
   @param t - tensor components
   @param ssst - strain/stress state
   @param ss - strain/stress identifier
   
   JK, 25.3.2002
*/
void vector_tensor (vector &v,matrix &t,strastrestate ssst,strastre ss)
{
  fillm (0.0,t);
  
  if (ss==strain){
    switch (ssst){
    case bar:{
      t[0][0]=v[0];
      break;
    }
    case plbeam:{
      break;
    }
    case planestress:{
      t[0][0]=v[0];      t[0][1]=v[2]/2.0;  t[0][2]=0.0;
      t[1][0]=v[2]/2.0;  t[1][1]=v[1];      t[1][2]=0.0;
      t[2][0]=0.0;       t[2][1]=0.0;       t[2][2]=v[3];
      break;
    }
    case planestrain:{
      t[0][0]=v[0];      t[0][1]=v[2]/2.0;  t[0][2]=0.0;
      t[1][0]=v[2]/2.0;  t[1][1]=v[1];      t[1][2]=0.0;
      t[2][0]=0.0;       t[2][1]=0.0;       t[2][2]=v[3];
      break;
    }
    case axisymm:{
      t[0][0]=v[0];      t[0][1]=v[3]/2.0;  t[0][2]=0.0;
      t[1][0]=v[3]/2.0;  t[1][1]=v[1];      t[1][2]=0.0;
      t[2][0]=0.0;       t[2][1]=0.0;       t[2][2]=v[2];
      break;
    }
    case spacestress:{
      t[0][0]=v[0];      t[0][1]=v[5]/2.0;  t[0][2]=v[4]/2.0;
      t[1][0]=v[5]/2.0;  t[1][1]=v[1];      t[1][2]=v[3]/2.0;
      t[2][0]=v[4]/2.0;  t[2][1]=v[3]/2.0;  t[2][2]=v[2];
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown strain state is required in function vector_tensor (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  if (ss==stress){
    switch (ssst){
    case bar:{
      t[0][0]=v[0];
      break;
    }
    case plbeam:{
      break;
    }
    case planestress:{
      t[0][0]=v[0];  t[0][1]=v[2];  t[0][2]=0.0;
      t[1][0]=v[2];  t[1][1]=v[1];  t[1][2]=0.0;
      t[2][0]=0.0;   t[2][1]=0.0;   t[2][2]=v[3];
      break;
    }
    case planestrain:{
      t[0][0]=v[0];  t[0][1]=v[2];  t[0][2]=0.0;
      t[1][0]=v[2];  t[1][1]=v[1];  t[1][2]=0.0;
      t[2][0]=0.0;   t[2][1]=0.0;   t[2][2]=v[3];
      break;
    }
    case axisymm:{
      t[0][0]=v[0];  t[0][1]=v[3];  t[0][2]=0.0;
      t[1][0]=v[3];  t[1][1]=v[1];  t[1][2]=0.0;
      t[2][0]=0.0;   t[2][1]=0.0;   t[2][2]=v[2];
      break;
    }
    case spacestress:{
      t[0][0]=v[0];  t[0][1]=v[5];  t[0][2]=v[4];
      t[1][0]=v[5];  t[1][1]=v[1];  t[1][2]=v[3];
      t[2][0]=v[4];  t[2][1]=v[3];  t[2][2]=v[2];
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown stress state is required in function vector_tensor (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }

}

/**
   function creates %vector components from tensor components
   
   @param v - %vector components
   @param t - tensor components
   @param ssst - strain/stress state
   @param ss - strain/stress identifier
   
   JK, 25.3.2002
*/
void tensor_vector (vector &v,matrix &t,strastrestate ssst,strastre ss)
{
  fillv (0.0,v);
  
  if (ss==strain){
    switch (ssst){
    case bar:{
      v[0]=t[0][0];
      break;
    }
    case plbeam:{
      break;
    }
    case planestress:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[0][1]*2.0;  v[3]=t[2][2];
      break;
    }
    case planestrain:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[0][1]*2.0;  v[3]=t[2][2];
      break;
    }
    case axisymm:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[2][2];  v[3]=t[0][1]*2.0;
      break;
    }
    case spacestress:{
      v[0]=t[0][0];      v[1]=t[1][1];      v[2]=t[2][2];
      v[3]=t[1][2]*2.0;  v[4]=t[0][2]*2.0;  v[5]=t[0][1]*2.0;
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown strain state is required in function tensor_vector (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  if (ss==stress){
    switch (ssst){
    case bar:{
      v[0]=t[0][0];
      break;
    }
    case plbeam:{
      break;
    }
    case planestress:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[0][1];  v[3]=t[2][2];
      break;
    }
    case planestrain:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[0][1];  v[3]=t[2][2];
      break;
    }
    case axisymm:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[2][2];  v[3]=t[0][1];
      break;
    }
    case spacestress:{
      v[0]=t[0][0];  v[1]=t[1][1];  v[2]=t[2][2];
      v[3]=t[1][2];  v[4]=t[0][2];  v[5]=t[0][1];
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown stress state is required in function tensor_vector (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }

}

/**
   function creates auxiliary %matrix m
   (useful for thermal strains, pore pressure, etc.)

   @param ssst - strain/stress state
   @param m - auxiliary %matrix
   
   JK, 24.10.2004
*/
void tensor_vector_matrix (strastrestate ssst,matrix &m)
{
  fillm (0.0,m);
  
  switch (ssst){
  case bar:{
    m[0][0]=1.0;
    break;
  }
  case planestress:{
    m[0][0]=1.0;
    m[1][0]=1.0;
    m[2][0]=0.0;
    break;
  }
  case axisymm:{
    m[0][0]=1.0;
    m[1][0]=1.0;
    m[2][0]=1.0;
    m[3][0]=0.0;
    break;
  }
  case spacestress:{
    m[0][0]=1.0;
    m[1][0]=1.0;
    m[2][0]=1.0;
    m[3][0]=0.0;
    m[4][0]=0.0;
    m[5][0]=0.0;
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown stress state is required in function tensor_vector_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function assembles %matrix corresponding to the fourth order tnesor
   
   @param m - %matrix form
   @param t - tensor form
   @param ssst - strain/stress state
   
   18.3.2005, JK
*/
void tensor4_matrix (matrix &m,matrix &t,strastrestate ssst)
{
  fillm (0.0,m);
  
  switch (ssst){
  case bar:{
    m[0][0]=t[0][0];
    break;
  }
  case plbeam:{
    break;
  }
  case planestress:{
    m[0][0]=t[0][0];  m[0][1]=t[0][1];  m[0][2]=t[0][5];  m[0][3]=0.0;
    m[1][0]=t[1][0];  m[1][1]=t[1][1];  m[1][2]=t[1][5];  m[1][3]=0.0;
    m[2][0]=t[5][0];  m[2][1]=t[5][1];  m[2][2]=t[5][5];  m[2][3]=0.0;
    m[3][0]=0.0;      m[3][1]=0.0;      m[3][2]=0.0;      m[3][3]=0.0;
    break;
  }
  case planestrain:{
    m[0][0]=t[0][0];  m[0][1]=t[0][1];  m[0][2]=t[0][2];  m[0][3]=t[0][5];
    m[1][0]=t[1][0];  m[1][1]=t[1][1];  m[1][2]=t[1][2];  m[1][3]=t[1][5];
    m[2][0]=t[2][0];  m[2][1]=t[2][1];  m[2][2]=t[2][2];  m[2][3]=t[2][5];
    m[3][0]=t[5][0];  m[3][1]=t[5][1];  m[3][2]=t[5][2];  m[3][3]=t[5][5];
    break;
  }
  case axisymm:{
    m[0][0]=t[0][0];  m[0][1]=t[0][1];  m[0][2]=t[0][2];  m[0][3]=t[0][5];
    m[1][0]=t[1][0];  m[1][1]=t[1][1];  m[1][2]=t[1][2];  m[1][3]=t[1][5];
    m[2][0]=t[2][0];  m[2][1]=t[2][1];  m[2][2]=t[2][2];  m[2][3]=t[2][5];
    m[3][0]=t[5][0];  m[3][1]=t[5][1];  m[3][2]=t[5][2];  m[3][3]=t[5][5];
    break;
  }
  case spacestress:{
    copym(t,m);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown strain/stress state is required in function tensor4_vector (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

/**
   function assembles %matrix corresponding to the fourth order tensor
   
   m - matrix form
   t - tensor form
   ssst - strain/stress state
   
   18.3.2005, JK
*/
void gentensor_matrix (matrix &m,matrix &t,strastrestate ssst)
{
  long i,n;
  fillm (0.0,m);
  
  n=t.n;

  switch (ssst){
  case bar:{
    for (i=0;i<n;i++){
      m[0][i]=t[0][i];
    }
    break;
  }
  case plbeam:{
    break;
  }
  case planestress:{
    for (i=0;i<n;i++){
      m[0][i]=t[0][i];
      m[1][i]=t[1][i];
      m[2][i]=t[5][i];
    }
    break;
  }
  case planestrain:{
    for (i=0;i<n;i++){
      m[0][i]=t[0][i];
      m[1][i]=t[1][i];
      m[2][i]=t[2][i];
      m[3][i]=t[5][i];
    }
    break;
  }
  case axisymm:{
    for (i=0;i<n;i++){
      m[0][i]=t[0][i];
      m[1][i]=t[1][i];
      m[2][i]=t[2][i];
      m[3][i]=t[5][i];
    }
    break;
  }
  case spacestress:{
    copym(t,m);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown strain/stress state is required in function tensor_vector (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

/**
   function converts long vectors (with 6 components) to short vectors (with appropriate number of components)
   e.g. in the case of plane stress, vectors with 6 components are converted to vectors with 3 components
   
   @param lv - long %vector
   @param sv - short %vector
   @param ssst - strain/stress state
   
   JK, 17.2.2007
*/
/*
void longvect_shortvect (double *lv,double *sv,strastrestate ssst)
{
  fillm (0.0,sv);
  
  switch (ssst){
  case bar:{
    sv[0]=lv[0];
    break;
  }
  case plbeam:{
    break;
  }
  case planestress:{
    sv[0]=lv[0];
    sv[1]=lv[1];
    sv[2]=lv[5];
    break;
  }
  case planestrain:{
    sv[0]=lv[0];
    sv[1]=lv[1];
    sv[2]=lv[5];
    break;
  }
  case axisymm:{
    sv[0]=lv[0];
    sv[1]=lv[1];
    sv[2]=lv[2];
    sv[3]=lv[5]
    break;
  }
  case spacestress:{
    sv[0]=lv[0];
    sv[1]=lv[1];
    sv[2]=lv[2];
    sv[3]=lv[3];
    sv[4]=lv[4];
    sv[5]=lv[5];
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown strain state is required in function longvect_shortvect (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

*/