slopesol.cpp

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

            Mm->updateipval();
          break; 
	}
      }

      modif=0;
      if (stop==0){
	//  modification of the arc length
	dl/=2.0;
	if (dl<dlmin){
          dl=dlmin;  
          if (nodiv_dl)     
            break; 
          else
            nodiv_dl = 1;
        }
	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
        copyv (r, ra, n);
	//  restoring of lambda parameter
	lambda=blambda;
      }
    }
  }
  
  // ------------------------------------
  //  finish of main iteration loop  ----
  // ------------------------------------
  
  if (Mp->nlman.hdbackupal==1)
    arclsave (i,n,lambda,dl,ra,fp);

  // termit
  if (Mp->adaptivity)
    if (i < ni)
      arclsave_adapt (i,lambda,0,dl);
    else
      if (Mp->nlman.hdbackupal%10 == 1)
	arclsave (i,n,lambda,dl,ra,fp);
  
  if ((Mp->adaptflag & 16) && (Mp->nlman.nial==100 || Mp->nlman.nial==1))
    Ada->run (0,0,0);
  
  
  delete [] r;		    
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
  return (i);
}



void newton_raphsonrv1 (long lcid, double rlambda, double rlerr)
  //  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;
  double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
  double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
  long back_incr = 1;
  long check_rv = 1;
  double binc;


  //  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;

  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;
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  print_init(-1, "wt");
  print_step(lcid, 0, 0.0, f);
  print_flush();
  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;
    if (back_incr)
      binc = dlambda;
   
    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++){
/*      if (check_rv)
      {
        f[j]=lambda*(fp[j]+fc[j]);
        fb[j]=dlambda*(fp[j]+fc[j]);
      }
      else
      {
        f[j]=rlambda*fc[j]+lambda*fp[j];
        fb[j]=dlambda*fp[j];
      }*/
      if (check_rv)
      {
        f[j]=lambda*(fc[j]);
        fb[j]=dlambda*(fc[j]);
      }
      else
      {
        f[j]=rlambda*fc[j]+(lambda-rlambda)*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++){
/*      if (check_rv)
        fb[j]=f[j]+dlambda*(fp[j]+fc[j]);
      else
        fb[j]=f[j]+dlambda*fp[j];*/
      if (check_rv)
        fb[j]=f[j]+dlambda*(fc[j]);
      else
        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;
    
    if (norf<err){
      lambda+=dlambda;
      if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
      {
        if (back_incr)
	{
          binc = dlambda;
          back_incr = 0;
	}
        //  modification of the newton-rhapson
        dlambda/=2.0;
        if (dlambda<=dlambdamin)
        {
          if (dlambda == dlambdamin)
            break;
          dlambda=dlambdamin;
        }
        if (Mespr==1)
          fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
        //  restoring of left hand side vector
        for (j=0;j<n;j++)
          r[j]=rb[j];
        //  restoring of lambda parameter
        lambda=blambda;
      }
      if (fabs(lambda-rlambda) <= rlerr)
      {
        dlambda = binc;
        check_rv = 0;
      }
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
      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++){
/*        if (check_rv)
          fb[k]=f[k]+dlambda*(fp[k]+fc[k])-fi[k];
        else
          fb[k]=f[k]+dlambda*fp[k]-fi[k];*/
        if (check_rv)
          fb[k]=f[k]+dlambda*(fc[k])-fi[k];
        else
          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);
      
      if (norf<err)  break;
    }
    
    if (j==ini || norf>err){
      if (dlambda == dlambdamin)
      {
	fprintf (stdout,"\n increment of lambda cannot be decreased, iteration will be stopped\n");
	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;
      if (check_rv && ((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
      {
        if (back_incr)
	{
          binc = dlambda;
          back_incr = 0;
	}
        //  modification of the newton-rhapson
        dlambda/=2.0;
        if (dlambda<=dlambdamin)
        {
          if (dlambda == dlambdamin)
            break;
          dlambda=dlambdamin;
        }
        if (Mespr==1)
          fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
        //  restoring of left hand side vector
        for (j=0;j<n;j++)
          r[j]=rb[j];
        //  restoring of lambda parameter
        lambda=blambda;
      }
      if (fabs(lambda-rlambda) <= rlerr)
      {
        if (back_incr == 0)
          dlambda = binc;
        check_rv = 0;
      }
      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();
}