output.cpp

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

 @retval 0 - if succes
 @retval 1 - if nodes haven't assigned any boundary condition
*/
long wr_bocon(FILE *out, long nprop)
{
  long i, j, k, n, nl,nbc, tcn;

  nbc = 0;
  tcn = 0;
  Nspd = 0;
  Ndpd = 0;
  Boccn = new long *[nprop];
  memset(Boccn, 0, nprop*sizeof(*Boccn));
  nl = Nlay;
  if (nl == 0)
    nl = 1;
  for (i = 0; i < Top->nodes; i++)
  {
    for (j = 0; j < nprop; j++)
    // loop over all properties from input
    {
      if (Top->N[i].property == Nnp[j])
      // test if number of property matchs number from input
      {
        if (Boc[j])
        // test if any boundary condition is prescribed
        {
          nbc++;
          if (Boccn[j] == NULL)
          {
            Boccn[j] = new long[Tndofn[i]];
            memset(Boccn[j], 0, Tndofn[i] * sizeof(long));
            for (k = 0; k < Tndofn[i]; k++)
            // loop for assigning bc code numbers
            // if direction has bc, increase cn
            {
              if ((Boc[j]->dir[k] >= 0) && (Boc[j]->con[k] != 0.0))
              {
                tcn++;
                Boccn[j][k] = tcn;
                if (Boc[j]->expr[k] == NULL)
                  Nspd++;
                else
                  Ndpd++;
              }
            }
          }
        }
        if (Nccn[j])
        // test if any common code numbers are prescribed
          nbc++;

      }
    }
  }
/*  if (nbc == 0)
  {
    fprintf(stderr, "\nError - problem hasn't assigned any boundary condition\n");
    return(1);
  }*/
  fprintf(out, "%ld\n", nbc*nl);
  for (i = 0; i < Top->nodes; i++)
  {
    for (j = 0; j < nprop; j++)
    // loop over all properties from input
    {
      if (Top->N[i].property == Nnp[j])
      // test if number of property matchs number from input
      {
        if ((Boc[j]) || (Nccn[j]))
        // test if any boundary condition is prescribed or any common code numbers are prescribed
        {
          for (n = 0; n < nl; n++)
          {
            fprintf(out, "%ld", i*nl+n+1);
            for (k = 0; k < Tndofn[i]; k++)
            // loop for writing bc
            // if direction has bc, write 0 else write 1
            {
              if (Boc[j] && Nccn[j])
              {
                if ((Boc[j]->dir[k] > 0) && (Nccn[j][k]))
                {
                  fprintf(stderr, "\n\nError - node %ld has prescribed boundary condition and common code number\n", i+1);
                  fprintf(stderr, " in the same direction\n.");
                  fprintf(stderr, " function wr_bocon (input.cpp)\n.");
                  return(2);
                }
                if ((Boc[j]->dir[k] < 0) && (Nccn[j][k] == 0))
                // no boundary condition nor common code number
                {
                  fprintf(out, " 1");
                  continue;
                }
              }
              if (Boc[j])
              {
                if (Boc[j]->dir[k] < 0)
                // no boundary condition
                {
                  fprintf(out, " 1");
                  continue;
                }
                if (Boc[j]->con[k] == 0.0)
                // rigid support
                {
                  fprintf(out, " 0");
                  continue;
                }
                if (Boc[j]->con[k] != 0.0)
                // prescribed displacement
                {
                  fprintf(out, " %ld", -Boccn[j][k]);
                  continue;
                }
              }
              if (Nccn[j])
              // common code number
              {
                if (Nccn[j][k])
                {
                  fprintf(out, " %ld", Nccn[j][k]);
                  continue;
                }
                if (Nccn[j][k] == 0)
                {
                  fprintf(out, " 1");
                  continue;
                }
              }
            }
            fprintf(out, "\n");
          }
          break;
        }
      }
    }
  }
  fprintf(out, "\n");
  return(0);
}



