creepb-f.cpp

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

#include "creepb.h"
#include "global.h"
#include "math.h"
#include "globmat.h"
#include "genfile.h"
#include "adaptivity.h"
#include "intpoints.h"
#include "stdlib.h"
#include "elastisomat.h"
#include "therisomat.h"

creepb::creepb (void)
{
  timeMax=Mp->timecon.endtime ();
  double m = log10(timeMax);
  nRetTime=7;  
  allocv (nRetTime,retTime);
  allocv (nRetTime,ert);
  retTime[0]=1.0e-9;
  retTime[1]=1.0e-3;
  retTime[2]=1.0;
                  //  retTime[3]=14.0;  //retTime[4]=100.0; //retTime[5]=800.0;
  retTime[3]=pow(10.,0.25*m);
  retTime[4]=pow(10.,0.5*m);
  retTime[5]=pow(10.,0.75*m);
  retTime[6]=timeMax;
  type_h=1;
  type_temp=1;
  h_s = 0.4;   // end humidity
  desht=0.0;
  k_s = 1.0;	   //shape factor
  k_d = 1.0;      // thicknes
  timemat=0.0;
  alfa=0.0;
  allocm(nRetTime, nRetTime, apom); 
  ccTime = -1.0e30;
}

creepb::~creepb (void)
{
}

void creepb::creepinit (long ipp,double val,nonmechquant nmq)
{
  long nc=Mm->ip[ipp].ncompstr;
  
  if (nmq==moist)
    Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=val;
  if (nmq==temperature)
    Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=val;
}

void creepb::read (FILE *in)
{
 fscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
	          &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&k_d,&type_h,&type_temp);
  
  if (type_h==1){fscanf (in,"%lf ",&h_s); }
  if (type_temp==1){fscanf (in,"%lf ",&temp_s); }

}
/**
   function approximates function defined by nodal values

   @param areacoord - vector containing area coordinates
   @param nodval - nodal values

*/
double creepb::approx (vector &areacoord,vector &nodval)
{
  double f;
  scprd (areacoord,nodval,f);
  return f;
}

/**
   function inversion sym. matrix

   @param a - matrix

*/
void creepb::inv_sym (matrix &a)
{
  long i,j,k;
  double diag,an;
    for (k=0;k<a.n;k++){
       diag = a[k][k];
       if (diag<=0.0) { abort();}
       for (i=0;i<a.n;i++){
         an = a[i][k]/diag;
         if (i!=k) { 
           for (j=i;j<a.n;j++){
             if(j!=k) {
               a[i][j] = a[i][j] - an*a[k][j];
               a[j][i] = a[i][j];
			 }
		   }
		 }
         a[i][k] = an;
         a[k][i] = an;
       }
       a[k][k] = - 1.0/diag;
    }

	for (i=0;i<a.n;i++){
       for (j=0;j<a.n;j++){
         a[i][j] = - a[i][j];
	   }
	}
}

void creepb::updateval()
{
  long i,j,k;
  double et,pomt,pomt1,qq,delYi,delYj,m0,m;
  
  if (ccTime != Mp->timecon.actualtime ())
  {
    //ddTime=Mp->timefun.getval(Mp->time);
    ddTime=Mp->timecon.actualforwtimeincr ();
    //ccTime=Mp->time;
    ccTime=Mp->timecon.actualtime ();
    pomt=ccTime+ddTime/2.;
    qq=2./3.;
    //   from concrete starts   tb 
    fillm(0.0, apom); 
    if (timemat<=pomt) {
      pomt1=pomt+0.0001;   
      m0 = log10(pomt1);
      m = log10(2.*timeMax)*0.999;
      k=0;
      while (pomt1<=2.*timeMax)
      {
        for (i=0;i<nRetTime;i++)
        {  
	  delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
	  for (j=0;j<nRetTime;j++)
	  {
	    delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
	    apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
	  }
        }
        k=k+1;
        pomt1=pow(10.,m0+(m-m0)/40*k);
      }
      inv_sym(apom);
    }
    fprintf(stdout, "\n ### inversion\n");
  }
}


/**
   function returns eps from history
   
   @param screep - vector of history
   @param epsscr - vector deformation of history
   
   10.10.2002
*/
void creepb::nlstresses (long ipp)
{
  //  number of components od strain/stress arrays
  nc=Mm->ip[ipp].ncompstr;
  
  //  type of stress state
  //bar=1,plbeam=2,spacebeam=5,
  //planestress=10,planestrain=11,plate=15,
  //axisymm=20,shell=25,spacestress=30
  ss=Mm->ip[ipp].ssst;  
  
  //timeMax=Mp->end_time;
  //ddTime=Mp->timefun.getval(Mp->time);
  ddTime=Mp->timecon.actualforwtimeincr ();
  //ccTime=Mp->time;
  ccTime=Mp->timecon.actualtime ();
  
  //  id of elastic material
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  
  //  actual moisture
  if(type_h!=1){
    h_s=Mm->moist[ipp];
  }
  
  //  actual temperature
  if(type_temp!=1){
    temp_s=Mm->tempr[ipp];
  }
  
  get_h (ipp);     //initial  moist
  get_temp (ipp);     //initial  tempr
  
  if (Mp->phase==1){
    phase1(ipp);
  }
  
  if (Mp->phase==2){
    phase2(ipp);
  }
  
}

void creepb::phase1 (long ipp)
{
  //  right hand side from time dependent value computation
  //  creep function for shrink from pomt to cctime
  long i,j,ii;
  double qq,dj,endTime;
  
  vector epscr(nc),depsshtem(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);
  
  //    Mm->matstiff(d,ipp);  //compute desht
  //    e0=Mm->eliso[imat].e;
  //    mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  endTime=ccTime+ddTime;
  get_desht (desht, ccTime, endTime);
  for (i=0; i<nc; i++){
    epscr[i]=0.0;
  }
  qq=2.0/3.0;
  for (i=0;i<nRetTime;i++){
    ii=nc*(i+1);
    dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
    for (j=0; j<nc; j++){
      screep[j][i]=Mm->ip[ipp].eqother[ii+j];	
      epscr[j]=epscr[j]+(1.-exp(-dj))*screep[j][i];
    }
  }
  //        fprintf (Out,"\n\n v case t  prirustek dt faze 1");
  //        fprintf (Out," %e\t%e", ccTime, ddTime);
  //        fprintf (Out,"\n\n eps faze1");
  //        fprintf (Out," %le %le %le %le %le %le",epscr[0],epscr[1],epscr[2],epscr[3],epscr[4],epscr[5]);
  //        fprintf (Out,"\n\n tau faze1");
  //        fprintf (Out," %le %le %le %le %le %le",screep[0][0],screep[1][0],screep[2][0],screep[3][0],screep[4][0],screep[5][0]);
  //  increment of shrink and temperature
  
  //  jen kvuli debuggovani
  //desht=0.0;

  epscr[0]=epscr[0]+desht;
  epscr[1]=epscr[1]+desht;
  if(ss==planestrain) epscr[2]=epscr[2]+desht;
  else if(ss==axisymm) epscr[2]=epscr[2]+desht;
  else if(ss==spacestress) epscr[2]=epscr[2]+desht;

  //	fprintf (Out,"\n\n desht prir faze 1"); fprintf (Out," %le",desht);
  //	fprintf (Out,"\n h_s prir faze 1"); fprintf (Out," %le %le",h_s, h_slast);
  //	fprintf (Out,"\n temp_s prir faze 1"); fprintf (Out," %le %le",temp_s, temp_slast);
  
  //		fprintf (Out,"\n\n phase 1 eps1, eps2 prir");
  //		fprintf (Out," %e %e",epscr[0], epscr[1]);
  
  Mm->stiff_deps_creep (ipp, epscr, nc,nc);  // save sig
}



