therisomattime.cpp

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

#include "therisomattime.h"
#include "global.h"
#include "vector.h"
#include "matrix.h"
#include "intpoints.h"

therisomattime::therisomattime (void)
{
  alpha=0.0;
}
therisomattime::~therisomattime (void)
{

}

void therisomattime::read (XFILE *in)
{
  xfscanf (in,"%lf",&alpha);
}

void therisomattime::matdilat (matrix &d,strastrestate ssst)
{
  switch (ssst){
  case bar:{
    matdilat_bar (d);
    break;
  }
  case plbeam:{
    matdilat_plbeam (d);
    break;
  }
  case planestress:{
    matdilat_plstress (d);
    break;
  }
  case planestrain:{
    matdilat_plstrain (d);
    break;
  }
  case plate:{
    matdilat_plate (d);
    break;
  }
  case axisymm:{
    matdilat_axi (d);
    break;
  }
  case spacestress:{
    matdilat_spacestr (d);
    break;
  }
  default:{
    fprintf (stderr,"\n unknown number of components of stress tensor is required");
    fprintf (stderr,"\n in function therisomattime::matdilat (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   function creates thermal dilatancy matrix of the elastic
   isotropic material for bar elements
   
   d - thermal dilatancy matrix of the material
   
   7.6.2005 TKo
*/
void therisomattime::matdilat_bar (matrix &d)
{
  d[0][0] = alpha;
}

/**
   function creates thermal dilatancy matrix of the elastic
   isotropic material for plane beam elements
   
   @param d - thermal dilatancy matrix of the material
   
   7.6.2005 TKo
*/
void therisomattime::matdilat_plbeam (matrix &d)
{
  d[0][0] = alpha;
  d[1][1] = 0.0;
  d[2][2] = 0.0;//tady vlozeno
  //d[2][2] = alpha;
}

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

   @param d - thermal dilatancy matrix of the material

   7.6.2005 TKo
*/
void therisomattime::matdilat_plstress (matrix &d)
{
  fillm(0.0,d);
  
  d[0][0] = alpha;  d[0][1] = 0.0;    d[0][2] = 0.0;
  d[1][0] = 0.0;    d[1][1] = alpha;  d[1][2] = 0.0;
  d[2][0] = 0.0;    d[2][1] = 0.0;    d[2][2] = 0.0;
}

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

   @param d - thermal dilatancy matrix of the material

   7.6.2005 TKo
*/
void therisomattime::matdilat_plstrain (matrix &d)
{
  fillm(0.0,d);
  
  d[0][0] = alpha;  d[0][1] = 0.0;     d[0][2] = 0.0;
  d[1][0] = 0.0;    d[1][1] = alpha;   d[1][2] = 0.0;
  d[2][0] = 0.0;    d[2][1] = 0.0;     d[2][2] = 0.0;
}

void therisomattime::matdilat_axi (matrix &d)
{
  fillm(0.0,d);
  
  d[0][0]=alpha;  d[0][1]=0.0;    d[0][2]=0.0;
  d[1][0]=0.0;    d[1][1]=alpha;  d[1][2]=0.0;
  d[2][0]=0.0;    d[2][1]=0.0;    d[2][2]=alpha;

  d[3][3]=0.0;
}

/**
   function creates thermal dilatancy matrix of the elastic
   isotropic material for plate elements
   
   @param d - thermal dilatancy matrix
   
   7.6.2005 TKo
*/
void therisomattime::matdilat_plate (matrix &d)
{
  fillm(0.0,d);
  
  d[0][0]=alpha;  d[1][1]=alpha;
}

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

   7.6.2005, TKo
*/
void therisomattime::matdilat_spacestr (matrix &d)
{
  fillm(0.0,d);
  
  d[0][0]=alpha;  d[1][1]=alpha;  d[2][2]=alpha;
}

/**
   Function initializes eqother array with initial temperature.
   Actual value of array Mm->tempr[ipp] is taken as initial temperature.
   
   @param ipp - integration point pointer
   @param ido - index of internal variables for given material in the ipp other array

   7.6.2005, TKo
*/
void therisomattime::initvalues (long ipp, long ido)
{
  Mm->ip[ipp].eqother[ido] = Mm->tempr[ipp];
  Mm->ip[ipp].eqother[ido+1] = Mm->tempr[ipp];
}

/**
   function computes eigenstrains caused by temperature changes
   function assumes temperature changes in the array Mm->tempr

   @param ipp - integration point id
   @param ido - index of internal variables in the ip's ipp other array
   
   7.6.2005, TKo
*/
void therisomattime::temprstrains (long ipp, long ido)
{
  double dt;
  strastrestate sss;
  
  //  type of strain/stress state
  sss = Mm->ip[ipp].ssst;
  
  //  change of temperature
  dt = Mm->tempr[ipp]-Mm->ip[ipp].eqother[ido];

  switch (sss){
  case bar:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    break;
  }
  case plbeam:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    Mm->tempstrains[ipp][1]=0.0;
    Mm->tempstrains[ipp][2]=0.0;//tady oprava
    //Mm->tempstrains[ipp][2]=alpha*dt;
    break;
  }
  case planestress:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    Mm->tempstrains[ipp][1]=alpha*dt;
    Mm->tempstrains[ipp][2]=0.0;
    Mm->tempstrains[ipp][3]=alpha*dt;
    break;
  }
  case planestrain:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    Mm->tempstrains[ipp][1]=alpha*dt;
    Mm->tempstrains[ipp][2]=0.0;
    Mm->tempstrains[ipp][3]=alpha*dt;
    break;
  }
  case axisymm:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    Mm->tempstrains[ipp][1]=alpha*dt;
    Mm->tempstrains[ipp][2]=alpha*dt;
    break;
  }
  case spacestress:{
    Mm->tempstrains[ipp][0]=alpha*dt;
    Mm->tempstrains[ipp][1]=alpha*dt;
    Mm->tempstrains[ipp][2]=alpha*dt;
    break;
  }
  default:{
    fprintf (stderr,"\n unknown strain/stress state is required in function therisomattime::temprstrains (file %s, line %d)\n",
	     __FILE__,__LINE__);
  }
  }
  
}



