anisodam.cpp

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

  @param lfc - %vector of actual load function values for compression
  @param pdamc - %vector of principal damage parameters for tension
  @param pdamt - %vector of principal damage parameters for compression
  
  TKo 9.2006
*/
void anisodam::pdam_dev(long ipp, vector &pyc, vector &pyt, vector &dyt, vector &dyc, 
                           double aat, double aac, vector &lft, vector &lfc, vector &pdamt, vector &pdamc)
{
  long i;
  double tmp, ltmp;

  fillv(0.0, pdamt);
  fillv(0.0, pdamc);

  for(i=0; i<3; i++)
  {    
    // tension 
    if ((dyt[i] > 0.0) && (lft[i] > 0.0))
    {
      ltmp = log(pyt(i) - y0t);
      pdamt(i) = aat*exp(ltmp*bt);
      tmp = 1.0/(1.0+pdamt(i));
      pdamt(i) = 1.0 - tmp;
    }
    // compression 
    if ((dyc[i] > 0.0) && (lfc[i] > 0.0))
    {
      ltmp = log(pyc(i) - y0c);
      pdamc(i) = aac*exp(ltmp*bc);
      tmp = 1.0/(1.0+pdamc(i));
      pdamc(i) = 1.0-tmp;
    }
  }
}



/**
  Function computes material parameters A, At and Ac. 
  If the correction of dissipated energy is required, they are computed with respect 
  to size of element for given ipp and variable softening modulus technique was used in 
  order to avoid mesh dependence. Otherwise, values read from the input file are used.

  @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
  @param aa -  actual value of parameter A 
  @param aat - actual value of parameter Ac 
  @param aac - actual value of parameter Ac 

  @retval The function returns computed 
  
  TKo 9.2006
*/
void anisodam::give_actual_param_a(long ipp, long ido, double &aa, double &aac, double &aat)
{
  long eid;
  double e, h, tmp;

  aa  = Mm->ip[ipp].eqother[ido];
  aat = Mm->ip[ipp].eqother[ido+1];
  aac = Mm->ip[ipp].eqother[ido+2];
  if ((aa == 0.0) || (aac == 0.0) || (aat == 0.0))
  {
    if (cde == corr_off)
    // no correction of dissipated energy is required
    {
      aa  = a;
      aat = at;
      aac = ac;
      Mm->ip[ipp].eqother[ido]   = aa;
      Mm->ip[ipp].eqother[ido+1] = aat;
      Mm->ip[ipp].eqother[ido+2] = aac;
      return;
    }
    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, "\nError determining dimension of element");
        fprintf(stderr, " in function anisodam::give_actual_param_a\n(file %s, line %d, elem : %ld)\n", __FILE__, __LINE__, eid+1);
    }
    e = Mm->give_actual_ym(ipp);
    // variable softening for volumetric damage
    tmp = (gf/h - y0/e)*b*e*sin(M_PI/b)/M_PI; 
    if (tmp < 0.0)
    {
      fprintf(stderr, "\nError while computing variable softening for volumetric damage");
      fprintf(stderr, " in function anisodam::give_actual_param_a\n(file %s, line %d, elem : %ld)\n", __FILE__, __LINE__, eid+1);     
    }
    aa = pow(tmp, -b);
    // variable softening for tension
    tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI; 
    if (tmp < 0.0)
    {
      fprintf(stderr, "\nError while computing variable softening for tension");
      fprintf(stderr, " in function anisodam::give_actual_param_a\n(file %s, line %d, elem : %ld)\n", __FILE__, __LINE__, eid+1);     
    }
    aat = pow(tmp, -bt);
    // variable softening for compression
    tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI; 
    if (tmp < 0.0)
    {
      fprintf(stderr, "\nError while computing variable softening for compression");
      fprintf(stderr, " in function anisodam::give_actual_param_a\n(file %s, line : %d, elem : %ld)\n", __FILE__, __LINE__, eid+1);     
    }
    aac = pow(tmp, -bc);
    Mm->ip[ipp].eqother[ido]   = aa;
    Mm->ip[ipp].eqother[ido+1] = aat;
    Mm->ip[ipp].eqother[ido+2] = aac;
  }
}



/**
  Function computes initializes material parameters A, At and Ac. 

  @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

  @retval The function returns computed 
  
  TKo 9.2006
*/
void anisodam::initvalues(long ipp, long ido)
{
  double aa, aat, aac;
  give_actual_param_a(ipp, ido, aa, aac, aat);
}



