nfssolver.cpp

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

      print_init(-1, "wt");
      print_step(lcid, li, lambda, fa);
    }
  
  
  
  // ***************************
  //  main iteration loop   ****
  // ***************************
  for (i=li;i<ni;i++){
    
    fprintf (Out,"\n\n arc-length  prirustek %ld",i);
    

    stop=1; // termit
    //  backup of left hand side vector
    copyv (ra, r, nlm);
    //  backup of reached lambda parameter
    blambda=lambda;

    //  backup of the fp, in ldl_sky the right hand side will be destroyed
    copyv(fp, f, n);

    //  solution of K(r).v=F
    Fsd->pmcg (f,Out);
    copyv (Fsd->w,v,nlm);

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

    //  compute new dlambda increment
    dlambda = dl/sqrt(norv+psi*psi*norfp);
    
    for (j=0;j<nlm;j++){
      ddr[j]=dlambda*v[j];
      ra[j]+=ddr[j];
    }
    for (j=0;j<n;j++){
      fa[j]=fc[j]+(lambda+dlambda)*fp[j];
    }
    ddlambda=dlambda;
    
    Fsd->add_mult (dlambda);
    Fsd->mult_correction ();

    Fsd->nonlinlagrmultdispl (d,fa);

    
    //  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     norv %e      norf %e",i,norv,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",dl);
//	  fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl);
	}
      }
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, fa);
      print_flush();
    }
    else{
      // ****************************
      //  inner iteration loop  ****
      // ****************************
      stop=0;
      for (j=0;j<ini;j++){
	
	fprintf (Out,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	
	//  solution of K(r).v=F
	Fsd->pmcg (f,Out);
	copyv (Fsd->w,u,nlm);
	
	
        //copyv(ddr, f, nlm);
        addv(ddr, u, nlm);
	//  coefficient of quadratic equation
	//quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
	quadeqcoeff (ddr,v,nlm,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;
	}

	if (fabs(l1)>fabs(l2))  dlambda=l2;
	else dlambda=l1;

	for (k=0;k<nlm;k++){
	  ddr[k]+=dlambda*v[k];
	  ra[k]+=u[k]+dlambda*v[k];
	}
	for (k=0;k<n;k++){
	  fa[k]+=dlambda*fp[k];
	}
	ddlambda+=dlambda;
	
	Fsd->add_mult (dlambda);
	Fsd->mult_correction ();
	
	Fsd->nonlinlagrmultdispl (d,fa);
	
	//  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    norf=%e  ierr=%e",norf,ierr);
	  //fprintf (Out,"\n increment  %ld     inner loop  %ld    norf=%e",i,j,norf);
	}

	if (norf<ierr){
	  lambda+=ddlambda;
	  Mm->updateipval ();
	  stop=1; 
	  compute_req_val (lcid);
          print_step(lcid, i+1, lambda, fa);
          print_flush();
	  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",dl);
	  //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl);
	}
	//  restoring of left hand side vector
        copyv(r, ra, nlm);
	//  restoring of lambda parameter
	lambda=blambda;
      }
    }
    
    fprintf (stdout,"\n increment %ld    total lambda %e",i,lambda);
    
    if (stop==0)
      continue;
    
  }
  
  // ------------------------------------
  //  finish of main iteration loop  ----
  // ------------------------------------
  
  
  print_close();
  delete [] r;		    
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
}
*/

