shefplast.cpp

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

    mxi[6] = -1.0* hf ; 
    scprd(tmpv,mxi,tmps) ;
    cmulv(1.0/tmps,tmpv,tmpv);
    vxv(mxi,tmpv,tmp);
    mxm(ainv,tmp,td7);
    subm(ainv,td7,td7);

    // We only keep the 6*6 (in 3D) first components for the tangent stiffness matrix
    for (i=0;i<d.m;i++) {
      for (j=0;j<d.n;j++) {
	td[i][j] = td7[i][j] ;
      }
    }
    destrm(tmp);destrm(dinv);destrm(sig);destrm(dfdstens);destrm(dfdsds);destrm(am);
    destrm(a);destrm(ainv);destrm(td7);

    destrv(dfds);destrv(str);destrv(dfdsdk);destrv(dhds);destrv(dfdk);
    destrv(av);destrv(mxi);destrv(tmpv);
  }
}

/**
   This function computes stresses at given integration point ipp,
   depending on the reached strains.
   The cutting plane algorithm is used. The stress and the other attribute of
   given integration point is actualized.

   @param ipp - integration point number in the mechmat ip array.
   @param ido - index of internal variables for given material in the ipp other array
*/
void shefplast::nlstresses (long ipp, long ido)
//
{
  long i,n=Mm->ip[ipp].ncompstr;
  double gamma;
  vector epsn(n),epsp(n),q(1);

  //  initial values
  for (i=0; i<n; i++) {
    epsn[i]=Mm->ip[ipp].strain[i];
    epsp[i]=Mm->ip[ipp].eqother[ido+i];
  }
  gamma=Mm->ip[ipp].eqother[ido+n];
  q[0] = Mm->ip[ipp].eqother[ido+n+1];
  if (q[0] == 0.0)
    q[0] = kh0;

  //  stress return algorithm
  //Mm->cutting_plane (ipp,gamma,epsn,epsp,q);
  stress_return (ipp,gamma,q,epsn,epsp,ido);
  
  //  new data storage
  for (i=0; i<n; i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+n]=gamma;
  Mm->ip[ipp].other[ido+n+1]=q[0];
}

/**
   This 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 ido - index of internal variables for given material in the ipp other array

*/
void shefplast::updateval (long ipp,long ido)
{
  long i,n = Mm->ip[ipp].ncompstr;

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

  Mm->ip[ipp].eqother[ido+n]=Mm->ip[ipp].other[ido+n];
  Mm->ip[ipp].eqother[ido+n+1]=Mm->ip[ipp].other[ido+n+1];
}

