nssolver.cpp

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

double displ_controlrv (long lcid, double inilambda, double dmplambda, double rlambda, double errl, long incrdispl)
/**  
Function solves system of nonlinear algebraic
equations by Newton-Raphson method solved system does not contain time variable.
Displacement control is assumed and only required value of lambda is reached 
(initial lambda value is assumed for global Mp->lambda and rlambda is  ).
@param lcid      - required load case id
@param inilambda - initial value of lambda for displacements
@param rlambda   - required value of lambda
@param errl      - required error for reaching rlambda
@incrdispl       - indicator for increasing Mp->lambda(1=yes/0=no)

    9.11.2005
    TKo.
*/
{
  long i,j,k,n,ni,ini;
  double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
  double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp, *fir2;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nlman.ninr;
  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilnr;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  err = Mp->nlman.errnr;
  //  increment size
  dlambda=Mp->nlman.incrnr;
  //  minimum increment size
  dlambdamin=Mp->nlman.minincrnr;

  rb = new double [n];
  f  = new double [n];
  fb = new double [n];
  fi = new double [n];
  fir2 = new double[n];
  dr = new double [n];
  memset (rb,0,n*sizeof(double));
  memset (f,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (dr,0,n*sizeof(double));
  
  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  lambda=0.0;
  Mp->lambda = inilambda;
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  print_init(-1, "wt");
  print_step(lcid, 0, 0.0, f);
  print_flush();
  for (i=0;i<ni;i++){
    
    // checking required value of lambda
    if (fabs(lambda-rlambda) < errl)
      break;
    if ((lambda + dlambda) > rlambda)
      dlambda = rlambda - lambda;

    //  backup of left hand side vector
    for (j=0;j<n;j++){
      rb[j]=r[j];
    }
    //  backup of reached lambda parameter
    blambda=lambda;
    if (incrdispl)
      Mp->lambda = inilambda+dmplambda*lambda;
   
    fprintf (stdout,"\n increment %ld   lambda %e,    dlambda %e",i,lambda,dlambda);
    
    //  vector of maximum load and vector of load increment
    if (incrdispl)
      Mp->lambda += dmplambda*dlambda;
    internal_forces(lcid, fir2);

    for (j=0;j<n;j++){
      f[j]=fc[j]+lambda*fp[j];
      fb[j]= f[j] + dlambda*fp[j] - fir2[j];
    }
    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == 0))
      stiffness_matrix (lcid);
    
    //  solution of K(r).v=F
    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    
    //  computation of internal forces
    internal_forces (lcid,fi);
    
    //  vector of unbalanced forces
    for (j=0;j<n;j++){
      fb[j]=f[j]+dlambda*fp[j];
    }
    norfa=ss(fb,fb,n);
    for (j=0;j<n;j++){
      fb[j] -= fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    if (norfa != 0.0)
      norf /= norfa;
    
    if (norf<err){
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
      continue;
    }
    
    //  internal iteration loop
    for (j=0;j<ini;j++){
      
      //  solution of K(r).v=F
      //Smat->solve_system (Gtm,dr,fb);
      Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

      for (k=0;k<n;k++){
	r[k]+=dr[k];
      }
      
      //  computation of internal forces
      internal_forces (lcid,fi);
      
      //  vector of unbalanced forces
      for (k=0;k<n;k++){
	fb[k]=f[k]+dlambda*fp[k]-fi[k];
      }
      
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,n);
      if (norfa != 0.0)
        norf /= norfa;
      
      fprintf (stdout,"\n increment %ld   inner loop j %ld     norf=%e  dlambda %e",i,j,norf,dlambda);
      
      if (norf<err)  break;
    }
    
    if (j==ini || norf>err)
    {
      if (dlambda == dlambdamin)
      {
        if (Mespr==1)  fprintf (stdout,"\n increment cannot be decreased dlambda=%e    norf=%e",dlambda,norf);
        break;
      }
      dlambda/=2.0;
      if (dlambda<dlambdamin)  dlambda=dlambdamin;
      if (Mespr==1)  fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
      
      for (j=0;j<n;j++)
      	r[j]=rb[j];

      lambda=blambda;     
    }
    else
    {
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
    }
  }
  print_close();
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f;
  return (lambda);
}

