probdesc.cpp

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

#include <string.h>
#include "probdesc.h"
#include "global.h"
#include "auxdatout.h"
#include "gmatrix.h"
#include "intools.h"
#include "iotools.h"


probdesc::probdesc ()
{
  niilnr=0;  errnr=0.0;

  tstorsm = (storagetype) 0;
  tstormm = (storagetype) 0;
  stmat = initial_stiff;
  tepsol = (epsolvertype) 0;

  //  initialization of damping description
  damp = nodamp;
  dmass=0.0;  dstiff=0.0;
  
  //  compute stresses during evaluation of internal forces
  strcomp = 1;

  //  for open dx
  detnodstrain = 0;
  
  stochasticcalc=0;
  eigstrains=0;
  temperature=0;
  
  adaptivity=0;  diagonalization = 0;   adaptflag=0;  adapt_accuracy=0.0;
  
  limit=0.0;  zero=1.0e-20;

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

  alphafvn=0.0;  deltafvn=0.0;

  time = 0.0;
  phase = 0;
  //  internal forces are computed at nodes
  nodeintfor = 1;
  lambda = 0.0;
  
  niep =0;  stepep = 0.0; errep = 0.0;  rlerrep = 0.0;  nselnodep = 0;  selnodep=NULL; loadcoef = NULL;

  matmodel = (materialmodel) 0;
  
  straincomp = 0;  strainaver = 0;  strainstate=0;  strainpos=0;
  stresscomp = 0;  stressaver = 0;  stressstate=0;  stresspos=0;
  othercomp = 0;   otheraver = 0;   otherstate=0;   otherpos=0;
  reactcomp = 0;

  
  memset(name, 0, sizeof(*name)*1001);

  tlm=0;
  
  homog=0;
  
  mcnne=0;  dim=0;
  
  ense=0;
  
  nivsm=0;
  path = filename = suffix = NULL;
  selnodep = NULL;
  loadcoef = NULL;

  //  type of time printing
  tpr = seconds;
}

