stochdriver.cpp

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

#include "stochdriver.h"
#include "matrix.h"
#include "vector.h"
#include "global.h"
#include "elastisomat.h"
#include "splas1d.h"
#include "j2flow.h"
#include "microM4.h"
#include "chen.h"
#include "crsec2dbar.h"
#include "crsec2dbeam.h"
#include "crsec3dbeam.h"
#include "crsecplstr.h"
#include "crsec3d.h"
#include "loadcase.h"
#include "loadn.h"
#include "intpoints.h"
#include "element.h"
#include "probdesc.h"

stochdriver::stochdriver ()
{
  //  number of stochastic materials
  nsmt=0;
  //  number of stochastic cross-sections
  nscs=0;
  //  number of stochastic nodal loads
  nsnl=0;
  
  //  number of samples
  nsampl=0;
  //  number of stochastic variables
  nstochvar=0;
  //  number of printed output variables
  nprunknowns=0;
  
  //  materials
  mt=NULL;  idm=NULL;  atm=NULL;
  //  cross-sections
  cst=NULL;  idcs=NULL;  atcs=NULL;
  //  nodal loads
  idln=NULL;  atln=NULL;
  
  //  nodes with output
  nna=NULL;  pnd=NULL;
  //  elements with output
  ena=NULL;  ev=NULL;
  
  neigv=0;

  avi = NULL;  avo = NULL;

  datin = NULL;
  datout = NULL;
  
  //  fuzzy numbers
  fn = NULL;
}

stochdriver::~stochdriver ()
{
  delete [] mt;  delete [] idm;  delete [] atm;
  delete [] cst;  delete [] idcs;  delete [] atcs;
  delete [] idln;  delete [] atln;
  delete [] nna;  delete [] pnd;
  delete [] ena;  delete [] ev;

  delete [] avi;  delete [] avo;
  
  delete [] fn;
  
  if (datin!=NULL)
    fclose (datin);
  if (datout!=NULL)
    fclose (datout);
}

/**
   function reads data from the input file
   
   @param in - input stream
   
   JK
*/
void stochdriver::read (XFILE *in)
{
  long i;
  
  if (Mp->stochasticcalc==1 || Mp->stochasticcalc==2){
    //  samples are read from files
    xfscanf (in,"%a",auxfilein);
    datin = fopen (auxfilein,"r");
  }
  if (Mp->stochasticcalc==3){
    //  samples are generated by generator
    fg.read (in);
  }
  
  //  output file for fuzzy results
  xfscanf (in,"%a",auxfileout);
  datout = fopen (auxfileout,"w");
  
  
  // ************************
  // **  input parameters  **
  // ************************

  //  number of stochastic materials
  xfscanf (in,"%ld",&nsmt);
  if (Mespr==1)  fprintf (stdout,"\n number of stochastic materials  %ld",nsmt);
  if (nsmt>0){
    mt = new mattype [nsmt];
    idm = new long [nsmt];
    atm = new atsel[nsmt];
    for (i=0;i<nsmt;i++){
      xfscanf (in,"%d %ld",(int*)mt+i,idm+i);
      idm[i]--;
      atm[i].read(in);
    }
  }
  
  //  number of stochastic cross-sections
  xfscanf (in,"%ld",&nscs);
  if (Mespr==1)  fprintf (stdout,"\n number of stochastic cross sections  %ld",nscs);
  if (nscs>0){
    cst = new crsectype [nscs];
    idcs = new long [nscs];
    atcs = new atsel [nscs];
    for (i=0;i<nscs;i++){
      xfscanf (in,"%d %ld",(int*)cst+i,idcs+i);
      idcs[i]--;
      atcs[i].read(in);
    }
  }
  
  //  number of stochastic nodal loads
  xfscanf (in,"%ld",&nsnl);
  if (Mespr==1)  fprintf (stdout,"\n number of stochastic loaded nodes  %ld",nsnl);
  if (nsnl>0){
    idln = new long [nsnl];
    atln = new atsel [nsnl];
    for (i=0;i<nsnl;i++){
      xfscanf (in,"%ld",idln+i);
      idln[i]--;
      atln[i].read(in);
    }
  }
  
  // *************************
  // **  output parameters  **
  // *************************

  //  number of printed nodal variables (displacements)
  xfscanf (in,"%ld",&npnd);
  if (Mespr==1)  fprintf (stdout,"\n number of printed nodal displacements  %ld",npnd);
  if (npnd>0){
    nna = new long [npnd];
    pnd = new atsel [npnd];
    for (i=0;i<npnd;i++){
      xfscanf (in,"%ld",nna+i);
      nna[i]--;
      pnd[i].read (in);
    }
  }
  
  
  //  number of elements with required output
  xfscanf (in,"%ld",&npev);
  if (Mespr==1)  fprintf (stdout,"\n number of elements with required output  %ld",npev);
  if (npev>0){
    ena = new long [npev];
    for (i=0;i<npev;i++){
      xfscanf (in,"%ld",ena+i);
      ena[i]--;
    }
  }
  
  
  
  // ***********************************************
  // **  initialization of the class stochdriver  **
  // ***********************************************
  
  //  number of samples and number of stochastic variables
  if (Mp->stochasticcalc==1 || Mp->stochasticcalc==2){
    fscanf (datin,"%ld %ld",&nsampl,&nstochvar);
  }
  if (Mp->stochasticcalc==3){
    nsampl=fg.tncomb;
    nstochvar=fg.nfv;
    
    fprintf (stdout,"\n number of samples  %ld",nsampl);
  }
  
  //  array of actual input variables
  avi = new double [nstochvar];
  
  //  data are read all at once
  if (Mp->stochasticcalc==1)
    readtable (datin);
  
  if (Mp->tprob==eigen_dynamics){
    neigv=Mp->eigsol.neigv;
  }
  
  
  //  computes number of output unknowns
  compute_nprunknowns ();
  fprintf (datout,"%ld\n",nprunknowns);
  
  
  //  array of actual output variables
  avo = new double [nprunknowns];
  
  //  data are stored in RAM and then printed out all at once
  if (Mp->stochasticcalc==1)
    allocm (nsampl,nprunknowns,stochtabout);

  //  output data are saved as fuzzy numbers
  if (Mp->stochasticcalc==3){
    fn = new fuzzynum [nprunknowns];
    for (i=0;i<nprunknowns;i++){
      fn[i].initiate (fg.nalph,fg.alpha);
    }
  }
  
}