/**
  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 anisodam::matstiff (matrix &d,long ipp,long ido)
{
//  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 anisodam::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 anisodam::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 anisodam::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, 9.2006
*/
void anisodam::nlstresses (long ipp, long im, long ido)
{
  long ime, idem, i;
  long ncompstr=Mm->ip[ipp].ncompstr;

  // actual material parameters a 
  double aa, aat, aac;
  // old values of driving forces
  double yo;
  vector pyto(3), pyco(3);
  // old values of damage parameters
  double dvo;
  vector pdto(3), pdco(3);
  // increments of driving forces;
  double dy; 
  vector pdyt(3), pdyc(3);
  // new values of driving forces
  double yn;
  vector pytn(3), pycn(3);
  // damage parameters
  double dv;
  vector pdt(3), pdc(3);
  // strains and principal strains
  vector epsn(ncompstr);
  // tensor form of strains and transformation matrix for principal direction of strains
  matrix epst(3,3), t(3,3);
  // stress tensor and stress vector
  matrix sigt(3,3);
  vector sigma(ncompstr);
  // bulk modulus*3, shear modulus, volumetric strain, Young modulus, Poissons ratio
  double k0, g0, epsv, e, nu;
  // actual values of load fuctions
  double lf;
  vector peps(3), psig(3), lft(3), lfc(3);
  double tmp;

  // initialize values of material parameter a
  give_actual_param_a(ipp, ido, aa, aat, aac);

  //  total new strains
  for (i=0;i<ncompstr;i++)
    epsn[i] = Mm->ip[ipp].strain[i];

  // previous principal driving forces and 
  // corresponding previous principal damage parameters
  yo      = Mm->ip[ipp].eqother[ido+3];       // volumetric damage driving forces
  dvo     = Mm->ip[ipp].eqother[ido+3+1+3+3]; // volumetric damage parameter
  for (i=0; i<3; i++)
  {
    // damage driving forces
    pyto[i] = Mm->ip[ipp].eqother[ido+3+1+i];   // deviatoric - tension
    pyco[i] = Mm->ip[ipp].eqother[ido+3+1+3+i]; // deviatoric - compression
    // damage parameters
    pdto[i] = Mm->ip[ipp].eqother[ido+3+1+3+3+1+i];   // deviatoric - tension
    pdco[i] = Mm->ip[ipp].eqother[ido+3+1+3+3+1+3+i]; // deviatoric - compression    
  }
  
  // new driving forces
  vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
  princ_val (epst,peps,t,nijac,limit,Mp->zero,3);
  yn = damdrvforce_vol(ipp, peps);
  damdrvforce_dev(ipp, peps, t, pytn, pycn);

  // incremements of damage driving forces (only positive are assumed)
  dy = yn - yo;
  if (dy < 0.0)
    dy = 0.0;
  for (i=0; i<3; i++)
  {
    // for tension
    pdyt[i] = pytn[i] - pyto[i];
    if (pdyt[i] < 0.0)
      pdyt[i] = 0.0;
    // for compression
    pdyc[i] = pycn[i] - pyco[i];
    if (pdyc[i] < 0.0)
      pdyc[i] = 0.0;
  }
  
  
  // new damage parameters
  lf = loadfuncvol(ipp, peps, dvo, aa);
  dv = dam_vol(ipp, yn, dy, aa, lf);
  if (dv < dvo)
    dv = dvo;
  loadfuncdev(ipp, peps, t, pdto, pdco, aat, aac, lft, lfc);
  pdam_dev(ipp, pycn, pytn, pdyt, pdyc, aat, aac, lft, lfc, pdt, pdc);
  for (i=0; i<3; i++)
  {
    if (pdt(i) < pdto(i))
      pdt(i) = pdto(i);
    if (pdc(i) < pdco(i))
      pdc(i) = pdco(i);
  }

  // new principal stresses
  e  = Mm->give_actual_ym(ipp);
  idem = Mm->ip[ipp].gemid();
  ime = Mm->ip[ipp].idm[idem];
  nu = Mm->eliso[ime].nu;
  k0 = e;
  k0 /= (1-2.0*nu);
  g0 = e;
  g0 /= 2.0*(1.0+nu);
  epsv = peps[0]+peps[1]+peps[2];
  epsv /= 3.0;
  psig[0] = k0 - 2.0*g0;
  if (epsv > 0.0)
    psig[0] -= k0*(dv);
  psig[0] *= epsv;
  psig[1] = psig[2] = psig[0];
  for (i=0; i<3; i++)
  {
    tmp = 1.0;
    if (peps[i] > 0.0)
      tmp -= pdt[i];
    if (peps[i] < 0.0)
      tmp -= pdc[i];
    psig[i] += 2.0*g0*tmp*peps[i];
  }

  // transformation of principal stresses to the global coordinate system
  fillm(0.0, sigt);
  for (i = 0; i < 3; i++)
    sigt[i][i] = psig[i];
  lgmatrixtransf(sigt, t);
  tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);

  //  new data storage
  for (i=0;i<ncompstr;i++)
    Mm->ip[ipp].stress[i]=sigma[i];
  Mm->ip[ipp].other[ido+3] = yn;       // volumetric damage driving forces
  Mm->ip[ipp].other[ido+3+1+3+3] = dv; // volumetric damage parameter
  for (i=0; i<3; i++)
  {
    // damage driving forces
    Mm->ip[ipp].other[ido+3+1+i]   = pytn[i]; // deviatoric - tension
    Mm->ip[ipp].other[ido+3+1+3+i] = pycn[i]; // deviatoric - compression
    // damage parameters
    Mm->ip[ipp].other[ido+3+1+3+3+1+i]   = pdt[i]; // deviatoric - tension
    Mm->ip[ipp].other[ido+3+1+3+3+1+3+i] = pdc[i]; // deviatoric - compression    
  }
}

/**
  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 anisodam::updateval (long ipp,long im,long ido)
{
  long i;
  long n = Mm->givencompeqother(ipp,im);

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