cpsolver.cpp

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

#include <string.h>
#include "cpsolver.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"
#include "node.h"
#include "element.h"


/**
   Function solves non-linear MEFEL equations by Newton-Raphson method for time-dependent problems
   for growing structure
   
   @param lcid - load case id
   
   15.6.2007, TKr
*/
void solve_prob_constr_phases ()
{
  long lcid,ii,j,k,nm,mncd,mnce,ini;
  double dt,prev_timem,end_time;
  double *dr,*fb,*r,*f,*fp,*fi;
  double norf, norfa, err;  
  char fname[1020];
  char *path;
  char *name;
  char *suffix;
  FILE *aux;
  
  //  load case id must be equal to zero
  lcid=0;

  Mp->filename_decomposition (Outdm->outfn,path,name,suffix);
  Mp->path = path;
  Mp->filename=name;
  Mp->suffix=suffix;
  sprintf (fname,"%s%s.bac",path, name);
  
  aux = fopen (fname,"w");
  
  //  maximum number of iterations in inner loop
  ini = Mp->niilnr;
  //  required norm of vector of unbalanced forces
  err = Mp->errnr;
  //  number of mechanical degrees of freedom
  nm=Ndofm;
  
  /* **************************** */
  /*  vectors of mechanical part  */
  /* **************************** */
  fp = new double [nm];
  fb = new double [nm];
  fi = new double [nm];
  dr = new double [nm];
  
  //  initial values
  nullv (fp,nm);
  nullv (fb,nm);
  nullv (fi,nm);
  nullv (dr,nm);
  
  r = Lsrs->give_lhs (lcid);
  f = Lsrs->give_rhs (lcid);
  
  //  provizorni umisteni, bude v mefelinit
  Gtm->comp_domain_sizes ();
  Mb->alloc_sumcomp ();
  //  konec provizoria??!!

  //  initial time
  Mp->time = Mp->timecon.starttime ();
  prev_timem = Mp->time;
  Mp->timecon.take_values (&Mp->timecon);
  //  end time
  end_time = Mp->timecon.endtime ();

  //  initiation of mechanical material models
  initmaterialmodels();

  Gtm->nm = new long [Gtm->nn];
  
  // **************************
  //  main iteration loop  ***
  // **************************
  ii=0;
  print_init(-1, "wt");
  print_step(lcid,ii,Mp->time,f);
  print_flush();

  //saving values for the first iteration
  Mt->lhs_save (r,1);

  do{
    
    //if (Mespr==1)  
    fprintf (stdout,"\n\n ------------------------------------------------------------------------------");
    //if (Mespr==1)  
    fprintf (stdout,"\n\n provadi se MEFEL,  mefel time %e",Mp->time);
    //if (Mespr==1)  
    //fprintf (stdout,"\n time increment %lf",dt);
    //if (Mespr==1)  
    fprintf (stdout,"\n -----------------------------------------------------------------------------\n");
    
            
      //  searching for new elements
      //  only new elements are switched on
      //  former elements are temporarily switched off
      mnce = Gtm->search_newelem (Mp->time,prev_timem);
    
      if (mnce!=0){
	//  attained displacements have to be saved on new elements
	//  otherwise stresses will not be correct
	Mt->initial_displ ();
      }
      
      //  searching for new DOFs
      //  only new DOFs are switched on
      //  former DOFs are temporarily switched off
      mncd = Gtm->search_newdofs (Mp->time,prev_timem);
      
      //  marking of elements switched on
      Gtm->update_elem (Mp->time);
      
      //  marking of nodes switched on
      Gtm->update_nodes ();
      
      //  marking of DOFs switched on
      Gtm->update_dofs (Mp->time);
      
      //  definition of DOFs
      Ndofm = Gtm->codenum_generation (Out);
      nm=Ndofm;
      
      //  update of nodes and DOFs
      //if (mncd>0){
      if (Smat != NULL){
	delete Smat;
	Smat=NULL;
      }
      //}
      
      //  assembling of lhs array for new code numbers
      //if (mnce!=0 || mncd!=0){
      Mt->lhs_restore (r,nm);
      Mt->rhs_restore (fp,nm);
      //}
            
      Mm->updateipval();

      //fprintf (stdout,"\nright hand side\n\n");
      //  assembling of the right hand side (prescribed loading)
      mefel_right_hand_side (lcid,f);
      
      //  computation of sums of force components in particular directions
      //  it is used in special creep and consolidation models
      Mb->comp_sum (f);
      
      Mp->phase=1;
      //  computation of unbalanced forces which will be on the right hand side
      internal_forces (lcid,fi);
      
      for (j=0;j<nm;j++){
	fb[j]=f[j]-fp[j]+fi[j];
	fp[j]=f[j];
      }    
      
      //  assembling of tangent stiffness matrix
      //fprintf (stdout,"\nphase1 = stiffness_matrix\n\n");
      stiffness_matrix (lcid);
      
      //  solution of K(r).v=F
      Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

      //  new displacement vector
      for (j=0;j<nm;j++){
	r[j]+=dr[j];
      }
      
      Mm->tempstr_eigstr();
      Mp->phase=2;

      //  update of computed values
      internal_forces (lcid,fi);

      norfa=ss(f,f,nm);
      for (j=0;j<nm;j++){
	fb[j] = f[j] - fi[j];
      }
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,nm);
      if (norfa > 0.0)
	norf /= norfa;
      
      if (norf<err){

	//  backup of nodal values because code numbers may be changed
	Mt->lhs_save (r,1);
	Mt->rhs_save (f);

	// computation of values required for output    
	compute_req_val (lcid);
	//  time increment
	prev_timem=Mp->time;
	Mp->time = Mp->timecon.newtime ();
	Mp->timecon.take_values (&Mp->timecon);

	if (Mp->timecon.isitimptime ()==1){
	  mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0);
	}
	
	print_step(lcid, ii, Mp->time, f);
	print_flush();
	ii++;	
      }
      
      // iteration of unbalanced forces caused by time independent models
      //  internal iteration loop
      else{
	for (j=0;j<ini;j++)
	  {
	    
	    if(Mp->nrsolv == fullnewton){
	      //cleaning of stiffness matrix
	      if (Smat != NULL){
		delete Smat;
		Smat=NULL;
	      }
	      //assembling of tangent stiffness matrix for fullnewton
	      stiffness_matrix (lcid);
	    }
	    
	    //  solution of K(r).v=F
	    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

	    for (k=0;k<nm;k++)
	      r[k]+=dr[k];
	    
	    //  computation of internal forces
	    internal_forces (lcid,fi);
	    
	    //  vector of unbalanced forces
	    for (k=0;k<nm;k++)
	      fb[k]=f[k]-fi[k];
	    
	    //  norm of vector of unbalanced forces
	    norf=ss(fb,fb,nm);
	    if (norfa > 0.0)
	      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;

	    //  backup of nodal values because code numbers may be changed
	    Mt->lhs_save (r,1);
	    Mt->rhs_save (f);
	    
	    // equilibrium state was reached not
	    compute_req_val (lcid);
	    //  time increment
	    prev_timem = Mp->time;
	    Mp->time = Mp->timecon.newtime ();
	    Mp->timecon.take_values (&Mp->timecon);
	    
	    if (Mp->timecon.isitimptime ()==1)
	      mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0);
	    print_step(lcid, ii, Mp->time, f);
	    print_flush();
	    ii++;
	  }
	else
	  {
	    //  backup of nodal values because code numbers may be changed
	    Mt->lhs_save (r,1);
	    Mt->rhs_save (f);
	    
	    // equilibrium state was reached
	    compute_req_val (lcid);
	    //  time increment
	    prev_timem = Mp->time;
	    Mp->time = Mp->timecon.newtime ();
	    Mp->timecon.take_values (&Mp->timecon);

	    if (Mp->timecon.isitimptime ()==1)
	      mtsolver_save (aux,r,fp,ii,Mp->time,nm,0,0);

	    print_step(lcid, ii, Mp->time, f);
	    print_flush();
	    ii++;
	  }
      }
      
      if (Mespr==1){
	//cleaning once again for more free memmory
	if (Smat != NULL){
	  delete Smat;
	  Smat=NULL;
	}
      }

  }while(Mp->time<=end_time);
  
  print_close ();
  
  delete [] fi;
  delete [] fb;
  delete [] fp;
  delete [] dr;
}