void shefplast::stress_return (long ipp,double &lambda,vector &k,vector &eps,vector &epsp,long ido)
{
  long i,ii,jj,stop,ni;
  double f,hf,nlambda,dlambda,dhdk,nory,err,tmp1,tmp2;
  vector dfds,dfdk,dfdsdk,dhds,nsig,dsig,mxi,tmp,x,y,av,sigtrial,nk,tmp6;
  matrix d,dinv,dfdsds,a,ainv,am,sigtens,nsigtens,dfdstens;
  
  //ni = Mp->nicp;
  //err = Mp->errcp;

  ni=sra.give_ni ();
  err=sra.give_err ();  
  
  allocv (6,dfds);
  allocv (6,dfdsdk);
  allocv (6,dhds);
  allocv (6,nsig);
  allocv (6,dsig);
  allocv (6,sigtrial);
  allocv (7,mxi);
  allocv (7,tmp);
  allocv (6,tmp6);
  allocv (7,x);
  allocv (7,y);
  allocv (6,av);
  allocv (1,nk);
  allocv (1,dfdk);

  allocm (6,6,d);
  allocm (6,6,dinv);
  allocm (6,6,dfdsds);
  allocm (7,7,a);
  allocm (7,7,ainv);
  allocm (6,6,am);
  allocm (3,3,sigtens);
  allocm (3,3,nsigtens);
  allocm (3,3,dfdstens);
  
  nlambda = 0.0;
  nk[0] = k[0];

  //  elastic stiffness matrix
  Mm->elmatstiff (d,ipp);
  // We try to get the same set up as in ASTER but does not change anything
  /*  for (i=3;i<6;i++) {
      d[i][i] *= 2.0 ;
      eps[i] /= 2.0 ;
      epsp[i] /= 2.0 ;
      }*/
  invm(d,dinv,Mp->zero);
  
  //  elastic strain
  subv(eps,epsp,av);

  //  trial stress
  mxv(d,av,sigtrial);
  
  vector_tensor(sigtrial,nsigtens,Mm->ip[ipp].ssst,stress);
  copyv(sigtrial,nsig);
  // ***********************
  //  main iteration loop
  // ***********************
  stop=0;
  for (i=0;i<ni;i++) {
    //  f computing
    f=yieldfunction (nsigtens,nk);
    // if f is negative at the first step material is still elastic
    if (f<err && i==0) {
      stop =1 ;
      break ;
    }

    //  dfds assembling
    deryieldfsigma (ipp,nsigtens,nk,dfdstens,ido);
    tensor_vector (dfds,dfdstens,Mm->ip[ipp].ssst,stress);
    
    //  dfdk computing
    deryieldfq (nsigtens,nk,dfdk);
    
    //  dfdsds assembling
    numdiff_dfdsdsc (ipp,nsigtens,nk,dfdsds,ido);
    
    //  dfdsdk assembling
    numdiff_dfdsdkc (ipp,nsigtens,nk,dfdsdk,ido);
    
    //  h computing
    hf = hardening(ipp,nsigtens,nk,ido);
    
    //  dhds assembling
    numdiff_dhdsc (ipp,nsigtens,nk,dhds,ido);

    //  dhdk computing
    numdiff_dhdkc (ipp,nsigtens,nk,dhdk,ido);
    
    //
    //  Jacobian assembling
    //

    //  block 1,1 (6*6)
    cmulm(nlambda,dfdsds,am);
    addm(am,dinv,am);
    
    for (ii=0;ii<6;ii++) {
      for (jj=0;jj<6;jj++) {
	a[ii][jj]=am[ii][jj];
      }
    }

    //  block 1,2
    cmulv(nlambda,dfdsdk,av);
    
    for (ii=0;ii<6;ii++) {
      a[ii][6]=av[ii];
      mxi[ii] = dfds[ii] ;
    }
    mxi[6] = dfdk[0] ;
    
    //  block 2,1
    cmulv(-1.0*nlambda,dhds,av);
    
    for (ii=0;ii<6;ii++) {
      a[6][ii]=av[ii];
    }
    
    //  block 2,2
    a[6][6]=1.0-nlambda*dhdk;
    // Compute the vector of residuals (7*1)

    //  block 1 (6*1)
    cmulv(nlambda,dfds,av);
    mxv(dinv,nsig,tmp6);
    addv(tmp6,av,av);

    mxv(dinv,sigtrial,tmp6);
    subv(av,tmp6,av);

    for (ii=0;ii<6;ii++) {
      y[ii]=av[ii];
    }
    
    //  block 2
    if(nk[0]==1.0)
      y[6] = 0.0 ;
    else
      y[6]=nk[0]-nlambda*hf-k[0];

    // Compute the Lambda increment dlambda
    invm(a,ainv,Mp->zero);
    vxm(mxi,ainv,tmp);
    scprd(tmp,y,tmp1) ;
    // We change the last value of vector mxi to get the vector (m,-h)
    mxi[6] = -1.0* hf ; 
    scprd(tmp,mxi,tmp2) ;
    dlambda = (f-tmp1)/tmp2 ;
    // Compute the sigma and kappa increment in the x vector
    cmulv(dlambda,mxi,mxi);
    addv(y,mxi,tmp);
    cmulv(-1.0,tmp,tmp);
    mxv(ainv,tmp,x) ;
    
    //
    //  new values
    //
    for (ii=0;ii<6;ii++) {
      dsig[ii]=x[ii];
    }
    
    //  new stresses
    addv(nsig,dsig,nsig);
    vector_tensor (nsig,nsigtens,Mm->ip[ipp].ssst,stress);

    //  new k
    nk[0]+=x[6];
    if(nk[0] < 0.0) {
      nk[0] = 0.0 ;
    }
    if(nk[0] > 1.0) {
      nk[0] = 1.0 ;
    }
    //  new lambda
    nlambda+=dlambda;

      //  norm of the vector of residuals
      nory = fabs(y[0]) ;
      for (ii=1;ii<7;ii++) {
        nory = maxim(nory,fabs(y[ii]));
     }
      nory = maxim(nory,f);

    if(ipp==0)
      fprintf(stderr,"local error %e/%e \n",nory,err);
    if (nory<err) {
      stop=1;
      break;
    }
  }

  if (stop==1) {
    copyv (nk,k);
    Mm->storestress (0,ipp,nsig);
    lambda=nlambda;
    mxv(dinv,nsig,tmp6); 
    subv(eps,tmp6,epsp);
  }
  if(i==ni) {
    fprintf (stderr,"\n\n stress return algo in %s is not successfull after %ld steps for ipp = %ld.\n",__FILE__,ni,ipp);
  }
  
  destrm (dfdstens);  destrm (nsigtens);  destrm (sigtens);
  destrm (am);  destrm (a);  destrm (dfdsds);
  destrm (d); destrm (dinv);
  
  destrv (dfdk);  destrv (nk);  destrv (av);
  destrv (y);  destrv (x);  destrv (sigtrial);
  destrv (nsig);  destrv (dhds);  destrv (dfdsdk);
  destrv (dfds);
}

