epsolver.cpp

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

  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilal;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  ierr = Mp->nlman.erral;
  //  length of arc
  dl = Mp->nlman.dlal;
  //  maximum length of arc
  dlmax = Mp->nlman.dlmaxal;
  //  minimum length of arc
  dlmin = Mp->nlman.dlminal;
  //  displacement-loading driving switch
  psi = Mp->nlman.psial;
  
  //  allocation of auxiliary arrays
  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];
  
  //  initialization phase
  ra = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  
  if (Mp->nlman.hdbackupal==2){
    arclopen (li,n,lambda,dl,ra,fp);    
  }
  else if (Mp->nlman.hdbackupal==11){
    arclopen_adapt (li,lambda,blambda,dl);
  }
  else{
    lambda=0.0;
    lambdao=0.0;
    li=0;
  }
  
  if (li == 0)
    mode = "wt";
  else
    mode = "at";  
  
  sprintf(file, "%s%s.epg", Mp->path, Mp->filename);
  grout = fopen (file, "wt");  
  
  for (j=0;j<n;j++)
    fa[j]=fc[j]+(lambda)*fp[j];

  //  norm of proportionality vector
  norfp = ss(fp,fp,n);
  modif=0;

  // pouze docasne pro jeden typ ulohy
  long kk;
  matrix sm(6,6);

  // ***************************
  //  main iteration loop   ****
  // ***************************
  for (i=li;i<ni;i++) {

    fprintf(grout, "Lambda %e\n", lambda);
    fprintf(grout, "Posuny ux --------\n");
    for (kk = 0; kk < Mt->nn; kk++)
    {
      cn = Mt->give_dof(kk,0)-1;
      fprintf(grout, "   %2ld %e\n", kk+1, ra[cn]);
    }
    fprintf(grout, "Zatizeni fx --------\n");
    for (kk = 0; kk < Mt->nn; kk++)
    {
      cn = Mt->give_dof(kk,0)-1;
      double fpv = fp[cn];
      fprintf(grout, "   %2ld %le\n", kk+1, fpv);
    }
    fprintf(grout, "Tuhosti a sily--------\n");
    for (kk = 0; kk < Mt->ne; kk++)
    {
      if (Mt->give_elem_type(kk) == spring_1)
      {
        Spring->res_stiffness_matrix (kk,sm);
        Mt->give_code_numbers (kk, ecn);
        if (ecn[0] > 0)
          fprintf(grout, "   %2ld %e %e\n", kk+1, sm[0][0], sm[0][0]*ra[ecn[0]-1]);
      }
    }
    fprintf(grout, "\n");
    fflush(grout);
    //  backup of left hand side vector
    copyv (ra, r, n);
    //  backup of reached lambda parameter
    blambda=lambda;

    fprintf (stdout,"\n arc-length: increment %ld   lambda %e  dl %e",i,lambda,dl);
/*    if (Mespr==1){
      fprintf (Out,"\n\n *******************************************************************");
      fprintf (Out,"\n arc-length: increment %ld,   lambda=%e   dl=%e",i,lambda,dl);
    }*/
    
    //  assembling of tangent stiffness matrix
    stiffness_matrix (lcid);
    
    //  backup of the fp, in ldl_sky the right hand side will be destroyed
    copyv (fp, f, n);
    //  solution of K(r).v=F
    //Smat->solve_system (Gtm,v,f);
    Mp->ssle.solve_system (Gtm,Smat,v,f,Out);
    //  generalized norm of displacement increments
    norv = displincr (v,n);
    //  compute new dlambda increment
    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);
    subv(fa, fi, f, n);    
    norf=sizev(f,n);
    norfa = sizev(fa, n);
    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++;

      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;
        }
        modif = 0;
        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 (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);
	}
      }
      if (fabs(lambda-rlambda) <= rlerr)
        break;
      Mm->updateipval ();
