epsolver.cpp

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

      fprintf(grout, "   %e\n", ra[3*kk]);
    fprintf(grout, "Tuhosti --------\n");
    for (kk = 15; kk < Mt->ne; kk++)
    {
      Spring->res_stiffness_matrix (kk,sm);
      fprintf(grout, "   %e\n", sm[0][0]);
    }
    fprintf(grout, "\n");
    fflush(grout);
    //  backup of left hand side vector
    for (j=0;j<n;j++)
      r[j]=ra[j];
    //  backup of reached lambda parameter
    blambda=lambda;
    //  backup of state variables
//    Mm->backup_state_var ();
//    aux_nonlin_print (gr,ra,lambda);

    fprintf (stdout,"\n arc-length: increment %ld   lambda %e  dl %e",i,lambda,dl);
//    if (Mespr==1)
//    {
//      if (Mp->adp==1)
//        Mp->ado->print (Out,i,lambda,0);
//    }

    //  assembling of tangent stiffness matrix
    stiffness_matrix (lcid);

    //  backup of the fp, in ldl_sky the right hand side will be destroyed
    for (j=0;j<n;j++)
      f[j]=fp[j];
    //  solution of K(r).v=F
    Smat->solve_system (v,f);

    //  generalized norm of displacement increments
    norv = displincr (v,n);

    dlambda = dl/sqrt(norv+psi*psi*norfp);

    for (j=0;j<n;j++)
    {
      ddr[j]=dlambda*v[j];
      ra[j]+=ddr[j];
      fa[j]=fc[j]+(lambda+dlambda)*fp[j];
    }

    ddlambda=dlambda;

    //  computation of internal forces
    internal_forces (lcid,fi);

    for (k=0;k<n;k++)
      f[k]=fa[k]-fi[k];

    norf=ss(f,f,n);
    norf=sqrt(norf);
    norfa = ss(fa,fa,n);
    norfa = sqrt(norfa);
    norf /= norfa;

    if (Mespr==1)  fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);

    if (norf<ierr)
    {
      // ******************************************
      //   no inner iteration loop is required  ***
      // ******************************************
      modif++;
      if (modif>1)
      {
        //  arc length modification
        dl*=2.0;
        if (dl>dlmax)  dl=dlmax;
        modif=0;

        if (Mespr==1)
        {
          fprintf (stdout,"\n arc length must be modified (dl increase) dl=%e    norf=%e",dl,norf);
          fprintf (Out,"\n arc length must be modified (dl increase) dl=%e    norf=%e",dl,norf);
        }
      }

      lambda+=dlambda;
      if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
      {
        //  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;
      }
      if (fabs(lambda-rlambda) <= rlerr)
        break;
      Mm->updateipval();
      continue;
    }

    // ****************************
    //   inner iteration loop  ****
    // ****************************
    stop=0;
    for (j=0;j<ini;j++)
    {
//Smat->solve_system (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;

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

      //  zacatek novinky

      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
//      fprintf(Out, "Dlambda %g, i = %ld, j = %ld\n", dlambda, i, j);

      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;

      //  computation of internal forces
      internal_forces (lcid,fi);

      for (k=0;k<n;k++)
        f[k]=fa[k]-fi[k];

      norf=ss(f,f,n);
      norf=sqrt(norf);
      norfa = ss(fa,fa,n);
      norfa=sqrt(norfa);
      norf /= norfa;

      if (Mespr==1)
      {
        fprintf (stdout,"\n increment  %ld     inner loop  %ld    norf=%e",i,j,norf);
        fprintf (Out,"\n increment  %ld     inner loop  %ld    norf=%e",i,j,norf);
      }

      if (norf<ierr)
      {
        lambda+=ddlambda;

        stop=1;
        if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
          stop=0;
        if (fabs(lambda-rlambda) <= rlerr)
	  stop=2;
        if (stop > 0)
          Mm->updateipval();
        break;
      }
    }

    // limit value of lambda was reached
    if (stop == 2)
      break;

    modif=0;

    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;
      //  restoring of state variables
//      Mm->restore_state_var ();
    }
  }
  fclose(grout);
  delete [] r;
  delete [] fi;
  delete [] fa;
  delete [] fc;
  delete [] f;
  delete [] v;
  delete [] u;
  delete [] ddr;
  return lambda;
}
*/

double newton_raphsonep (long lcid, double *pfp)
  //  function solves system of nonlinear algebraic
  //  equations by Newton-Raphson method
  //  solved system does not contain time variable, it takes into account
  //  previous reached displacement/stress state and starts with nonzero 
  //  displacement vector. pfp is proportional load vector of the previous 
  //  reached initial displacements.
  //
  //  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;

  //  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;
  
  i=0;
/*  char fn[1001];
  FILE *femcad;
  sprintf(fn,"%s%s.jk0", Mp->path, Mp->filename);
  femcad = fopen(fn, "wt");
  export_femcad(femcad, i);
  write_displ(femcad, 0, 0);
  Mm->compute_nodeothers(0);
  for (j=0; j<Mm->ip[0].ncompother; j++)
    write_nodscalar(femcad, other, 0, j, 0);
  Mm->compute_nodestresses(0);
  for (j=0; j<Mm->ip[0].ncompstr; j++)
    write_nodscalar(femcad, stress, 0, j, 0);*/
  // ***************************
  //  main iteration loop  ****
  // ***************************
  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;
    //  backup of state variables
//    Mm->backup_state_var ();
   
    fprintf (stdout,"\n lambda %e,    dlambda %e",lambda,dlambda);
//    fprintf (Out,"\n\n i=%ld     lambda %f",i,lambda);
    //print_displacements (Out,0);

    
    //  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
    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]=pfp[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 ();
/*      write_displ(femcad, 0, i+1);
      Mm->compute_nodeothers(0);
      for (j=0; j<Mm->ip[0].ncompother; j++)
        write_nodscalar(femcad, other, 0, j, i+1);
      Mm->compute_nodestresses(0);
      for (j=0; j<Mm->ip[0].ncompstr; j++)
        write_nodscalar(femcad, stress, 0, j, i+1);*/
      print_step(lcid, i, 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++){
	fb[k]=pfp[k]+f[k]+dlambda*fp[k]-fi[k];
      }
      
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,n);
      norf /= norfa;
      
      fprintf (stdout,"\n j=%ld     norf=%e",j,norf);
//      fprintf (Out,"\n j=%ld     norf=%e",j,norf);
      
      if (norf<err)  break;
    }
    
    if (j==ini || norf>err){
      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);
//	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;
      //      Mm->restore_state_var ();
    }
    else{
      lambda+=dlambda;
      Mm->updateipval ();
      print_step(lcid, i, lambda, f);
      print_flush();
/*      write_displ(femcad, 0, i+1);
      Mm->compute_nodeothers(0);
      for (j=0; j<Mm->ip[0].ncompother; j++)
        write_nodscalar(femcad, other, 0, j, i+1);
      Mm->compute_nodestresses(0);
      for (j=0; j<Mm->ip[0].ncompstr; j++)
        write_nodscalar(femcad, stress, 0, j, i+1);*/
    }
    
  }
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f;

  
//  fclose (femcad);
  return lambda;
}