/**
   function computes zeta
   
   14.1.2004
*/
void shefplast::compzeta()
{
  if (xi>0.0)
    /*    zeta = -ah + sqrt(ah*ah+ch);*/
    zeta = 2.0e-5;
  else
    zeta = -ah + sqrt(ah*ah-bh*xi+ch);
}

double shefplast::hardening(long ipp,matrix &sigtens,vector &q,long ido)
{
  double h;
  matrix dfds;

  if (q[0]<1.0) {
    allocm (3,3,dfds);
    deryieldfsigma (ipp,sigtens,q,dfds,ido);
    compzeta ();
    // The tensor norm already computes the square root !!!
    h = sqrt(2.0/3.0)*tensornorm (dfds)/zeta;
    destrm (dfds);
  }
  else {
    h=0.0;
  }
  
  return h;
}

/**
   function computes numerically the second derivatives of yield function
   with respect to stress tensor
   the first derivatives are expressed explicitly, the second derivatives
   are computed numerically using a centred integration scheme
   
   14.1.2004
*/

void shefplast::numdiff_dfdsdsc(long ipp,matrix &sigtens,vector &q,matrix &dfdsds,long ido)
{
  long i,j;
  double dh;
  long ncomp = Mm->ip[ipp].ncompstr;
  vector sig(ncomp),sigdh(ncomp),sigmdh(ncomp);
  matrix dfds(3,3),dfdsdh(3,3),sigtensdh(3,3),sigtensmdh(3,3);
  
  tensor_vector (sig,sigtens,Mm->ip[ipp].ssst,stress);
  
  for (i=0;i<6;i++) {
    copyv (sig,sigdh);
    copyv (sig,sigmdh);
    sizev(sig, dh);
    dh *= sqrt(DBL_EPSILON);

    sigdh[i]+=dh;
    sigmdh[i]-=dh;
    vector_tensor (sigdh,sigtensdh,Mm->ip[ipp].ssst,stress);
    vector_tensor (sigmdh,sigtensmdh,Mm->ip[ipp].ssst,stress);
    
    deryieldfsigma (ipp,sigtensmdh,q,dfds,ido);
    deryieldfsigma (ipp,sigtensdh,q,dfdsdh,ido);
    
    subm(dfdsdh,dfds,dfds);
    cmulm(0.5/dh,dfds,dfds);
    
    tensor_vector (sigdh,dfds,Mm->ip[ipp].ssst,stress);
    
    for (j=0;j<6;j++) {
      dfdsds[j][i]=sigdh[j];
    }
  }
}

