mtsolver.cpp

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

#include <string.h>
#include "mtsolver.h"
#include "global.h"
#include "globmat.h"
#include "mechprint.h"
#include "matrix.h"
#include "vector.h"
#include "gmatrix.h"
#include "gtopology.h"
#include "dloadcase.h"



void solve_time_dep_prob ()
{
  long i;
  
  //for (i=0;i<Lsrs->nlc;i++){
  i=0;
  visco_solver (i);
  
  //print_stresses_old (Out,i);
  //    print_other (Out,i);
  
}

/**
   function solves time dependent problem (viscoplasticity)
   
   @param lcid - load case id
   
   27.10.2001
*/
void visco_solver (long lcid)
{
  long i,j,k,n, ini;
  double end_time,zero,*f,*fi,*fb,*fp,*r,*dr;
  double norf, norfa, err;
  char fname[1020];
  char *path;
  char *name;
  char *suffix;
  FILE *aux;

  Mp->filename_decomposition (Outdm->outfn,path,name,suffix);
  sprintf (fname,"%s%s.bac",path, name);

  delete [] path; delete [] name; delete [] suffix;

  if (Mp->timecon.nit > 0)
    aux = fopen (fname,"w");


  ini = 100;
  err = 1.0e-5;  
  //  number of rows of the matrix
  n = Ndofm;
  //  computer zero
  zero = Mp->zero;

  fp = new double [n];
  fb = new double [n];
  fi = new double [n];
  dr = new double [n];
  memset (fp,0,n*sizeof(double));
  memset (fb,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (dr,0,n*sizeof(double));
  
  //  initialization phase
  r = Lsrs->give_lhs (lcid);
  f = Lsrs->give_rhs (lcid);
  
  //  initial time
  Mp->time=Mp->timecon.starttime ();
  //  end time
  end_time = Mp->timecon.endtime ();
  
  //  allocation of backup arrays for state variables in integration points
  //Mm->alloc_backup ();
  
  //  provizorni umisteni, bude v mefelinit
  Gtm->comp_domain_sizes ();
  Mb->alloc_sumcomp ();
  //  konec provizoria

  //  initiation of mechanical material models
  initmaterialmodels();

  nullv (fp,n);
  i=0;
  
  print_init(-1, "wt");
  print_step(lcid, i, Mp->time, f);
  print_flush();
  //print_close();
  
  // ***************************
  //   main iteration loop  ****
  // ***************************
  do{

    //  indicator of strain computation
    //  it means, strains at integration points have not been computed
    Mp->strainstate=0;
    //  indicator of stress computation
    //  it means, stresses at integration points have not been computed
    Mp->stressstate=0;
    //  indicator of computation of other array
    //  it means, stresses at integration points have not been computed
    Mp->otherstate=0;

    
    Mm->updateipval();

    //if (Mespr==1)
    fprintf (stdout,"\n time %e",Mp->time);

    //    fprintf (Out,"\n\n\n\n time %e",Mp->time);
    
    //  prescribed nodal forces
    mefel_right_hand_side (lcid,f);

    // initialization of material models, initial temperatures are
    // computed during right hand side assembling -> right hand side
    // have to be reassembled due to correct temprstrains
    /* if (i==0)
       {
       initmaterialmodels();
       mefel_right_hand_side (lcid,f);
       }
    */
    Mb->comp_sum (f);
    
    
    //  unbalanced forces
    Mp->phase=1;
    //Mm->est=eigstrain;
    internal_forces (lcid,fi);
    
    //  right hand side assembling
    for (j=0;j<n;j++){
      fb[j]=f[j]-fp[j]+fi[j];
      //fb[j]=f[j]-fp[j];
      fp[j]=f[j];
    }
    
    //  assembling of tangent stiffness matrix
    stiffness_matrix (lcid);

    //  solution of K(r).v=F
    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);
    
    //  new displacement vector
    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    
    
    Mm->tempstr_eigstr();
    //  computation of internal forces
    Mp->phase=2;
    internal_forces (lcid,fi);

    norfa=ss(f,f,n);
    for (j=0;j<n;j++){
      fb[j] = f[j] - fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    if (norfa > 0.0)
      norf /= norfa;
    
    if (norf<err){
      // computation of values required for output    
      compute_req_val (lcid);
      //  time increment
      Mp->time = Mp->timecon.newtime ();
      i++;
      if (Mp->timecon.isitimptime ()==1){
        mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
      }
      //print_init(i, "wt");
      print_step(lcid, i, Mp->time, f);
      print_flush();
      
    }
    else{
      // iteration of unbalanced forces caused by time independent models
      //  internal iteration loop
      for (j=0;j<ini;j++)
	{
	  //  solution of K(r).v=F
	  //Smat->solve_system (Gtm,dr,fb);
	  Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);
	  
	  for (k=0;k<n;k++)
	    r[k]+=dr[k];
	  
	  //  computation of internal forces
	  internal_forces (lcid,fi);
	  
	  //  vector of unbalanced forces
	  for (k=0;k<n;k++)
	    fb[k]=f[k]-fi[k];
	  
	  //  norm of vector of unbalanced forces
	  norf=ss(fb,fb,n);
	  norf /= norfa;
	  
	  //if (Mespr==1)
	  fprintf (stdout,"\n Time increment %lf   inner loop j %ld     norf=%e",Mp->time,j,norf);
	  if (norf<err)  break;
	}
      if (j==ini || norf>err)
	{
	  //if (Mespr==1)
	  fprintf (stdout,"\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
	  
	  //backup of norf in case of unbalanced internal forces
	  err = norf;
	  
	  // equilibrium state was not reached: added by TKr
	  compute_req_val (lcid);
	  //  time increment
	  Mp->time = Mp->timecon.newtime ();
	  i++;
	  if (Mp->timecon.isitimptime ()==1)
	    mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
	  print_step(lcid, i, Mp->time, f);
	  print_flush();
	}
      else
	{
	  // equilibrium state was reached
	  compute_req_val (lcid);
	  //  time increment
	  Mp->time = Mp->timecon.newtime ();
	  i++;
	  if (Mp->timecon.isitimptime ()==1)
	    mtsolver_save (aux,r,fp,i,Mp->time,n,0,0);
	  print_step(lcid, i, Mp->time, f);
	  print_flush();
	}
    }
  }while(Mp->time<end_time);
  
  print_close ();
  
  if (Mp->timecon.nit > 0)
    fclose (aux);

  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f, delete [] fp;
}

/**
   function executes one step of time dependent mechanical algorithm
   (visco-elasticity, visco-plasticity, etc.)
   
   @param lcid - load case id
   @param f - %vector of prescribed nodal forces
   @param fi - %vector of computed nodal forces
   @param fp - %vector of previous nodal forces
   @param dr - %vector of displacement increments
   @param r - %vector of total displacements
   
   JK, 19.9.2004
*/
void one_step (long lcid,long n,double *f,double *fi,double *fp,double *fb,double *dr,double *r)
{
  long i;
  
  //  evaluation of the prescribed nodal forces at actual time
  mefel_right_hand_side (lcid,f);
  
  //  estimation of total forces in each direction
  //  used in creep
  Mb->comp_sum (f);
  
  //  nodal forces from attained stresses
  Mp->phase=1;
  internal_forces (lcid,fi);
  
  //  evaluation of the right hand side
  for (i=0;i<n;i++){
    fb[i]=f[i]-fp[i]+fi[i];
    fp[i]=f[i];
  }
  
  //  assembling of tangent stiffness matrix
  stiffness_matrix (lcid);
  
  //  solution of K(r).v=F  => increments of displacements
  //Smat->solve_system (Gtm,dr,fb);
  Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);
  
  //  new displacement vector
  for (i=0;i<n;i++){
    r[i]+=dr[i];
  }
  
  //  update of values with respect to attained state
  Mp->phase=2;
  internal_forces (lcid,fi);
  
}


