microm4.cpp

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

      projN[mPlane][i] = n(ii)*n(jj);
      projM[mPlane][i] = 0.5*(m(ii)*n(jj) + m(jj)*n(ii)) ;
      projL[mPlane][i] = 0.5*(l(ii)*n(jj) + l(jj)*n(ii)) ;
    }
  }
  return;
}

void microM4::matstiff (matrix &d,long ipp,long ido)
/**  function returns elastic stiffness matrix
  @param d - elastic stiffness matrix 
  @param ipp - integration point pointer
  @param ido - index of internal variables for given material in the ipp other array
*/
{
  Mm->elmatstiff (d,ipp);
}

void microM4::nlstresses (long ipp, long ido)
{
  long i, i_other,mPlaneIndex;
  double previousSigmaV,previousSigmaL,previousSigmaM,previousSigmaN,previousSigmaD;
  double sel,sem,sev,sed,f;
  double cv,cd,meanN=0.;
  double sigmaN,sigmaL,sigmaM,sigmaV,sigmaD;
  double epsV,depsV, epsD, epsN, depsD, depsM, depsL,depsN;
  vector macroStrain(6), macroStrainIncrement(6), macroSigma(6);

  //  initial values of strain vector and strain increment
  for (i=0;i<6;i++){
    macroStrain[i] = Mm->ip[ipp].strain[i];
    macroStrainIncrement[i] = Mm->ip[ipp].strain[i] - Mm->ip[ipp].eqother[ido+i+numberOfMicroplanes*3+1];

    //UPDATE history variables= store total macroStrain
    Mm->ip[ipp].other[ido+i+numberOfMicroplanes*3+1]=Mm->ip[ipp].strain[i];
  }

  //volumetric microstrain and microstrain increment
  epsV = (macroStrain[0]+macroStrain[1]+macroStrain[2])/3.0;
  depsV= (macroStrainIncrement[0]+macroStrainIncrement[1]+macroStrainIncrement[2])/3.0;

  //ask for previous volumetric stress
  previousSigmaV=Mm->ip[ipp].eqother[ido+0];

  //loop over all microplanes
  for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {

    //compute total microstrains on microplane and its increments
    epsN=0., depsL=0., depsM=0.; depsN=0.;
    for (i=0; i<6; i++) {
      epsN += projN [mPlaneIndex][i]*macroStrain[i];
      depsL += projL [mPlaneIndex][i]*macroStrainIncrement[i];
      depsM += projM [mPlaneIndex][i]*macroStrainIncrement[i];
      depsN += projN [mPlaneIndex][i]*macroStrainIncrement[i];
    }
    epsD = epsN - epsV;
    depsD = depsN - depsV;

    //ask for previous stress components on microplane
    i_other=ido+(mPlaneIndex+1)*3-2;//poradi komponenty v poli history parametru L,M,N
    previousSigmaL=Mm->ip[ipp].eqother[i_other];
    previousSigmaM=Mm->ip[ipp].eqother[i_other+1];
    previousSigmaN=Mm->ip[ipp].eqother[i_other+2];
    previousSigmaD=previousSigmaN-previousSigmaV;

    //evaluation of new microplane stresses

    //secant moduli
    cv = ev;
    cd = ed;

    //volumetric microstress
    sev=previousSigmaV + cv*depsV;
    sigmaV=minim( maxim(sev,FVminus(epsV)), FVplus(epsV));
    //deviatoric microstress
    sed=previousSigmaD + cd*depsD;
    sigmaD=minim(maxim(sed,FDminus(epsD)),FDplus(epsD));
    //normal microstress
    sigmaN=minim( (sigmaV+sigmaD) ,FN(epsN,previousSigmaV));
    //shear microstresses
    sel=previousSigmaL + et*depsL;
    sem=previousSigmaM + et*depsM;
    f=FT(sigmaN,epsV);
    if (sel>f) sigmaL=f; else if (sel<-f) sigmaL=-f; else sigmaL=sel;
    if (sem>f) sigmaM=f; else if (sem<-f) sigmaM=-f; else sigmaM=sem;

    //mean normal microstress (computed after sweeping throgh all microplanes)
    meanN += sigmaN* microplaneWeights[mPlaneIndex]*6.0;

    //UPDATE history variables, i.e. store microstresses
    Mm->ip[ipp].other[i_other]=sigmaL;
    Mm->ip[ipp].other[i_other+1]=sigmaM;
    Mm->ip[ipp].other[i_other+2]=sigmaN;

  }//end 1st loop over nmp

  //volumetric stress is the same for all  mplanes
  //and does not need to be homogenized .
  //Only updating accordinging to mean normal stress must be done.
     if (sigmaV > (meanN/3.))  {sigmaV = meanN/3.;}

  //UPDATE history variables= store posibly new sigmaV
  Mm->ip[ipp].other[ido+0]=sigmaV;

  fillv(0,macroSigma);//zero macroSigma

  //compute macrostresstensor
  for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {
    i_other=ido+(mPlaneIndex+1)*3-2;//poradi komponenty v poli history parametru L,M,N
    sigmaL=Mm->ip[ipp].other[i_other];
    sigmaM=Mm->ip[ipp].other[i_other+1];
    sigmaN=Mm->ip[ipp].other[i_other+2];
    //recalculation of sigmaD according to possibly new sigmaV
    sigmaD = sigmaN - sigmaV;


    for (i=0; i<6; i++) {
                      macroSigma[i] +=((projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD +
		       projL[mPlaneIndex][i]*sigmaL +
		       projM[mPlaneIndex][i]*sigmaM)* microplaneWeights[mPlaneIndex]*6.0;

		      /*		   macroSigma[i]+=(projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD*microplaneWeights[mPlaneIndex]*6.0;  */
    }

  }//end 2nd loop over nmp

  //2nd constraint, addition of volumetric part
  macroSigma[0]+=sigmaV;
  macroSigma[1]+=sigmaV;
  macroSigma[2]+=sigmaV;

  //send new macrostress vector to the solver
  Mm->storestress(0,ipp,macroSigma);

}


void microM4::nonloc_nlstresses (long ipp, long ido)
{
  long i, j;
  matrix e(6,6);//elastic stiffness matrix
  double sigma_elast;
  vector sigma_nonloc(6);

  //compute elas. stiff. matrix for ipp
  Mm->elmatstiff(e,ipp);

  for (i=0;i<6;i++){
    sigma_elast=0;
    for (j=0;j<6;j++){
      sigma_elast += e[i][j] * (Mm->ip[ipp].strain[j]);
    }
    //nonlocal stess
    sigma_nonloc[i]=Mm->ip[ipp].nonloc[i] + sigma_elast;
  }

  //send nonlocal stress vector to the solver
  Mm->storestress(0,ipp,sigma_nonloc);
}

void microM4::updateval(long ipp, long im, long ido)
{
  long i,n = Mm->givencompother(ipp, im);

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

void microM4::changeparam (atsel &atm,vector &val)
{
  long i;
  
  for (i=0;i<atm.num;i++){
    switch (atm.atrib[i]){
    case 0:{
      k1=val[i];
      break;
    }
    case 1:{
      k2=val[i];
      break;
    }
    case 2:{
      k3=val[i];
      break;
    }
    case 3:{
      k4=val[i];
      break;
    }
    case 4:{
      c3=val[i];
      break;
    }
    case 5:{
      c20=val[i];
      break;
    }
    case 6:{
      e=val[i];
      break;
    }
    case 7:{
      nu=val[i];
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong number of atribute in function changeparam (%s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }

  ev = e/(1-2*nu);
  ed = 5*e/(2+3*mu)/(1+nu);
  et = mu*ed;
}

inline double
microM4  :: macbra(double x) /* Macauley bracket = positive part of x */
{
  return(maxim(x,0.0));
}

inline double
microM4 :: FVplus (double epsV)
/*positive volumetric boundary */
{
  return(ev*k1*c13/(1+(c14/k1)*macbra(epsV-c13*c15*k1)));
}

inline double
microM4 :: FVminus (double epsV)
/*negative volumetric boundary */
{
return(-e*k1*k3*exp(-epsV/(k1*k4)));
}

inline double
microM4 :: FDminus(double epsD)
 /*negative deviatoric boundary */
 {
   double a;
   a=macbra(-epsD-c8*c9*k1)/(k1*c7);
   return(-e*k1*c8/(1+a*a));
 }

inline double
microM4  :: FDplus(double epsD)
/*positive deviatoric bondary */
{
  double a;
  a=(macbra(epsD-c5*c6*k1)/(k1*c20*c7));
  return (e*k1*c5/(1+a*a));
}

inline double
microM4 :: FN(double epsN,double sigmaV)
 /*normal boundary */
{
  return(e*k1*c1*exp(-macbra(epsN-c1*c2*k1)/(k1*c3+macbra(-c4*(sigmaV/ev)))));
}

inline double
microM4 ::  FT (double sigmaN,double epsV)
/*shear boundary */
{
  double a, sn0;
  sn0= et*k1*c11/(1+c12*fabs(epsV));
  a=macbra(-sigmaN+sn0);
  return(et*k1*k2*c10*a/(et*k1*k2+c10*a));
}

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

inline double
microM4 :: minim (double a,double b)
{
  return(a<b ? a:b);
}