/**
   function computes number of output unknowns
   
   12.2.2005, JK
*/
void stochdriver::compute_nprunknowns ()
{
  long i;
  
  nprunknowns=0;
  
  //  contribution from nodal displacements
  for (i=0;i<npnd;i++){
    nprunknowns+=pnd[i].num;
  }
  ndispl=nprunknowns;
  
  //  nastaveno jen na 2D pruty a fuzzy ram
  //  je treba to rozmyslet a dodelat
  for (i=0;i<npev;i++){
    nprunknowns+=6;
    abort ();
  }
  
  nelem=nprunknowns-ndispl;

  if (Mp->tprob==eigen_dynamics){
    nprunknowns*=neigv;
  }
  
}

/**
   function reads table of stochastic values
   number of rows = number of samples
   number of columns = number of stochastic variables in the problem
   
   13.3.2003
*/
void stochdriver::readtable (FILE *in)
{
  allocm (nsampl,nstochvar,stochtabin);
  readm(in, stochtabin);
}

/**
   function reads table of stochastic values
   number of rows = number of samples
   number of columns = number of stochastic variables in the problem
   
   13.3.2003
*/
void stochdriver::writetable ()
{
  printm(stochtabout,datout,10,20);
}



/**
   function changes stochastic values
   
   @param sampleid - sample id
   
   13.3.2003
*/
void stochdriver::changevalues (long sampleid)
{
  assemble_new_values (sampleid);
  replace_values ();
}




/**
   function assembles new input data for the next sample
   
   JK
*/
void stochdriver::assemble_new_values (long sampleid)
{
  long i;
  
  //  data are obtained from the list in RAM
  if (Mp->stochasticcalc==1){
    for (i=0;i<nstochvar;i++){
      avi[i]=stochtabin[sampleid][i];
    }
  }
  //  data are obtained from the file
  if (Mp->stochasticcalc==2){
    for (i=0;i<nstochvar;i++){
      fscanf (datin,"%le",&avi[i]);
    }
  }
  //  data are generated by fuzzy generator
  if (Mp->stochasticcalc==3){
    fg.give_new_values (sampleid,avi,Out);
  }
  
}

/**
   function replaces nondeterministic variables by new values
   
   JK
*/
void stochdriver::replace_values ()
{
  long i,j,k,l,m;
  vector val;
  
  //  stochastic materials
  j=0;
  for (i=0;i<nsmt;i++){
    k=atm[i].num;  m=0;
    allocv (k,val);
    for (l=j;l<j+k;l++){
      val[m]=avi[l];
      m++;
    }
    j+=k;
    changematerials (i,val);
    destrv (val);
  }
  
  //  stochastic cross sections
  for (i=0;i<nscs;i++){
    k=atcs[i].num;  m=0;
    allocv (k,val);
    for (l=j;l<j+k;l++){
      val[m]=avi[l];
      m++;
    }
    j+=k;
    changecrsections (i,val);
    destrv (val);
  }
  
  //  stochastic nodal loads
  for (i=0;i<nsnl;i++){
    k=atln[i].num;  m=0;
    allocv (k,val);
    for (l=j;l<j+k;l++){
      val[m]=avi[l];
      m++;
    }
    j+=k;
    changenodloads (i,val);
    destrv (val);
  }

}