void creepb::phase2 (long ipp)
{
//  new total strain and stress
  long i,j,ii;
  double qq,dj, endTime;
  vector epscr(nc),deps(nc),depsshtem(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);
  
  endTime=ccTime+ddTime;
  //  fprintf (Out,"\n\n phase 2 v case t   prirustek dt");
  //  fprintf (Out," %e\t%e", ccTime, ddTime);
  // arr. strain = is total strain (including e from screep) from t=0
  for (i=0; i<nc; i++){
    deps[i]=Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[i];
  } //increment deps
  qq=2.0/3.0;
  for (i=0;i<nRetTime;i++){
    ii=nc*(i+1);
    dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
    for (j=0; j<nc; j++){
      screep[j][i]=Mm->ip[ipp].eqother[ii+j];
      deps[j]=deps[j]-(1.-exp(-dj))*screep[j][i];     // deps=deps-e(screep)
    }
  }
  get_desht (desht, ccTime, endTime);
  //  fprintf (Out,"\n\n desht prir faze 1"); fprintf (Out," %le",desht);
  
  //  debuggovani
  //desht=0.0;

  deps[0]=deps[0]-desht;
  deps[1]=deps[1]-desht;
  if(ss==planestrain) deps[2]=deps[2]-desht;
  else if(ss==axisymm) deps[2]=deps[2]-desht;
  else if(ss==spacestress) deps[2]=deps[2]-desht;
  //  increment of stress components
  //        fprintf (Out,"\n\n eps celkove,   eps celkove minule,  delta eps");
  //        fprintf (Out,"\n %le %le %le %le %le %le",epscr[0],epscr[1],epscr[2],epscr[3],epscr[4],epscr[5]);
  //        fprintf (Out,"\n %le %le %le ",Mm->ip[ipp].eqother[0],Mm->ip[ipp].eqother[1],Mm->ip[ipp].eqother[2]);
  //        fprintf (Out,"\n %le %le %le ",Mm->ip[ipp].eqother[3],Mm->ip[ipp].eqother[4],Mm->ip[ipp].eqother[5]);
  //        fprintf (Out,"\n\n delta eps");
  //        fprintf (Out,"\n %e %e %e",deps[0],deps[1],deps[2]);
  //  stiffness matrix of material
  Mm->matstiff(d,ipp);  
  e0=Mm->eliso[imat].e;
  mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  mxv (d,deps,sig);
  // save total strain (from t=0)
  for (i=0; i<nc; i++){Mm->ip[ipp].eqother[i]=Mm->ip[ipp].strain[i];}// it is last strain in arr. strain
  //  time history
  seps_time (screep,sig);
  for (i=0; i<nRetTime; i++){
    ii=nc*(i+1);
    for (j=0; j<nc; j++){
      Mm->ip[ipp].eqother[ii+j]=screep[j][i];
    }
  }
  //       fprintf (Out,"\n\n tau faze2");
  //       fprintf (Out," %le %le %le %le %le %le",screep[0][0],screep[1][0],screep[2][0],screep[3][0],screep[4][0],screep[5][0]);
  //  Stress data storage first position, Hide (strain) seconde position
  //    fprintf (Out,"\n\n *****prir sigma phase 2");
  //    fprintf (Out,"\n %e %e %e ",sig[0],sig[1],sig[2]);
  //    fprintf (Out,"\n %e %e %e %e %e %e",sig[0],sig[1],sig[2],sig[3],sig[4],sig[5]);
  Mm->ip[ipp].eqother[nc*(nRetTime+1)+0]+=desht;
  Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=h_s;
  Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=temp_s;
}

void creepb::get_h (long ipp)
{
  h_slast=1.0;
  if(type_h !=1 ){h_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]; }
}


void creepb::get_temp (long ipp)
{
  temp_slast=283.0;   //10 C
  if(type_temp != 1 ){temp_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]; }
}



/**
   function apdates material parameters
   
   @param E - function of creep

   10.10.2002
*/
void creepb::matstiff (matrix &d, long ipp)
//  function returns elastic stiffness matrix