visplast.cpp

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

#include "visplast.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.
*/
visplast::visplast (void)
{

}

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

}

/**
  This function computes material stiffnes %matrix.

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

*/
void visplast::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 visplast::matstiff (%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 visplast::nlstresses (long ipp,long im,long ido)
{
  long i,ncomp,ncompo,nhp,nm;
  double dt;
  
  //  number of strain/stress components in the problem
  ncomp = Mm->ip[ipp].ncompstr;
  //  number of components of array eqother belonging to viscous material
  ncompo = Mm->givencompother(ipp,im+1);
  
  vector sig(ncomp),epsn(ncomp),epso(ncomp),epsp(ncomp),q;
  matrix d(ncomp,ncomp);
  
  if (Mp->phase==1){
    /*********************************/
    //  right hand side computation  //
    /*********************************/
    
    //  stress components
    for (i=0;i<ncomp;i++){
      sig[i]=Mm->ip[ipp].eqother[ido+i];
    }
    
    //  number of hardening parameters
    nhp = Mm->give_num_interparam (ipp,im+2);
    allocv (nhp,q);
    
    //  hardening parameters
    Mm->give_interparam (ipp,im+2,ido+ncompo,q);
    
    //  time increment
    dt=Mp->timecon.actualforwtimeincr ();
    
    //  computation of trial stress increments
    Mm->stiff_deps_vispl (ipp,im,ido,sig,q,dt);
    
    destrv (q);

    // 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);
      // new total stress increment storage
      for (i=0; i < ncomp; i++)
        Mm->ip[ipp].stress[i] += sig[i];
    }
  }    
  
  if (Mp->phase==2){
    
    for (i=0;i<ncomp;i++){
      //  new total strain
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strain
      epso[i] = Mm->ip[ipp].eqother[ido+ncomp+i];
      //  increment of irreversible strain
      epsp[i] = Mm->ip[ipp].eqother[ido+ncompo+i];
      //epsp[i] = Mm->ip[ipp].eqother[ido+2*ncomp+i];
    }
    
    //  total strain increment
    subv (epsn,epso,epso);
    //  elastic strain increment
    subv (epso,epsp,epsp);
    //  stiffness matrix of material
    Mm->matstiff(d,ipp);
    //  stress increment
    mxv (d,epsp,sig);
    
    //  new data storage
    for (i=0;i<ncomp;i++){
      //  new stress
      Mm->ip[ipp].eqother[ido+i] += sig[i];
      //  total strains
      Mm->ip[ipp].eqother[ido+ncomp+i] = epsn[i];
      //  stress increment changed 27.7.2005 by T.Ko.
      Mm->ip[ipp].stress[i]=Mm->ip[ipp].eqother[ido+i];
    }
    if ((im == 0) && (Mm->ip[ipp].hmt & 1))
    {
      nm=Mm->ip[ipp].nm-1;
      ncompo = Mm->givencompeqother(ipp, 0);
      ncompo -= Mm->givencompeqother(ipp, nm);
      // actualization of previous temperature
      Mm->computenlstresses(ipp, nm, ncompo);
    }

  }
}