void newton_raphson_restart (long lcid)
  //  function solves system of nonlinear algebraic
  //  equations by Newton-Raphson method
  //  solved system does not contain time variable
  //
  //  16.8.2001
{
  long i,j,k,n,ni,ini,nii;
  double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
  double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;


  //  number of rows of the matrix
  n = Ndofm;
  //  maximum number of increments
  ni = Mp->nlman.ninr;
  // number of initial increments
  nii = Mp->nlman.nienr;
  //  maximum number of iterations in one increment
  ini = Mp->nlman.niilnr;
  //  computer zero
  zero = Mp->zero;
  //  required error in inner loop
  err = Mp->nlman.errnr;
  //  increment size
  dlambda=Mp->nlman.incrnr;
  //  minimum increment size
  dlambdamin=Mp->nlman.minincrnr;

  rb = new double [n];
  f  = new double [n];
  fb = new double [n];
  fi = new double [n];
  dr = new double [n];
  memset (rb,0,n*sizeof(double));
  memset (f,0,n*sizeof(double));
  memset (fi,0,n*sizeof(double));
  memset (dr,0,n*sizeof(double));
  
  //  initialization phase
  r  = Lsrs->give_lhs (lcid);
  fp = Lsrs->give_rhs (lcid*2);
  fc = Lsrs->give_rhs (lcid*2+1);
  lambda=0.0;




  long ncompstr;
  double tmp;
  FILE *aux;

  // if restart from backup file is required
  if (Mp->nlman.hdbr)
  {
    aux = fopen (Mp->nlman.backupfname, "rt");
    if (aux == NULL)
    {
      fprintf(stderr, "\n\nError - unable to open backup file %s\n", Mp->nlman.backupfname);
      fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
      fprintf(stderr, "Program will be terminated\n");
      abort();
    }
    // searching required backup
    fprintf(stdout, "\n");
    for (k=0; k<Mp->nlman.hdbid; k++)
    {
      fprintf(stdout, "\r reading backup record %ld", k+1);
      fflush(stdout);
      mtsolver_restore(aux, r, fc, i, tmp, j, 39, 8);
    }
    fprintf(stdout, "\n backup time : %e\n", tmp);
    if (j != Ndofm)
    {
      fprintf(stderr, "\n\nError - incorrect number of DOFs is stored in backup file\n");
      fprintf(stderr, "file %s, line %d\n", __FILE__, __LINE__);
    }
/*
  Begining of Hinkley modification
*/
    for (i=0; i<Mm->tnip; i++)
    {
      if (Mm->ip[i].ncompeqother >= 8)
      {
        ncompstr = Mm->ip[i].ncompstr;
        for (j=0; j<ncompstr; j++)
          Mm->ip[i].stress[j] = Mm->ip[i].eqother[ncompstr+j];
      }
    }
  }

  Mp->strcomp = 0;
  internal_forces(lcid, fc);
  for (i=0; i < Ndofm; i++)
    fc[i] *= 1.0e+6;
  Mp->strcomp = 1;
/*
 End of Hinkley modification
*/


  // compute and print values initial values
  compute_req_val (lcid);
  print_init(-1, "wt");
  print_step(lcid, 0, 0, fc);
  print_flush();


  //  assembling of tangent stiffness matrix
  stiffness_matrix (lcid);
  // norm of load vector
  norfa=ss(fc,fc,n);
  //  internal iteration loop
  for (j=0;j<nii;j++)
  {
    //  computation of internal forces
    internal_forces (lcid,fi);
    //  vector of unbalanced forces
    for (k=0;k<n;k++)
      fb[k]=fc[k]-fi[k];
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    norf /= norfa;
    
    if (norf<err)
    {
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, j+1, -0.0000001, fc);
      print_flush();
      break;
    }
    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

/*
    double tmp1, tmp2, tmp3, tmp4;
    long ncomp;
    tmp1=tmp2=tmp3=tmp4=0.0;
    fprintf(Out, "\n\n*Krok %ld\n", j);
    fprintf(Out,"\n\nn\tfc\tfi\tfc-fi\tr\tr+dr\tddr\n");
    for(k=0; k<n; k++)
    {
      tmp1+=fc[k]*fc[k];
      tmp2+=fi[k]*fi[k];
      tmp3+=r[k]*r[k];
      r[k]+=dr[k];
      tmp4+=r[k]*r[k];
      if (k%100==0)
        fprintf(Out, "%ld\t%e\t%e\t%e\t%e\t%e\t%e\n",k, tmp1, tmp2, tmp1-tmp2, tmp3, tmp4, tmp3-tmp4);
    }
    fprintf(Out, "%ld\t%e\t%e\t%e\t%e\t%e\t%e\n",k, tmp1, tmp2, tmp1-tmp2, tmp3, tmp4, tmp3-tmp4);
    fprintf(Out,"\n\nIpp\t omega\n");
   
    for (k=0; k < Mm->tnip; k++)
    {
      if (Mm->ip[k].ncompother > 8)
      {
        ncomp = 2*Mm->ip[k].ncompstr;
        fprintf(Out, "%ld\t%e\t%e\n", k, Mm->ip[k].other[ncomp],Mm->ip[k].other[ncomp+1]);
      }
    }
    fflush(Out);*/
    for(k=0; k<n; k++)
      r[k]+=dr[k];
/*    swap_other();
    compute_req_val (lcid);
    print_step(lcid, j+1, j+1, fc);
    print_flush();
    swap_other();*/
    fprintf (stdout,"\n Equilibrium interation loop j %ld     norf=%e", j, norf);
  }
  if ((nii > 0) && (norf > err))
  {
    fprintf(stdout, "\n Initial equilibrium could not be reached\n");
    fprintf(stdout, " Program will be terminated\n");
    print_close();
    return;
  }
  else
    fprintf(stdout, "\n\n Starting incremental loading . . .\n");


  // ***************************
  //  main iteration loop  ****
  // ***************************
  for (i=0;i<ni;i++){
    
    //  backup of left hand side vector
    for (j=0;j<n;j++){
      rb[j]=r[j];
    }
    //  backup of reached lambda parameter
    blambda=lambda;
   
    fprintf (stdout,"\n increment %ld   lambda %e,    dlambda %e",i,lambda,dlambda);
    
    //  vector of maximum load and vector of load increment
    for (j=0;j<n;j++){
      f[j]=fc[j]+lambda*fp[j];
      fb[j]=dlambda*fp[j];
    }
    
    //  assembling of tangent stiffness matrix
    if ((Mp->stmat==tangent_stiff) || (i == 0))
      stiffness_matrix (lcid);

    //  solution of K(r).v=F

    //Smat->solve_system (Gtm,dr,fb);
    Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

    for (j=0;j<n;j++){
      r[j]+=dr[j];
    }
    //  computation of internal forces
    internal_forces (lcid,fi);
    //  vector of unbalanced forces
    for (j=0;j<n;j++){
      fb[j]=f[j]+dlambda*fp[j];
    }
    norfa=ss(fb,fb,n);
    for (j=0;j<n;j++){
      fb[j] -= fi[j];
    }
    //  norm of vector of unbalanced forces
    norf=ss(fb,fb,n);
    norf /= norfa;
    
    
    
//    fprintf (Out,"\n norf=%e",norf);

    if (norf<err){
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
      continue;
    }
    
    //  internal iteration loop
    for (j=0;j<ini;j++){
      
      //  solution of K(r).v=F
      //Smat->solve_system (Gtm,dr,fb);
      Mp->ssle.solve_system (Gtm,Smat,dr,fb,Out);

      for (k=0;k<n;k++){
	r[k]+=dr[k];
      }
      
      //  computation of internal forces
      internal_forces (lcid,fi);
      
      //  vector of unbalanced forces
      for (k=0;k<n;k++){
	fb[k]=f[k]+dlambda*fp[k]-fi[k];
      }
      
      //  norm of vector of unbalanced forces
      norf=ss(fb,fb,n);
      norf /= norfa;
      
      fprintf (stdout,"\n increment %ld   inner loop j %ld     norf=%e  dlambda %e",i,j,norf,dlambda);
//      fprintf (Out,"\n j=%ld     norf=%e",j,norf);
      
      if (norf<err)  break;
    }
    
    if (j==ini || norf>err){
      dlambda/=2.0;
      if (dlambda<dlambdamin)  dlambda=dlambdamin;

      if (Mespr==1){
	fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
//	fprintf (Out,"\n increment must be modified (dlambda decrease) dlambda=%e    norf=%e",dlambda,norf);
      }
      
      for (j=0;j<n;j++){
      	r[j]=rb[j];
      }
      lambda=blambda;
    }
    else{
      lambda+=dlambda;
      Mm->updateipval ();
      compute_req_val (lcid);
      print_step(lcid, i+1, lambda, f);
      print_flush();
    }
    
  }
  
  delete [] dr;  delete [] fi;  delete [] fb;  delete [] f;
  print_close();
}

