nfssolver.cpp

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

   d stands for delta
   dd stands for capital delta
   
   fc - nonproportional %vector
   fp - proportional %vector
   n - order of the system
   
   @param lcid - load case id
   
   25.7.2001
*/
/*
void solve_nonlinear_floating_subdomains ()
{
  long i,j,k,n,nlm,ni,ini,li,numr,stop,modif,lcid;
  double dl,dlmax,dlmin,psi,ierr,zero,lambda,dlambda,ddlambda,blambda,norf,norfa;
  double a0,a1,a2,l1,l2,norv;
  double *r,*ddr,*fc,*fp,*f,*dv,*du,*dnu,*rhs;
  double *br,*bmu,*av,*av1,*av2,*fa,*fi,*aa1,*aa2,*a,*b,*c;

  
  //long i,j,k,n,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;

  
  
  //  number of rows of the matrix
  n = Ndofm;
  //  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;
  
  //  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);


  //  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
  r  = 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;
  
  
  //  allocation of auxiliary arrays
  //  array for backup of nodal displacements
  br = new double [n];
  //  array for backup of Lagrange multipliers
  bmu = new double [nlm];

  du = new double [n];
  dv = new double [n];
  f = new double [n];
  fa = new double [n];
  fi = new double [n];
  av = new double [n];
  av1 = new double [n];
  av2 = new double [n];
  aa1 = new double [n];
  aa2 = new double [n];

  dnu = new double [nlm];
  
  a = new double [n];
  b = new double [n];
  c = new double [n];
  

  //  proportional factor
  lambda = 0.0;
  dlambda = 0.0;
  ddlambda = 0.0;
  
  //  attained level in arclength method
  //  it is 0 by default
  //  otherwise it is equal to performed number of steps in previous computations
  li=0;
  
  modif=0;
  
  
  
  //  initiation of load vector
  for (j=0;j<n;j++){
    fa[j]=fc[j]+lambda*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 increment %ld",i);
    fprintf (stdout,"\n\n arc-length increment   %ld",i);


    //stop=1; // termit



    //  backup of nodal displacements
    copyv (r, br, n);
    //  backup of Lagrange multipliers
    copyv (Fsd->muo,bmu,nlm);
    //  backup of reached lambda parameter
    blambda=lambda;



    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == li)){
      stiffness_matrix (lcid);
      Fsd->factorize ();
    }
    
    //  solution of K(r).dv = fp-B.dnu
    Fsd->pmcg (fp,Out);

    for (k=0;k<nlm;k++){
      fprintf (Out,"\n krok %ld                 lambda %ld   %lf",i,k,Fsd->w[k]);
    }
    

    Fsd->lagrmultdispl (dv,fp);
    copyv (Fsd->w,dnu,nlm);
    
    
    //  generalized norm of displacement increments
    norv = displincr (dv,n);
    
    Gtm->flsub.coarse_local (dnu,av);
    
    norf=0.0;
    for (k=0;k<n;k++){
      norf+=(fp[k]-av[k])*(fp[k]-av[k]);
    }
    
    //  compute new dlambda increment
    dlambda = dl/sqrt(norv+psi*psi*norf);
    
    
    for (k=0;k<n;k++){
      ddr[k]=dlambda*dv[k];
      r[k]+=ddr[k];
      fa[k]=fc[k]+(lambda+dlambda)*fp[k];
    }
    ddlambda=dlambda;
    for (k=0;k<nlm;k++){
      Fsd->ddmu[k]=dlambda*dnu[k];
      Fsd->mu[k]+=Fsd->ddmu[k];
    }
    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 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;
      
      for (k=0;k<nlm;k++){
	Fsd->muo[k]+=Fsd->ddmu[k];
      }
      
      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++){
	if (Mp->stmat==tangent_stiff)
	  stiffness_matrix (lcid);
	
	Fsd->factorize ();
	
	//  back substitution
	Fsd->pmcg (f,Out);
	Fsd->lagrmultdispl (du,f);

	for (k=0;k<nlm;k++){
	  fprintf (Out,"\n krok %ld iterace %ld     lambda %ld   %lf",i,j,k,Fsd->w[k]);
	}
	
	fprintf (Out,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	
	
	Gtm->flsub.coarse_local (Fsd->ddmu,av1);
	Gtm->flsub.coarse_local (Fsd->w,av2);
	
	for (k=0;k<n;k++){
	  a[k]=ddr[k]+du[k];
	}
	
	for (k=0;k<n;k++){
	  b[k]=ddlambda*fp[k]-av1[k]-av2[k];
	}
	
	
	Gtm->flsub.coarse_local (dnu,av);
	for (k=0;k<n;k++){
	  c[k]=fp[k]-av[k];
	}
	
	a2=0.0;  a1=0.0;  a0=0.0;
	for (k=0;k<n;k++){
	  a2+=c[k]*c[k]*psi*psi+dv[k]*dv[k];
	  a1+=2.0*b[k]*c[k]*psi*psi+2.0*dv[k]*a[k];
	  a0+=a[k]*a[k]+b[k]*b[k]*psi*psi;
	}
	a0-=dl*dl;
	
	//  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<n;k++){
	  ddr[k]+=du[k]+dlambda*dv[k];
	  r[k]+=du[k]+dlambda*dv[k];
	  fa[k]+=dlambda*fp[k];
	}
	ddlambda+=dlambda;
	for (k=0;k<nlm;k++){
	  Fsd->ddmu[k]+=Fsd->w[k]+dlambda*dnu[k];
	  Fsd->mu[k]+=Fsd->w[k]+dlambda*dnu[k];
	}
	Fsd->mult_correction ();
	
	//fprintf (stdout,"   ddlambda %e",ddlambda);
	
	//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    norf=%e  ierr=%e",norf,ierr);
	  //fprintf (Out,"\n increment  %ld     inner loop  %ld    norf=%e",i,j,norf);
	}

	if (norf<ierr){
	  lambda+=ddlambda;
	  
	  for (k=0;k<nlm;k++){
	    Fsd->muo[k]+=Fsd->ddmu[k];
	  }
	  
	  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 nodal displacements
        copyv(br, r, n);
	//  restoring of Lagrange multipliers
        copyv(bmu,Fsd->muo,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  ----
  // ------------------------------------
  
  
  
  
  //if (Mp->nlman.hdbackupal==1)
  //arclsave (i,n,lambda,dl,ra,fp);

  print_close();
  

  delete [] br;
  delete [] bmu;
  delete [] dv;
  delete [] f;
  delete [] fa;  
  delete [] fi;
  delete [] av;
  delete [] av1;
  delete [] av2;
  delete [] aa1;
  delete [] aa2;
  delete [] dnu;

}
*/




