scaldamtime.cpp

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

  @param kappa - %vector of the parameters of damage function
  
  TKo
*/
double scaldamtime::der_damfunction(long ipp, vector &kappa)
{
  double domega = 0.0, e, tmp,dtmp,ft, uf_a;
  long idem = Mm->ip[ipp].gemid();
  
  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 scaldamtime::der_damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }
  // Young modulus
  e = Mm->give_actual_ym(ipp);
  ft = Mm->give_actual_ft(ipp);
  uf_a = uf/(ft/st);
  tmp = exp(-(kappa[0]-ft/e)/(uf_a-ft/e)); 
  dtmp = -1.0/(uf_a-ft/e)*tmp;
  domega = -ft/e*(kappa[0]*dtmp-tmp)/(kappa[0]*kappa[0]);
  return domega;
}

/**
  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 scaldamtime::matstiff (matrix &d,long ipp,long ido)
{
  double dp;
  long ncompstr = Mm->ip[ipp].ncompstr;

  switch (Mp->stmat)
  {
    case initial_stiff :
      elmatstiff(d, ipp);
      break;
    case tangent_stiff :
      elmatstiff(d, ipp);
      dp=Mm->ip[ipp].eqother[ido+2*ncompstr+2];
      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 scaldamtime::matstiff (file %s, line %d)\n", __FILE__, __LINE__);
  }
}


/**
  This function computes elastic material stiffnes %matrix from actual Young modulus.

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

  TKo
*/
void scaldamtime::elmatstiff (matrix &d,long ipp)
{
  double e, e0;
  long idem = Mm->ip[ipp].gemid();

  Mm->elmatstiff (d,ipp);
  e = Mm->give_actual_ym(ipp);
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\n Invalid type of elastic material is required");
    fprintf(stderr, "\n  in function scaldamtime::elmatstiff, (file %s, line %d)\n", __FILE__, __LINE__);
  }
  e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  cmulm(e/e0, d);
}


/**
  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 scaldamtime::nlstresses (long ipp, long im, long ido)
{
  long i,ncompstr=Mm->ip[ipp].ncompstr, ncompo;
  double epse, kappao, omegao, omega, domega, tmp;
  vector epsn(ncompstr),epso(ncompstr), deps(ncompstr), depseq(ncompstr), sigma(ncompstr), kappa(1);
  vector deps_dam(ncompstr), auxv(ncompstr);
  matrix d(ncompstr, ncompstr);
  long nm;

  if (Mp->phase == 1)
  {
    for (i=0;i<ncompstr;i++){
      //  total new strains
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strains
      epso[i] = Mm->ip[ipp].eqother[ido+i];
      // strain increment
      deps[i] = Mm->ip[ipp].eqother[ido+ncompstr+i];
    }
    
    // previous equivalent strain
    kappao = Mm->ip[ipp].eqother[ido+2*ncompstr+0];
    // previous damage
    omegao = Mm->ip[ipp].eqother[ido+2*ncompstr+1];
    // new equivalent strain
    damfuncpar(ipp, epsn, kappa);
    epse = epsefunction (ipp);
    if (kappa[0] < epse)
      omega = omegao;
    else
      // new damage 
      omega = damfunction(ipp,kappa);
    if ((omega <= omegao))
    {
      omega = omegao;
      domega = 0.0;    
    }
    else
    {
      der_damfuncpar(ipp, epsn, kappa, depseq);
      domega = der_damfunction(ipp, kappa);
      scprd(depseq, deps, tmp);
      domega *= tmp;
    }
    cmulv(domega, epsn, deps_dam);
    cmulv(omega, deps, auxv);
    addv(auxv, deps_dam, deps_dam);
    elmatstiff (d,ipp);
    mxv(d, deps_dam, sigma);

    // stress increment storage
    for (i=0;i<ncompstr;i++)
      Mm->ip[ipp].stress[i]=sigma[i];
    // new equivalent strain storage 
    Mm->ip[ipp].eqother[ido+2*ncompstr+0] = kappa[0];
    // previous damage parameter storage
    Mm->ip[ipp].eqother[ido+2*ncompstr+1] = omega;

    // 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 < ncompstr; i++)
        Mm->ip[ipp].stress[i] += sigma[i];
    }
  }   
  if (Mp->phase==2)
  {
    for (i=0;i<ncompstr;i++)
    {
      //  total new strains
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strains
      epso[i] = Mm->ip[ipp].eqother[ido+i];
      // total strain increment
      deps[i] = epsn[i] - epso[i];
    }
    // previous equivalent strain
    kappao = Mm->ip[ipp].eqother[ido+2*ncompstr+0];
    // previous damage
    omegao = Mm->ip[ipp].eqother[ido+2*ncompstr+1];
    // new equivalent strain
    damfuncpar(ipp, epsn, kappa);
    epse = epsefunction (ipp);
    if (kappa[0] < epse)
      omega = omegao;
    else
      // new damage 
      omega = damfunction(ipp,kappa);
    if ((omega <= omegao))
      omega = omegao;
    elmatstiff (d,ipp);
    cmulm(1.0-omega, d);
    mxv(d, epsn, sigma);
    for (i=0;i<ncompstr;i++)
      Mm->ip[ipp].stress[i]=sigma[i];

    for (i=0;i<ncompstr;i++)
    {
      // total strain storage
      Mm->ip[ipp].eqother[ido+i] = epsn[i];
      // previous increment of total strain
      Mm->ip[ipp].eqother[ido+ncompstr+i] = deps[i];
    }
    // actual omega storage
    Mm->ip[ipp].eqother[ido+2*ncompstr+2] = omega;

    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);
    }
  }
}
*/