//      aux_mech_nonlin_print (gr,ra,lambda);
//      print_output(lcid);
        print_step(lcid, i, lambda, fa);
        print_flush();
    }
    else{
      // ****************************
      //  inner iteration loop  ****
      // ****************************
      stop=0;
      for (j=0;j<ini;j++){
	//  back substitution
	//Smat->solve_system (Gtm,u,f);
	Mp->ssle.solve_system (Gtm,Smat,u,f,Out);
        copyv (ddr, f, n);
        addv(ddr, u, n);
	//  coefficient of quadratic equation
	quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
	/*
	  fprintf (Out,"\n\n\n Kontrola kvadraticke rovnice");
	  fprintf (Out,"\n norfp     %15.10le",norfp);
	  fprintf (Out,"\n norv      %15.10le",norv);
	  fprintf (Out,"\n (ddr,v)   %15.10le",ss(ddr,v,n));
	  fprintf (Out,"\n (ddr,ddr) %15.10le",ss(ddr,ddr,n));
	  fprintf (Out,"\n dl        %15.10le",dl);
	  fprintf (Out,"\n ddlambda  %15.10le",ddlambda);
	  fprintf (Out,"\n psi       %15.10le",psi);
	  fprintf (Out,"\n a2        %15.10le",a2);
	  fprintf (Out,"\n a1        %15.10le",a1);
	  fprintf (Out,"\n a0        %15.10le",a0);
	  fprintf (Out,"\n discrim   %15.10le",a1*a1-4.0*a2*a0);
	  fprintf (Out,"\n\n");
	*/
	
	//  solution of quadratic equation
        numr = solv_polynom_2(a2, a1, a0, l1, l2);
        switch (numr)
	{
	  case -1 :
            fprintf (stderr,"\n\n infinite number of solution of constrained condition in function arclength");
            break;
	  case 0 :
            fprintf (stderr,"\n\n nonsolvable constrained condition in function arclength");
            break;
    	  case 1 :
            dlambda = l1;
            break;
	  default :
            break;
	}
	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;

	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;
	
//	fprintf (Out,"\n ddlambda %e     dlambda %e",ddlambda,dlambda);
	//  computation of internal forces
	internal_forces (lcid,fi);
        subv(fa, fi, f, n);	
	norf=sizev(f,n);
	norfa = sizev(fa, n);
	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;
//	  aux_mech_nonlin_print (gr,ra,lambda);
          print_step(lcid, i, lambda, fa);
          print_flush();
	  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
        copyv (r, ra, n);
	//  restoring of lambda parameter
	lambda=blambda;
      }
    }
    if (stop==0)
      continue;

    // termit
    if (i%10 && Mm->plast)
      continue;
    newmesh = 0;
    if (Mp->nlman.nial > (i+1) && Mp->adaptivity && !( (Mp->adaptflag & 16) && (Mp->nlman.nial==100 || Mp->nlman.nial==1))){
      // print of fxlambda
      FILE *ff; sprintf (file,"%sfxlambda.dat",Mp->path);
      ff = fopen (file,"a");  if (ff==NULL) fprintf (stderr,"\n ff file has not been opened.\n\n");
      fprintf (ff,"  %20.10f",lambda);  fclose (ff);
      long width = 1 + (long)log10((double)Mp->nlman.nial);
      newmesh = Ada->run (1,width,i);
    }
    if (newmesh){
      i++; // toto se da odstranit
      break;
    }
    // timret
    print_output(lcid);
  }
  
  // ------------------------------------
  //  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);
  
  //  if (i==ni){   sprintf (file,"%s%s.dx.%ld.finito",Mp->path,Mp->filename,i);   print_default_2_dx (Gtm,Mp,Mt,0,file);  }
  // timret
  
  delete [] r;		    
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
  //fclose (gr); 
  return (lambda);
}














/* -------------------------------------------------------------------------------------- */







/**
  This function solves system of nonlinear algebraic equations
  by the arc-length method. This method was modified to stop iteration
  on the given value rlambda of the lambda parameter. Rlambda value is reached with
  required error rlerr. Only one right hand side vector is supported with respect
  to nonlinearity and absence of superposition.


  @param lcid    - load case id
  @param rlambda - required value of lambda parameter
  @param rlerr   - required error between lambda and rlambda

  @retval The function returns reached lambda parameter.

double arclengthrv (long lcid, double rlambda, double rlerr)
  //
  //  d odpovida delte
  //  dd odpovida DELTE (trojuhelniku)
  //
  //  fc - konstantni vektor
  //  fp - vektor proporcionality
  //
  //  n - pocet radku matice
{
  long i,j,k,n,ni,ini,stop,modif,li, kk;
  double a0,a1,a2,d,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norfa,norv,zero,ierr;
  double lambdao,ss1,ss2,ss3,ss4,ss5;
  double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
  FILE *grout;

  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nial;
  //  maximum number of iterations in one increment
  ini = Mp->niilal;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  ierr = Mp->erral;
  //  length of arc
  dl = Mp->dlal;
  //  maximum length of arc
  dlmax = Mp->dlmaxal;
  //  minimum length of arc
  dlmin = Mp->dlminal;
  //  displacement-loading driving switch
  psi = Mp->psial;

  //  allocation of auxiliary arrays
  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];
  fc = new double [n];
  memset(fc, 0, sizeof(*fc)*n);

  //  allocation of backup arrays for state variables in integration points
  //  Mm->alloc_backup ();

  //  initialization phase
  ra = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid);
//  fc = Lsrs->give_rhs (lcid*2+1);

  //  array containing selected DOFs
  seldofinit ();

  lambda=0.0;
  lambdao=0.0;
  li=0;

  //  norm of proportionality vector
  norfp = ss (fp,fp,n);
  matrix sm(6,6);
  modif=0;
  grout = fopen("gr.dat", "wt");
  // ***************************
  //   main iteration loop  ****
  // ***************************
  for (i=li;i<ni;i++){

    fprintf(grout, "%e\n", lambda);
    for (kk = 0; kk < Mt->nn; kk++)