microfiber.cpp

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

      if(fabs(aux) > 1.0e-12){
	m[0] = n[2]/aux;
	m[1] = 0.0;
	m[2] = -n[0]/aux;
      }
      else {
	m[0] = 0.0;
	m[1] = 0.0;
	m[2] = 1.0;
      }
    }
    
    //l.beVectorProductOf (m,n);
    crprd(m,n,l);
    
    for (i=0; i<6; i++) {
      ii = ij[i][0]-1;
      jj = ij[i][1]-1;
      
      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 microfiber::matstiff (matrix &d,long ipp,long ido)
  //  function returns elastic stiffness matrix
  //  d - elastic stiffness matrix
{
  Mm->elmatstiff (d,ipp);
}

void microfiber::nlstresses (long ipp,long ido)
{
  long i, i_other,mPlaneIndex;
  double previousSigmaV,previousSigmaL,previousSigmaM,previousSigmaN,previousSigmaD,
    previousSigmaN_fiber,previousEpsN;
  double sel,sem,sev,sed,f,senf;
  double cv,cd,meanN=0., fiber_module=0;
  //double sigmaN,sigmaL,sigmaM,sigmaV,sigmaD;
  double sigmaL,sigmaM,sigmaV,sigmaD;
  double sigmaN_concrete, sigmaD_concrete, sigmaN_fiber;
  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*5+1];

    //UPDATE history variables= store total macroStrain
    Mm->ip[ipp].other[ido+i+numberOfMicroplanes*5+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)*5-4;//poradi komponenty v poli history parametru L,M,N,sigN_fib,EpsN_fib
    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;

    previousSigmaN_fiber=Mm->ip[ipp].eqother[i_other+3];
    previousEpsN=Mm->ip[ipp].eqother[i_other+4];

    //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_concrete=minim( (sigmaV+sigmaD) ,FN(epsN,previousSigmaV));

    //napeti ve vlakne:
    fiber_module=(previousEpsN!=0.0 ? (previousSigmaN_fiber/previousEpsN) : 0.0 );
    senf=previousSigmaN_fiber+fiber_module*depsN;

   
    /*if(depsN<0) {
if (ipp==1) printf("mp=%d : depsN=%ef : senf=%lf : FF=%lf\n",mPlaneIndex,depsN,senf,FFiber(macbra(epsN)));
}*/

    sigmaN_fiber=(depsN<0.0 ? senf : FFiber(macbra(epsN)));
    
    //// sigmaN_fiber=macbra( minim(senf, FFiber(macbra(epsN))));
    //sigmaN_fiber= FFiber(macbra(epsN));
    // sigmaN=sigmaN_concrete+ sigmaN_fiber;

    //shear microstresses
    sel=previousSigmaL + et*depsL;
    sem=previousSigmaM + et*depsM;    
    f=FT(sigmaN_concrete,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_concrete* 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_concrete;
    Mm->ip[ipp].other[i_other+3]=sigmaN_fiber;
    Mm->ip[ipp].other[i_other+4]=epsN;
    
  }//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)*5-4;//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_concrete=Mm->ip[ipp].other[i_other+2];
    //recalculation of sigmaD according to possibly new sigmaV 
    sigmaD_concrete = sigmaN_concrete - sigmaV;
    sigmaN_fiber=Mm->ip[ipp].other[i_other+3];
    
    
    
      for (i=0; i<6; i++) {
    
	  macroSigma[i] +=((projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD_concrete + 
		       projN[mPlaneIndex][i]*sigmaN_fiber +
		       projL[mPlaneIndex][i]*sigmaL +
		       projM[mPlaneIndex][i]*sigmaM)* microplaneWeights[mPlaneIndex]*6.0;
		       }
	

      /*  macroSigma[i] += ((projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD_concrete + 
	  projN[mPlaneIndex][i]*sigmaN_fiber)* 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 microfiber::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];
  }
}

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

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

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

inline double
microfiber :: 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
microfiber  :: 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
microfiber :: 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
microfiber ::  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
microfiber ::  FFiber (double epsN)
/*fiber boundary */
{
  double f1, f2, a, x;
  x=epsN-k10;
  f1=e*k13*(2*sqrt(x/k11) - x/k11);
  a=1-(epsN-k10-k11)/k12;
  f2=e*k13*a*a;

  if ((epsN<k10) || (epsN>(k10+k11+k12)) ) return(0.0);
  else { 
    if( epsN<=(k10+k11) ) return (f1);
    else return (f2);
  }
}

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

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