/**
  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 scaldamtime::nlstresses (long ipp, long im, long ido)
{
  long i,ncompstr=Mm->ip[ipp].ncompstr, ncompo, nm;
  double epse, kappao, omegao, omega, domega;
  vector epsn(ncompstr),epso(ncompstr), deps(ncompstr), depseq(ncompstr), sigma(ncompstr), kappa(1);
  vector deps_dam(ncompstr), auxv(ncompstr);
  matrix d(ncompstr, ncompstr);

  if (Mp->phase == 1)
  {
    for (i=0;i<ncompstr;i++){
      //  total new strains
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strains
      epso[i] = Mm->ip[ipp].eqother[ido+i];
      // strain increment
      deps[i] = Mm->ip[ipp].eqother[ido+ncompstr+i];
    }
    
    // previous equivalent strain
    kappao = Mm->ip[ipp].eqother[ido+2*ncompstr+0];
    // previous damage
    omegao = Mm->ip[ipp].eqother[ido+2*ncompstr+1];
    // new equivalent strain
    damfuncpar(ipp, epsn, kappa);
    epse = epsefunction (ipp);
    if (kappa[0] < epse)
      omega = omegao;
    else
      // new damage 
      omega = damfunction(ipp,kappa);
    if ((omega <= omegao))
    {
      omega = omegao;
      domega = 0.0;    
    }
    else
      domega = omega - omegao;
    cmulv(domega, epsn, deps_dam);
    cmulv(omega, deps, auxv);
    addv(auxv, deps_dam, deps_dam);
    elmatstiff (d,ipp);
    mxv(d, deps_dam, sigma);
    // stress increment storage
    for (i=0;i<ncompstr;i++)
      Mm->ip[ipp].stress[i]=sigma[i];
    // new equivalent strain storage 
    Mm->ip[ipp].eqother[ido+2*ncompstr+0] = kappa[0];
    // previous damage parameter storage
    Mm->ip[ipp].eqother[ido+2*ncompstr+1] = omega;
    // 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 < ncompstr; i++)
        Mm->ip[ipp].stress[i] += sigma[i];
    }
    
  }   
  if (Mp->phase==2)
  {
    for (i=0;i<ncompstr;i++)
    {
      //  total new strains
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strains
      epso[i] = Mm->ip[ipp].eqother[ido+i];
      // total strain increment
      deps[i] = epsn[i] - epso[i];
    }
    // previous equivalent strain
    kappao = Mm->ip[ipp].eqother[ido+2*ncompstr+0];
    // previous damage
    omegao = Mm->ip[ipp].eqother[ido+2*ncompstr+1];
    // new equivalent strain
    damfuncpar(ipp, epsn, kappa);
    epse = epsefunction (ipp);
    if (kappa[0] < epse)
      omega = omegao;
    else
      // new damage 
      omega = damfunction(ipp,kappa);
    if ((omega <= omegao))
      omega = omegao;
    elmatstiff (d,ipp);
    cmulm(1.0-omega, d);
    mxv(d, epsn, sigma);
    for (i=0;i<ncompstr;i++)
      Mm->ip[ipp].stress[i]=sigma[i];

    for (i=0;i<ncompstr;i++)
    {
      // total strain storage
      Mm->ip[ipp].eqother[ido+i] = epsn[i];
      // previous increment of total strain
      Mm->ip[ipp].eqother[ido+ncompstr+i] = deps[i];
    }
    // actual omega storage
    Mm->ip[ipp].eqother[ido+2*ncompstr+2] = omega;

    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);
    }
  }
}


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

  @param ipp - integration point number in the mechmat ip array.
  
  TKo
*/
double scaldamtime::epsefunction (long ipp)
{
  double ft;

  double e = Mm->give_actual_ym(ipp);
  ft = Mm->give_actual_ft(ipp);

  return (ft /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 scaldamtime::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__);
    }
    }
  }
}



/**
  This function returns the value of damage parameter

  @param ipp - integration point number in the mechmat ip array.
  @param ido - index of internal variables for given material in the ipp other array
  
  TKo
*/
double scaldamtime::givedamage (long ipp, long ido)
{
  long ncompstr = Mm->ip[ipp].ncompstr;
  return Mm->ip[ipp].eqother[ido+2*ncompstr+2];
}



/**
  This function returns the value of tensile strength

  @param ipp - integration point number in the mechmat ip array.
  @param im - material index
  @param ido - index of internal variables for given material in the ipp other array
  
  TKo
*/
double scaldamtime::give_actual_ft (long ipp, long im, long ido)
{
  return st;
}