probdesc.jk

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

#include "probdesc.h"
#include "global.h"
#include "auxdatout.h"
#include <string.h>

probdesc::probdesc (void)
{
  layers=0;

  tstorsm = (storagetype) 0;
  tstormm = (storagetype) 0;

  tlinsol = (linsolvertype) 0;
  tprec = (precondtype) 0;
  tnlinsol = (nonlinsolvertype) 0;
  displnorm = (displacementnorm) 0;
  tsra = (stressretalgtype) 0;
  tepsol = (epsolvertype) 0;

  straincomp = 0;  stresscomp = 0;  reactcomp = 0;
  strainaver = 0;  stressaver = 0;

  adaptivity=0;  diagonalization=0;  adapt_accuracy=0.0;

  nicg=0;  anicg=0;  errcg=0.0;  aerrcg=0.0;
  ssoromega=0.0;  indegamma=0.0;
  limit=0.0;  zero=1.0e-20;

  npun = 0;  requn = NULL;  stmat = 0;

  neigv=0;  nev=0;  nies=0;  anies=0;  erres=0.0;  aerres=0.0;
  nijmr=0;  njacthr=0;  jacthr=NULL;

  timeincr=0.0;  nifv=0;  alphafvn=0.0;  deltafvn=0.0;  betafvh=0.0;  gammafvh=0.0;

  nial = 0;  niilal = 0;  erral = 0.0;  dlal = 0.0;  dlmaxal = 0.0;  psial = 0.0;
  nsnal=0;  selnodal=NULL;  hdbackupal=0;

  ninr = 0;  niilnr = 0;  errnr = 0.0;  incrnr = 0.0;  minincrnr = 0.0;
  nivs = 0;  incrvs = 0.0;  time = 0.0;  phase = 0;

  nicp = 0;  errcp = 0.0;
  nicppvp = 0;  errcppvp = 0.0;  ddtcppvp = 0.0;

  niep =0; niiep = 0; errep = 0.0;
  
}

probdesc::~probdesc (void)
{

}