/**
 Function writes section with description of elements to the text file given by out.
  @param out - pointer to the opened text file, where the data will be written
  @param nprop - number of entries in the input property file

  @retval 0 : on succes
  @retval 1 : if element hasn't assigned type
  @retval 2 : if cross-section was not found in database
  @retval 3 : if element hasn't assigned material
  @retval 4 : if material was not found in database
  @retval 5 : if element hasn't assigned cross-section
  @retval 6 : if element hasn't assigned time function
*/
long wr_elements(FILE *out, long nprop, descrip *d)
{
  long i, j, k, l, tip, tfid;
  crsectype *ce;
  long *ci, *ic;
  strastrestate ssst;
  long rindex;
  elemtype et;

  mattype **me;
  long **mi, **im, *nm;

  me   = new mattype* [Top->elements];
  memset(me, 0, sizeof(*me)*Top->elements);
  mi   = new long* [Top->elements];
  memset(mi, 0, sizeof(*mi)*Top->elements);
  im   = new long* [Top->elements];
  memset(im, 0, sizeof(*im)*Top->elements);
  nm   = new long  [Top->elements];

  Tet = new elemtype[Top->elements];
  memset(Tet, 0, sizeof(*Tet)*Top->elements);
  ce   = new crsectype [Top->elements];
  memset(ce, 0, sizeof(*ce)*Top->elements);
  ci   = new long  [Top->elements];
  memset(ci, 0, sizeof(*ci)*Top->elements);
  ic   = new long  [Top->elements];

  fprintf(out, "%ld\n", Top->elements);
  for (i = 0; i < Top->elements; i++)
  {
    for (j = 0, et = elemtype(-1); j < nprop; j++)
    // loop over all properties from input
    {
      // test if number of property matchs number from input
      if ((Top->E[i].property == Nep[j]) && (Sel[j] == Top->E[i].type))
      {
        if (Tel[j] >= 0)
        // assigning element type
        {
          et = Tel[j];
          ssst = Ssse[j];
          break;
        }
      }
    }
    if (et < 0)
    {
      fprintf(stderr, "\nError - element number %ld hasn't assigned type\n", i+1);
      for (k = 0; k < Top->elements; k++)
      {
        delete [] me[k]; delete [] mi[k]; delete [] im[k];
      }
      delete [] me; delete [] mi; delete [] nm;
      delete [] ce; delete [] ci; delete [] ic;
      return(1);
    }
    for (j = 0; j < nprop; j++)
    // loop over all properties from input
    {
      if ((Top->E[i].property == Nep[j]) && (Ecs[j]))
      // assigning cross-section to the element with number i
      {
        ce[i] = Ecs[j];
        ci[i] = Ecsi[j];
      }
    }
    if (ce[i] > 0)
    {
      if ((ic[i] = Dbcrs->search_crs(ce[i], ci[i]-1)) < 0)
      {
        fprintf(stderr, "\nError - cross-section with number %d and index %ld not found\n", ce[i], ci[i]);
        delete [] ce; delete [] ci; delete [] ic;
        return(2);
      }
      else
      {
        if (Dbcrs->crs[ic[i]].propu[ci[i]-1] == 0)
          Dbcrs->crs[ic[i]].npropu++; // indicator whether given cs type is used
        if (Ncu[ic[i]] == 0)
          Ncu[ic[i]]++; // indicator whether cs type from cs database is used in the given task
        Dbcrs->crs[ic[i]].propu[ci[i]-1]++;  // indicator whether given cs index is used
      }
    }
    for (j = 0; j < Dbcrs->numt; j++)
    {
      if (Dbcrs->crs[j].npropu)
      {
        rindex = 1;
        for (k = 0; k < Dbcrs->crs[j].nprop; k++)
        {
          if (Dbcrs->crs[j].propu[k])
          {
            if (Dbcrs->crs[j].ridx[k] != 0)
              rindex++;
          }
        }
        for (k = 0; k < Dbcrs->crs[j].nprop; k++)
        {
          if (Dbcrs->crs[j].propu[k])
          {
            if (Dbcrs->crs[j].ridx[k] == 0)
	    {
              Dbcrs->crs[j].ridx[k] = rindex;
              rindex++;
	    }
          }
        }
      }
    }
    nm[i] = 0;
    for (j = 0; j < nprop; j++)
    // loop over all properties from input
    {
      // test if number of property matchs number from input
      if (Top->E[i].property == Nep[j])
      {
        if (Mel[j])
        // assigning matrial to the element with number i
        {
          nm[i] = Nmel[j];
          me[i] = new mattype[Nmel[j]];
          mi[i] = new long[Nmel[j]];
          im[i] = new long[Nmel[j]];
          for (k = 0; k < Nmel[j]; k++)
          {
            me[i][k] = Mel[j][k];
            mi[i][k] = Melid[j][k];
          }
        }
      }
    }
    if (me[i] ==  NULL) // unassigned element material
    {
      fprintf(stderr, "\nError - element %ld hasn't assigned material\n", i+1);
      for (k = 0; k < Top->elements; k++)
      {
        delete [] me[k]; delete [] mi[k]; delete [] im[k];
      }
      delete [] me; delete [] mi; delete [] nm;
      return(3);
    }
    for (j = 0; j < nm[i]; j++)
    {
      // searching element material in the material database
      if ((im[i][j] = Dbmat->search_mat(me[i][j], mi[i][j]-1)) < 0)
      {
        fprintf(stderr,
                "\nError - material with number %d and index %ld not found\n", me[i][j], mi[i][j]);
        for (k = 0; k < Top->elements; k++)
        {
          delete [] me[k]; delete [] mi[k]; delete [] im[k];
        }
        delete [] me; delete [] mi; delete [] nm;
        return(4);
      }
      else
      {
        if (Dbmat->mat[im[i][j]].propu[mi[i][j]-1] == 0)
          Dbmat->mat[im[i][j]].npropu++; // indicator whether given material type is used
        if (Nmu[im[i][j]] == 0)
          Nmu[im[i][j]]++;  // indicator whether material type from material database is used in the given task
        Dbmat->mat[im[i][j]].propu[mi[i][j]-1]++; // indicator whether given material index is used
      }
    }
    for (j = 0; j < Dbmat->numt; j++)
    {
      if (Dbmat->mat[j].npropu)
      {
        rindex = 1;
        for (k = 0; k < Dbmat->mat[j].nprop; k++)
        {
          if (Dbmat->mat[j].propu[k])
          {
            if (Dbmat->mat[j].ridx[k] != 0)
              rindex++;
          }
        }
        for (k = 0; k < Dbmat->mat[j].nprop; k++)
        {
          if (Dbmat->mat[j].propu[k])
          {
            if (Dbmat->mat[j].ridx[k] == 0)
	    {
              Dbmat->mat[j].ridx[k] = rindex;
              rindex++;
	    }
          }
        }
      }
    }
    tip = 0;
    switch (et)  // determining number of integration points on the element and
                 //  whether ip needs cross-section
    {
      case beam2d:
      {
        tip = 1;
        break;
      }
      default :
        break;
    }
    if ((tip) && (ce[i] == 0))
    {
      fprintf(stderr, "\nError - element %ld hasn't assigned cross-section\n", i+1);
      delete [] ce; delete [] ci; delete [] ic;
      return(5);
    }
    if (Mp->tprob == growing_mech_structure)
    {
      for (j = 0, tfid = -1; j < nprop; j++)
      // loop over all properties from input
      {
        // test if number of property matchs number from input
        if ((Top->E[i].property == Nep[j]) && (Tfel[j]))         
        {
          // assigning element time function
          {
            tfid = Tfel[j];
            break;
          }
        }
      }
      if (tfid < 0)
      {
        fprintf(stderr, "\nError - element number %ld hasn't assigned time function\n", i+1);
        for (k = 0; k < Top->elements; k++)
        {
          delete [] me[k]; delete [] mi[k]; delete [] im[k];
        }
        delete [] me; delete [] mi; delete [] nm;
        delete [] ce; delete [] ci; delete [] ic;
        return(6);
      }
    }
    Mt->elements[i].te = et;
//    Top->E[i].type = et;

    // writing type of element, element number
    fprintf(out, "%6ld ", i+1);
    if ((et >= planeelementlt) && (et <= planeelementsubqt))
      fprintf(out, "%3d %3d", et, ssst);
    else
      fprintf(out, "%3d", et);
    Tet[i] = et;
    // writing element nodes
    for (j = 0; j < Top->E[i].size; j++)
    {
/*      if ((et == planeelementqt) && (d->t3d) && (j == 3))
      // nodes in T3d file has in this case folowing order 1 2 3 5 6 4 in view of
      // SIFEL order 1 2 3 4 5 6
        continue;*/
      fprintf(out, " %6ld", Top->E[i].node[j]);
    }
/*    if ((et == planeelementqt) && (d->t3d))
    // nodes in T3d file has in this case folowing order 1 2 3 5 6 4 in view of
    // SIFEL order 1 2 3 4 5 6
    // so it is need to write fourth node
      fprintf(out, " %6ld", Top->E[i].node[3]);*/

    // no code numbers on element (code numbers not supported at this time)
    fprintf(out, " 0");
  
    // time function for growing structure problem
    if (Mp->tprob == growing_mech_structure)
      fprintf(out, " %ld", tfid);
      
    if (ce[i] > 0)
      fprintf(out, " %d %ld", ce[i], Dbcrs->crs[ic[i]].ridx[ci[i]-1]); // if element has cs, write cs type and index
    else
      fprintf(out, " 0"); // writing element cross-section indicator
    fprintf(out, " %ld", nm[i]);  
    for (l = 0; l < nm[i]; l++)
      fprintf(out, " %d %ld", me[i][l], Dbmat->mat[im[i][l]].ridx[mi[i][l]-1]); // writing element material type and index
    fprintf(out, "\n");
  }
  fprintf(out, "\n");
  for (k = 0; k < Top->elements; k++)
  {
    delete [] me[k];
    delete [] mi[k];
    delete [] im[k];
  }
  delete [] me; delete [] mi; delete [] nm;
  delete [] ce; delete [] ci; delete [] ic;
  return(0);
}



