nfssolver.cpp

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

   JK, 20.6.2006 (25.7.2001)
*/
/*
void solve_nonlinear_floating_subdomains_24_6 ()
{
  long i,j,k,n,l,ni,ini,stop,modif,li,numr,lcid,nlm;
  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,*blm;
  
  Fsd = new flsubdom;
  
  lcid=0;

  rhs=Lsrs->give_rhs (0);
  mefel_right_hand_side (lcid,rhs);
  
  //  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;
  
  dlambda=0.0;
  
  //  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
  //  vector of nodal displacements
  ra = Lsrs->give_lhs (lcid);
  //  vector of proportional load
  fp = Lsrs->give_rhs (lcid*2);
  //  vector of nonproportional load
  fc = Lsrs->give_rhs (lcid*2+1);
  
  
  if (Mp->nlman.hdbackupal == 2){            // termit
    arclopen (li,n,lambda,dl,ra,fp);
  }
  else{
    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);
    }
  
  
  
  stiffness_matrix (lcid);
  Fsd->initialization (Ndofm,Mp->ense,f);
  //  number of Lagrange multipliers
  nlm = Fsd->nlm;
  
  blm = new double [nlm];

  double *buddr;
  buddr = new double [n];
  double *buf;
  buf = new double [n];
  double *bulm;
  bulm = new double [nlm];

  // ***************************
  //  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 reached Lagrange multipliers
    copyv (Fsd->tw,blm,nlm);

    //fprintf (stdout,"\n\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
    if ((Mp->stmat==tangent_stiff) || (i == li)){
      //stiffness_matrix (lcid);
      //Fsd->initialization (Ndofm,Mp->ense,f);
    }
    
    stiffness_matrix (lcid);
    

    //  backup of the fp, in ldl_sky the right hand side will be destroyed
    copyv(fp, f, n);
    
    Fsd->initialization (Ndofm,Mp->ense,f);
    
    Fsd->solve_lin_alg_system (v,f);
    
    //  solution of K(r).v=F
    //Smat->solve_system (v,f);
    
    //  generalized norm of displacement increments
    norv = displincr (v,n);
    
    //  compute new dlambda increment
    dlambda = dl/sqrt(norv+psi*psi*norfp);
    
    Fsd->add_mult (dlambda);
    Fsd->mult_correction ();
    
    //  novinka
    if (Fsd->nch!=0){
      for (j=0;j<ini;j++){

	fprintf (stdout,"\n cyklime v casti jedna  %ld",j);

	//  restoring of Lagrange multipliers
	copyv (blm,Fsd->tw,nlm);
	
	Fsd->solve_lin_alg_system (v,f);
	
	//  generalized norm of displacement increments
	norv = displincr (v,n);
	
	//  compute new dlambda increment
	dlambda = dl/sqrt(norv+psi*psi*norfp);
	
	Fsd->add_mult (dlambda);
	Fsd->mult_correction ();
	
	if (Fsd->nch==0)
	  break;
      }
    }
    //  konec novinky
    
    
    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);
    //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++){
	if (Mp->stmat==tangent_stiff){
	  stiffness_matrix (lcid);
	  Fsd->initialization (Ndofm,Mp->ense,f);
	}
	
	fprintf (Out,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	
	
        copyv(f, buf, n);
	
	Fsd->solve_lin_alg_system (u,f);
	
	//  back substitution
	//Smat->solve_system (u,f);
	
	//  backup of the vector ddr
        copyv(ddr, buddr, n);
	
	
        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;
	}
	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 (fabs(l1)>fabs(l2))  dlambda=l2;
	else dlambda=l1;
	
	copyv (Fsd->tw,bulm,nlm);

	Fsd->add_mult (dlambda);
	Fsd->mult_correction ();
	
	//  novinka
	if (Fsd->nch!=0){
	  for (l=0;l<ini;l++){
	    
	    fprintf (stdout,"\n cyklime v casti dve  %ld",j);
	    
	    
	    //  restoring of Lagrange multipliers
	    copyv (bulm,Fsd->tw,nlm);
	    copyv (buddr,ddr,n);
	    copyv (buf,f,n);
	    
	    Fsd->solve_lin_alg_system (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;
	      }
	    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 (fabs(l1)>fabs(l2))  dlambda=l2;
	    else dlambda=l1;
	    
	    
	    Fsd->add_mult (dlambda);
	    Fsd->mult_correction ();
	    
	    if (Fsd->nch==0)
	      break;
	  }
	}
	//  konec novinky




	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 (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 (norfa<1.0e-8){}
	//else 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;
	//  restoring of Lagrange multipliers
	copyv (blm,Fsd->tw,nlm);
      }
    }
    
    fprintf (stdout,"\n increment %ld    total lambda %e",i,lambda);

    if (stop==0)
      continue;
    
  }
  
  // ------------------------------------
  //  finish of main iteration loop  ----
  // ------------------------------------
  
  
  fprintf (stdout,"\n\n multiplikatory   %lf",Fsd->tw[0]);
  
  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;
  
}
*/


/**
   JK, 26.6.2006
*/
/*
void solve_nonlinear_floating_subdomains_29_6 ()
{
  long i,j,k,n,lcid,ni,ini,stop,modif,li,numr,nlm;
  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,*d,*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
  d   = 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 [nlm];
  ra = new double [nlm];
  ddr = new double [nlm];
  u = new double [nlm];
  v = new double [nlm];
  f = new double [n];
  fa = new double [n];
  fi = new double [n];
  
  
  lambda=0.0;
  lambdao=0.0;
  dlambda=0.0;
  li=0;
  modif=0;  
  
  
  
  
  
  //  norm of proportionality vector
  norfp = ss(fp,fp,n);
  
  
  
  // 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
    {