void seldofinit ()
{
  long i,j,k,l,ndofn;
  double x1,x2,y1,y2,z1,z2,length;
  
  switch (Mp->nlman.displnorm){
  case alldofs:{  break; }
  case seldofs:
  case seldofscoord:{
    for (i=0;i<Mp->nlman.nsdofal;i++){
      j=Mp->nlman.seldofal[i];
      Mp->nlman.seldofal[i]=Mt->give_dof (Mp->nlman.selnodal[i],j)-1;
    }
    break;
  }
  case selecnodes:{
    Mp->nlman.nsdofal=0;
    for (i=0;i<Mp->nlman.nsnal;i++){
      Mp->nlman.nsdofal+=Mt->give_ndofn (Mp->nlman.selnodal[i]);
    }
    Mp->nlman.seldofal = new long [Mp->nlman.nsdofal];
    k=0;
    for (i=0;i<Mp->nlman.nsnal;i++){
      ndofn=Mt->give_ndofn (Mp->nlman.selnodal[i]);
      for (j=0;j<ndofn;j++){
	l=Mt->give_dof (Mp->nlman.selnodal[i],j);
	if (l>0){
	  Mp->nlman.seldofal[k]=l-1;
	  k++;
	}
      }
    }
    Mp->nlman.nsdofal=k;
    break;
  }
  case nodesdistincr:{
    Mp->nlman.nsdofal=0;
    for (i=0;i<Mp->nlman.nsnal;i++){
      Mp->nlman.nsdofal+=Mt->give_ndofn (Mp->nlman.selnodal[i]);
    }
    Mp->nlman.seldofal = new long [Mp->nlman.nsdofal];
    k=0;
    for (i=0;i<Mp->nlman.nsnal;i++){
      ndofn=Mt->give_ndofn (Mp->nlman.selnodal[i]);
      for (j=0;j<ndofn;j++){
	Mp->nlman.seldofal[k]=Mt->give_dof (Mp->nlman.selnodal[i],j);
	k++;
      }
    }
    x1=Gtm->gnodes[Mp->nlman.selnodal[0]].x;    x2=Gtm->gnodes[Mp->nlman.selnodal[1]].x;
    y1=Gtm->gnodes[Mp->nlman.selnodal[0]].y;    y2=Gtm->gnodes[Mp->nlman.selnodal[1]].y;
    z1=Gtm->gnodes[Mp->nlman.selnodal[0]].z;    z2=Gtm->gnodes[Mp->nlman.selnodal[1]].z;
    length=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
    Mp->nlman.nxal=(x2-x1)/length;  Mp->nlman.nyal=(y2-y1)/length;  Mp->nlman.nzal=(z2-z1)/length;
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown displacement norm is required in function seldofinit (%s, line %d)",__FILE__,__LINE__);
  }
  }
}

