mechmat.cpp

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

#include <math.h>
#include <string.h>
#include "mechmat.h"
#include "global.h"
#include "globmat.h"
#include "elemhead.h"
#include "vecttens.h"
#include "intpoints.h"
#include "element.h"
#include "node.h"
#include "elastisomat.h"
#include "elastgmat3d.h"
#include "splas1d.h"
#include "j2flow.h"
#include "j2flow2.h"
#include "microM4.h"
#include "microSIM.h"
#include "microfiber.h"
#include "mohrc.h"
#include "mohrcparab.h"
#include "boermat.h"
#include "drprag.h"
#include "camclay.h"
#include "shefplast.h"
#include "glasgmech.h"
#include "glasgowdam.h"
#include "creepdam.h"
#include "timeswmat.h"
#include "simviscous.h"
#include "svispl.h"
#include "lemaitre.h"
#include "anisodam.h"
#include "scaldam.h"
#include "scaldamcc.h"
#include "scaldamtime2.h"
#include "aci78.h"
#include "cebfip78.h"
#include "graphmat.h"
#include "graphmattime.h"
#include "geoelast.h"
#include "creep_b3.h"
#include "creep_rspec.h"
#include "creepb.h"
#include "creep_dpl.h"
#include "consol.h"
#include "winpast.h"
#include "nonlocmicroM4.h"
#include "therisomat.h"
#include "therisomattime.h"
#include "nonlocplast.h"
#include "nonlocdamg.h"
#include "damplast.h"
#include "tensor.h"
#include "visplast.h"
#include "elasttime.h"
#include "normmat.h"
#include "chen.h"
#include "lenjonesmat.h"
#include "contactmat.h"
#include "elemswitch.h"


mechmat::mechmat ()
{
  nmt=0;  tnip=0;
  plast=0;

  max_nncompo=0; max_encompo=0;
  max_ncompstrn=0; max_ncompstre=0;

  ip=NULL;
  ic=NULL; 
  elip = NULL;
  eliso=NULL;  elgm3d=NULL;
  spl1d=NULL;  j2f=NULL;  j2f2=NULL;  mcoul=NULL;  mcpar=NULL;  boerm=NULL;  drprm=NULL;  cclay=NULL;  shpl=NULL;
  mpM4=NULL;  mpSIM=NULL;  mpfib=NULL;  nmpM4=NULL;
  svis=NULL;  svipl=NULL;
  scdam=NULL; scdamtime=NULL; anidam=NULL;
  aci78mod=NULL;  cebfip78mod=NULL;
  grmat=NULL;  grmattime=NULL; geoel=NULL;
  crb3=NULL;
  crrs=NULL;
  crdpl=NULL;
  crbaz=NULL;  csol=NULL;
  tidilat=NULL;
  tidilattime=NULL;
  shpl = NULL; 
  glasgmat = NULL;
  glasgdam = NULL;
  crdam = NULL;

  dampl = NULL;  visplas = NULL;  eltimemat = NULL;  normm = NULL;  chplast = NULL;
  lenjon = NULL;
  conmat=NULL;

  ipv = NULL;  eigstrains = NULL;  eigstresses = NULL;
  tempstrains = NULL;  tempr = NULL;  moist = NULL;
}

mechmat::~mechmat ()
{
  long i;

  delete [] ip;
  delete [] elip;
  delete [] eliso;  delete [] elgm3d;
  delete [] spl1d;  delete [] j2f;  delete [] j2f2;  delete [] mcoul;  delete [] mcpar;
  delete [] boerm;  delete [] drprm;  delete [] cclay;  delete [] shpl;
  delete [] mpM4;  delete [] mpSIM;  delete [] mpfib;   delete [] nmpM4;
  delete [] svis;  delete [] svipl;
  delete [] scdam; delete [] glasgmat;    delete [] glasgdam; delete [] crdam; delete [] scdamtime; delete [] anidam;
  delete [] shpl;
  delete [] aci78mod;  delete [] cebfip78mod;
  delete [] grmat; delete [] grmattime; delete [] geoel;
  delete [] crb3;
  delete [] crrs;
  delete [] crdpl;
  delete [] crbaz;  delete [] csol;
  delete [] tidilat;
  delete [] tidilattime;

  delete [] dampl;  delete [] visplas;  delete [] eltimemat;  delete [] normm;  delete [] chplast;
  delete [] lenjon;
  delete [] conmat;

  if (ic != NULL)
  {
    for (i = 0; i < tnip; i++)
      delete [] ic;
  }
  delete [] ic;
  
  delete [] ipv;
  delete [] tempr;
  delete [] moist;
  
  if (eigstrains != NULL){
    for (i=0;i<tnip;i++){
      delete [] eigstrains[i];
    }
    delete [] eigstrains;
  }
  if (eigstresses != NULL){
    for (i=0;i<tnip;i++){
      delete [] eigstresses[i];
    }
    delete [] eigstrains;
  }
  if (tempstrains != NULL){
    for (i=0;i<tnip;i++){
      delete [] tempstrains[i];
    }
    delete [] tempstrains;
  }
  
}