probdesc::~probdesc ()
{
  delete [] selnodep;
  delete [] loadcoef;

  delete [] path; 
  delete [] filename; 
  delete [] suffix;
}

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

  //  problem name
  //inputln (in,name,1000);
  xfscanf (in,"%a",name);
  
  //  detail messages
  Mespr=0;
  xfscanf (in,"%k%ld","mespr",&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
  xfscanf (in,"%k%m","problemtype",&problemtype_kwdset,(int*)&tprob);
  
  // *************************
  //  computation of strains
  // *************************
  xfscanf (in,"%ld",&straincomp);
  if (straincomp!=0 && straincomp!=1)
    fprintf (stderr,"\n\n wrong definition of strain computation in function probdesc::read (file %s, line %d).\n",
	     __FILE__,__LINE__);
  if (straincomp==1){
    xfscanf (in,"%ld %ld",&strainpos,&strainaver);

    if (strainpos!=1 && strainpos!=2 && strainpos!=3){
      fprintf (stderr,"\n\n wrong definition of position where strains are computed");
      fprintf (stderr,"\n in function probdesc::read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    if (strainaver!=0 && strainaver!=1)
      fprintf (stderr,"\n\n wrong definition of strain averaging in function probdesc::read (file %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 (strainpos==1)   fprintf (stdout,"\n strains will be computed at integration points");
      if (strainpos==2)   fprintf (stdout,"\n strains will be copied to nodes from the closest integration point");
      if (strainpos==3)   fprintf (stdout,"\n strains will be computed at nodes");
      
      if (strainaver==0)  fprintf (stdout,"\n strains will not be averaged");
      if (strainaver==1)  fprintf (stdout,"\n strains will be averaged at nodes");
    }
  }
  
  
  // **************************
  //  computation of stresses
  // **************************
  xfscanf (in,"%ld",&stresscomp);
  if (stresscomp!=0 && stresscomp!=1)
    fprintf (stderr,"\n\n wrong definition of stress computation in function probdesc::read (file %s, line %d).\n",
	     __FILE__,__LINE__);
  if (stresscomp==1){
    xfscanf (in,"%ld %ld",&stresspos,&stressaver);

    if (strainpos!=1 && strainpos!=2 && strainpos!=3){
      fprintf (stderr,"\n\n wrong definition of position where stresses are computed");
      fprintf (stderr,"\n in function probdesc::read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    if (stressaver!=0 && stressaver!=1)
      fprintf (stderr,"\n\n wrong definition of stress averaging in function probdesc::read (file %s, line %d).\n",
	       __FILE__,__LINE__);
  }
  
  if (Mespr==1){
    if (stresscomp==0)  fprintf (stdout,"\n stresses will not be computed and stored");
    else{
      fprintf (stdout,"\n stresses will be computed and stored");
      if (stresspos==1)   fprintf (stdout,"\n stresses will be computed at integration points");
      if (stresspos==2)   fprintf (stdout,"\n stresses will be copied to nodes from the closest integration point");
      if (stresspos==3)   fprintf (stdout,"\n stresses will be computed at nodes");
      
      if (stressaver==0)  fprintf (stdout,"\n stresses will not be averaged");
      if (stressaver==1)  fprintf (stdout,"\n stresses will be averaged at nodes");
    }
  }
  
  
  // ***************************
  //  computation of reactions
  // ***************************
  xfscanf (in,"%ld",&reactcomp);
  if (reactcomp!=0 && reactcomp!=1)
    fprintf (stderr,"\n\n wrong definition of reaction computation in function probdesc::read (file %s, line %d).\n",
	     __FILE__,__LINE__);
  if (Mespr==1){
    if (reactcomp==0)  fprintf (stdout,"\n reactions will not be computed and stored");
    else  fprintf (stdout,"\n reactions will be computed and stored");
  }
  

  // **************************
  //   adaptivity computation
  // **************************
  xfscanf (in,"%k%ld","adaptivity",&adaptivity);
  if (adaptivity<0 || adaptivity>50)
    fprintf (stderr,"\n\n wrong adaptivity requirements in function probdesc::read (file %s, line %d).\n",__FILE__,__LINE__);
  
  if (adaptivity>0){
    c = getc (in->file);
    if (c == 'd')
      diagonalization = 1;
    else
      ungetc (c,in->file);
    
    switch (tprob){
    case linear_statics:   { xfscanf (in,"%ld %lf %lf",&adaptflag,&adapt_accuracy,&corr);                      break; }
    case mat_nonlinear_statics:{ xfscanf (in,"%ld %lf %lf %lf %lf %lf",&adaptflag,&adapt_accuracy,&corr,&limf,&clos,&plra);  break; }
    default:{
      fprintf (stderr,"\n\n unknown problem type is required in function probdesc::read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  
  // **********************************
  //  stochastic or fuzzy computation
  // **********************************
  //  stochasticcalc=0 - deterministic computation
  //  stochasticcalc=1 - stochastic/fuzzy computation, data are read all at once
  //  stochasticcalc=2 - stochastic/fuzzy computation, data are read sequentially
  //  stochasticcalc=3 - stochastic/fuzzy computation, data are generated in the code
  xfscanf (in,"%k%ld","stochasticcalc",&stochasticcalc);
  if (Mespr==1){
    if (stochasticcalc==0)
      fprintf (stdout,"\n deterministic computation will be computed");
    if (stochasticcalc>0)
      fprintf (stdout,"\n stochastic computation will be computed");
    
  }
  
  // ****************************
  //  homogenization procedures
  // ****************************
  xfscanf (in,"%k%ld","homogenization",&homog);
  if (Mespr==1){
    if (homog==0)
      fprintf (stdout,"\n homogenization will not be computed");
    if (homog!=0)//changed 1.2.2006, [e-mail]
      fprintf (stdout,"\n homogenization will be computed");
  }
  
  // ***************************
  //  renumbering of the nodes
  // ***************************
  xfscanf (in,"%k%m","noderenumber",&noderenumb_kwdset,(int*)&Gtm->nodren);
  if (Mespr==1){
    switch (Gtm->nodren){
    case no_renumbering:{
      fprintf (stdout,"\n nodes will not be renumbered");
      break;
    }
    case rev_cuthill_mckee:{
      fprintf (stdout,"\n nodes will be renumbered by reverse Cuthill-McKee algorithm");
      break;
    }
    case sloan:{
      fprintf (stdout,"\n nodes will be renumbered by Sloan algorithm");
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown type of node reordering is required in function read (file %s, line %d)\n",__FILE__,__LINE__);
    }
    }
  }

  
  switch (tprob){
    // *****************
    //  LINEAR STATICS
    // *****************
  case linear_statics:{
    if (Mespr==1)  fprintf (stdout,"\n linear static problem");
    
    //  tstorsm - type of storage of stiffness matrix
    xfscanf (in,"%k%m","stiffmatstor",&storagetype_kwdset,(int*)&tstorsm);
    
    //  data about solver of system of linear equations
    ssle.read (in,Mespr);
    
    break;
  }
    // *********************************************************************************
    //  EIGENDYNAMICS - EIGENVALUES AND EIGENVECTORS (EIGENFREQUENCIES AND EIGENMODES)
    // *********************************************************************************
  case eigen_dynamics:{
    if (Mespr==1)  fprintf (stdout,"\n eigenvalue vibration problem");
    
    //  type of storage of stiffness matrix
    //  type of storage of mass matrix
    xfscanf (in,"%k%m%k%m","stiffmatstor",&storagetype_kwdset,(int*)&tstorsm,"massmatstor",&storagetype_kwdset,(int*)&tstormm);
    
    //  data about solver of eigenvalues and eigenvectors
    eigsol.read (in);
    
    //  data about solver of system of linear equations
    ssle.read (in,Mespr);

    break;
  }
    // ******************
    //  FORCED DYNAMICS
    // ******************
  case forced_dynamics:{
    if (Mespr==1)  fprintf (stdout,"\n forced vibration problem");
    
    //  type of solver of forced vibration problem
    xfscanf (in,"%k%m","forceddyn",&forcedsolver_kwdset,(int*)&tforvib);
    
    switch (tforvib){
    case newmark:{
      xfscanf (in,"%lf %lf",&alphafvn,&deltafvn);
      break;
    }
    case findiff:{
      break;
    }
    case modal_analysis:{
      //  data about eigenmode/eigenvectors
      eigsol.read (in);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown solver of forced vibration is required");
      fprintf (stderr,"\n in function probdesc::read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
    
    //  type of damping
    xfscanf (in,"%k%m","damping",&damping_kwdset,(int*)&damp);
    switch (damp){
    case nodamp:{
      if (Mespr==1)  fprintf (stdout,"\n no damping is prescribed");
      break;
    }
    case massdamp:{
      xfscanf (in,"%lf",&dmass);
      if (Mespr==1)  fprintf (stdout,"\n damping proportional to the mass matrix is prescribed");
      break;
    }
    case stiffdamp:{
      xfscanf (in,"%lf",&dstiff);
      if (Mespr==1)  fprintf (stdout,"\n damping proportional to the stiffness matrix is prescribed");
      break;
    }
    case rayleighdamp:{
      xfscanf (in,"%lf %lf",&dstiff,&dmass);
      if (Mespr==1)  fprintf (stdout,"\n Rayleigh damping (proportional to the mass and stiffness matrices) is prescribed");
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown type of damping is required in function probdesc::read (file %s, line %d),\n",__FILE__,__LINE__);
    }
    }

    //  type of storage of stiffness matrix
    //  type of storage of mass matrix
    //  type of storage of damping matrix is defined by type of damping
    xfscanf (in,"%k%m%k%m","stiffmatstor",&storagetype_kwdset,(int*)&tstorsm,"massmatstor",&storagetype_kwdset,(int*)&tstormm);

    //  time controller
    timecon.read (in);

    //  data about solver of system of linear equations
    ssle.read (in,Mespr);

    break;
  }
    // *******************************
    //  MATERIALLY NONLINEAR STATICS
    // *******************************
  case mat_nonlinear_statics:{
    if (Mespr==1)  fprintf (stdout,"\n nonlinear static problem");
    
    //  data about nonlinear solver
    nlman.read (in,Mespr);
    
    //  type of storage of stiffness matrix
    //  type of stiffness matrix (initial, tangent)
    xfscanf (in,"%d %d",(int*)&tstorsm,(int*)&stmat);
    
    //  data about solver of system of linear equations
    ssle.read (in,Mespr);
    
    break;
  }
    // *********************************************************************
    //  TIME-DEPENDENT MECHANICAL PROBLEMS WITH NEGLIGIBLE INERTIAL FORCES
    // *********************************************************************
  case mech_timedependent_prob:{
    if (Mespr==1)  fprintf (stdout,"\n system of nonlinear equations will be solved by visco-solver");
    
    //  time controller
    timecon.read (in);

    //  type of storage of stiffness matrix
    //  type of stiffness matrix (initial, tangent)
    xfscanf (in,"%d %d",(int*)&tstorsm,(int*)&stmat);
    
    //  data about solver of system of linear equations
    ssle.read (in,Mespr);
    
    break;
  }
    // **********************************************************************************************
    //  TIME-DEPENDENT MECHANICAL PROBLEMS WITH NEGLIGIBLE INERTIAL FORCES AND GRADUAL CONSTRUCTION
    // **********************************************************************************************
  case growing_mech_structure:{
    if (Mespr==1)  fprintf (stdout,"\n system with gradual construction will be solved by visco-solver");
    
    //  time controller
    timecon.read (in);
    
    //  type of storage of stiffness matrix
    //  type of stiffness matrix (initial, tangent)
    xfscanf (in,"%d %d",(int*)&tstorsm,(int*)&stmat);
    
    //  data about solver of system of linear equations
    ssle.read (in,Mespr);
    
    xfscanf (in,"%ld %lf",&niilnr,&errnr);
    
    //  type of nonlinear solver
    xfscanf (in,"%d",(int*)&nrsolv);
    
    break;
  }
    // **************************
    //  EARTH PRESSURE PROBLEMS
    // **************************
  case earth_pressure:
  {
    if (Mespr==1)  fprintf (stdout,"\n system of nonlinear equations will be solved by earth pressure solver");
    
    xfscanf(in, "%d", (int*)&tepsol);
    switch (tepsol)
    {
      
      case gep_sol :
        xfscanf (in,"%ld %lf %lf %lf",&niep,&errep,&stepep,&rlerrep);

	//  manager of solver of nonlinear equations
	nlman.read (in,Mespr);
	
        //  type of storage of stiffness matrix
        xfscanf (in,"%d",(int*)&tstorsm);
        //  data about solver of system of linear equations
        ssle.read (in,Mespr);
        break;  
      case gepvarsup_sol :
        xfscanf (in,"%ld %lf %lf",&niep,&errep, &rlerrep);

	//  manager of solver of nonlinear equations
	nlman.read (in,Mespr);

        xfscanf (in, "%ld", &nselnodep);
        if (niep < nselnodep)
        {
          fprintf(stderr, "\n\n Number of iteration steps is less then number of selected nodes");
          fprintf(stderr, "\n in file %s, line %d", __FILE__, __LINE__);
        }
        selnodep = new long [nselnodep];
        memset(selnodep, 0, sizeof(*selnodep)*nselnodep);
        for (i = 0; i < nselnodep; i++)
        {
          xfscanf (in, "%ld", selnodep+i);
          selnodep[i]--;
        }
        //  type of storage of stiffness matrix
        xfscanf (in,"%d",(int*)&tstorsm);
        //  data about solver of system of linear equations
        ssle.read (in,Mespr);
        break;  
      case epplast_sol :