graphmattime.cpp

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

#include "graphmattime.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.
*/
graphmattime::graphmattime(void)
{
   gt = glinear;
   numf = 0;
   k = 0.0;
   eq = NULL;
   tab = NULL;
}

/**
  This destructor deallocates used memory.
*/
graphmattime::~graphmattime(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 graphmattime::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 graphmattime::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 graphmattime::read() (%s, line %d).\n", __FILE__, __LINE__);
        return(2);
      }
      deq = new Equation* [1];
      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 graphmattime::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 graphmattime::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 graphmattime::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;
        }
        d[i][i] = tab->getval2(Mm->ip[ipp].strain[i],val);
        d[i][i] = val;
      }
      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 graphmattime::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 graphmattime::nlstresses (long ipp, long im, long ido)
  //
{
  long i, id, nc=Mm->ip[ipp].ncompstr;
  double kn, ko,val;

  if (Mp->phase == 1)
  {
    switch (gt)
    {
      case glinear :
        for (i = 0; i < nc; i++)
          Mm->ip[ipp].stress[i] = 0;
        break;
      case gtable :
        for (i = 0; i < nc; i++)
        {
          kn = tab->getval2(Mm->ip[ipp].strain[i],val);
          kn = val;
          ko = tab->getval2(Mm->ip[ipp].eqother[ido+i],val);
          ko = val;
          Mm->ip[ipp].stress[i] = (kn-ko) * (Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[ido+i]);
        }
        break;
      case gfunc :
        for (i = 0; i < nc; i++)
        {
          if (deq[0]->Variables.at(0))
            deq[0]->Variables.at(0)->Value = Mm->ip[ipp].eqother[ido+i];
          ko = deq[0]->Evaluate();
          if (deq[0]->Variables.at(0))
            deq[0]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
          kn = deq[0]->Evaluate();
          Mm->ip[ipp].stress[i] = (kn-ko) * (Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[ido+i]);
        }
        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 (deq[id]->Variables.at(0))
                deq[id]->Variables.at(0)->Value = Mm->ip[ipp].eqother[ido+i];
              ko = deq[id]->Evaluate();
              if (deq[id]->Variables.at(0))
                deq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
              kn = deq[id]->Evaluate();
              Mm->ip[ipp].stress[i] = (kn-ko) * (Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[ido+i]);
              break;
            }
            if (id == numf-1)
            {
              if (deq[id]->Variables.at(0))
                deq[id]->Variables.at(0)->Value = Mm->ip[ipp].eqother[ido+i];
              ko = deq[id]->Evaluate();
              if (deq[id]->Variables.at(0))
                deq[id]->Variables.at(0)->Value = Mm->ip[ipp].strain[i];
              kn = deq[id]->Evaluate();
              Mm->ip[ipp].stress[i] = (kn-ko) * (Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[ido+i]);
            }
          }
        }
        break;
      }
      default:
        fprintf (stderr,"\n\n Unknown graph type is required in function graphmattime::nlstresses() (%s, line %d).\n\n",
                 __FILE__, __LINE__);
    }
    for(i=0; i<Mm->ip[ipp].ncompstr; i++)
      Mm->ip[ipp].eqother[ido+i] = Mm->ip[ipp].strain[i];
  }
  if (Mp->phase == 2)
  {
    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++)
        {
          k = tab->getval2(Mm->ip[ipp].strain[i],val);
          k = val;
          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 graphmattime::nlstresses() (%s, line %d).\n\n",
                 __FILE__, __LINE__);
    }
  }
}