scaldamtime2.cpp

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

  }
  // average size of element
  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__);
  }
  ft = Mm->give_actual_ft(ipp);
  uf_a = uf/(ft/st);
  tmp = -omegao*kappao[0]*h/uf_a;
  dsigma = -ft*h/eo*exp(tmp)*(dkappa[0]*eo*eo+(e-eo)*(1.0-omegao)*eo*kappao[0]);
  dsigma /= uf_a*eo-(ft*h)*exp(tmp);
//  dsigma = -ft*h/uf_a*dkappa[0]*exp(tmp);
//  dsigma /= 1 - (ft/eo)*(h/uf_a)*exp(tmp);

  tmp = (1.0-omegao)*(eo*dkappa[0]+(e-eo)*kappao[0])-dsigma;
  domega = tmp/(eo*kappao[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+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 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, idem;
  double epse, omegao, omega, domega, e, nu;
  vector epsn(ncompstr),epso(ncompstr), deps(ncompstr), depseq(ncompstr);
  vector sigma(ncompstr), kappa(1), dkappa(1), kappao(1);
  vector deps_dam(ncompstr), auxv(ncompstr), dsigma(ncompstr), sigmao(ncompstr);
  vector dsc(ncompstr);
  matrix d(ncompstr, ncompstr);

  
  e = Mm->give_actual_ym(ipp);
  if (Mm->ip[ipp].ssst == planestress)
  {
    idem = Mm->ip[ipp].gemid();
    nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
    Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
  }
  if (Mp->phase == 1)
  {
/*    for (i=0;i<ncompstr;i++){
      // strain increment
      deps[i] = Mm->ip[ipp].eqother[ido+ncompstr+i];
      // previous stresses increment
      dsigma[i] = Mm->ip[ipp].eqother[ido+3*ncompstr+3+i];
    }
    
    // previous Young modulus
    eo = Mm->ip[ipp].eqother[ido+2*ncompstr+2];
    
    // computation of dsigma-*E_i*deps=dsigma_dam
    elmatstiff (d,ipp);
    cmulm(eo/e, d, d); // elastic siffness matrix for old Young modulus
    mxv(d, deps, auxv); 
    subv(dsigma, auxv, dsigma); 

    for (i=0; i < ncompstr; i++)
      Mm->ip[ipp].stress[i] = dsigma[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);
      // new total stress increment storage
      for (i=0; i < ncompstr; i++)
        Mm->ip[ipp].stress[i] += sigma[i];
    }*/
    for (i=0; i < ncompstr; i++)
      Mm->ip[ipp].stress[i] = 0.0;
  }   
  if (Mp->phase==2)
  {
    for (i=0;i<ncompstr;i++){
      //  total new strains
      epsn[i] = Mm->ip[ipp].strain[i];
    }
    
    // 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;
      domega = 0.0;
    }
    else
    {
      omega = damfunction(ipp, kappa);
      if ((omega <= omegao))
      {
        omega = omegao;
        domega = 0.0;    
      }
      else
        domega = omega - omegao;
    }
    
    elmatstiff (d,ipp);

    /*//old
      cmulm(eo/e, d, d); // elastic siffness matrix for old Young modulus
      
      // computation of -domega*E_i*eps_i component of dsigma
      fillv(0.0, dsigma);
      cmulv(domega, epso, auxv);
      mxv(d, auxv, dsigma); 
      cmulv(-1.0, dsigma); 
      
      // computation of (1-omega_i)*dE*eps_i component of dsigma
      cmulv(1.0-omegao, epso, auxv);
      mxv(d, auxv, dsc);
      cmulv(e-eo, dsc);
      addv(dsigma, dsc, dsigma);
      
      // computation of (1-omega_i)*E_i*deps component of dsigma
      cmulv(1.0-omegao, deps, auxv);
      mxv(d, auxv, dsc);
      addv(dsigma, dsc, dsigma);
    */
    
    //new 7.11.2006
    // computation of (1-omega_i)*E_i*eps component of dsigma
    fillv(0.0, sigma);
    cmulv(1.0-omega, epsn, auxv);
    mxv(d, auxv, sigma);
    
    // new total stresses
    for (i=0; i < ncompstr; i++)
      Mm->ip[ipp].stress[i] = sigma[i];

    // previous total omega storage
    Mm->ip[ipp].eqother[ido+2*ncompstr+1] = 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+1];
}



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

void scaldamtime::initvalues (long ipp, long im, long ido)
{
  //long ncompstr = Mm->ip[ipp].ncompstr;
  // initial value of Young modulus
  //Mm->ip[ipp].eqother[ido+2*ncompstr+2] = Mm->give_actual_ym(ipp);
}