/**
   function reads basic data about solved problem
   
   @param in - input stream
   
   20.7.2001
*/
void probdesc::read (FILE *in)
{
  long i;

  //  problem name
  inputln (in,name,1000);
  
  //  output file
  inputln (in,outfile,1000);
  
  //  detail messages
  Mespr=0;
  fscanf (in,"%ld",&Mespr);
  if (Mespr==1)  fprintf (stdout,"\n detail information will be printed\n");
  else fprintf (stdout,"\n only important messages will be printed\n");
  
  //  problem type
  fscanf (in,"%d",(int*)&tprob);
  
  //  computation of strains, stresses and reactions
  fscanf (in,"%ld",&straincomp);
  if (straincomp!=0 && straincomp!=1)
    fprintf (stderr,"\n\n wrong identification of strain computation in function probdesc::read (%s, line %d).\n",__FILE__,__LINE__);
  if (straincomp==1)  fscanf (in,"%ld",&strainaver);
  if (strainaver!=0 && strainaver!=1)
    fprintf (stderr,"\n\n wrong identification of strain averaging in function probdesc::read (%s, line %d).\n",__FILE__,__LINE__);

  fscanf (in,"%ld",&stresscomp);
  if (stresscomp!=0 && stresscomp!=1)
    fprintf (stderr,"\n\n wrong identification of stress computation in function probdesc::read (%s, line %d).\n",__FILE__,__LINE__);
  if (stresscomp==1)  fscanf (in,"%ld",&stressaver);
  if (stressaver!=0 && stressaver!=1)
    fprintf (stderr,"\n\n wrong identification of stress averaging in function probdesc::read (%s, line %d).\n",__FILE__,__LINE__);

  fscanf (in,"%ld",&reactcomp);
  if (reactcomp!=0 && reactcomp!=1)
    fprintf (stderr,"\n\n wrong identification of reaction computation in function probdesc::read (%s, line %d).\n",__FILE__,__LINE__);

  if (Mespr==1){
    if (straincomp==0)  fprintf (stdout,"\n strains will not be computed and stored");
    else{
      fprintf (stdout,"\n strains will be computed and stored");
      if (strainaver==0)  fprintf (stdout,"\n strains will not be averaged");
      if (strainaver==1)  fprintf (stdout,"\n strains will be averaged at nodes");
    }

    if (stresscomp==0)  fprintf (stdout,"\n stresses will not be computed and stored");
    else{
      fprintf (stdout,"\n stresses will be computed and stored");
      if (stressaver==0)  fprintf (stdout,"\n stresses will not be averaged");
      if (stressaver==1)  fprintf (stdout,"\n stresses will be averaged at nodes");
    }
    
    if (reactcomp==0)  fprintf (stdout,"\n reactions will not be computed and stored");
    else  fprintf (stdout,"\n reactions will be computed and stored");
  }
  
  //  presence of adaptivity
  fscanf (in,"%ld",&adaptivity);
  if (adaptivity<0 || adaptivity>50)
    fprintf (stderr,"\n\n wrong adaptivity requirements in function probdesc::read (probdesc.cpp).\n");
  
  if (adaptivity>0){
    fscanf (in,"%ld %lf",&diagonalization,&adapt_accuracy);
  }



  switch (tprob){
  case linear_statics:{
    if (Mespr==1)  fprintf (stdout,"\n linear static problem");

    //  tstorsm - type of storage of stiffness matrix
    //  tlinsol - type of solver of linear equation system
    fscanf (in,"%d %d",(int*)&tstorsm,(int*)&tlinsol);

    switch (tlinsol){
    case gauss_elim:{
      if (Mespr==1)  fprintf (stdout,"\n system of linear equations will be solved by gaussian elimination method");
      break;
    }
    case ldl:{
      if (Mespr==1)  fprintf (stdout,"\n system of linear equations will be solved by LDL method");
      break;
    }
    case lu:{
      if (Mespr==1)  fprintf (stdout,"\n system of linear equations will be solved by LU method");
      break;
    }
    case cg:{
      if (Mespr==1)  fprintf (stdout,"\n system of linear equations will be solved by CG method");
      fscanf (in,"%ld %lf",&nicg,&errcg);
      break;
    }
    case bicg:{
      if (Mespr==1)  fprintf (stdout,"\n system of linear equationes will be solved by bi-conjugate gradient method");
      fscanf (in,"%ld %lf",&nicg,&errcg);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown type of solver of linear equations is required");
      fprintf (stderr,"\n in function probdesc::read (probdesc.cpp).\n");
    }
    }
    
    if (tlinsol==cg){
      fscanf (in,"%d",(int*)&tprec);
      switch (tprec){
      case noprecond:{
        if (Mespr==1)  fprintf (stdout,"\n no preconditioner will be used");
        break;
      }
      case diagprec:{
        if (Mespr==1)  fprintf (stdout,"\n diagonal (Jacobi) preconditioner will be used");
        break;
      }
      case ssorprec:{
        if (Mespr==1)  fprintf (stdout,"\n SSOR preconditioner will be used");
        fscanf (in,"%lf",&ssoromega);
        break;
      }
      case incomdec:{
        if (Mespr==1)  fprintf (stdout,"\n preconditioner based on incomplete decomposition of matrix will be used");
        fscanf (in,"%lf",&indegamma);
        break;
      }
      default:{
        fprintf (stderr,"\n\n unknown preconditioner type is required in function probdesc::read (probdesc.cpp).\n");
      }
      }
    }
    
    break;
  }
  case eigen_dynamics:{
    if (Mespr==1)  fprintf (stdout,"\n eigenvalue vibration problem");
    
    //  type of eigenvector solver
    //  type of storage of stiffness matrix
    //  type of storage of mass matrix
    fscanf (in,"%d %d %d %d",(int*)&teigsol,(int*)&tlinsol,(int*)&tstorsm,(int*)&tstormm);
    
    switch (teigsol){
    case inv_iteration:{
      if (Mespr==1)  fprintf (stdout,"\n eigenvalue problem will be solved by inverse iteration method");
      
      //  number of eigenvectors
      nev=1;  neigv=1;
      
      fscanf (in,"%ld %lf",&nies,&erres);
      break;
    }
    case subspace_it_jacobi:{
      if (Mespr==1){
        fprintf (stdout,"\n eigenvalue problem will be solved by subspace");
        fprintf (stdout,"\n iteration method with Jacobi method of rotations");
      }
      
      fscanf (in,"%ld %ld %ld %lf",&neigv,&nev,&nies,&erres);
      if (neigv>nev){
        fprintf (stderr,"\n\n wrong number of vectors used in subspace iteration method");
        fprintf (stderr,"\n nev must be greater or equal to neigv. \n");
      }
      fscanf (in,"%ld %ld",&nijmr,&njacthr);
      jacthr = new double [njacthr];
      for (i=0;i<njacthr;i++){
        fscanf (in,"%lf",jacthr+i);
      }
      
      break;
    }
    case subspace_it_gsortho:{
      if (Mespr==1){
        fprintf (stdout,"\n eigenvalue problem will be solved by subspace");
        fprintf (stdout,"\n iteration method with Gram-Schmidt orthonormalization technique");
      }
      
      fscanf (in,"%ld %ld %ld %le",&neigv,&nev,&nies,&erres);
      if (neigv>nev){
        fprintf (stdout,"\n\n wrong number of vectors used in subspace iteration method");
        fprintf (stdout,"\n nev must be greater or equal to neigv. \n");
      }
      
      break;
    }
    default:{
      fprintf (stderr,"\n\n wrong type of eigenvalue problem solver is required");
      fprintf (stderr,"\n in function probdesc::read (probdesc.cpp).\n");
    }
    }

    break;
  }
  case forced_dynamics:{
    if (Mespr==1)  fprintf (stdout,"\n forced vibration problem");
    
    //  auxiliary output file
    getnextln(in);
    inputln (in,auxfile,1000);
    
    fscanf (in,"%d %d %d %d",(int*)&tforvib,(int*)&tlinsol,(int*)&tstorsm,(int*)&tstormm);
    switch (tforvib){
    case newmark:{
      fscanf (in,"%ld %lf %lf %lf",&nifv,&timeincr,&alphafvn,&deltafvn);
      break;
    }
    case hughes:{
      fscanf (in,"%ld %lf %lf %lf",&nifv,&timeincr,&betafvh,&gammafvh);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown solver of forced vibration is required");
      fprintf (stderr,"\n in function probdesc::read (probdesc.cpp).\n");
    }
    }
    

    //  number of printed unknowns in non-linear solvers - especially for graphics
    fscanf (in,"%ld",&npun);
    requn = new long [2*npun];
    memset (requn,0,2*npun*sizeof(long));
    //  printed unknowns in non-linear solvers
    read_prunk (in,requn,npun);
    
    //  auxiliary data output
    fscanf (in,"%ld",&adp);
    if (adp==1){
      ado = new auxdatout;
      ado->read (in);
    }

    break;
  }
  case nonlinear_statics:{
    if (Mespr==1)  fprintf (stdout,"\n nonlinear static problem");

    //  auxiliary output file
    getnextln(in);
    inputln (in,auxfile,1000);

    //  type of non-linear solver, type of stress retrun algorithm
    fscanf (in,"%d %d",(int*)&tnlinsol,(int*)&tsra);
    fscanf (in,"%d %d",(int*)&tlinsol,(int*)&tstorsm);
    
    switch (tnlinsol){
    case arcl:{
      if (Mespr==1)  fprintf (stdout,"\n system of nonlinear equations will be solved by arc-length method");
      fscanf (in,"%ld %ld %lf %lf %lf %lf %lf",&nial,&niilal,&erral,&dlal,&dlminal,&dlmaxal,&psial);
      
      fscanf (in,"%ld %d",&hdbackupal,(int*)&displnorm);
      
      switch (displnorm){
      case alldofs:{  break; }
      case seldofs:{
        fscanf (in,"%ld",&nsdofal);
        selnodal = new long [nsdofal];
        seldofal = new long [nsdofal];
        read_seldof (in,selnodal,seldofal,nsdofal);
        break;
      }
      case selnodes:{
        fscanf (in,"%ld",&nsnal);
        selnodal = new long [nsnal];
        read_selnod (in,selnodal,nsnal);
        break;
      }
      case nodesdistincr:{
        nsnal=2;
        fscanf (in,"%ld",&probdimal);
        selnodal = new long [nsnal];
        read_selnod (in,selnodal,nsnal);
        break;
      }
      default:{
        fprintf (stderr,"\n\n unknown norm measurement of displacement increment (%s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
      break;
    }
    case newton:{
      if (Mespr==1)  fprintf (stdout,"\n system of nonlinear equations will be solved by Newton-Raphson method");
      fscanf (in,"%ld %ld %lf %lf %lf",&ninr,&niilnr,&errnr,&incrnr,&minincrnr);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown solver of nonlinear system of equations");
      fprintf (stderr,"\n is required in function probdesc::read (probdesc.cpp).\n");
    }
    }
    
    switch (tsra){
    case 0:{
      fprintf (stdout,"\n\n no stress return algorithm is required.\n");
      break;
    }
    case cp:{
      fscanf (in,"%ld %lf",&nicp,&errcp);
      break;
    }
    case cppvp:{
      fscanf (in,"%ld %lf %lf",&nicppvp,&errcppvp,&ddtcppvp);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown stress return algorithm is required in function");
      fprintf (stderr,"\n probdesc::read (probdesc.cpp).\n");
    }
    }
    
    //  number of printed unknowns in non-linear solvers - especially for graphics
    fscanf (in,"%ld",&npun);
    requn = new long [2*npun];
    memset (requn,0,2*npun*sizeof(long));
    //  printed unknowns in non-linear solvers
    read_prunk (in,requn,npun);