graphmat.cpp

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

#include "graphmat.h"
#include "global.h"
#include "intpoints.h"
#include "../../GEFEL/tablefunct.h"
#include "../../GEFEL/matrix.h"
#include "../../GEFEL/PARSER/parser.h"

/**
  This constructor inializes attributes to zero or default values.
*/
graphmat::graphmat(void)
{
   gt = glinear;
   numf = 0;
   k = 0.0;
   eq = NULL;
   tab = NULL;
}

/**
  This destructor deallocates used memory.
*/
graphmat::~graphmat(void)
{
   for (long i = 0; i < numf; i++)
   {
     delete deq[i];
     delete eq[i];
   }
   delete [] deq;
   delete [] eq;
   delete tab;
}

/**
  This function reads material parameters from the opened text file given
  by the parameter in. First it reads the graph type attribute which
  denotes type of used input of the diagram. Then the appropriate
  input function is called.

  @param in - pointer to the opned text file

  @retval 0 - on success
  @retval 1 - unable to read graph type
  @retval 2 - error parsing single expression
  @retval 3 - error parsing multiple expression
  @retval 4 - unknown graph type
*/
long graphmat::read (XFILE *in)
{
  long i;
  Parser par;
  char func[256];

  xfscanf(in, "%d", (int *)&gt);
  if ((gt < glinear) || (gt > gfunc_ser))
  {
    fprintf (stderr,"\n\n Wrong type of material graph graphmat::read() (%s, line %d).\n", __FILE__, __LINE__);
    return(1);
  }

  switch (gt)
  {
    case glinear :
    {
      xfscanf (in,"%lf", &k);
      break;
    }
    case gtable:
    {
      tab = new tablefunct;
      tab->read(in);
      break;
    }
    case gfunc:
    {
      eq = new Equation* [1];
      fgets(func, 255, in->file);
      eq[0] = par.TextToTree(func);
      if (eq[0] == NULL)
      {
        fprintf(stderr, "\nError parsing expression in function graphmat::read() (%s, line %d).\n", __FILE__, __LINE__);
        return(2);
      }
      deq[0] = new Equation;
      eq[0]->Differentiate(deq[0], 1, &par);
      break;
    }
    case gfunc_ser:
    {
      xfscanf(in, "%ld", &numf);
      eq  = new Equation* [numf];
      deq = new Equation* [numf];
      limval = new double [numf];
      for (i = 0; i < numf; i++)
      {
        xfscanf(in, "%le", limval+i);
        fgets(func, 255, in->file);
        eq[i]  = par.TextToTree(func);
        if (eq[i] == NULL)
        {
          fprintf(stderr, "\nError parsing expression in function graphmat::read() (%s, line %d).\n", __FILE__, __LINE__);
          return(3);
        }
        deq[i] = new Equation;
        eq[i]->Differentiate(deq[i], 1, &par);
      }
      break;
    }
    default:
    {
      fprintf (stderr,"\n\n Unknown type of graph is required in function graphmat::read() (%s, line %d).\n", __FILE__, __LINE__);
      return (4);
    }
  }
  return(0);
}

/**
  This function computes material stiffnes matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number

*/
void graphmat::matstiff (matrix &d, long ipp)
{
  long i, id;
  double val;

  switch (gt)
  {
    case glinear :
      for (i = 0; i < d.m; i++)
        d[i][i] = k;
      break;
    case gtable :
      for (i = 0; i < d.m; i++)
      {
        if (i >= Mm->ip[ipp].ncompstr)
        {
          d[i][i] = 0.0;
          break;
        }
        tab->lininterpol(Mm->ip[ipp].strain[i], k);
        d[i][i] = k;
      }
      break;
    case gfunc :
      for (i = 0; i < d.m; i++)
      {
        if (i >= Mm->ip[ipp].ncompstr)
        {
          d[i][i] = 0.0;
          break;
        }
        if (deq[0]->Variables.at(0))
          deq[0]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
        d[i][i] = deq[0]->Evaluate();
      }
      break;
    case gfunc_ser :
    {
      for (i = 0; i < d.m; i++)
      {
        for (id = 0; id < numf; id++)
        {
          val = Mm->ip[ipp].strain[i];
          if (limval[id] >= val)
          {
            if (i >= Mm->ip[ipp].ncompstr)
            {
              d[i][i] = 0.0;
              break;
            }
            if (deq[id]->Variables.at(0))
              deq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
            d[i][i] = deq[id]->Evaluate();
            break;
          }
          if (id == numf-1)
          {
            if (deq[id]->Variables.at(0))
              deq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
            d[i][i] = deq[id]->Evaluate();
          }
        }
      }
      break;
    }
    default:
      fprintf (stderr,"\n\n Unknown graph type is required in function graphmat::matstiff() (%s, line %d).\n\n",
               __FILE__, __LINE__);
  }
  return;
}


/**
  This function computes stresses at given integration point ipp,
  depending on the reached strains.
  The stress of the given integration point is actualized.

  @param ipp - integration point number in the mechmat ip array.
*/
void graphmat::nlstresses (long ipp)
  //
{
  long i, id, nc=Mm->ip[ipp].ncompstr;

  switch (gt)
  {
    case glinear :
      for (i = 0; i < nc; i++)
        Mm->ip[ipp].stress[i] = k * Mm->ip[ipp].strain[i];
      break;
    case gtable :
      for (i = 0; i < nc; i++)
      {
        tab->lininterpol(Mm->ip[ipp].strain[i], k);
        Mm->ip[ipp].stress[i] = k * Mm->ip[ipp].strain[i];
      }
      break;
    case gfunc :
      for (i = 0; i < nc; i++)
      {
        if (eq[0]->Variables.at(0))
          eq[0]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
        Mm->ip[ipp].stress[i] = eq[0]->Evaluate();
      }
      break;
    case gfunc_ser :
    {
      for (i = 0; i < nc; i++)
      {
        for (id = 0; id < numf; id++)
        {
          if (limval[id] >= Mm->ip[ipp].strain[i])
          {
            if (eq[id]->Variables.at(0))
              eq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
            Mm->ip[ipp].stress[i] = eq[id]->Evaluate();
            break;
          }
          if (id == numf-1)
          {
            if (eq[id]->Variables.at(0))
              eq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
            Mm->ip[ipp].stress[i] = eq[id]->Evaluate();
          }
        }
      }
      break;
    }
    default:
      fprintf (stderr,"\n\n Unknown graph type is required in function graphmat::nlstresses() (%s, line %d).\n\n",
               __FILE__, __LINE__);
  }
}