nssolver_old.cpp

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

#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "nssolver.h"
#include "slopesol.h"
#include "global.h"
#include "globmat.h"
#include "mechprint.h"
#include "gmatrix.h"
#include "gtopology.h"
#include "loadcase.h"
#include "dloadcase.h"
#include "intpoints.h"
#include "auxdatout.h"
#include "mechprint.h"
#include "j2flow.h" // pak vyhodit
#include "chen.h" // pak vyhodit
#include "hardsoft.h" // pak vyhodit
#include "mathem.h"
#include "vector.h"
#include "mtsolver.h"

#include "adaptivity.h"
#include "meshtransfer.h"

/**
   function solves non-linear statics problems
   
   JK
*/
void solve_nonlinear_statics ()
{
  long i;
  double *rhs;

  Mp->lambda = 1.0;  
  rhs=Lsrs->give_rhs (0);
  
  for (i=0;i<Lsrs->nlc;i++){
    
    //  vectors of loading
    //Mb->lc[i*2+0].assemble (0,rhs+(i*2+0)*Ndofm);
    //Mb->lc[i*2+1].assemble (0,rhs+(i*2+1)*Ndofm);
    mefel_right_hand_side (i,rhs);
    
    if (!Mp->adaptivity)
      nonlinear_solver (i);
    else
      nonlinear_solver_adaptiv (i);

 }
  
}

void print_output(long lcid)
{
/*
  //  strains
  if (Mp->straincomp==1){
    Mm->computestrains (lcid);
    Mm->stra.transformvalues (1);
  }
    
  //  stresses
  if (Mp->stresscomp==1){
    Mm->computenlstresses (lcid);
    Mm->stre.transformvalues(0);
  }
  
  //  reactions
  if (Mp->reactcomp==1)
    Mb->lc[lcid].compute_reactions (lcid);
  
    
  //  output data
  Lsrs->output (Out,lcid);    
  // other at the integration points  
  if ((Mp->straincomp==1) || (Mp->stresscomp==1))
    print_other (Out,lcid);
*/
}

void nonlinear_solver (long lcid)
{
  switch (Mp->nlman.tnlinsol){
  case arcl:{
    arclength (lcid,0);
    break;
  }
  case newton:{
    newton_raphson (lcid);
    break;
  }
  case arclrv1:{
    arclengthrv1 (lcid,1.1,1.0e-3);
    break;
  }
  case newtonrv1:{
    newton_raphsonrv1 (lcid,Mp->nlman.rvlam,1.0e-3);
    break;
  }
  case newtonrestart:{
    newton_raphson_restart (lcid);
    break;
  }
  case displctrl:{
    displ_control (lcid);
    break;
  }
  case adaptram:{
    arclengthadapt (lcid,0);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
    fprintf (stderr,"\n in function nonlinear_solver (%s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

void refresh_plast ()
{
  for (long i=0;i<Mm->tnip;i++)
    if ( Mm->j2f[Mm->ip[i].idm[0]].give_consparam (i,0) ){
      Mm->plast = 1;
      return;
    }
}

/**
   if adaptcontrol =
   1 - arclength nyni bezi poprve
   2 - arclength nyni nebezi poprve
   
   created  1.12.2002, Ladislav Svoboda, [e-mail]
 */
void nonlinear_solver_adaptiv (long lcid)
{
  int width;
  long i,adaptcontrol;
  char file[255];
  
  switch (Mp->nlman.tnlinsol){
  case arcl:{
    if (Mespr==1){
      fprintf (stdout,"\n\n ----------------------------");
      fprintf (stdout,  "\n ***  Arc-length - start  ***\n");
    }
    
    width = 1 + (long)log10((double)Mp->nlman.nial);
    adaptcontrol = 1;
    
    for (i=0; i < Mp->nlman.nial ;){
      // i = number of performed steps
      i = arclength (lcid,adaptcontrol);
      
      if ( i<Mp->nlman.nial || Mp->nlman.nial==1 && (Mp->adaptflag & 16) ){
	if (Mespr==1)  fprintf (stdout,"\n\n NEW MESH will be generated\n");
	
	newmeshgen (width,i);
	sprintf (file,"%s%s.in.%0*ld",Mp->path,Mp->filename,width,i);
	newmeshread (file,lcid);
	
	if ( i<Mp->nlman.nial ){
	  if (Mm->plast == 0)
	    refresh_plast ();
	  
	  arclength15  (lcid);
	}
	else
	  Ada->run (4,width,i); // only for testing
      }
      
      //if (!(Mp->adaptflag & 16)){  sprintf (file,"%s%s.dx.%ld.1",Mp->path,Mp->filename,i);  print_default_2_dx (Gtm,Mp,Mt,0,file);  }
      
      adaptcontrol++;
    }
    
    if (Mespr==1){
      fprintf (stdout,"\n\n ***   Arc-length is finished  ***");
      fprintf (stdout,  "\n *********************************");
    }
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
    fprintf (stderr,"\n in function nonlinear_solver_adaptiv (%s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   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
   
   25.7.2001
*/
long arclength (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];
  
  //  pokus
  //double *zaloha,*zaloha2;
  //zaloha = new double [n];
  //zaloha2 = new double [n];
  //for (i=0;i<n;i++){
  //zaloha[i]=0.0;
  //zaloha2[i]=0.0;
  //}
  // konec pokusu
  
  //  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;
  }
  
  //refresh_plast ();  // termit
  
  //  Mp->ns;

  //  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);
  }
  
  //fprintf (Out,"\n\n osel   epspuc %lf   sc %lf\n",Mm->chplast[0].hs.epspuc,Mm->chplast[0].hs.sc);
  
  // ***************************
  //  main iteration loop  ****
  // ***************************
  for (i=li;i<ni;i++){

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

    //if (i>250)
    //Mp->phase=1;

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

    //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);
      Smat->printmat (Out);
    }
    
    //  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 (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;
    
    //  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);

	fprintf (Out,"\n arc-length  vnitrni smycka  %ld %ld",i,j);
	
	//  back substitution
	Smat->solve_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);
	/*
	  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;

	
	/*
	k=n-1;
	if (ra[k]+u[k]+l1*v[k] > ra[k]+u[k]+l2*v[k])  dlambda=l1;
	else  dlambda=l2;
	*/
	/*
	if (Mp->phase==1){
	  k=n-1;
	  if (ra[k]+u[k]+l1*v[k] > ra[k]+u[k]+l2*v[k])  dlambda=l1;
	  else  dlambda=l2;
	}
	else{
	  if (fabs(l1)>fabs(l2))  dlambda=l2;
	  else dlambda=l1;
	}
	fprintf (Out,"\n\n arc-length   %le  %le     %le %le   %le\n",ra[k]+u[k]+l1*v[k],ra[k]+u[k]+l2*v[k],l1,l2,dlambda);
	*/
	
	//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 (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;
      }
    }
    
    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;