nssolver.cpp

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

  delete [] ddr;
  
  return (i);
}

long arclength15 (long lcid)
{
  long j,k,n,ini,stop,li;
  double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,bdlambda,dlambda,ddlambda,norf,norfp,norv,zero,ierr;
  double ss1,ss2,ss3,ss4,ss5;
  double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
  //  FILE *gr;
  //  char file[255];
  // double d;


  //  gr = fopen (Mp->auxfile,"a");
  
  n = Ndofm;                   
  ini = Mp->nlman.niilal;           
  zero = Mp->zero;          
  ierr = Mp->nlman.erral;           
  dl = Mp->nlman.dlal;              
  dlmax = Mp->nlman.dlmaxal;        
  dlmin = Mp->nlman.dlminal;        
  psi = Mp->nlman.psial;          
  
  r = new double [n];     ddr = new double [n];   u = new double [n];     v = new double [n];   
  f = new double [n];     fa = new double [n];    fi = new double [n];  
  
  //  Mm->alloc_backup (); 
  ra = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  
  arclopen_adapt (li,lambda,dlambda,dl);
  
  //  aux_mech_nonlin_print (gr,ra,lambda);
  // *******************************************************************
  //  backup of left hand side vector
  fprintf (stdout,"\n arc-length15: lambda %e  dl %e",lambda,dl);
  
  stiffness_matrix (lcid);
  for (j=0;j<n;j++)    f[j]=fp[j];

  //Smat->solve_system (Gtm,v,f);
  Mp->ssle.solve_system (Gtm,Smat,v,f,Out);
  norv = displincr (v,n);
  norfp = ss (fp,fp,n);
    
  if (0 && !Mm->plast){
    for (j=0;j<n;j++)
      ra[j]=(lambda+dlambda)*v[j];
    
    //    aux_mech_nonlin_print (gr,ra,lambda);
  }
  else{  
    stop = 0;

    for (j=0;j<n;j++){      r[j]=ra[j];    }
    blambda=lambda;
    bdlambda=dlambda;
    //    Mm->backup_state_var ();
    
    for (j=0;j<n;j++){
      ddr[j]=dlambda*v[j];
      fa[j]=fc[j]+(lambda+dlambda)*fp[j];
    }
    ddlambda=dlambda;
    
    internal_forces (lcid,fi);
    
    for (k=0;k<n;k++)
      f[k]=fa[k]-fi[k];
    
    norf=ss(f,f,n);
    norf=sqrt(norf);    if (Mespr==1)  fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
    
    if (norf<ierr){
      lambda+=dlambda;
      Mm->updateipval ();
      stop = 1;
    }

    while (!stop){


      
      // ****************************
      //  inner iteration loop  ****   arclength15
      // ****************************
      for (j=0;j<10*ini;j++){
	
      //      aux_mech_nonlin_print (gr,ra,lambda);
      //  back substitution
	//Smat->solve_system (Gtm,u,f);
	Mp->ssle.solve_system (Gtm,Smat,u,f,Out);

	for (k=0;k<n;k++){
	  f[k]=ddr[k];
	  ddr[k]+=u[k];
	}
	
	//  coefficient of quadratic equation
	//a2=norv+psi*psi*norfp;
	//a1=2.0*ss(ddr,v,n)+2.0*ddlambda*psi*psi*norfp;
	//a0=ss(ddr,ddr,n)+ddlambda*ddlambda*psi*psi*norfp-dl*dl;
	
	quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
	
	
	//  solution of quadratic equation
	/*
	if (fabs(a2)<zero){
	  if (fabs(a1)<zero){
	    if (fabs(a0)>zero)  fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
	  else  fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
	  }
	  else{
	    dlambda=(0.0-a0)/a1;
	  }
	}
	else{
	  d=a1*a1-4.0*a2*a0;
	  
	  if (d<0.0){
	    fprintf (stderr,"\n negative discriminant in function arclength");
	    break;
	  }
	  
	  
	  l1=(0.0-a1+sqrt(d))/2.0/a2;
	  l2=(0.0-a1-sqrt(d))/2.0/a2;
	  
	  if (l1>l2)  dlambda=l1;
	  else  dlambda=l2;
	  \\if (fabs(l1)<fabs(l2))  dlambda=l1;
	  \\else  dlambda=l2;
	}
	*/
	//  zacatek novinky   arclength15
        /*
	ss1=0.0;  ss2=0.0;
	ss3=0.0;  ss4=0.0;  ss5=0.0;
	for (k=0;k<n;k++){
	  ss1+=(ddr[k]+l1*v[k])*f[k];
	  ss2+=(ddr[k]+l2*v[k])*f[k];
	  ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
	  ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
	  ss5+=f[k]*f[k];
	}
	*/
	if (ss1/ss3/ss5>ss2/ss4/ss5)  dlambda=l1;
	else          dlambda=l2;
	//  konec novinky
	
        dlambda = 0;
	
	for (k=0;k<n;k++){
	  ddr[k]+=dlambda*v[k];
	  ra[k]+=u[k]+dlambda*v[k];
	  fa[k]+=dlambda*fp[k];
	}
	ddlambda+=dlambda;
	
	internal_forces (lcid,fi);
	
	for (k=0;k<n;k++)
	  f[k]=fa[k]-fi[k];
	
	norf=ss(f,f,n);   norf=sqrt(norf);
	
	if (Mespr==1){ 
          fprintf (stdout,"\n increment  %ld     inner loop  %ld    norf=%e",li,j,norf); 
//          fprintf (Out,"\n increment  %ld     inner loop  %ld    norf=%e",li,j,norf);
        }
	
	if (norf<0.10*ierr){
	  lambda+=ddlambda;
	  //	  aux_mech_nonlin_print (gr,ra,lambda);
	  Mm->updateipval ();
	  stop=1;  break;
	}
      }   // arclength15
      
      if (stop==0){
	//  modification of the arc length
	dl/=2.0;
	if (dl<dlmin){  dl=dlmin;  break; }
	
	if (Mespr==1){
	  fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e    norf=%e",dl,norf);
//	  fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e    norf=%e",dl,norf);
	}
	//  restoring of left hand side vector
	for (j=0;j<n;j++){
	  ra[j]=r[j];
	}
	//  restoring of lambda parameter
	lambda=blambda;
	ddlambda=bdlambda;
	//  restoring of state variables
	//	Mm->restore_state_var ();

	//****************
	stop=1;
	//****************

      }
    }
  }
  
  // ------------------------------------
  //  finish of main iteration loop  ----  arclength15
  // ------------------------------------
  
  arclsave_adapt (li,lambda,0.0,dl);
  
  delete [] r;
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
  //  fclose (gr); 
  
  return (li);
} //arclength15

