elasttime.cpp

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

#include "elasttime.h"
#include "matrix.h"
#include "vector.h"
#include "elastisomat.h"
#include "global.h"
#include "intpoints.h"
#include "vecttens.h"
#include <math.h>

/**
  This constructor initializes attributes to zero values.
*/
elasttime::elasttime (void)
{

}

/**
  This destructor is only for the formal purposes.
*/
elasttime::~elasttime (void)
{

}

/**
  This function computes material stiffnes %matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number

*/
void elasttime::matstiff (matrix &d,long ipp,long im,long ido)
{
  switch (Mp->stmat){
  case initial_stiff:{
    Mm->elmatstiff (d,ipp);
    break;
  }
  default:{
    fprintf(stderr, "\n\nError - unknown type of stifness matrix");
    fprintf(stderr, "\n in function elasttime::matstiff (file %s, line %d)\n",__FILE__,__LINE__);
  }
  }
  
}


/**
  This function computes correct stresses in the integration point and stores
  them into ip stress array.

  @param ipp - integration point pointer
  @param im - index of material
  @param ido - index of the viscous material in the array eqother
  
  25.6.2004
*/
void elasttime::nlstresses (long ipp,long im,long ido)
{
  long i,ncomp;
  //  number of strain/stress components in the problem
  ncomp = Mm->ip[ipp].ncompstr;
  
  vector sig(ncomp),epsn(ncomp);
  matrix d(ncomp,ncomp);
  long nm, ncompo;
  double e, nu;
  long idem;
  
  if (Mp->phase==1)
  {
    /*********************************/
    //  right hand side computation  //
    /*********************************/

    fillv(0.0, sig);
    Mm->storestress (0,ipp,sig);
    
    // computation of increments of stresses due to the temperature strain increment
    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      Mm->computenlstresses(ipp, nm, ncompo);
    }
  }
  if (Mp->phase==2){
    //  new total strain
    for (i=0;i<ncomp;i++)
      epsn[i] = Mm->ip[ipp].strain[i];
  
    idem = Mm->ip[ipp].gemid();
    nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
    e = Mm->give_actual_ym(ipp);
    if (Mm->ip[ipp].ssst == planestress)
      Mm->ip[ipp].strain[3] = epsn[3] = -nu / (1.0 - nu) * (epsn[0]+epsn[1]);

    //  stiffness matrix of material
    Mm->matstiff(d,ipp);
    //  stress increment
    mxv (d,epsn,sig);
    
    //  new data storage
    for (i=0;i<ncomp;i++)
      //  stress increment
      Mm->ip[ipp].stress[i] = sig[i];
 
    // computation of increments of stresses due to the temperature strain increment
    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      Mm->computenlstresses(ipp, nm, ncompo);
    }

  }

}