/*
void solve_nonlinear_floating_subdomains ()
{
  long i,j,k,n,lcid,nlm,ni,ini,stop,modif,li,newmesh, numr;
  double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
  double lambdao,ss1,ss2,ss3,ss4,ss5;
  double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi,*rhs;
  
  
  //  load case id must be equal to zero
  //  only one load case can be solved
  //  one load case may contain several proportional vectors
  lcid = 0;
  
  //  assembling of stiffness matrix
  stiffness_matrix (lcid);
  
  //  number of unknown displacements
  n = Ndofm;
  
  //  auxiliary assigment of pointer
  rhs=Lsrs->give_rhs (lcid);
  
  //  assembling of right hand side vector (vector of nodal forces)
  mefel_right_hand_side (lcid,rhs);
  
  //  vector of nodal displacements
  ra = Lsrs->give_lhs (lcid);
  //  vector of proportional load
  fp = Lsrs->give_rhs (lcid*2);
  //  vector of constant load
  fc = Lsrs->give_rhs (lcid*2+1);
  
  
  
  //  allocation of object floating subdomains
  Fsd = new flsubdom;
  
  //  computation of rigid body motions
  //  list of dependent equations
  //  determination of number of Lagrange multipliers
  Fsd->initialization (Ndofm,Mp->ense,fp);
  
  //  number of Lagrange multipliers
  nlm = Fsd->nlm;
  
  
  
  // *******************************
  //  data about arc-length solver
  // *******************************

  //  maximum number of increments
  ni = Mp->nlman.nial;
  //  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];
  
  
  dlambda=0.0;
  lambda=0.0;
  lambdao=0.0;
  li=0;

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


  // assembling reached load vector
  for (j=0;j<n;j++){
    fa[j]=fc[j]+(lambda+dlambda)*fp[j];
  }
  if (li)
    print_init(-1, "at");
  else
  {
    print_init(-1, "wt");
    print_step(lcid, li, lambda, fa);
  }
  

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

    fprintf (Out,"\n\n arc-length  prirustek %ld",i);


    stop=1; // termit
    //  backup of left hand side vector
    copyv (ra, r, n);
    //  backup of reached lambda parameter
    blambda=lambda;

    
    //  backup of the fp, in ldl_sky the right hand side will be destroyed
    copyv(fp, f, n);
    
    
    stiffness_matrix (lcid);
    Fsd->factorize ();


    //  solution of K(r).v=F
    Fsd->pmcg (f,Out);
    
    Fsd->lagrmultdispl (v,f);
    
    
    //  generalized norm of displacement increments
    norv = displincr (v,n);

    //  compute new dlambda increment
    dlambda = dl/sqrt(norv+psi*psi*norfp);
    
    //if (norv+psi*psi*norfp < 1.0e-8)  dlambda=0.0;
    //else 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;
    
    Fsd->add_mult (dlambda);
    Fsd->mult_correction ();
    
    //  computation of internal forces
    internal_forces (lcid,fi);
    subv(fa, fi, f, n);    
    norf=sizev(f,n);
    norfa = sizev(fa, n);
    //if (norfa<1.0e-8){}
    //else norf /= norfa;
    norf /= norfa;

    
    //if (Mespr==1)  fprintf (stdout,"\n ddlambda %e    dlambda %e   norf %e  ierr %e",ddlambda,dlambda,norf,ierr);
    if (Mespr==1)  fprintf (stdout,"\n increment %ld     norv %e      norf %e",i,norv,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",dl);
//	  fprintf (Out,"\n arc-length must be modified (dl increase) dl=%e",dl);
	}
      }
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, fa);
      print_flush();
    }
    else{
      // ****************************
      //  inner iteration loop  ****
      // ****************************
      stop=0;
      for (j=0;j<ini;j++){


	fprintf (Out,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	fprintf (stdout,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	
	
	stiffness_matrix (lcid);
	Fsd->factorize ();
	
	
	//  back substitution
	Fsd->pmcg (f,Out);
	Fsd->lagrmultdispl (u,f);
	
        copyv(ddr, f, n);
        addv(ddr, u, n);
	//  coefficient of quadratic equation
	quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2);
	
	//  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;
	}

	if (fabs(l1)>fabs(l2))  dlambda=l2;
	else dlambda=l1;

	//fprintf (stdout,"\n increment  %ld     inner loop  %ld   x1=%e x2=%e   dlambda %e",i,j,l1,l2,dlambda);
	//fprintf (stdout,"\n increment  %ld     inner loop  %ld   x1=%e x2=%e   ddlambda %e",i,j,l1,l2,ddlambda);

	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;
	
	Fsd->add_mult (dlambda);
	Fsd->mult_correction ();
	
	
	//  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    norf=%e  ierr=%e",norf,ierr);
	  //fprintf (Out,"\n increment  %ld     inner loop  %ld    norf=%e",i,j,norf);
	}

	if (norf<ierr){
	  lambda+=ddlambda;
	  Mm->updateipval ();
	  stop=1; 
	  compute_req_val (lcid);
          print_step(lcid, i+1, lambda, fa);
          print_flush();
	  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",dl);
	  //fprintf (Out,"\n arc length must be modified (dl decrease) dl=%e",dl);
	}
	//  restoring of left hand side vector
        copyv(r, ra, n);
	//  restoring of lambda parameter
	lambda=blambda;
      }
    }
    
    fprintf (stdout,"\n increment %ld    total lambda %e",i,lambda);

    if (stop==0)
      continue;
    
  }
  
  // ------------------------------------
  //  finish of main iteration loop  ----
  // ------------------------------------
  
  
  if (Mp->nlman.hdbackupal==1)
    arclsave (i,n,lambda,dl,ra,fp);

  print_close();
  delete [] r;		    
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
}
*/


/**
   function solves system of nonlinear algebraic equation
   by arc-length method (continuation method)
   only one right hand side vector is supported with respect
   to nonlinearity and absence of superposition