mechprint.cpp

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

        if (print_header)
	{
          fprintf(out, "MESH Quads4 dimension 3  Elemtype Quadrilateral Nnode 4\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 planeelementqq :
      case axisymmqq :
        if (print_header)
	{
          fprintf(out, "MESH Quads8 dimension 3  Elemtype Quadrilateral Nnode 8\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 lineartet:
        if (print_header)
	{
          fprintf(out, "MESH Tetras4 dimension 3  Elemtype Tetrahedra Nnode 4\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 linearhex :
        if (print_header)
	{
          fprintf(out, "MESH Brick8 dimension 3  Elemtype Hexahedra Nnode 8\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 quadrhex:
        if (print_header)
	{
          fprintf(out, "MESH Brick20 dimension 3  Elemtype Hexahedra Nnode 20\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");
}


/**
  The function exports sets of gausspoints 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

  Return :
  The function does not return anything
*/
void export_gid_gauss_pt(FILE *out)
{
  long i, j, k, ii, jj, brk;
  vector gp1, gp2, gp3, gp, w;


  brk = 0;
  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 :
        fprintf(out, "GaussPoints \"Lin_1D\" Elemtype Linear\n");
        fprintf(out, "Number Of Gauss Points: 3\n");
        fprintf(out, "Nodes included\n");
        fprintf(out, "Natural coordinates: internal\n");
        fprintf(out, "end GaussPoints\n\n");
        brk = 1;
        break;    
      default:
        break;
    }
    if (brk)
      break;
  }

  brk = 0;
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    switch(Mt->elements[i].te)
    {
      case barq2d :
        fprintf(out, "GaussPoints \"Quad_1D\" Elemtype Linear\n");
        fprintf(out, "Number Of Gauss Points: 3\n");
        fprintf(out, "Nodes not included\n");
        fprintf(out, "Natural coordinates: internal\n");
        fprintf(out, "end GaussPoints\n\n");
        brk = 1;
        break;    
      default:
        break;
    }
    if (brk)
      break;
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == axisymmlt)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Triangle\n", axisymmlt);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Asymlt->nb;ii++)
      {
        for (jj=0;jj<Asymlt->nb;jj++)
        {
          if (Asymlt->intordsm[ii][jj]==0)  
            continue;
          allocv (Asymlt->intordsm[ii][jj],w);
          allocv (Asymlt->intordsm[ii][jj],gp1);
          allocv (Asymlt->intordsm[ii][jj],gp2);
          gauss_points_tr (gp1.a,gp2.a,w.a,Asymlt->intordsm[ii][jj]);
          for (i=0;i<Asymlt->intordsm[ii][jj];i++)
            fprintf(out, "%le %le\n", gp1[i], gp2[i]);
          destrv(gp1); destrv(gp2); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == cctel)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Triangle\n", cctel);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Cct->nb;ii++)
      {
        for (jj=0;jj<Cct->nb;jj++)
        {
          if (Cct->intordsm[ii][jj]==0)  
            continue;
          allocv (Cct->intordsm[ii][jj],w);
          allocv (Cct->intordsm[ii][jj],gp1);
          allocv (Cct->intordsm[ii][jj],gp2);
          gauss_points_tr (gp1.a,gp2.a,w.a,Cct->intordsm[ii][jj]);
          for (i=0;i<Cct->intordsm[ii][jj];i++)
            fprintf(out, "%le %le\n", gp1[i], gp2[i]);
          destrv(gp1); destrv(gp2); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == planeelementlt)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Triangle\n", planeelementlt);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: internal\n");
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == planeelementrotlt)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Triangle\n", planeelementrotlt);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Perlt->nb;ii++)
      {
        for (jj=0;jj<Perlt->nb;jj++)
        {
          if (Perlt->intordsm[ii][jj]==0)  
            continue;
          allocv (Perlt->intordsm[ii][jj],w);
          allocv (Perlt->intordsm[ii][jj],gp1);
          allocv (Perlt->intordsm[ii][jj],gp2);
          gauss_points_tr (gp1.a,gp2.a,w.a,Perlt->intordsm[ii][jj]);
          for (i=0;i<Perlt->intordsm[ii][jj];i++)
            fprintf(out, "%le %le\n", gp1[i], gp2[i]);
          destrv(gp1); destrv(gp2); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == planeelementqt)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Triangle\n", planeelementqt);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Peqt->nb;ii++)
      {
        for (jj=0;jj<Peqt->nb;jj++)
        {
          if (Peqt->intordsm[ii][jj]==0)  
            continue;
          allocv (Peqt->intordsm[ii][jj],w);
          allocv (Peqt->intordsm[ii][jj],gp1);
          allocv (Peqt->intordsm[ii][jj],gp2);
          gauss_points_tr (gp1.a,gp2.a,w.a,Peqt->intordsm[ii][jj]);
          for (i=0;i<Peqt->intordsm[ii][jj];i++)
            fprintf(out, "%le %le\n", gp1[i], gp2[i]);
          destrv(gp1); destrv(gp2); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == axisymmlq)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", axisymmlq);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Asymlq->nb;ii++)
      {
        for (jj=0;jj<Asymlq->nb;jj++)
        {
          if (Asymlq->intordsm[ii][jj]==0)  
            continue;
          allocv (Asymlq->intordsm[ii][jj],w);
          allocv (Asymlq->intordsm[ii][jj],gp);
          gauss_points (gp.a,w.a,Asymlq->intordsm[ii][jj]);
          for (i=0;i<Asymlq->intordsm[ii][jj];i++)
            for (j=0; j<Asymlq->intordsm[ii][jj]; j++)
              fprintf(out, "%le %le\n", gp[i], gp[j]);
          destrv(gp); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }

  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    /*********************************************/
    /* old=original */
    /* if (Mt->elements[i].te == planeelementlq)
       {
       fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", planeelementlq);
       fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
       fprintf(out, "Nodes not included\n");
       fprintf(out, "Natural coordinates: given\n");
       for (ii=0;ii<Pelq->nb;ii++)
       {
       for (jj=0;jj<Pelq->nb;jj++)
       {
       if (Pelq->intordsm[ii][jj]==0)  
       continue;
       allocv (Pelq->intordsm[ii][jj],w);
       allocv (Pelq->intordsm[ii][jj],gp);
       gauss_points (gp.a,w.a,Pelq->intordsm[ii][jj]);
       for (i=0;i<Pelq->intordsm[ii][jj];i++)
       for (j=0; j<Pelq->intordsm[ii][jj]; j++)
       fprintf(out, "%le %le\n", -gp[i], -gp[j]);
       destrv(gp); destrv(w);
       }
       }
       fprintf(out, "end GaussPoints\n\n");
       break;
       }
    */
    /*********************************************/
    /* small correction for GiD printing */
    if (Mt->elements[i].te == planeelementlq)
      {
	fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", planeelementlq);
	fprintf(out, "Number Of Gauss Points: %ld\n", 4);
	fprintf(out, "Nodes not included\n");
	fprintf(out, "Natural coordinates: given\n");
	for (ii=0;ii<1;ii++)
	  {
	    for (jj=0;jj<1;jj++)
	      {
		if (Pelq->intordsm[ii][jj]==0)  
		  continue;
		allocv (Pelq->intordsm[ii][jj],w);
		allocv (Pelq->intordsm[ii][jj],gp);
		gauss_points (gp.a,w.a,Pelq->intordsm[ii][jj]);
		for (i=0;i<Pelq->intordsm[ii][jj];i++)
		  for (j=0; j<Pelq->intordsm[ii][jj]; j++)
		    fprintf(out, "%le %le\n", -gp[i], -gp[j]);
		destrv(gp); destrv(w);
	      }
	  }
	fprintf(out, "end GaussPoints\n\n");
	break;
      }
  }
  
  for (i=0; i < Mt->ne; i++)
  {
    if (Gtm->leso[i]==0)
      continue;

    if (Mt->elements[i].te == planeelementrotlq)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", planeelementrotlq);
      fprintf(out, "Number Of Gauss Points: %ld\n", Mt->give_totnip(i));
      fprintf(out, "Nodes not included\n");
      fprintf(out, "Natural coordinates: given\n");
      for (ii=0;ii<Perlq->nb; ii++)
      {
        for (jj=0;jj<Perlq->nb;jj++)
        {
          if (Perlq->intordsm[ii][jj]==0)  
            continue;
          allocv (Perlq->intordsm[ii][jj],w);
          allocv (Perlq->intordsm[ii][jj],gp);
          gauss_points (gp.a,w.a,Perlq->intordsm[ii][jj]);
          for (i=0;i<Perlq->intordsm[ii][jj];i++)
            for (j=0; j<Perlq->intordsm[ii][jj]; j++)
              fprintf(out, "%le %le\n", -gp[i], -gp[j]);
          destrv(gp); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");