mechmat.cpp

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

      if (Mespr==1)  fprintf (stdout,"\n number of artificial material damage-plasticity  %ld",numtype);
      dampl = new damplast [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of artificial material damge-plasticity");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
      }
      break;
    }

    case viscoplasticity:{
      if (Mespr==1)  fprintf (stdout,"\n number of artificial material of viscoplasticity  %ld",numtype);
      visplas = new visplast [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of artificial material of viscoplasticity");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
      }
      break;
    }

    case elasttimemat:{
      if (Mespr==1)  fprintf (stdout,"\n number of artificial elastic isotropic material for time dependent problems %ld",numtype);
      eltimemat = new elasttime [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of artificial elastic isotropic material for time dependent problems");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
      }
      break;
    }
      
    case rcmatmodnorm:{
      if (Mespr==1)  fprintf (stdout,"\n number of artificial material of viscoplasticity  %ld",numtype);
      normm = new normmat [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of artificial material of viscoplasticity");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
	normm[k-1].read (in);
      }
      break;
    }
    
    case lenjonespot:{
      if (Mespr==1)  fprintf (stdout,"\n number of Lennard-Jones interatomic potentials   %ld",numtype);
      lenjon = new lenjonesmat [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of Lennard-Jones interatomic potentials");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
	lenjon[k-1].read (in);
      }
      break;
    }
    
    case contmat:{
      if (Mespr==1)  fprintf (stdout,"\n number of materials for contact problems   %ld",numtype);
      conmat = new contactmat [numtype];
      for (j=0;j<numtype;j++){
        k=numtype+1;
        xfscanf (in,"%ld",&k);
        if (k>numtype || k<1){
          fprintf (stderr,"\n\n wrong number of materials for contact problems");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
	conmat[k-1].read (in);
      }
      break;
    }
      
    default:{
      fprintf (stderr,"\n\n unknown material type is required in function readmatchar (file %s, line %d).\n",
	       __FILE__,__LINE__);
    }
    }

  }
}

/**
   function reads integration points, material characteristics,
   auxiliary points for strain and stress computation
   
   @param in - input stream
   
   21.7.2001
*/
void mechmat::read (XFILE *in)
{

  //  reading of material characteristics
  readmatchar (in);

  //  computation of number of all integration points
  tnip = intpnum ();
  if (Mespr==1){
    fprintf (stdout,"\n number of integration points  %ld",tnip);
  }

  //  allocation of integration points
  intpointalloc ();

  //  reading of integration points
  //readip (in);
  intpointinit ();

/*
  //  reading of strain and stress points
  if (Mp->straincomp==1){
    stra.read(in);
  }
  if (Mp->stresscomp==1){
    stre.read(in);
  }
  */
}



/**
   function returns actual value of the Young modulus for given material and
   in the required integration point
   
   @param ipp - number of integration point
   @param im  - index of the material in the ipp tm and idm arrays. The standard value is zero.
   @param ido - index of internal variables in the ip's ipp eqother array. The standard value is zero.
   
   TKo, 8.8.2005
*/
double mechmat::give_actual_ym(long ipp, long im, long ido)
{
  long idem = ip[ipp].gemid();
  double e = 0.0;

  switch (ip[ipp].tm[im])
  {
    case elisomat:
    case winklerpasternak:
    case rcmatmodnorm:
    case simplas1d:
    case jflow:
    case jflow2:
    case microplaneM4:
    case microsimp:
    case microfibro:
    case mohrcoul:
    case mohrcoulparab:
    case boermaterial:
    case druckerprager:
    case chenplast:
    case modcamclaymat:
    case shefpl:
    case glasgowmechmat:
    case glasgowdamage:
    case simvisc:
    case simvisplas:
    case lemaitr:
    case scaldamage:
    case scaldamagecc:
    case scaldamage_time:
    case anisodamage:
    case creepdpl:
    case creepbaz:
    case consolidation:
    case nonlocplastmat:
    case nonlocdamgmat:
    case nonlocalmicroM4:
    case damage_plasticity:
    case viscoplasticity:
    case elasttimemat:
      e = eliso[ip[ipp].idm[idem]].e;
      break;
    case creep_damage:
      e = give_actual_ym(ipp, im+1, ido);
      break;
    case creepb3:
      e = crb3[ip[ipp].idm[im]].give_actual_ym(ipp, im, ido);
      break;
    case creeprs:
      e = crrs[ip[ipp].idm[im]].give_actual_ym(ipp, im, ido);
      break;
    case time_switchmat:
      e = tswmat[ip[ipp].idm[im]].give_actual_ym(ipp, im, ido);
      break;
    default:
      fprintf (stderr,"\n\n Unknown material type is required in function");
      fprintf (stderr,"\n give_actual_ym (file %s, line %d).\n",__FILE__,__LINE__);
  }
  return e;
}