/**
 Function writes section with description of elements for layered problems to the text
  file given by out.
  @param out - pointer to the opened text file, where the data will be written
  @param nprop - number of entries in the input property file

  @retval 0 : on succes
  @retval 1 : if element hasn't assigned type
*/
long wr_lay_elements(FILE *out, long nprop, descrip *d)
{
  long i, j, k, l, tip, ln;
  crsectype *ce;
  long *ci, *ic;
  strastrestate ssst;
  long rindex;
  elemtype et;

  mattype **me;
  long **mi, **im, *nm;

  me   = new mattype* [Top->elements];
  memset(me, 0, sizeof(*me)*Top->elements);
  mi   = new long* [Top->elements];
  memset(mi, 0, sizeof(*mi)*Top->elements);
  im   = new long* [Top->elements];
  memset(im, 0, sizeof(*im)*Top->elements);
  nm   = new long  [Top->elements];

  Tet = new elemtype[Top->elements];
  memset(Tet, 0, sizeof(*Tet)*Top->elements);
  ce   = new crsectype [Top->elements];
  memset(ce, 0, sizeof(*ce)*Top->elements);
  ci   = new long  [Top->elements];
  memset(ci, 0, sizeof(*ci)*Top->elements);
  ic   = new long  [Top->elements];

  fprintf(out, "%ld\n", Top->elements*Nlay);
  for (i = 0; i < Top->elements; i++)
  {