mechprint.cpp

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

#include "mechprint.h"
#include "global.h"
#include "globmat.h"
#include "alias.h"
#include "gtopology.h"
#include "meshtransfer.h"
#include "loadcase.h"
#include "element.h"
#include "node.h"
#include "intpoints.h"
#include "plelemqq.h"
#include "mathem.h"
#include "genfile.h"
#include "siftop_element_types.h"
#include "probdesc.h"
#include "siftop.h"
#include "elemhead.h"
#include "intp.h"
#include "vector.h"
//#include "gmatrix.h"


#include <stdlib.h>
#include <math.h>
#include <string.h>



/**
  Function opens data files for each type of outputfile

  @param istep - step id. The parameter enables to open file with name enhanced
                 by the step id - istep >= 0 or it leaves required filename untouched for instance that
                 istep < 0. In case of stochastic calculations and in case istep >= 0 the istep precedes the 
                 stochastic step id.
  @param mode - string with control sequence for the file opening. It enables to open 
                new file (mode = "wt") or to append existing ones (mode = "at").
  @param idn1 - id of the first node for GiD mesh (default is idn1 = 1 -> nodes are numbered from 1)
  @param ide1 - id of the first element for GiD mesh (default is ide1 = 1 -> elements are numbered from 1)
  Return :
  The function does not return anything
*/
void print_init(long istep, char *mode, long idn1, long ide1)
{
  char fname[FNAMELEN+20];
  char *path;
  char *name;
  char *suffix;
  long i;
  long sl;


  if ((Outdm->outf == NULL) && (Outdm->textout==1))
  {
    Mp->filename_decomposition (Outdm->outfn,path,name,suffix);
    Mp->path = path;
    Mp->filename=name;
    Mp->suffix=suffix;
    if (St == NULL)
    {
      if (istep < 0)
        sprintf(fname, "%s%s%s", path, name, suffix);
      else
        sprintf(fname, "%s%s.%ld%s", path, name, istep+1, suffix);
    }
    else
    {
      if (istep < 0)
        sprintf(fname, "%s%s.%ld%s", path, name, Mp->ns+1, suffix);
      else
        sprintf(fname, "%s%s.%ld.%ld%s", path, name, istep+1, Mp->ns+1, suffix);
    }
    Outdm->outf = fopen(fname, mode);
    if (Outdm->outf == NULL)
    {
      fprintf(stderr, "\n\nUnable to open output text file '%s' in function print_init()\n", fname);
      fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
      abort();
    }
    Outdm->print_header(Outdm->outf);
  }

  if ((Outdm->gf != grfmt_no) && (Outdm->outgr == NULL))
  {
    Mp->filename_decomposition (Outdm->outgrfn,path,name,suffix);
    if (St == NULL)
    {
      if (istep < 0)
        sprintf(fname, "%s%s%s", path, name, suffix);
      else
        sprintf(fname, "%s%s.%ld%s", path, name, istep+1, suffix);
    }
    else
    {
      if (istep < 0)
        sprintf(fname, "%s%s.%ld%s", path, name, Mp->ns+1, suffix);
      else
        sprintf(fname, "%s%s.%ld.%ld%s", path, name, istep+1, Mp->ns+1, suffix);
    }
    if ((Outdm->gf == grfmt_gid) || (Outdm->gf == grfmt_gid_sep))
    {
      sl = strlen(fname);
      sprintf(fname+sl, ".flavia.msh");
    }
    Outdm->outgr = fopen(fname, mode);
    if (Outdm->outgr == NULL)
    {
      fprintf(stderr, "\n\nUnable to open graphics file '%s' in function print_init()\n", fname);
      fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
      abort();
    }
    switch(Outdm->gf)
    {
      case grfmt_no:
        break;   
      case grfmt_open_dx:
        break;
      case grfmt_femcad:
        export_femcad(Outdm->outgr);
        break;
      case grfmt_gid:
        export_gid_mesh(Outdm->outgr, idn1, ide1);
        fclose(Outdm->outgr);
        if (Outdm->ncut > 0)
	{
          sprintf(fname+sl, "2d.flavia.msh");
          Outdm->outgr = fopen(fname, mode);
          export_gid_2dmesh(Outdm->outgr, Outdm->ncut, idn1, ide1);
          fclose(Outdm->outgr);
	}
        sprintf(fname+sl, ".flavia.res");
        Outdm->outgr = fopen(fname, mode);
        if (Outdm->outgr == NULL)
        {
          fprintf(stderr, "\n\nUnable to open graphics file '%s' in function print_init()\n", fname);
          fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
          abort();
        }
        if (strchr(mode, 'w') != NULL)
	{
          fprintf(Outdm->outgr, "GiD Post Results File 1.0\n");
          export_gid_gauss_pt(Outdm->outgr);
	}
        break;
      case grfmt_gid_sep:
        export_gid_mesh(Outdm->outgr, idn1, ide1);
        fclose(Outdm->outgr);
        if (Outdm->ncut > 0)
	{
          sprintf(fname+sl, "2d.flavia.msh");
          Outdm->outgr = fopen(fname, mode);
          export_gid_2dmesh(Outdm->outgr, Outdm->ncut, idn1, ide1);
          fclose(Outdm->outgr);
	}
        strncpy(Outdm->outgrfngs, fname, sl);
        Outdm->outgr = NULL;
        Outdm->create_files_gidsp(mode);
        break;
      default:
        fprintf(stderr, "\n\nUnknown type of graphics format is required in function print_init\n");
        fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
    }
    delete [] path;
    delete [] name;
    delete [] suffix;
  }
  if ((Outdm->ndiag > 0) && (Outdm->outdiagf[0] == NULL))
  {
    Mp->filename_decomposition (Outdm->outdiagfn,path,name,suffix);
    for (i=0; i<Outdm->ndiag; i++)
    {
      if (St == NULL)
      {
        if (Outdm->ndiag > 1)
          sprintf(fname, "%s%s.%ld%s", path, name, i+1, suffix);
        else
          sprintf(fname, "%s%s%s", path, name, suffix);
      }
      else
      {
        if (Outdm->ndiag > 1)
          sprintf(fname, "%s%s.%ld.%ld%s", path, name, i+1, Mp->ns+1, suffix);
        else
          sprintf(fname, "%s%s.%ld%s", path, name, Mp->ns+1, suffix);
      }
      if (Outdm->outdiagf[i] == NULL)
      {
        Outdm->outdiagf[i] = fopen(fname, mode);
        if (Outdm->outdiagf[i] == NULL)
        {
          fprintf(stderr, "\n\nUnable to open diagram file '%s' in function print_init()\n", fname);
          fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
          abort();
        }
      }
    }
    delete [] path;
    delete [] name;
    delete [] suffix;
  }
  Outdm->idn1 = idn1;
  Outdm->ide1 = ide1;


}