/**
   function returns actual value of the tensile strenght for given material and
   in the required integration point
   
   @param ipp - number of integration point
   @param im  - index of the material in the ipp tm and idm arrays. The standard value is zero.
   @param ido - index of internal variables in the ip's ipp eqother array. The standard value is zero.
   
   TKo, 8.8.2005
*/
double mechmat::give_actual_ft(long ipp, long im, long ido)
{
  double ft = 0.0;

  switch (ip[ipp].tm[im])
  {
    case scaldamage_time:
      ft = scdamtime[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
      break;
    case creep_damage:
      ft = crdam[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
      break;
    case creepb3:
      ft = crb3[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
      break;
    case creeprs:
      ft = crrs[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
      break;
    case time_switchmat:
      ft = tswmat[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
      break;
    default:
      fprintf (stderr,"\n\n Unknown material type is required in function");
      fprintf (stderr,"\n give_actual_ym (file %s, line %d).\n",__FILE__,__LINE__);
  }
  return ft;
}



/**
   function computes stiffness %matrix of the material
   in the required integration point
   
   @param d - stiffness %matrix
   @param ipp - number of integration point
   @param im  - index of the material in the ipp tm and idm arrays. The standard value is zero.
   @param ido - index of internal variables in the ip's ipp other array
   
   JK, 17.7.2001
*/
void mechmat::matstiff (matrix &d,long ipp,long im,long ido)
{
  long i;
  switch (ip[ipp].tm[im]){
  case elisomat:{
    i=ip[ipp].idm[im];
    eliso[i].matstiff (d,ip[ipp].ssst);
    break;
  }
  case simplas1d:{
    i=ip[ipp].idm[im];
    spl1d[i].matstiff (d,ipp,ido);
    break;
  }
  case jflow:{
    i=ip[ipp].idm[im];
    j2f[i].matstiff (d,ipp,ido);
    break;
  }
  case jflow2:{
    i=ip[ipp].idm[im];
    j2f2[i].matstiff (d,ipp,ido);
    break;
  }
  case microplaneM4:{
    i=ip[ipp].idm[im];
    mpM4[i].matstiff (d,ipp,ido);
    break;
  }
  case microsimp:{
    i=ip[ipp].idm[im];
    mpSIM[i].matstiff (d,ipp,ido);
    break;
  }
  case microfibro:{
    i=ip[ipp].idm[im];
    mpfib[i].matstiff (d,ipp,ido);
    break;
  }
  case mohrcoul:{
    i=ip[ipp].idm[im];
    mcoul[i].matstiff (d,ipp,ido);
    break;
  }
  case mohrcoulparab:{
    i=ip[ipp].idm[im];
    mcpar[i].matstiff (d,ipp,ido);
    break;
  }
  case boermaterial:{
    i=ip[ipp].idm[im];
    boerm[i].matstiff (d,ipp,ido);
    break;
  }
  case druckerprager:{
    i=ip[ipp].idm[im];
    drprm[i].matstiff (d,ipp,ido);
    break;
  }
  case chenplast:{
    i=ip[ipp].idm[im];
    chplast[i].matstiff (d,ipp,ido);
    break;
  }
  case modcamclaymat:{
    i=ip[ipp].idm[im];
    cclay[i].matstiff (d,ipp,ido);
    break;
  }
  case shefpl:{
    i=ip[ipp].idm[im];
    shpl[i].matstiff (d,ipp,ido);
    break;
  }
  case glasgowmechmat:{
    i=ip[ipp].idm[im];
    glasgmat[i].matstiff (d,ipp,ido);
    break;
  }
  case glasgowdamage:{
    i=ip[ipp].idm[im];
    glasgdam[i].matstiff (d,ipp,ido);
    break;
  }
  case creep_damage:{
    i=ip[ipp].idm[im];
    crdam[i].matstiff (d,ipp,im,ido);
    break;
  }
  case time_switchmat:{
    i=ip[ipp].idm[im];
    tswmat[i].matstiff (d,ipp,im,ido);
    break;
  }
  case simvisplas:{
    i=ip[ipp].idm[im];
    svipl[i].matstiff (d,ipp,ido);
    break;
  }
  case lemaitr:{
    i=ip[ipp].idm[im];
    lmtr[i].matstiff (d,ipp,ido);
    break;
  }
  case scaldamage:{
    i=ip[ipp].idm[im];
    scdam[i].matstiff (d,ipp,ido);
    break;
  }
  case scaldamagecc:{
    i=ip[ipp].idm[im];
    scdamcc[i].matstiff (d,ipp,ido);
    break;
  }
  case scaldamage_time:{
    i=ip[ipp].idm[im];
    scdamtime[i].matstiff (d,ipp,ido);
    break;
  }
  case anisodamage:{
    i=ip[ipp].idm[im];
    anidam[i].matstiff (d,ipp,ido);
    break;
  }
  case graphm:{
    i=ip[ipp].idm[im];
    grmat[i].matstiff (d,ipp);
    break;
  }
  case graphmtime:{
    i=ip[ipp].idm[im];
    grmattime[i].matstiff (d,ipp);
    break;
  }
  case geoelast:{
    i=ip[ipp].idm[im];
    geoel[i].matstiff (d,ipp,ido);
    break;
  }
  case creeprs:
  case creepb3:{
    creep_matstiff (d,ipp,im,ido);
    break;
  }
  case creepdpl:{
    creep_matstiff (d,ipp,im,ido);
    break;
  }
  case creepbaz:{
    i=ip[ipp].idm[im];
    crbaz[i].matstiff (d,ipp);
    break;
  }
  case consolidation:{
    i=ip[ipp].idm[im];
    csol[i].matstiff (d,ipp);
    break;
  }
  case winklerpasternak:{
    i=ip[ipp].idm[im];
    wpast[i].matstiff (d,ipp);
    break;
  }

  case nonlocdamgmat:{
    matstiff (d,ipp,im+1,ido);
    break;
  }
  case nonlocplastmat:{
    matstiff (d,ipp,im+1,ido);
    break;
  }
  case nonlocalmicroM4:{
    matstiff (d,ipp,im+1,ido);
    break;
  }

  case damage_plasticity:{
    i=ip[ipp].idm[im];
    dampl[i].matstiff (d,ipp);
    break;
  }

  case viscoplasticity:{
    i=ip[ipp].idm[im];
    visplas[i].matstiff (d,ipp,im,ido);
    break;
  }

  case elasttimemat:{
    i=ip[ipp].idm[im];
    eltimemat[i].matstiff (d,ipp,im,ido);
    break;
  }
  case rcmatmodnorm:{
    i=ip[ipp].idm[im+1];
    eliso[i].matstiff (d,ip[ipp].ssst);
    break;
  }

  case contmat:{
    i=ip[ipp].idm[im+1];
    conmat[i].matstiff (d,ipp);
    break;
  }

  default:{
    fprintf (stderr,"\n\n unknown material type is required in function matstiff (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

/**
   function computes elastic stiffness %matrix of the material
   in the required integration point