/**
   function returns number of integration points and determines pointers
   to integration points on elements
   
   8.7.2001
*/
long mechmat::intpnum (void)
{
  long i,j,k,n,nb;
  
  n=0;
  for (i=0;i<Mt->ne;i++){
    nb=Mt->give_nb (i);
    Mt->elements[i].ipp = new long* [nb];
    for (j=0;j<nb;j++){
      Mt->elements[i].ipp[j] = new long [nb];
    }
    for (j=0;j<nb;j++){
      for (k=0;k<nb;k++){
        Mt->elements[i].ipp[j][k]=n;
        n+=Mt->give_nip (i,j,k);
      }
    }
  }
  
  return n;
}

/**
   function allocates integration points
   
   JK, 28.11.2006
*/
void mechmat::intpointalloc ()
{
  ip = new intpoints [tnip];
}

/**
   function reads material types and material id

   @param in - input stream
   
   8.7.2001
*/
void mechmat::readip (FILE *in)
{
  long i;
  for (i=0;i<tnip;i++){
    ip[i].read (in);
  }
}

/**
   function initiates integration point variables
   
   21.8.2001
*/
void mechmat::init_ip_1 (void)
{
  long i,j,ne,ii,jj,ipp,nb,nip,ncomp;
  strastrestate ssst;
  
  ne=Mt->ne;

  for (i=0;i<ne;i++){
    nb=Mt->give_nb (i);
    ncomp=Mt->give_ncomp (i);
    for (ii=0;ii<nb;ii++){
      ssst=Mt->give_ssst (i,ii);
      for (jj=0;jj<nb;jj++){
        ipp=Mt->elements[i].ipp[ii][jj];
        nip=Mt->give_nip (i,ii,jj);
        for (j=0;j<nip;j++){
          ip[ipp].ncompstr=ncomp;
          ip[ipp].ssst=ssst;
          
          if ((ssst == planestrain) || (ssst == planestress))
            ip[ipp].ncompstr=4;
          
          ipp++;
        }
      }
    }
  }
}

/**
   function initiates integration point variables
   
   JK, 21.8.2001
*/
void mechmat::init_ip_2 (void)
{
  long i,j,ne,ii;
  
  ne=Mt->ne;

  for (i=0;i<tnip;i++){
    ip[i].alloc(Mb->nlc,i);

    ii = ip[i].hmt & 2; // the material is nonlocal
    if (ii > 0)
      // Additional allocation for the nonlocal material models is neccessary
      {
	switch (ip[i].tm[0])
	  {
	  case nonlocplastmat:
	    nlplast->allocation(i);
	    break;
	  case nonlocdamgmat:
	    nldamg->allocation(i);
	    break;
	  case nonlocalmicroM4:
	    nmpM4->allocation(i);
	    break;
	  default:
	    fprintf (stderr,"\n\n unknown nonlocal material type is required in function");
	    fprintf (stderr,"\n alloc_ip (file %s, line %d).\n",__FILE__,__LINE__);
	  }
      }
    
    ii = ip[i].hmt & 1; // presence of thermal dilatation material model
    if (ii > 0){
      if (tempr==NULL){
	tempr = new double [tnip];
	for (j=0;j<tnip;j++){
	  tempr[j]=0.0;
	}
	
      }
      
      
    }

  }
  

}


/**
   function defines material models on integration points
   material models are read on elements
   
   JK, 28.11.2006
 */
