ndsolver.cpp

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

#include "ndsolver.h"
#include "global.h"
#include "globmat.h"
#include "gmatrix.h"
#include "loadcase.h"
#include "dloadcase.h"
#include "auxdatout.h"
#include "mechprint.h"
#include "elemswitch.h"

void solve_nonlinear_dynamics ()
{
  switch (Mp->tforvib){
  case newmark:{
    nonlin_newmark_method ();
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown solver of nonlinear dynamics is required");
    fprintf (stderr,"\n in function solve_nonlinear_dynamics () (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

/**
   Nonlinear Newmark method for solution of nonlinear dynamic problems
   the method is implicit, equilibrium is solved by the Newton-Raphson method
   
   tam co byla alpha je nyni beta
   tam co byla delta je nyni gamma
   
   5.3.2007, JK
*/
void nonlin_newmark_method ()
{
  long i,j,n;
  double beta,gamma,zero,time,dt,end_time;
  double *a,*v,*p,*dd,*vv,*lhs,*rhs,*fs,*brhs;
  
  //  number of degrees of freedom
  n = Ndofm;
  //  computer zero
  zero=Mp->zero;
  //  coefficient of the method
  beta=Mp->alphafvn;
  //  coefficient of the method
  gamma=Mp->deltafvn;
  
  //  initial time
  time=Mp->timecon.starttime ();
  //  end time
  end_time = Mp->timecon.endtime ();
  //  first time increment
  dt = Mp->timecon.initialtimeincr ();
  
  //  load case id
  //  this solver enables only one load case with respect to nonlinearity
  lcid=0;
  
  //  nodal values - displacements
  //  (assigment of initial values)
  lhs = Lsrs->give_lhs (lcid);
  //  time derivatives of nodal values - velocities
  //  (assigment of initial values)
  v = Lsrs->give_tdlhs (lcid);
  //  vector of dynamic load
  rhs = Lsrs->give_rhs (2*lcid);
  //  vector of static load
  fs = Lsrs->give_rhs (2*lcid+1);
  //  backup of the right hand side
  brhs = new double [n];

  //  vector of acceleration
  a = Lsrs->give_stdlhs (lcid);
  nullv (a,n);
  //  predictor of displacements
  dd = new double [n];
  //  predictor of velocities
  vv = new double [n];
  //  auxiliary vector
  p = new double [n];
  
  //  assembling of static loads  
  nullv (fs+lcid*n, n);
  Mb->lc[lcid].assemble (lcid,fs,1.0);
  
  //  number of step
  i=0;
  print_init(-1, "wt");  
  print_step(lcid, i, Mp->time, rhs);
  print_flush();
  
  
  //  time initialization
  Mp->time = Mp->timecon.newtime ();
  //  time increment
  dt = Mp->timecon.actualbacktimeincr ();
  
  
  //  loop over time interval
  do{
    
    //  stiffness matrix
    stiffness_matrix (lcid);
    //  mass matrix
    mass_matrix (lcid);
    //  damping matrix
    damping_matrix ();
    
    if (Mespr==1)  fprintf (stdout,"\n                                time  %f",Mp->time);
    
    //  assembling of dynamic loads
    nullv (rhs+lcid*n, n);
    Mb->dlc[lcid].assemble (lcid,rhs,n,Mp->time);
    
    //  total load
    addv (rhs,fs,n);
    
    //  modification of vector of prescribed nodal forces
    cmulv (dt*dt*beta,rhs,n);
    //  copy of the actual prescribed forces
    copyv (rhs,brhs,n);

    
    //  predictor vectors
    for (j=0;j<n;j++){
      dd[j] = lhs[j] + dt*v[j] + dt*dt*(0.5-beta)*a[j];
      vv[j] = v[j] + dt*(1.0-gamma)*a[j];
    }
    
    //  M.dd
    Mmat->gmxv (dd,p);
    
    //  contribution to the right hand side
    addv (rhs,p,n);
    
    //  C.dd
    Dmat->gmxv (dd,p);
    cmulv (dt*gamma,p,n);
    addv (rhs,p,n);
    
    //  C.vv
    Dmat->gmxv (vv,p);
    cmulv (-1.0*dt*dt*beta,p,n);
    addv (rhs,p,n);
    
    //  dynamic stiffness matrix
    //  (stored in array for stiffness matrix)
    //  M + gamma.dt.C + beta.dt.dt.K
    Smat->scalgm (dt*dt*beta);
    Smat->addgm (1.0,*Mmat);
    Smat->addgm (dt*gamma,*Dmat);
    
    
    //  solution of system of algebraic equations
    Smat->solve_system (Gtm,lhs,rhs);
    
    for (j=0;j<n;j++){
      a[j] = 1.0/dt/dt/beta*(lhs[j]-dd[j]);
      v[j] = vv[j] + gamma*dt*a[j];
    }
    
    
    //  restoration of the prescribed forces
    copyv (brhs,rhs,n);
    
    
    //  inertial forces
    //  M.a
    Mmat->gmxv (a,p);
    cmulv (-1.0*dt*dt*beta,p,n);
    addv (rhs,p,n);
    
    //  damping forces
    //  C.v
    Dmat->gmxv (v,p);
    cmulv (-1.0*dt*dt*beta,p,n);
    addv (rhs,p,n);
    
    //  internal forces
    internal_forces (lcid,p);
    cmulv (-1.0*dt*dt*beta,p,n);
    addv (rhs,p,n);
    
    //  norm of the residual vector
    norres = ss (rhs,rhs,n);
    fprintf (stdout,"\n norm of residual vector  %le",norres);
    if (norres>err){
      //  norm of residual vector is greater than required tolerance
      //  the Newton-Raphson method has to be used for reaching equilibrium
      for (j=0;j<ni;j++){
	
	//  solution of system of algebraic equations
	Smat->solve_system (Gtm,p,rhs);
	
	for (k=0;k<n;k++){
	  lhs[k]+=p[k];
	  a[k] = 1.0/dt/dt/beta*(lhs[k]-dd[k]);
	  v[k] = vv[j] + gamma*dt*a[k];
	}

	
	
      }
      
      
    }
    
    
    
    compute_req_val (lcid);
    print_step(lcid, i, Mp->time, rhs);
    print_flush();
    
    //  time increment
    Mp->time = Mp->timecon.newtime ();
    dt = Mp->timecon.actualbacktimeincr ();
    //  step number increment
    i++;

  }while(Mp->time<end_time);

  print_close();  
  delete [] p;
  delete [] dd;
  delete [] vv;
}


/**
   difference method for solution of forced dynamic problems
   the method is explicit
   
   @param lcid - load case id
   
   JK, 26.7.2005
*/
void difference_method (long lcid)
{
  long i,j,n;
  double zero,time,dt,end_time;
  double *a,*v,*p,*da,*dp,*lhs,*rhs,*fs;
  
  //  number of degrees of freedom
  n = Ndofm;
  //  computer zero
  zero=Mp->zero;
  
  //  initial time
  time=Mp->timecon.starttime ();
  //  end time
  end_time = Mp->timecon.endtime ();
  //  first time increment
  dt = Mp->timecon.initialtimeincr ();
  
  //  nodal values - displacements at next step
  //  (assigment of initial values)
  lhs = Lsrs->give_lhs (lcid);
  //  time derivatives of nodal values - velocities
  //  (assigment of initial values)
  v = Lsrs->give_tdlhs (lcid);
  //  second time derivatives of nodal values - accelerations
  a = Lsrs->give_stdlhs (lcid);
  //  vector of dynamic load
  rhs = Lsrs->give_rhs (2*lcid);
  //  vector of static load
  fs = Lsrs->give_rhs (2*lcid+1);
  

  //  vector of displacements at previous step
  dp = new double [n];
  //  vector of displacements at actual step
  da = new double [n];
  //  auxiliary vector
  p = new double [n];
  
  //  assembling of static loads  
  nullv (fs+lcid*n, n);
  Mb->lc[lcid].assemble (lcid,fs,1.0);
  

  i=0;
  print_init(-1, "wt");  
  print_step(lcid, i, Mp->time, rhs);
  print_flush();
  
  
  //  initial values
  for (j=0;j<n;j++){
    da[j] = lhs[j];
    dp[j] = da[j] - dt*v[j];
  }

  
  //  loop over time interval
  do{
    
    //  stiffness matrix
    stiffness_matrix (lcid);
    //  mass matrix
    mass_matrix (lcid);
    //  damping matrix
    damping_matrix ();
    
    if (Mespr==1)  fprintf (stdout,"\n                                time  %f",Mp->time);
    
    //  assembling of dynamic loads
    nullv (rhs+lcid*n, n);
    Mb->dlc[lcid].assemble (lcid,rhs,Ndofm,Mp->time);
    
    //  total load
    addv (rhs,fs,n);
    
    //  K.da
    Smat->gmxv (da,p);
    cmulv (-1.0,p,n);
    addv (rhs,p,n);
    
    //  M.da.2/dt/dt
    Mmat->gmxv (da,p);
    cmulv (2.0/dt/dt,p,n);
    addv (rhs,p,n);
    
    //  M.dp/dt/dt
    Mmat->gmxv (dp,p);
    cmulv (-1.0/dt/dt,p,n);
    addv (rhs,p,n);
    
    //  C.vv
    Dmat->gmxv (dp,p);
    cmulv (1.0/2.0/dt,p,n);
    addv (rhs,p,n);
    
    //  dynamic stiffness matrix
    //  (stored in array for mass matrix)
    //  M/dt/dt + C/2/dt
    Mmat->scalgm (1.0/dt/dt);
    Mmat->addgm (1.0/2.0/dt,*Dmat);
    
    //  solution of system of algebraic equations
    Mmat->solve_system (Gtm,lhs,rhs);
    
    for (j=0;j<n;j++){
      v[j]=(lhs[j]-dp[j])/2.0/dt;
      a[j]=(lhs[j]-2.0*da[j]+dp[j])/dt/dt;
      
      dp[j]=da[j];
      da[j]=lhs[j];
    }
    
    //  time increment
    Mp->time = Mp->timecon.newtime ();
    dt = Mp->timecon.actualbacktimeincr ();
    i++;

    compute_req_val (lcid);
    print_step(lcid, i, Mp->time, rhs);
    print_flush();
    

  }while(Mp->time<end_time);

  print_close();  
  delete [] da;
  delete [] dp;
  delete [] p;
}




/**
   function solves system of ordinary differential equations of second order
   such equations are used e.g. in dynamics
   function is intended for molecular dynamics
   
   @param lcid - load case id
   
   JK, 24.7.2005
*/
void verlet_method (long lcid)
{
  long i,n,ni;
  double mm,zero,dt,end_time;
  double *m,*dp,*da,*dn,*fp,*fi;
  
  //  number of rows of the matrix
  n = Ndofm;
  //  number of iterations
  ni = 0;
  //  computer zero
  zero = Mp->zero;
  
  //  vector of displacements at previous step
  dp = new double [n];
  //  vector of displacements at new step
  dn = new double [n];
  //  vector of displacements at actual step
  da = Lsrs->give_lhs (lcid);
  //  vector of prescribed forces
  fp = Lsrs->give_rhs (lcid);
  //  vector of internal forces
  fi = new double [n];
  //  vector of reciprocal values of weights
  m = new double [n];
  
  //  initial time
  Mp->time=Mp->timecon.starttime ();
  //  end time
  end_time = Mp->timecon.endtime ();
  
  //  extraction of weights
  for (i=0;i<n;i++){
    mm=Mmat->give_entry (i,i);
    if (mm<zero){

    }
    else{
      m[i]=1.0/mm;
    }
  }
  
  
  //  first step
  //  prescribed forces
  mefel_right_hand_side (i,fp);
  
  //  internal forces
  internal_forces (lcid,fi);
  
  for (i=0;i<n;i++){
    //  initial displacements
    da[i]=Lsrs->lhsi[i];
    //  initial velocities (temporarily stored in vector dn)
    dn[i]=Lsrs->tdlhsi[i];
    //  computation of vector of displacements at time t_init - dt
    dp[i]=da[i]-dt*dn[i]+dt*dt/2.0*m[i]*(fp[i]+fi[i]);
  }
  
  
  print_init(-1, "wt");
  //print_step(lcid, i, Mp->time, f);
  print_flush();
  
  
  do{
    fprintf (stdout,"\n time %e",Mp->time);
    
    //  prescribed forces
    mefel_right_hand_side (i,fp);
    
    //  internal forces
    internal_forces (lcid,fi);
    
    for (i=0;i<n;i++){
      dn[i]=2.0*da[i]-dp[i]+dt*dt*m[i]*(fi[i]+fp[i]);
      dp[i]=da[i];
      da[i]=dn[i];
    }
    
    
    //print_step(lcid, i, Mp->time, f);
    print_flush();
    
    //  time increment
    Mp->time = Mp->timecon.newtime ();
    ni++;
    
    
    
  }while(Mp->time<end_time);
  
  print_close ();

}

/**
   function solves seismic analysis using response spectrum
   
   @param lcid - load case id
   
   JK, 20.8.2005
*/
void response_spectrum_method (long lcid)
{
  long i,j,n,m;
  double coeff,per;
  double *v,*auxv;
  
  //  number of DOFs
  n=Ndofm;
  //  number of eigenvectors used in analysis
  m=Mp->eigsol.neigv;
  
  //  computation of eigenvalues and eigenmodes
  solve_eigen_dynamics (Lsrs->lhs+n,Lsrs->w);
  
  v = new double [n];
  auxv = new double [n];
  
  nullv (Lsrs->lhs,n);
  
  for (i=0;i<m;i++){
    //  period (computed from circular frequency)
    per=2.0*3.14159265358979/sqrt(Lsrs->w[i]);
    
    nullv (v,n);
    
    //  vector of directions multiplied by response spectrum values
    Mb->dlc[0].stool.assemble (v,per);
    
    //  mass matrix multiplication
    Mmat->gmxv (v,auxv);
    
    //  dot product with eigenmode
    coeff = ss (auxv,Lsrs->lhs+(i+1)*n,n)/(Lsrs->w[i]);
    
    //fprintf (Out,"\n response spectrum - circ. frequency  %le   period %le     contribution to the %ld-th mode    %le",Lsrs->w[i],per,i,coeff);

    for (j=0;j<n;j++){
      Lsrs->lhs[j]-=coeff*Lsrs->lhs[(i+1)*n+j];
    }
    
  }
  
  
  for (j=0;j<n;j++){
    Lsrs->lhs[j]=Lsrs->lhs[n+j];
  }
  
  /*
  fprintf (Out,"\n\n\n displacements");
  for (i=0;i<n;i++){
    fprintf (Out,"\n %5ld   %le",i,Lsrs->lhs[i]);
  }
  */
  
  compute_ipstresses (0);
  
  /*
  long ipp;
  fprintf (Out,"\n\n");
  for (i=0;i<Mt->ne;i++){
    ipp=Mt->elements[i].ipp[0][0];
    fprintf (Out,"\n %le %le %le",Mm->ip[ipp].stress[0],Mm->ip[ipp].stress[1],Mm->ip[ipp].stress[2]);
    ipp++;
    fprintf (Out,"     %le %le %le",Mm->ip[ipp].stress[0],Mm->ip[ipp].stress[1],Mm->ip[ipp].stress[2]);
  }
  */

  delete [] auxv;
  delete [] v;
}