/**
   function changes nondeterministic material parameters
   
   @param id - material id
   @param val - %vector of sample values
   
   JK
*/
void stochdriver::changematerials (long id,vector &val)
{
  switch (mt[id]){
  case elisomat:{
    Mm->eliso[idm[id]].changeparam (atm[id],val);
    break;
  }
  case simplas1d:{
    Mm->spl1d[idm[id]].changeparam (atm[id],val);
    break;
  }
  case jflow:{
    Mm->j2f[idm[id]].changeparam (atm[id],val);
    break;
  }
    
    
  case microplaneM4:{
    Mm->mpM4[idm[id]].changeparam (atm[id],val);
    break;
  }
    
  case chenplast:{
    Mm->chplast[idm[id]].changeparam (atm[id],val);
    break;
  }
    
/*
  case microsimp:{
    break;
  }
  case microfibro:{
    break;
  }
  case mohrcoul2d:{
    mc2d[i].matstiff (d,ipp);
    break;
  }
  case mohrcoul3d:{
    break;
  }
  case mohrcoulparab:{
    break;
  }
  case boermaterial:{
    boerm[i].matstiff (d,ipp);
    break;
  }
  case druckerprager:{
    drprm[i].matstiff (d,ipp);
    break;
  }
  case simvisplas:{
    svipl[i].matstiff (d,ipp);
    break;
  }
  case lemaitr:{
    lmtr[i].matstiff (d,ipp);
    break;
  }
  case scaldamage:{
    scdam[i].matstiff (d,ipp);
    break;
  }
  case nonlocalmod:{
    break;
  }
  case graphm:{
    break;
    }
    case creepbaz:{
    break;
    }
    case consolidation:{
    break;
    }
    */
  default:{
    fprintf (stderr,"\n\n unknown material type is required in function changematerials (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

void stochdriver::changecrsections (long id,vector &val)
{
  switch (cst[id]){
  case nocrosssection:{
    break;
  }
  case csbar2d:{
    Mc->cs2dbar[idcs[id]].changeparam (atcs[id],val);
    break;
  }
  case csbeam2d:{
    Mc->cs2dbeam[idcs[id]].changeparam (atcs[id],val);
    break;
  }
  case csbeam3d:{
    Mc->cs3dbeam[idcs[id]].changeparam (atcs[id],val);
    break;
  }
  case csplanestr:{
    Mc->csplstr[idcs[id]].changeparam (atcs[id],val);
    break;
  }
  case cs3dprob:{
    Mc->cs3d[idcs[id]].changeparam (atcs[id],val);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown material type is required in function changecrsections (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
}

void stochdriver::changenodloads (long id,vector &val)
{
  Mb->lc[0].lon[idln[id]].changeparam (atln[id],val);
}

/**
   function extracts required output values
   
   JK
*/
void stochdriver::extractor ()
{
  long i,j,k,ii,ipp,ci,nid,nloops;
  
  if (Mp->tprob==eigen_dynamics){
    nloops=neigv;
  }
  else{
    nloops=1;
  }
  
  ci=0;
  for (ii=0;ii<nloops;ii++){
    
    for (i=0;i<npnd;i++){
      nid=nna[i];
      for (j=0;j<pnd[i].num;j++){
	k=Mt->give_dof (nid,pnd[i].atrib[j])-1;
	if (k<0){
	  fprintf (stderr,"\n\n wrong number of DOF in function extractor (file %s, line %d).\n",__FILE__,__LINE__);
	}
	avo[ci]=Lsrs->lhs[ii*Ndofm+k];
	ci++;
      }
    }
    
    
    //  prvky jsou hotove jen pro pruty do fuzzy ramu
    //  jinak je to nepouzitelne a musi se to rozmyslet a dodelat
    /*
    for (i=0;i<npev;i++){
      j=ena[i];
      ipp=Mt->elements[j].ipp[0][0];

      avo[ci]=Mm->ip[ipp].stress[0];
      ci++;
      avo[ci]=Mm->ip[ipp].stress[1];
      ci++;
      avo[ci]=Mm->ip[ipp].stress[2];
      ci++;
      
      ipp++;
      avo[ci]=Mm->ip[ipp].stress[0];
      ci++;
      avo[ci]=Mm->ip[ipp].stress[1];
      ci++;
      avo[ci]=Mm->ip[ipp].stress[2];
      ci++;
      
    }
    */
  }
  
}

/**
   function stores extracted values
   required values are extracted to the array avo

   @param sampleid - sample id
   
   JK
*/
void stochdriver::save_results (long sampleid)
{  
  long i;
  
  if (Mp->stochasticcalc==1){
    for (i=0;i<nprunknowns;i++){
      stochtabout[sampleid][i]=avo[i];
    }
  }
  if (Mp->stochasticcalc==2){
    for (i=0;i<nprunknowns;i++){
      fprintf (datout,"%le  ",avo[i]);
    }
    fprintf (datout,"\n");
  }
  if (Mp->stochasticcalc==3){
    fg.save_values (fn,nprunknowns,avo);
  }
  
}

/**
   function makes postprocessing of diagrams obtained during stochastic
   or fuzzy computation of nonlinear problems
   
   JK, 12.6.2005
*/
void stochdriver::diagpostproc ()
{
  long i,j,nc,nr,minnr;
  double x,y,minx,dx;
  char *oname,fname[1000];
  char *path,*name,*suffix;
  FILE *in,*out;
  gfunct *gf;
  
  //  number of columns in diagrams
  nc = Outdm->odiag[0].npun;
  
  if (nc!=2){
    fprintf (stderr,"\n\n unsupported number of columns in function diagpostproc (file %s, line %d)\n",__FILE__,__LINE__);
    abort ();
  }

  //  name of diagram files
  oname = new char [1000];
  for (i=0;i<1000;i++){
    oname[i] = Outdm->outdiagfn[i];
  }
  
  gf = new gfunct [nsampl];

  //  loop over number of samples
  //  (over number of simulations)
  minx=1.0e8;
  minnr=10000000;
  for (i=0;i<nsampl;i++){
    
    Mp->filename_decomposition (oname,path,name,suffix);

    sprintf(fname, "%s.%ld.dat",name,i+1);
    in=fopen (fname,"r");
    
    delete [] path;
    delete [] name;
    delete [] suffix;

    //  computation of number of rows in particular diagrams
    nr=0;
    do{
      fscanf (in,"%le %le",&x,&x);
      nr++;
      
      //  indicator of end of file
      //  function feof returns 1 at the end of file
      j=feof(in);
    }while (j==0);
    
    nr--;
    if (minnr>nr)  minnr=nr;
    
    rewind (in);
    gf[i].init_tab (nr);

    for (j=0;j<nr;j++){
      fscanf (in,"%le %le",&x,&y);
      gf[i].tabf->x[j]=fabs(x);
      gf[i].tabf->y[j]=fabs(y);
    }
    
    if (fabs(x)<minx)  minx=fabs(x);
    
    fclose (in);
  }
  
  printf ("\n minnr %ld     minx %le\n",minnr,minx);
  



  out = fopen ("fuzzy.dat","w");
  dx=minx/minnr;
  x=1.0e-8;
  for (i=0;i<minnr;i++){
    fprintf (out,"%le ",x);
    for (j=0;j<nsampl;j++){
      fprintf (out,"%le ",gf[j].getval(x));
    }
    fprintf (out,"\n");
    x+=dx;
  }
  fclose (out);
  
  /*
  out = fopen ("fuzzy.dat","w");
  for (i=0;i<nsampl;i++){
    gf[i].print (out);
  }
  fclose (out);
  */
  
}

void stochdriver::update_auxparam ()
{
  if (Mp->stochasticcalc==3){

    //  update of auxiliary parameters
    fg.actcomb++;
    if (fg.actcomb==fg.ncomb){
      fg.actalph++;
      fg.actcomb=0;
    }
    
  }
  
}

/**
   function copies new input data to the buffer
   
   @param buff - buffer
   
   4.2.2006, JK
*/
/*
void stochdriver::give_new_invalues (double *buff)
{
  long i;
  for (i=0;i<nstochvar;i++){
    buff[i]=avi[i];
  }
}
*/

/**
   function copies new output data from the buffer
   
   @param buff - buffer
   
   4.2.2006, JK
*/
 /*
void stochdriver::save_new_invalues (double *buff)
{
  long i;
  for (i=0;i<nstochvar;i++){
    avi[i]=buff[i];
  }
}
 */

/**
   function copies new input data to the buffer
   
   @param buff - buffer
   
   4.2.2006, JK
*/
  /*
void stochdriver::give_new_outvalues (double *buff)
{
  long i;
  for (i=0;i<nprunknowns;i++){
    buff[i]=avo[i];
  }
}
  */
/**
   function copies new output data from the buffer
   
   @param buff - buffer
   
   4.2.2006, JK
*/
   /*
void stochdriver::save_new_outvalues (double *buff)
{
  long i;
  for (i=0;i<nprunknowns;i++){
    avo[i]=buff[i];
  }
}
   */