/**
   function solves system of nonlinear algebraic equations
   load is applied by prescribed displacements
   
   JK, 19.1.2007
*/
/*
void solve_nonlinear_floating_subdomains_19_1 ()
{
  long i,j,ni,ini,lcid;
  double lambda,dlambda,mindlambda,err,norfi;
  double *f,*lhs,*rhs,*arhs,*fi;
  
  //  only one load case is assumed
  lcid=0;
  //  number of iteration steps/increments
  ni=Mp->nlman.ninr;
  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilnr;
  //  total magnitude of the right hand side
  lambda=0.0;
  //  magnitude of the step increment
  dlambda=Mp->nlman.incrnr;
  //  minimum magnitude of the step increment
  mindlambda=Mp->nlman.minincrnr;
  //  required error
  err=Mp->nlman.errnr;
  
  //
  // fi = new double [?]
  
  Fsd = new flsubdom;
  
  //  nodal displacements
  lhs=Lsrs->give_lhs (lcid);
  //  proportional vector of the prescribed forces
  rhs=Lsrs->give_rhs (lcid);
  //  auxiliary array
  f = new double [Ndofm];
  //  auxiliary vector
  arhs = new double [Ndofm];
  
  
  //  proportional vector of the prescribed forces
  mefel_right_hand_side (lcid,rhs);
  
  //  stiffness matrix
  stiffness_matrix (lcid);
  
  Fsd->initialization (Ndofm,Mp->ense,rhs);
  
  //  loop over iteration steps/increments
  for (i=0;i<ni;i++){
    
    //  increment of the right hand side
    for (j=0;j<Ndofm;j++){
      f[j]=dlambda*rhs[j];
    }
    
    //  solution of the system of algebraic equations
    //  increment of the Lagrange multipliers are obtained
    Fsd->pmcg (f,Out);
    
    //  update of Lagrange multipliers
    Fsd->add_mult (1.0);
    
    //  correction of the Lagrange multipliers
    Fsd->mult_correction ();
    
    
    for (j=0;j<Ndofm;j++){
      arhs[j]=lambda*rhs[j];
    }
    
    //  computation of displacements from Lagrange multipliers
    Fsd->lagrmultdispl (lhs,arhs);
    
    
    //  actual internal forces
    internal_forces (lcid,fi);
    
    //  norm of the vector of internal forces
    norfi = ss (fi,fi,Ndofm);
    
    if (norfi<err){
      //  no inner loop is required
    }
    else{
      //  inner loop is required
      for (j=0;j<ini;j++){
	
	Fsd->solve_lin_alg_system (lhs,fi);
	
	
	internal_forces (lcid,fi);
	
      }
      
    }
    
    
    print_init(-1, "wt");    
    for (i=0;i<Lsrs->nlc;i++){
      //compute_ipstrains (i);
      //compute_ipstresses (i);
      //compute_req_val (i);
      print_step(i, 0, 0.0, NULL);    
    }

    
  }
  
  print_close();
  

}
*/