/**
   function solves problem of construction phases
   
   only one load case may be defined, one load case may contain
   several effects which are activated by time functions
   variable Mb->nlc must be equal to 1
   
   6.2.2006, JK
*/
void solve_prob_constr_phases_old ()
{
  long i,j,n,ncd,nce,lcid,idcp,ini;
  double prev_time,end_time,zero,*f,*fi,*fb,*fp,*r,*dr;
  double err;  
  char fname[1020];
  char *path;
  char *name;
  char *suffix;
  FILE *aux;
  
  //  load case id must be equal to zero
  lcid=0;
  
  //  definition of all output and auxiliary files
  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 ();
  prev_time=Mp->time;
  //  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
  
  //Gtm->cn_clean ();
  
  nullv (fp,n);
  i=0;
  idcp = 0;
  
  //BEGIN OF GRAPHICS:
  //print_init(-1, "wt");
  //print_step(lcid, i, Mp->time, f);
  //print_flush();
  //print_close();
  //END OF GRAPHICS

  //FILE *kon;
  //kon = fopen ("kontrolavondr","w");
  //kon = fopen ("caskon","w");
  
  //  initiation of mechanical material models
  initmaterialmodels();
  
  
  Gtm->nm = new long [Gtm->nn];
  
  
  // ***************************
  //   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;

    
    fprintf (stdout,"\n time %e",Mp->time);

    //  searching for new elements
    //  only new elements are switched on
    //  former elements are temporarily switched off
    nce = Gtm->search_newelem (Mp->time,prev_time);
    
    if (nce!=0){
      //  attained displacements have to be saved on new elements
      //  otherwise stresses will not be correct
      Mt->initial_displ ();
    }
    
    //  searching for new DOFs
    //  only new DOFs are switched on
    //  former DOFs are temporarily switched off
    ncd = Gtm->search_newdofs (Mp->time,prev_time);

    /***************************************************/
    //////BEGIN OF GRAPHICS
    /*
    long graphics=0;
    if ((nce!=0 || ncd!=0) && graphics==1){
      //if (nce!=0 || ncd!=0){
      // ************************
      //  correction of graphic
      // ************************
      
      
      //fprintf (kon,"\n cas %le  grafika",Mp->time);
      //for (j=0;j<Gtm->nn;j++){
      //fprintf (kon,"\n uzel %6ld",j+1);
      //long ndofn=Gtm->gnodes[j].ndofn;
      //for (long k=0;k<ndofn;k++){
      //fprintf (kon," %ld",Gtm->gnodes[j].cn[k]);
      //}
      //}

      //  generation of new code numbers
      Ndofm = Gtm->gencodnum ();
      n=Ndofm;
      
      //  cleaning of the stiffness matrix
      if (Smat != NULL){
	delete Smat;
	Smat=NULL;
      }
      

      //fprintf (kon,"\n cas %le  grafika",Mp->time);
      //for (j=0;j<Gtm->nn;j++){
      //fprintf (kon,"\n uzel %6ld",j+1);
      //long ndofn=Gtm->gnodes[j].ndofn;
      //for (long k=0;k<ndofn;k++){
      //fprintf (kon," %ld",Gtm->gnodes[j].cn[k]);
      //}
      //}
      
      //  allocation and assembling of the stiffness matrix
      stiffness_matrix (lcid);
      
      //  nodal forces caused by attained displacements
      rhs_graphic (fb,lcid=0,n);
      
      //  solution of system of equations
      Smat->solve_system (dr,fb);
      
      //  storage of auxiliary displacements to nodes
      Mt->lhs_save (dr,2);
    }
    */
    //////END OF GRAPHICS
    /***************************************************/


    //  marking of elements switched on
    Gtm->update_elem (Mp->time);
    
    //  marking of nodes switched on
    Gtm->update_nodes ();
    
    //  marking of DOFs switched on
    Gtm->update_dofs (Mp->time);
    
    /*
    fprintf (Out,"\n\n kontrola pole leso v case %le\n",Mp->time);
    for (j=0;j<Mt->ne;j++){
      fprintf (Out,"\n leso %4ld    %ld",j+1,Gtm->leso[j]);
    }
    
    fprintf (Out,"\n\n kontrola pole cn v case %le\n",Mp->time);
    for (j=0;j<Mt->nn;j++){
      fprintf (Out,"\n node %4ld    %4ld %4ld",j+1,Gtm->gnodes[j].cn[0],Gtm->gnodes[j].cn[1]);
    }
    */

    //  generation of new code numbers
    Ndofm = Gtm->codenum_generation (Out);
    n=Ndofm;
    
    /*
    fprintf (Out,"\n\n kontrola pole cn v case %le po generovani kodovych cisel\n",Mp->time);
    for (j=0;j<Mt->nn;j++){
      fprintf (Out,"\n node %4ld    %4ld %4ld",j+1,Gtm->gnodes[j].cn[0],Gtm->gnodes[j].cn[1]);
    }
    fprintf (Out,"\n\n  NDOF %ld",n);
    fprintf (stdout,"  NDOF %ld",n);
    */
    
    /*
    fprintf (kon,"\n cas %le  vypocet",Mp->time);
    for (j=0;j<Gtm->nn;j++){
      fprintf (kon,"\n uzel %6ld",j+1);
      long ndofn=Gtm->gnodes[j].ndofn;
      for (long k=0;k<ndofn;k++){
	fprintf (kon," %ld",Gtm->gnodes[j].cn[k]);
      }
    }