void newton_raphson (long lcid)
  //  function solves system of nonlinear algebraic
  //  equations by Newton-Raphson method
  //  solved system does not contain time variable
  //
  //  16.8.2001
{
  long i,j,k,n,ni,ini,modif;
  double lambda,blambda,dlambda,dlambdamin,dlambdamax,zero,err,norf,norfa;
  double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nlman.ninr;
  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilnr;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  err = Mp->nlman.errnr;
  //  increment size
  dlambda=Mp->nlman.incrnr;
  //  minimum increment size
  dlambdamin=Mp->nlman.minincrnr;
  //  maximum increment size
  dlambdamax=Mp->nlman.maxincrnr;

  rb = new double [n];
  f  = new double [n];
  fb = new double [n];
  fi = new double [n];
  dr = new double [n];
  memset (rb,0,n*sizeof(double));
  memset (f,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (dr,0,n*sizeof(double));
  
  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  lambda=0.0;
  
  modif=0;
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  print_init(-1, "wt");
  print_step(lcid, 0, 0.0, f);
  print_flush();


  for (i=0;i<ni;i++){
    
    //  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;


    //  backup of left hand side vector
    for (j=0;j<n;j++){
      rb[j]=r[j];
    }
    //  backup of reached lambda parameter
    blambda=lambda;
   
    fprintf (stdout,"\n increment %ld   lambda %e,    dlambda %e",i,lambda,dlambda);
    
    //  vector of maximum load and vector of load increment
    for (j=0;j<n;j++){
      f[j]=fc[j]+lambda*fp[j];
      fb[j]=dlambda*fp[j];
    }
    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == 0))
      stiffness_matrix (lcid);
    
    //  solution of K(r).v=F
    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    
    //  computation of internal forces
    internal_forces (lcid,fi);
    
    //  vector of unbalanced forces
    for (j=0;j<n;j++){
      fb[j]=f[j]+dlambda*fp[j];
    }
    norfa=ss(fb,fb,n);
    for (j=0;j<n;j++){
      fb[j] -= fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    norf /= norfa;
    
//    fprintf (Out,"\n norf=%e",norf);

    if (norf<err){

      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();

      
      modif++;
      if (modif>1){
	dlambda*=2.0;
	if (dlambda>dlambdamax)
	  dlambda=dlambdamax;

	if (Mespr==1){
	  fprintf (stdout,"\n increment must be modified (dlambda increase) dlambda=%e",dlambda);
	}
	modif=0;
      }
      
      continue;
    }
    
    //  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]+dlambda*fp[k]-fi[k];
      }
      
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,n);
      norf /= norfa;
      
      fprintf (stdout,"\n increment %ld   inner loop j %ld     norf=%e  dlambda %e",i,j,norf,dlambda);