void aux_nonlintime_print (FILE *aux,double *r,double l)
{
/*  long i,j,k,m,n;
  
  n=Mp->npun;
  
  fprintf (aux,"%e",l);
  for (i=0;i<n;i++){
    //  node number
    j=Mp->requn[2*i+0];
    //  unknown number
    k=Mp->requn[2*i+1];
    //  code number
    m=Mt->give_dof (j,k)-1;
    fprintf (aux,"  %e",r[m]);
  }
  fprintf (aux,"\n");
  fflush (aux);*/
}

void mtsolver_save (FILE *aux,double *r,double *fp,long ni,double time,long n,long ido1,long ncompo)
{
  long i;
  
  fprintf (aux,"\n\n\n\n");

  //  attained time
  fprintf (aux,"%le\n",time);
  //  attained number of steps
  fprintf (aux,"%ld\n",ni);
  //  number of degrees of freedom
  fprintf (aux,"%ld\n",n);

  //  left hand side vector
  for (i=0;i<n;i++){
    fprintf (aux,"%e\n",r[i]);
  }
  
  fprintf (aux,"\n");
  
  //  vector of the right hand side from previous step
  for (i=0;i<n;i++){
    fprintf (aux,"%e\n",fp[i]);
  }
  
  //  number of integration points in the problem
  fprintf (aux,"\n\n %ld\n",Mm->tnip);

  //  data on integration points
  Mm->save_intpoints (aux,ido1,ncompo);
  
  fflush (aux);
}

void mtsolver_restore (FILE *aux,double *r,double *fp,long &ni,double &time,long &n,long ido1,long ncompo)
{
  long i;
  
  //  attained time
  fscanf (aux,"%le",&time);
  //  attained number of steps
  fscanf (aux,"%ld\n",&ni);
  //  number of degrees of freedom
  fscanf (aux,"%ld\n",&n);

  if (n!=Ndofm){
    fprintf(stderr, "\n\nNumber of DOFs in backup file and in the actual problem is not same\n"); 
    fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
    fprintf(stderr, "Program will be terminated\n");
    abort();
  }

  //  left hand side vector
  for (i=0;i<n;i++){
    fscanf (aux,"%le",r+i);
  }
  //  vector of the right hand side from previous step
  for (i=0;i<n;i++){
    fscanf (aux,"%le",fp+i);
  }
  
  //  number of integration points in the problem
  fscanf (aux,"%ld",&i);
  
  if (i!=Mm->tnip){
    fprintf(stderr, "\n\nNumber of integration points in backup file and in the actual problem is not same\n"); 
    fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
    fprintf(stderr, "Program will be terminated\n");
    abort();
  }

  //  data on integration points
  Mm->restore_intpoints (aux,ido1,ncompo);

}