creepbbeam.cpp

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

#include "creepbbeam.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"

creepbbeam::creepbbeam (void)
{
  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[6]=3000.0;
  type_h=1;
  type_temp=1;
  timeMax=Mp->timecon.endtime ();
  h_s = 0.4;   // end humidity
  temp_s=20;    // temperature
  esht=0.0;
  nc=1;
  k_s = 1.0;    //shape factor
  r_s = 0.8;
  k_d = 0.12;  //effective cross area 2V/S
  timemat=0.0;
}

creepbbeam::~creepbbeam (void)
{
}

void creepbbeam::creepinit (long mie)
{
}

void creepbbeam::read (FILE *in)
{
  fscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
	          &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&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 creepbbeam::approx (vector &areacoord,vector &nodval)
{
  double f;
  scprd (areacoord,nodval,f);
  return f;
}

/**
   function inversion sym. matrix

   @param a - matrix

*/
void creepbbeam::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];
	   }
	}
}

/**
   function returns eps from history
   
   @param screep - vector of history
   @param epsscr - vector deformation of history
   
   10.10.2002
*/
void creepbbeam::nlstresses (long ipp)
{
  //long nc=Mm->ip[ipp].ncompstr;
  //timeMax=Mp->end_time;
  //ddTime=Mp->timefun.getval(Mp->time);
  ddTime=Mp->timecon.actualforwtimeincr ();
  //ccTime=Mp->time;
  ccTime=Mp->timecon.actualtime ();
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  //if(type_h!=1) {h_s=Mm->moist[ipp];}   //actual moist
  //if(type_temp!=1) {temp=Mm->tempr[ipp];}  //ctual tempr
  get_h (ipp);     //initial  moist
  get_temp (ipp);     //initial  tempr
  if (Mp->phase==1){
	phase1(ipp);
  }

  if (Mp->phase==2){
	phase2(ipp);
  }

}

void creepbbeam::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;
  vector epscr(nc),deps(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);

    Mm->matstiff(d,ipp);
    e0=Mm->eliso[imat].e;
    mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
	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((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);
	    for (j=0; j<nc; j++){
		 screep[j][i]=Mm->ip[ipp].other[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]);
//  shrink and temp
	epscr[0]=epscr[0]+deps0;
	epscr[1]=epscr[1]+deps0;
	if(epscr.n==6) epscr[2]=epscr[2]+deps0;
//		fprintf (Out,"\n\n time, esht prir faze 1");
  
//		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 creepbbeam::phase2 (long ipp)
{
//  new total strain and stress
  long i,j,ii;
  double qq,dj;
  vector epscr(nc),deps(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);

		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].other[i];} //increment deps
	qq=2.0/3.0;
	for (i=0;i<nRetTime;i++){
		ii=nc*(i+1);
		dj=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);
	    for (j=0; j<nc; j++){
		 screep[j][i]=Mm->ip[ipp].other[ii+j];
		 deps[j]=deps[j]-(1.-exp(-dj))*screep[j][i];     // deps=deps-e(screep)
		}
	}
//  shrink and temp
	deps[0]=deps[0]-deps0;
	deps[1]=deps[1]-deps0;
	if(deps.n==6) deps[2]=deps[2]-deps0;
//  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].other[0],Mm->ip[ipp].other[1],Mm->ip[ipp].other[2]);
//        fprintf (Out,"\n %le %le %le ",Mm->ip[ipp].other[3],Mm->ip[ipp].other[4],Mm->ip[ipp].other[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].other[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].other[ii+j]=screep[j][i];
		}
	}
//  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].other[nc*(nRetTime+1)+1]=esht;
	Mm->ip[ipp].other[nc*(nRetTime+1)+2]=h_s;
	Mm->ip[ipp].other[nc*(nRetTime+1)+3]=temp_s;
}


void creepbbeam::get_h (long ipp)
{
    h_slast=0.0;
	if(type_h != 1){h_slast=Mm->ip[ipp].other[nc*(nRetTime+1)+2]; }	
}


void creepbbeam::get_temp (long ipp)
{
    templast=0.0;
    if(type_temp != 1){templast=Mm->ip[ipp].other[nc*(nRetTime+1)+2]; }
}

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

   10.10.2002
*/
void creepbbeam::matstiff (matrix &d, long ipp)
//  function returns elastic stiffness matrix
//  d - elastic stiffness matrix
{
  long i,j,k;
  double jt,et,pomt,pomt1,qq,delYi,delYj;
  vector bpom(nRetTime);
  matrix apom(nRetTime,nRetTime);
//  matrix c(d.m,d.n);
  
  
  //  initial elastic stiffness matrix
  Mm->elmatstiff (d,ipp);
  //  initial elastic compliance matrix    Mm->elmatcompl (c,ipp);

  e0=Mm->eliso[imat].e;
  mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  
  //ddTime=Mp->timefun.getval(Mp->time);
  ddTime=Mp->timecon.actualforwtimeincr ();
  //ccTime=Mp->time;
  ccTime=Mp->timecon.actualtime ();
  pomt=ccTime+ddTime/2.;
  qq=2./3.;
  jt=0.0;
  //for (ii=0;ii<numMat;ii++){
  //  long ii;
  //     ii=imat;
  //   from concrete starts   tb 
  if(timemat<=pomt) {
    //        ert[nRetTime+1][imat]=0.;
    et=0.;
    for (i=0;i<nRetTime;i++){
      bpom[i]=0.;
      for (j=0;j<nRetTime;j++){
	apom[i][j]=0.;
      }
    }
    k=0;
    pomt1=pomt+0.0001;   
    while (pomt1<=timeMax*5.){
      //   creep function od casu pomt do casu pomt1
      b3_law(jt,esht, pomt, pomt1);
      //		fprintf (Out,"\n\n jt v case t   od t0 prirustek dt");
      //		fprintf (Out," %e\t%e\t%e\t%e",jt, pomt, ccTime, ddTime);
      for (i=0;i<nRetTime;i++){