//      fprintf (Out,"\n j=%ld     norf=%e",j,norf);
      
      if (norf<err)  break;
    }
    
    modif=0;

    if (j==ini || norf>err){
      dlambda/=2.0;
      if (dlambda<dlambdamin){
	dlambda=dlambdamin;
	break;
      }

      if (Mespr==1){
	fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
//	fprintf (Out,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
      }
      
      for (j=0;j<n;j++){
      	r[j]=rb[j];
      }
      lambda=blambda;
    }
    else{
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
    }
    
  }
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f;
  print_close();
}

void displ_control (long lcid)
  //  function solves system of nonlinear algebraic
  //  equations by Newton-Raphson method
  //  solved system does not contain time variable
  //  displacement contorl is assumed
  //
  //  2.11.2005
  //  TKo.
{
  long i,j,k,n,ni,ini,incs;
  double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa, dlambdamax;
  double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp, *fir2;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nlman.ninr;
  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilnr;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  err = Mp->nlman.errnr;
  //  increment size
  dlambda=Mp->nlman.incrnr;
  //  minimum increment size
  dlambdamin=Mp->nlman.minincrnr;
  // maximum increment size
  dlambdamax=Mp->nlman.incrnr;
 
  rb = new double [n];
  f  = new double [n];
  fb = new double [n];
  fi = new double [n];
  fir2 = new double[n];
  dr = new double [n];
  memset (rb,0,n*sizeof(double));
  memset (f,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (dr,0,n*sizeof(double));
  
  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  lambda=0.0;
  Mp->lambda = lambda;
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  print_init(-1, "wt");
  print_step(lcid, 0, 0.0, f);
  print_flush();
  incs = 0;
  for (i=0;i<ni;i++){

    //  backup of left hand side vector
    for (j=0;j<n;j++){
      rb[j]=r[j];
    }
    //  backup of reached lambda parameter
    blambda=lambda;
    Mp->lambda = lambda;
   
    fprintf (stdout,"\n increment %ld   lambda %e,    dlambda %e",i,lambda,dlambda);
    
    //  vector of maximum load and vector of load increment
    Mp->lambda += dlambda;
    internal_forces(lcid, fir2);

    for (j=0;j<n;j++){
      f[j]=fc[j]+lambda*fp[j];
      fb[j]= f[j] + dlambda*fp[j] - fir2[j];
    }
    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == 0))
      stiffness_matrix (lcid);
    
    //  solution of K(r).v=F
    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    
    //  computation of internal forces
    internal_forces (lcid,fi);
    
    //  vector of unbalanced forces
    for (j=0;j<n;j++){
      fb[j]=f[j]+dlambda*fp[j];
    }
    norfa=ss(fb,fb,n);
    for (j=0;j<n;j++){
      fb[j] -= fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    if (norfa != 0.0)
      norf /= norfa;
    
    if (norf<err){
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
      continue;
    }
    
    stiffness_matrix (lcid);
    //  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]+dlambda*fp[k]-fi[k];
      }
      
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,n);
      if (norfa != 0.0)
        norf /= norfa;
      
      fprintf (stdout,"\n increment %ld   inner loop j %ld     norf=%e  dlambda %e",i,j,norf,dlambda);
      
      if (norf<err)  
        break;
    }
    
    if (j==ini || norf>err)
    {
      if (dlambda == dlambdamin)
      {
        if (Mespr==1)  fprintf (stdout,"\n increment cannot be decreased dlambda=%e    norf=%e",dlambda,norf);
        break;
      }
      dlambda/=2.0;
      if (dlambda<dlambdamin)  dlambda=dlambdamin;
      if (Mespr==1)  fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
      
      for (j=0;j<n;j++)
      	r[j]=rb[j];

      lambda=blambda;     
    }
    else
    {
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
    }
  }
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f;
  print_close();
}