scaldam.cpp

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

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

#define nijac 20
#define limit 1.0e-8


/**
  This constructor inializes attributes to zero values.
  
  TKo
*/
scaldam::scaldam (void)
{
  indam=0.0;  st=0.0;  uf = 0.0; c = 0.0; k = 0.0; cde = corr_on;
}

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

}

/**
  This function reads material parameters from the opened text file given
  by the parameter in.

  @param in - pointer to the opened input text file
  
  TKo
*/
void scaldam::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf %d %d",&indam,&st, &uf, (int *)(&ft), (int *)(&cde));
  if (ft == norenergy)
    xfscanf(in, "%lf", &c);
  if (ft == vonmises)
    xfscanf(in, "%lf", &k);
  
  sra.read (in);
}

/**
  This function computes parameters for the damage function
  Different types of the parameter computing are given by the
  attribute ft.

  @param  ipp - integration point number
  @param  eps - %vector of the strains
  @param  kappa - %vector where the computed parameters are stored
  
  TKo
*/
void scaldam::damfuncpar(long ipp, vector &eps, vector &kappa)
{
  double epseq;
  vector poseps;
  vector princsig;
  vector princeps;
  vector posprincsig;
  vector tmp;
  vector sig;
  matrix pvect;
  matrix sigt;
  matrix epst;
  matrix dev;
  matrix d;
  double i1e, j2e, a, b, e, nu, temp; 
  long ncomp, idem = Mm->ip[ipp].gemid();
  long im;

  switch (ft)
  {
    case norstrain :
      // strain norm
      scprd(eps, eps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norenergy :
      // energy norm
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(ncomp, ncomp, d);
      Mm->elmatstiff (d,ipp);
      vxm(eps, d, tmp);
      scprd(tmp, eps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norposstrain :
      // positive strain norm
      allocv(eps.n,poseps);
      extractposv (eps, poseps);
      scprd(poseps, poseps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norposenergy :
      // positive energy norm
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(ncomp, ncomp, d);
      allocv(ncomp, tmp);
      Mm->elmatstiff (d,ipp);
      allocv(eps.n,poseps);
      extractposv (eps, poseps);
      vxm(poseps, d, tmp);
      scprd(tmp, poseps, epseq);
      kappa[0] = sqrt(epseq)*c;
      break;
    case norrankine :
    // Rankine norm
    {
      ncomp=Mm->ip[ipp].ncompstr;
      allocv(ncomp, sig);
      allocm(ncomp, ncomp, d);
      allocm(3, 3, sigt);
      allocm(3, 3, pvect);
      allocv(3, princsig);
      Mm->elmatstiff (d,ipp);
      mxv(d, eps, sig);
      vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
      princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3);
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function scaldam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      // principal values are sorted -> maximal value of principal stresses is at the last
      // position of vector princsig
      if (princsig[2] > 0.0)
        kappa[0] = princsig[2] / Mm->eliso[Mm->ip[ipp].idm[idem]].e;
      else
        kappa[0] = 0.0;
      break;
    }
    case norrankinesmooth :
    // smoothed Rankine norm
    {
      ncomp=Mm->ip[ipp].ncompstr;
      allocv(ncomp, sig);
      allocm(ncomp, ncomp, d);
      allocm(3, 3, sigt);
      allocm(3, 3, pvect);
      allocv(3, princsig);
      allocv(3, posprincsig);
      Mm->elmatstiff (d,ipp);
      mxv(d, eps, sig);
      vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
      princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3);
      extractposv (princsig, posprincsig);
      scprd(posprincsig, posprincsig, epseq);
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function scaldam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      kappa[0] = sqrt(epseq) / Mm->eliso[Mm->ip[ipp].idm[idem]].e;
      break;
    }
    case normazar :
    // Mazar's norm
    {
      allocm(3, 3, epst);
      allocm(3, 3, pvect);
      allocv(3, princeps);
      allocv(3,poseps);
      fillv(0.0, princeps);
      vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
      princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3);
      extractposv (princeps, poseps);
      scprd(poseps, poseps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    }
    case vonmises :
      // glasgow modified von Mises norm
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(3, 3, epst);
      allocm(3, 3, dev);
      vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
      i1e = first_invar(epst);
      deviator (epst,dev);
      j2e = second_invar (dev);
      idem = Mm->ip[ipp].gemid();
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function glasgowdam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      im = Mm->ip[ipp].idm[idem];
      e  = Mm->eliso[im].e;
      nu = Mm->eliso[im].nu;
      a   = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
      b   = 3.0/k/(1+nu)/(1.0+nu);
      temp = a*a*i1e*i1e + b*j2e;
      kappa[0] = a*i1e+sqrt(temp);
/*
      // original von Mises norm
      kappa[0]  = (k-1.0)/(2.0*k)/(1.0-2.0*nu)*i1e;
      tmp = (k-1.0)*(k-1.0)/(1.0-2.0*nu)/(1.0-2.0*nu)*i1e*i1e -12.0*k/(1.0+nu)/(1.0+nu)*j2e;
      kappa[0] += 1.0/(2.0*k)*sqrt(tmp);*/
      break;
    default :
    {
      fprintf(stderr, "\n\n Unknown damage parameter function type");
      fprintf(stderr, "\n  in function scaldam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
    }
  }
  return;
}

/**
  This function computes damage parameter omega which is the result of the
  damage function.

  @param ipp - integration point number
  @param kappa - %vector of the parameters of damage function
  
  TKo
*/
double scaldam::damfunction(long ipp, vector &kappa)
{
  double err,omega = 0.0, dtmp, tmp, h, e;
  long i,ni, eid, idem = Mm->ip[ipp].gemid();
  
  ni=sra.give_ni ();
  err=sra.give_err ();
  
  if (kappa[0] < indam)
  // elastic loading
    return 0.0;
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
    fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
    fprintf(stderr, "\n Function scaldam::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }
  // Young modulus
  e = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  if ((Mm->ip[ipp].hmt & 2) > 0)
  {
    omega = 1.0-(st/e/kappa[0])*exp(-(kappa[0]-st/e)/(uf-st/e));
    return omega;
  }
/*
  if (cde == corr_off)
  // no correction of disipated energy only simple scalar damge
  { 
    omega = 1.0-(st/e/kappa[0])*exp(-(kappa[0]-st/e)/(uf-st/e));
    return omega;
  }
*/
  // correction of dispiated energy - mesh independent scalar damage
  eid = Mm->elip[ipp];
  switch (Mt->give_dimension(eid))
  {
    case 1 :
      h = Mt->give_length(eid);
      break;
    case 2 :
      h = sqrt(Mt->give_area(eid));
      break;
    case 3 :
      h = pow(Mt->give_volume(eid), 1.0/3.0);
      break;
    default :
      fprintf(stderr, "\n\nError determining dimension of element");
      fprintf(stderr, "\n function damfunction, (file %s, line %d)\n", __FILE__, __LINE__);
  }
  if (cde == corr_off)
  // no correction of disipated energy only simple scalar damge
  { 
    omega = 1.0-(st/e/kappa[0])*exp(-h*(kappa[0]-st/e)/(uf-st/e));
    return omega;
  }
  for (i = 0; i < ni; i++)
  {
    dtmp = -e*kappa[0]+st/uf*h*kappa[0]*exp(-h*omega*kappa[0]/uf);
    tmp = (1.0-omega)*e*kappa[0]-st*exp(-h*omega*kappa[0]/uf);
    if (fabs(tmp) < err)
    {
      omega += - tmp/dtmp;
      break;
    }
    omega += - tmp/dtmp;
  }
  if (fabs(tmp) > err)
  {
    fprintf(stderr, "\n\nError - iteration of omega doesn't converge with required error");
    fprintf(stderr, "\n in function damfunction(), (file %s, line %d)\n", __FILE__, __LINE__);
  }
  if ((omega > 1.0) || (omega < 0.0))
  {
    fprintf(stderr, "\n\nError - iteration of omega doesn't converge to range <0,1>");
    fprintf(stderr, "\n in function damfunction(), (file %s, line %d)\n", __FILE__, __LINE__);
  }

  return omega;
}



/**
  This function computes material stiffnes %matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number
  @param ido - index of internal variables for given material in the ipp other array
  
  TKo
*/
void scaldam::matstiff (matrix &d,long ipp,long ido)
{
  double dp;

  switch (Mp->stmat)
  {
    case initial_stiff :
      Mm->elmatstiff (d,ipp);
      break;
    case tangent_stiff :
      Mm->elmatstiff (d,ipp);
      dp=Mm->ip[ipp].eqother[ido+1];
      if (dp > 0.999999)
        dp = 0.999999;
      cmulm (1.0-dp,d);
      break;
    default :
      fprintf(stderr, "\n\nError - unknown type of stifness matrix");
      fprintf(stderr, "\n in function scaldam::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 type for given ip
  @param ido - index of internal variables for given material in the ipp other array

  TKo, 7.10.2001
*/
void scaldam::nlstresses (long ipp, long im, long ido)
{
  long i,ncomp=Mm->ip[ipp].ncompstr;
  double dp;
  vector epsn(ncomp),sigma(ncomp), kappa(1);

  //  initial values
  for (i=0;i<ncomp;i++){
    epsn[i] = Mm->ip[ipp].strain[i];
  }
  kappa[0] = Mm->ip[ipp].eqother[ido+0];

  // damage stress solver
  dp = Mm->scal_dam_sol (ipp, im, ido, epsn, kappa, sigma);

  //  new data storage
  for (i=0;i<ncomp;i++){
    Mm->ip[ipp].stress[i]=sigma[i];
  }
  Mm->ip[ipp].other[ido+0]=kappa[0];
  Mm->ip[ipp].other[ido+1]=dp;
}

/**
  This function updates values in the other array reached in the previous equlibrium state to
  values reached in the new actual equilibrium state.

  @param ipp - integration point number in the mechmat ip array.
  @param im  - index of material type for given ip
  @param ido - index of internal variables for given material in the ipp other array
  
  TKo
*/
void scaldam::updateval (long ipp,long im,long ido)
{
  long i,n = Mm->givencompeqother(ipp,im);

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

}

/**
  This function returns the value of the limit elastic strain .

  @param ipp - integration point number in the mechmat ip array.
  
  TKo
*/
double scaldam::epsefunction (long ipp)
{
  long idem = Mm->ip[ipp].gemid();
  double e = Mm->eliso[Mm->ip[ipp].idm[idem]].e;

  return (st /e) ;
}

/**
   function changes material parameters for stochastic analysis
   
   @param atm - selected material parameters (parameters which are changed)
   @param val - array containing new values of parameters
   
   TKo
*/
void scaldam::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      indam=val[i];
      break;
    }
    case 1:{
      st=val[i];
      break;
    }
    case 2:{
      uf=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
}