double displincr (double *displincr,long n)
{
  long i;
  double norm,u1,u2,v1,v2,w1,w2;
  
  switch (Mp->nlman.displnorm){
  case alldofs:{
    norm = ss (displincr,displincr,n);
    break;
  }
  case seldofs:
  case seldofscoord:
  case selecnodes:{
    norm=0.0;
    for (i=0;i<Mp->nlman.nsdofal;i++){
      norm+=displincr[Mp->nlman.seldofal[i]]*displincr[Mp->nlman.seldofal[i]];
    }
    break;
  }
  case nodesdistincr:{
    if (Mp->nlman.probdimal==2){
      u1=0.0;  u2=0.0;  v1=0.0;  v2=0.0;
      if (Mp->nlman.seldofal[0]>0)  u1=displincr[Mp->nlman.seldofal[0]-1];
      if (Mp->nlman.seldofal[1]>0)  v1=displincr[Mp->nlman.seldofal[1]-1];
      if (Mp->nlman.seldofal[2]>0)  u2=displincr[Mp->nlman.seldofal[2]-1];
      if (Mp->nlman.seldofal[3]>0)  v2=displincr[Mp->nlman.seldofal[3]-1];
      
      norm=(u2-u1)*Mp->nlman.nxal + (v2-v1)*Mp->nlman.nyal;
    }
    if (Mp->nlman.probdimal==3){
      u1=0.0;  u2=0.0;  v1=0.0;  v2=0.0;  w1=0.0;  w2=0.0;
      if (Mp->nlman.seldofal[0]>0)  u1=displincr[Mp->nlman.seldofal[0]-1];
      if (Mp->nlman.seldofal[1]>0)  v1=displincr[Mp->nlman.seldofal[1]-1];
      if (Mp->nlman.seldofal[2]>0)  w1=displincr[Mp->nlman.seldofal[2]-1];
      if (Mp->nlman.seldofal[3]>0)  u2=displincr[Mp->nlman.seldofal[3]-1];
      if (Mp->nlman.seldofal[4]>0)  v2=displincr[Mp->nlman.seldofal[4]-1];
      if (Mp->nlman.seldofal[5]>0)  w2=displincr[Mp->nlman.seldofal[5]-1];
      
      norm=(u2-u1)*Mp->nlman.nxal + (v2-v1)*Mp->nlman.nyal + (w2-w1)*Mp->nlman.nzal;
    }
    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__);
  }
  }
  norm=sqrt (norm);
  return norm;
}


void quadeqcoeff (double *ddr,double *v,long n,double ddlambda,double psi,double norfp,double dl,