void mechmat::intpointinit ()
{
  long i, j, ii, jj, k, nb, nip, ipp;
  
  if ((scdam != NULL) || (scdamtime != NULL))
    elip = new long[tnip];
  
  for (i=0;i<Mt->ne;i++){
    nb=Mt->give_nb (i);
    for (ii=0;ii<nb;ii++){
      for (jj=0;jj<nb;jj++){
        nip=Mt->give_nip(i,ii,jj);
        ipp=Mt->elements[i].ipp[ii][jj];
        for (j=0;j<nip;j++){
          if (elip != NULL)
            elip[ipp] = i;
          ip[ipp].nm  = Mt->elements[i].nm;
          ip[ipp].tm  = new mattype[ip[ipp].nm];
          ip[ipp].idm = new long[ip[ipp].nm];
          for (k = 0; k < ip[ipp].nm; k++)
          {
            ip[ipp].tm[k]  = Mt->elements[i].tm[k];
            ip[ipp].idm[k] = Mt->elements[i].idm[k];
            switch (ip[ipp].tm[k])
            {
              case therisodilat :
                ip[ipp].hmt ^= 1;
                break;
	      case therisodilattime :
                ip[ipp].hmt ^= 1;
                break;
              case nonlocplastmat :
                ip[ipp].hmt ^= 2;
                break;
              case nonlocdamgmat :
                ip[ipp].hmt ^= 2;
                break;
              case nonlocalmicroM4 :
                ip[ipp].hmt ^= 2;
                break;
              default :
                break;
	    }
	  } 
          ipp++;
        }
      }
    }
  }
}

/**
   function reads material characteristics
   
   @param in - input stream
   
   8.7.2001
*/
void mechmat::readmatchar (XFILE *in)
{
  long i,j,k;
  long numtype;
  mattype mtype;

  xfscanf (in,"%k%ld","number_of_material_types",&nmt);
  if (nmt<1){
    fprintf (stderr,"\n\n wrong number of material types in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
  }
  if (Mespr==1)  fprintf (stdout,"\n number of different types of materials  %ld",nmt);

  for (i=0;i<nmt;i++){
    xfscanf (in,"%d %ld",(int*)&mtype,&numtype);
    if (numtype<1){
      fprintf (stderr,"\n\n wrong number of material characteristics in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
    }
    
    switch (mtype){
    case elisomat:{
      if (Mespr==1)  fprintf (stdout,"\n number of elastic isotropic materials  %ld",numtype);
      eliso = new elastisomat [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 elastic isotropic material");
          fprintf (stderr,"\n in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        eliso[k-1].read (in);
      }
      break;
    }
    case simplas1d:{
      if (Mespr==1)  fprintf (stdout,"\n number of simple plastic 1D materials   %ld",numtype);
      spl1d = new splas1d [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 simple plastic 1D material in function");
          fprintf (stderr,"\n mechmat::readmatchar (mechmat.cpp).\n");
        }
        spl1d[k-1].read (in);
      }
      break;
    }
    case jflow:{
      if (Mespr==1)  fprintf (stdout,"\n number of J2 flow materials   %ld",numtype);
      j2f = new j2flow [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 J2 flow material in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        j2f[k-1].read (in);
      }
      break;
    }
    case jflow2:{
      if (Mespr==1)  fprintf (stdout,"\n number of J2 flow materials   %ld",numtype);
      j2f2 = new j2flow2 [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 J2 flow material in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        j2f2[k-1].read (in);
      }
      break;
    }
    case microplaneM4:{
      if (Mespr==1)  fprintf (stdout,"\n number of microplane materials   %ld",numtype);
      mpM4 = new microM4 [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 microplane material in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        mpM4[k-1].read (in);
      }
      break;
    }
    case nonlocalmicroM4:{
      if (Mespr==1)  fprintf (stdout,"\n number of nonlocal microplane materials   %ld",numtype);
      nmpM4 = new nonlocmicroM4 [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 nonlocal microplane material in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        nmpM4[k-1].read (in);
      }
      break;
    }
    case microsimp:{
      if (Mespr==1)  fprintf (stdout,"\n number of microplane materials   %ld",numtype);
      mpSIM = new microSIM [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 microplane material in function mechmat::readmatchar (file %s, line %d).\n",__FILE__,__LINE__);
        }
        mpSIM[k-1].read (in);
      }
      break;
    }
    case microfibro:{
      if (Mespr==1)  fprintf (stdout,"\n number of microplane materials   %ld",numtype);