/**
   function computes increment of stresses caused by increment of temperature changes
   function assumes temperature changes in the array Mm->tempr

   @param ipp - integration point id
   @param ido - index of internal variables in the ip's ipp other array
   
   7.6.2005, TKo
*/
void therisomattime::nlstresses (long ipp, long ido)
{
  double dt;
  strastrestate sss;
  long ncompstr = Mm->ip[ipp].ncompstr;
  vector dstrain(ncompstr), dsig(ncompstr);
  matrix d(ncompstr, ncompstr);
  
  //  type of strain/stress state
  sss = Mm->ip[ipp].ssst;
  if (Mp->phase == 1)
  {
    //  change of temperature
    dt = Mm->tempr[ipp]-Mm->ip[ipp].eqother[ido+1];

    switch (sss)
    {
    case bar:
      dstrain[0]=alpha*dt;
      break;
    case plbeam:{
    dstrain[0]=alpha*dt;
    dstrain[1]=0.0;
    dstrain[2]=0.0;//tady oprava
    //dstrain[2]=alpha*dt;
    break;
    }
    case planestress:
      dstrain[0]=alpha*dt;
      dstrain[1]=alpha*dt;
      dstrain[2]=0.0;
      dstrain[3]=alpha*dt;
      break;
    case planestrain:
      dstrain[0]=alpha*dt;
      dstrain[1]=alpha*dt;
      dstrain[2]=0.0;
      dstrain[3]=alpha*dt;
      break;
    case axisymm:
      dstrain[0]=alpha*dt;
      dstrain[1]=alpha*dt;
      dstrain[2]=alpha*dt;
      break;
    case spacestress:
      dstrain[0]=alpha*dt;
      dstrain[1]=alpha*dt;
      dstrain[2]=alpha*dt;
      break;
    default:
      fprintf (stderr,"\n unknown strain/stress state is required in function therisomattime::nlstresses (file %s, line %d)\n",
	       __FILE__,__LINE__);
    }
    Mm->matstiff(d, ipp);
    mxv(d, dstrain, dsig);
    Mm->storestress(0, ipp, dsig);
  }
  if (Mp->phase == 2)
    Mm->ip[ipp].eqother[ido+1] = Mm->tempr[ipp];
  
}