nssolver.cpp

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

		  double &a0,double &a1,double &a2)
{
  long i;
  
  switch (Mp->nlman.displnorm){
  case alldofs:{
    a2=ss(v,v,n)+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;
    break;
  }
  case seldofs:
  case seldofscoord:
  case selecnodes:{
    a0=0.0;  a1=0.0;  a2=0.0;
    for (i=0;i<Mp->nlman.nsdofal;i++){
      a0+=ddr[Mp->nlman.seldofal[i]]*ddr[Mp->nlman.seldofal[i]];
      a1+=ddr[Mp->nlman.seldofal[i]]*v[Mp->nlman.seldofal[i]]*2.0;
      a2+=v[Mp->nlman.seldofal[i]]*v[Mp->nlman.seldofal[i]];
    }
    a0-=dl*dl;
    break;
  }
  case nodesdistincr:{
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown norm of displacement increment is required in function");
    fprintf (stderr,"\n displincr (%s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

void arclsave (long fi,long n,double blambda,double dl,double *r,double *fp)
{
  long i,j;
  char *file;
  FILE *aux;
  file = new char[strlen(Mp->path)+11+1];
  sprintf (file,"%sarcl.backup",Mp->path);
  aux = fopen (file,"w");

  //aux = fopen ("arcl.backup","w");
  
  //  attained number of steps
  fprintf (aux,"%ld\n",fi);
  //  number of rows
  fprintf (aux,"%ld\n",n);
  //  attained coefficient of proportionality
  fprintf (aux,"%e\n",blambda);
  //  arc-length
  fprintf (aux,"%e\n",dl);
  //  number of selected dofs
  fprintf (aux,"%ld\n",Mp->nlman.nsdofal);
  //  selected dofs
  for (i=0;i<Mp->nlman.nsdofal;i++){
    fprintf (aux,"%ld\n",Mp->nlman.seldofal[i]);
  }
  //  left hand side vector
  for (i=0;i<n;i++){
    fprintf (aux,"%e\n",r[i]);
  }
  //  proportionality vector
  for (i=0;i<n;i++){
    fprintf (aux,"%e\n",fp[i]);
  }
  //  number of integration points
  fprintf (aux,"%ld\n",Mm->tnip);
  //  integration points
  for (i=0;i<Mm->tnip;i++){
    fprintf (aux,"%ld %ld\n",Mm->ip[i].ncompstr,Mm->ip[i].ncompother);
    for (j=0;j<Mm->ip[i].ncompstr;j++){
      fprintf (aux,"%e\n",Mm->ip[i].strain[j]);
    }
    for (j=0;j<Mm->ip[i].ncompstr;j++){
      fprintf (aux,"%e\n",Mm->ip[i].stress[j]);
    }
    for (j=0;j<Mm->ip[i].ncompother;j++){
      fprintf (aux,"%e\n",Mm->ip[i].eqother[j]);
    }
  }
  
  fclose (aux);
}

void arclopen (long &fi,long &n,double &blambda,double &dl,double *r,double *fp)
{
  long i,j;
  char *file;
  FILE *aux;
  //  aux = fopen ("arcl.backup","r");
  // termit
  file = new char[strlen(Mp->path)+11+1];
  sprintf (file,"%sarcl.backup",Mp->path);
  aux = fopen (file,"r");
  
  //  attained number of steps
  fscanf (aux,"%ld",&fi);
  //  number of rows
  fscanf (aux,"%ld",&n);
  //  attained coefficient of proportionality
  fscanf (aux,"%le",&blambda);
  //  arc-length
  fscanf (aux,"%le",&dl);
  //  number of selected dofs
  fscanf (aux,"%ld",&Mp->nlman.nsdofal);
  //  selected dofs
  for (i=0;i<Mp->nlman.nsdofal;i++){
    fscanf (aux,"%ld",&Mp->nlman.seldofal[i]);
  }
  //  left hand side vector
  for (i=0;i<n;i++){
    fscanf (aux,"%le",&r[i]);
  }
  //  proportionality vector
  for (i=0;i<n;i++){
    fscanf (aux,"%le",&fp[i]);
  }
  //  number of integration points
  fscanf (aux,"%ld",&Mm->tnip);
  //  integration points
  for (i=0;i<Mm->tnip;i++){
    fscanf (aux,"%ld %ld",&Mm->ip[i].ncompstr,&Mm->ip[i].ncompother);
    for (j=0;j<Mm->ip[i].ncompstr;j++){
      fscanf (aux,"%le",&Mm->ip[i].strain[j]);
    }
    for (j=0;j<Mm->ip[i].ncompstr;j++){
      fscanf (aux,"%le",&Mm->ip[i].stress[j]);
    }
    for (j=0;j<Mm->ip[i].ncompother;j++){
      fscanf (aux,"%le",&Mm->ip[i].eqother[j]);
    }
  }
  
  fclose (aux);
}



/**

   created  1.12.2002, Ladislav Svoboda, [e-mail]
 */
void arclsave_adapt (long fi,double lambda,double dlambda,double dl)
{
  char *file;
  FILE *aux;
  
  file = new char[strlen(Mp->path)+19+1];
  //  sprintf (file,"%sarcl.backup.ada.%03ld",Mp->path,Mp->nial);
  sprintf (file,"%sarcl.backup.ada",Mp->path);
  aux = fopen (file,"w");
  
  //  attained number of steps
  fprintf (aux,"%ld\n",fi);
  
  //  attained coefficient of proportionality
  fprintf (aux,"%.16e\n",lambda);
  
  //  
  fprintf (aux,"%.16e\n",dlambda);
  
  //  arc-length
  fprintf (aux,"%.16e\n",dl);
  
  fclose (aux);
}
/**

   created  1.12.2002, Ladislav Svoboda, [e-mail]
 */
void arclopen_adapt (long &fi,double &lambda,double &dlambda,double &dl)
{
  char *file;
  FILE *aux;
  
  file = new char[strlen(Mp->path)+19+1];
  //  sprintf (file,"%sarcl.backup.ada.%03ld",Mp->path,Mp->nial);
  sprintf (file,"%sarcl.backup.ada",Mp->path);
  aux = fopen (file,"r");
  
  //  attained number of steps
  fscanf (aux,"%ld",&fi);
  
  //  attained coefficient of proportionality
  fscanf (aux,"%le",&lambda);
  
  //  
  fscanf (aux,"%le",&dlambda);
  
  //  arc-length
  fscanf (aux,"%le",&dl);

  fclose (aux);
}

/**

   created  1.12.2002, Ladislav Svoboda, [e-mail]
 */
void newmeshgen (int w,long i)
{
  char procname[255];
  
  if (Mespr==1)  fprintf (stdout,"\n\n ***  T3D generator - start  ***\n");
  
  if (Mespr) sprintf (procname,"/home/dr/Bin/T3d -i %sXmesht3d.in -o %smesh.out.%0*ld -m %smesh.bgm.%0*ld -d 100 -p 8 -r 1 -k %ld"
	                      ,Mp->path,Mp->path,w,i,Mp->path,w,i-1,Mt->give_degree(0));
  else       sprintf (procname,"/home/dr/Bin/T3d -i %sXmesht3d.in -o %smesh.out.%0*ld -m %smesh.bgm.%0*ld -d 100 -p 8 -r 1 -k %ld  > /dev/null 2>&1"
	                      ,Mp->path,Mp->path,w,i,Mp->path,w,i-1,Mt->give_degree(0));
  
  if (Mespr==1)  fprintf (stdout,"\n\n ***  T3D generator - end  *****\n");


  system (procname);
  
  
  sprintf (procname,"sed  -e s/out.i/out.%0*ld/g  -e s/dat.i/dat.%d/g  %sprep.pre > %sprep.%0*ld",w,i,0,Mp->path,Mp->path,w,i);
  system (procname);
  
  if (Mespr)
    sprintf (procname,"./mechprep %sprep.%0*ld %ssifel.in.%0*ld",Mp->path,w,i,Mp->path,w,i);
  else
    sprintf (procname,"./mechprep %sprep.%0*ld %ssifel.in.%0*ld  > /dev/null 2>&1",Mp->path,w,i,Mp->path,w,i);

  system (procname);
  

  sprintf (procname,"rm -f %sprep.%0*ld",Mp->path,w,i);
  system (procname);
}

/*
      FILE *ff;
      ff = fopen ("problemnl/aux.txt","a");      if (ff==NULL) fprintf (stderr,"\n ff file has not been opened.\n\n");
      for (long g=0;g<Ndofn;g++)
	fprintf (ff,"%40.20lf\n",Lsrs->lhs[i]);
      fclose (ff);
      
      


*/

/**
   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

   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
   
   23.12.2004
*/
long arclengthadapt (long lcid,long adaptcontrol)
{
  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;
  
  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
  ra = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  
  if (Mp->nlman.hdbackupal == 2){            // termit
    arclopen (li,n,lambda,dl,ra,fp);
  }
  else if (adaptcontrol >= 2){
    arclopen_adapt (li,lambda,blambda,dl);
  }
  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];
    f[j]=fp[j];
  }
  if (li)
    print_init(-1, "at");
  else
  {
    print_init(-1, "wt");
    print_step(lcid, li, lambda, fa);
  }
  
  
  //  elastic trial response of structure
  Mp->phase=1;
  //  stiffness matrix based on elastic models
  stiffness_matrix (lcid);
  //  solution of K.r=f
  //Smat->solve_system (Gtm,ra,f);
  Mp->ssle.solve_system (Gtm,Smat,ra,f,Out);

  computestresses (0);

  //  switch to the special material models
  Mp->phase=2;
  
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  for (i=li;i<ni;i++){
    stop=1; // termit
    //  backup of left hand side vector
    copyv (ra, r, n);
    //  backup of reached lambda parameter
    blambda=lambda;

    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == li))
      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
    //Mm->computestresses (0);
    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;
      Mm->updateipval ();
      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);
	
	//  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;
	*/
	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;
	
	//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;
	  Mm->updateipval ();
	  stop=1; 
          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;
    
    if (Mp->adaptflag & 16 && i%10 && Mm->plast)
      continue;
    
    newmesh = 0;
    if (Mp->adaptivity && Mp->nlman.nial > (i+1) && !( (Mp->adaptflag & 16) && (Mp->nlman.nial==100 || Mp->nlman.nial==1))){
      
      // print of fxlambda second part is in adaptivity.cpp
      //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 ((Mp->adaptflag & 16) && (Mp->nlman.nial==100 || Mp->nlman.nial==1))
      Ada->run (0,0,0);  // only for testing
    else
      if (i < ni)
	arclsave_adapt (i,lambda,0,dl);
  
  //if (Mp->adaptivity)   arclsave_adapt (i,blambda,ddlambda,dl);
  // timret
  print_close();
  delete [] r;		    
  delete [] fi;  
  delete [] fa;  
  delete [] f;
  delete [] v;  
  delete [] u;  
  delete [] ddr;
  
  return (i);
}