/**
   function computes numerically the second derivatives of yield function
   with respect to stress tensor and hardening parameters
   the first derivatives are expressed explicitly, the second derivatives
   are computed numerically
   
   14.1.2004
*/

void shefplast::numdiff_dfdsdkc(long ipp,matrix &sigtens,vector &q,vector &dfdsdk,long ido)
{
  double dh;
  vector qdh,qmdh;
  matrix dfds,dfdsdh,dfdsmdh;
  
  allocv (1,qdh);
  allocv (1,qmdh);
  allocm (3,3,dfds);
  allocm (3,3,dfdsdh);
  allocm (3,3,dfdsmdh);
  
  sizev(q, dh);
  dh *= sqrt(DBL_EPSILON);

  qdh[0]=q[0]+dh;
  qmdh[0]=q[0]-dh;
  
  deryieldfsigma (ipp,sigtens,qmdh,dfds,ido);
  deryieldfsigma (ipp,sigtens,qdh,dfdsdh,ido);
  
  subm(dfdsdh,dfds,dfds);
  cmulm(0.5/dh,dfds,dfds);
  
  tensor_vector (dfdsdk,dfds,Mm->ip[ipp].ssst,stress);
  
  destrm (dfdsdh); destrm (dfdsmdh);  destrm (dfds);
  destrv (qdh); destrv (qmdh);
}

/**
   function computes numerically derivatives of hardening function
   with respect to stress tensor
   
   14.1.2004
*/
void shefplast::numdiff_dhdsc(long ipp,matrix &sigtens,vector &q,vector &dhds,long ido)
{
  long i;
  long ncomp = Mm->ip[ipp].ncompstr;
  double h,hdh,dh;
  vector sig,sigdh,sigmdh;
  matrix sigtensdh;
  
  allocv (ncomp,sig);
  allocv (ncomp,sigdh);
  allocv (ncomp,sigmdh);
  allocm (3,3,sigtensdh);
  
  tensor_vector (sig,sigtens,Mm->ip[ipp].ssst,stress);
  
  for (i=0;i<6;i++) {
    copyv (sig,sigdh);
    copyv (sig,sigmdh);
    sizev(sig, dh);
    dh *= sqrt(DBL_EPSILON);

    sigdh[i]+=dh;
    sigmdh[i]-=dh;

    vector_tensor (sigdh,sigtensdh,Mm->ip[ipp].ssst,stress);
    vector_tensor (sigmdh,sigtens,Mm->ip[ipp].ssst,stress);
    
    hdh = hardening (ipp,sigtensdh,q,ido);
    h = hardening (ipp,sigtens,q,ido);

    dhds[i]=(hdh-h)*0.5/dh;
  }
  
  destrm (sigtensdh);
  destrv (sigdh);  destrv (sigmdh);  destrv (sig);
}

/**
   function computes numerically derivatives of hardening function
   with respect to hardening parameters
   
   14.1.2004
*/
void shefplast::numdiff_dhdkc(long ipp,matrix &sigtens,vector &q,double &dhdk,long ido)
{
  double h,hdh,dh;
  vector qdh,qmdh;
  
  allocv (1,qdh);
  allocv (1,qmdh);
  
  sizev(q, dh);
  dh *= sqrt(DBL_EPSILON);
  qdh[0]=q[0]+dh;
  qmdh[0]=q[0]-dh;

  hdh = hardening (ipp,sigtens,qdh,ido);
  h = hardening (ipp,sigtens,qmdh,ido);
  
  dhdk = (hdh-h)*0.5/dh;

  destrv (qdh);
  destrv (qmdh);
}


inline double
shefplast :: maxim (double a,double b)
{
  return(a>b ? a:b);
}