/**
  The function performs all prints requirements for given step.

  @param lcid - load case id
  @param istep - step id
  @param lambda - load coefficient or time (it depends on problem type and solver)
  @param fi - array of additional values at nodes which should be printed.
              It depends on problem and solver type. For example for nonlinear
              statics it contains the internal forces.
  
  Return :
  The function does not return anything.
*/
void print_step(long lcid, long istep, double lambda, double *fi)
{

  //compute_req_val (lcid);
  switch(Mp->tprob)
    {
    case linear_statics:
    case var_stiff_method:
    case lin_floating_subdomain:
    case layered_linear_statics:
      Outdm->print_out(Outdm->outf, lcid, istep, lambda);
      Outdm->print_graphics(Outdm->outgr, lcid, lambda, istep, fi);
      break;
    case mat_nonlinear_statics:
    case geom_nonlinear_statics:
    case earth_pressure:
    case forced_dynamics:
    case mech_timedependent_prob:
    case nonlin_floating_subdomain:
    case growing_mech_structure:
      Outdm->print_newstep(Outdm->outf, istep, lambda);
      Outdm->print_out(Outdm->outf, lcid, istep, lambda);
      Outdm->print_diags(lcid, lambda, istep, fi);
      Outdm->print_graphics(Outdm->outgr, lcid, lambda, istep, fi);
      break;
    case eigen_dynamics:
      Outdm->print_newstep(Outdm->outf, istep, lambda);
      Outdm->print_out(Outdm->outf, lcid, istep, lambda);
      Outdm->print_graphics(Outdm->outgr, lcid, lambda, istep, fi);
      break;
    default:
      fprintf(stderr, "\n\nUnsupported problem type is required in function print_step\n");
      fprintf(stderr, " (file %s, line %d)\n", __FILE__, __LINE__);
      break;
  }
}



/**
  The function performs buffer flush of all opened files managed by the outdriver.
  It is usefull for imadiate output of required values and should be called
  at solver after the function print_step is called. 
  
  Return :
  The function does not return anything.
*/
void print_flush()
{
  long i;
  if (Outdm->outf)
    fflush(Outdm->outf);
  if (Outdm->outgr)
    fflush(Outdm->outgr);
  for (i=0; i<Outdm->ndiag; i++)
  {
    if (Outdm->outdiagf[i])
      fflush(Outdm->outdiagf[i]);
  }
}



/**
  The function performs closing of all opened files managed by the outdriver.
  
  Return :
  The function does not return anything.
*/
void print_close()
{
  long i;
  if (Outdm->outf)
    fclose(Outdm->outf);
  Outdm->outf = NULL;
  if (Outdm->outgr)
    fclose(Outdm->outgr);
  Outdm->outgr = NULL;
  for (i=0; i<Outdm->ndiag; i++)
  {
    if (Outdm->outdiagf[i])
      fclose(Outdm->outdiagf[i]);
    Outdm->outdiagf[i] = NULL;
  }
}



/**
  The function exports sets of used elements to the file given by parameter out in GiD format. 

  Parameters :
  @param out  - pointer to the opened text file where the output will be produced
  @param idn1 - id of the first node for GiD mesh (default should be idn1 = 1 -> nodes are numbered from 1)
  @param ide1 - id of the first element for GiD mesh (default should be ide1 = 1 -> elements are numbered from 1)

  Return :
  The function does not return anything
*/
void export_gid_mesh(FILE *out, long idn1, long ide1)
{
  long i, print_header, print_coord = 1;

  print_header = 1;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;
    
    switch(Mt->elements[i].te)
    {
      case bar2d :
      case bar3d :
      case beam2d :
      case beam3d :
      case subsoilbeam :
        if (print_header)
	{
          fprintf(out, "MESH beams2 dimension 3  Elemtype Linear Nnode 2\n");
          print_header = 0;
          if (print_coord)
	  {
            write_gid_nodes(out, idn1);
            print_coord = 0;
	  }
          fprintf(out, "Elements\n");
	}
        write_gid_element(out, i, idn1, ide1);
        break;    
      default:
        break;
    }
  }
  if (print_header == 0)
    fprintf(out, "end Elements\n");
  print_header = 1;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    switch(Mt->elements[i].te)
    {
      case barq2d :
        if (print_header)
	{
          fprintf(out, "MESH beams3 dimension 3  Elemtype Linear Nnode 3\n");
          print_header = 0;
          if (print_coord)
	  {
            write_gid_nodes(out, idn1);
            print_coord = 0;
	  }
          fprintf(out, "Elements\n");
	}
        write_gid_element(out, i, idn1, ide1);
        break;    
      default:
        break;
    }
  }
  if (print_header == 0)
    fprintf(out, "end Elements\n");
  print_header = 1;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    switch(Mt->elements[i].te)
    {
      case cctel :
      case dktel:
      case dstel:
      case subsoilplatetr:
      case axisymmlt :
      case planeelementlt :
      case planeelementrotlt :
        if (print_header)
	{
          fprintf(out, "MESH trias3 dimension 3  Elemtype Triangle Nnode 3\n");
          print_header = 0;
          if (print_coord)
	  {
            write_gid_nodes(out, idn1);
            print_coord = 0;
	  }
          fprintf(out, "Elements\n");
	}
        write_gid_element(out, i, idn1, ide1);
        break;    
      default:
        break;
    }
  }
  if (print_header == 0)
    fprintf(out, "end Elements\n");
  print_header = 1;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    switch(Mt->elements[i].te)
    {
      case planeelementqt :
        if (print_header)
	{
          fprintf(out, "MESH trias6 dimension 3  Elemtype Triangle Nnode 6\n");
          print_header = 0;
          if (print_coord)
	  {
            write_gid_nodes(out, idn1);
            print_coord = 0;
	  }
          fprintf(out, "Elements\n");
	}
        write_gid_element(out, i, idn1, ide1);
        break;    
      default:
        break;
    }
  }
  if (print_header == 0)
    fprintf(out, "end Elements\n");
  print_header = 1;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    switch(Mt->elements[i].te)
    {
      case axisymmlq :
      case planeelementlq :
      case planeelementrotlq :
      case q4plateel:
      case subsoilplateq: