mechprint.cpp

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

      break;
    }
  }

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

    if (Mt->elements[i].te == planeelementqq)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", planeelementqq);
      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<Peqq->nb;ii++)
      {
        for (jj=0;jj<Peqq->nb;jj++)
        {
          if (Peqq->intordsm[ii][jj]==0)  
            continue;
          allocv (Peqq->intordsm[ii][jj],w);
          allocv (Peqq->intordsm[ii][jj],gp);
          gauss_points (gp.a,w.a,Peqq->intordsm[ii][jj]);
          for (i=0;i<Peqq->intordsm[ii][jj];i++)
            for (j=0; j<Peqq->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 == axisymmqq)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Quadrilateral\n", axisymmqq);
      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<Asymqq->nb;ii++)
      {
        for (jj=0;jj<Asymqq->nb;jj++)
        {
          if (Asymqq->intordsm[ii][jj]==0)  
            continue;
          allocv (Asymqq->intordsm[ii][jj],w);
          allocv (Asymqq->intordsm[ii][jj],gp);
          gauss_points (gp.a,w.a,Asymqq->intordsm[ii][jj]);
          for (i=0;i<Asymqq->intordsm[ii][jj];i++)
            for (j=0; j<Asymqq->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 == lineartet)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype TetraHedra\n", lineartet);
      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<Ltet->nb;ii++)
      {
        for (jj=0;jj<Ltet->nb;jj++)
        {
          if (Ltet->intordsm[ii][jj]==0)  
            continue;
          allocv (Ltet->intordsm[ii][jj],w);
          allocv (Ltet->intordsm[ii][jj],gp1);
          allocv (Ltet->intordsm[ii][jj],gp2);
          allocv (Ltet->intordsm[ii][jj],gp3);
          gauss_points_tet (gp1.a,gp2.a,gp3.a,w.a,Ltet->intordsm[ii][jj]);
          for (i=0;i<Ltet->intordsm[ii][jj];i++)
            fprintf(out, "%le %le %le\n", gp1[i], gp2[i], gp3[i]);
          destrv(gp1); destrv(gp2); destrv(gp3); 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 == linearhex)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Hexahedra\n", linearhex);
      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<Lhex->nb;ii++)
      {
        for (jj=0;jj<Lhex->nb;jj++)
        {
          if (Lhex->intordsm[ii][jj]==0)  
            continue;
          allocv (Lhex->intordsm[ii][jj],w);
          allocv (Lhex->intordsm[ii][jj],gp);
          gauss_points(gp.a,w.a,Lhex->intordsm[ii][jj]);
          for (i=0;i<Lhex->intordsm[ii][jj];i++)
          {
            for (j=0;j<Lhex->intordsm[ii][jj];j++)
            {
              for (k=0;k<Lhex->intordsm[ii][jj];k++)
                fprintf(out, "%le %le %le\n", -gp[i], -gp[j], -gp[k]);
	    }
	  }
          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 == quadrhex)
    {
      fprintf(out, "GaussPoints \"%d\" Elemtype Hexahedra\n", quadrhex);
      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<Qhex->nb;ii++)
      {
        for (jj=0;jj<Qhex->nb;jj++)
        {
          if (Qhex->intordsm[ii][jj]==0)  
            continue;
          allocv (Qhex->intordsm[ii][jj],w);
          allocv (Qhex->intordsm[ii][jj],gp);
          gauss_points(gp.a,w.a,Qhex->intordsm[ii][jj]);
          for (i=0;i<Qhex->intordsm[ii][jj];i++)
          {
            for (j=0;j<Qhex->intordsm[ii][jj];j++)
            {
              for (k=0;k<Qhex->intordsm[ii][jj];k++)
                fprintf(out, "%le %le %le\n", -gp[i], -gp[j], -gp[k]);
	    }
	  }
          destrv(gp); destrv(w);
	}
      }
      fprintf(out, "end GaussPoints\n\n");
      break;
    }
  }
}

/**
  The function exports from 3D brick element mesh the 2D plane element cuts 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 icut - index of the cut - will be used in the element property
  @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_2dmesh(FILE *out, long icut, long idn1, long ide1)
{
  long i, print_header, print_coord = 1;
  long range = Mt->ne/icut;

  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 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_2delement(out, i, 0, 4, i/range, 0, idn1, ide1);
        break;    
      default:
        break;
    }
  }
  for (i=Mt->ne-range; i < Mt->ne; i++)
  {
    write_gid_2delement(out, i, 4, 8, icut, range, idn1, ide1);
  }
  if (print_header == 0)
    fprintf(out, "end Elements\n");
}

/**
  The function exports nodes 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)

  Return :
  The function does not return anything
*/
void write_gid_nodes(FILE *out, long idn1)
{
  long i;
  fprintf(out, "Coordinates\n");
  for (i=0; i<Mt->nn; i++)
    fprintf(out, "%ld %e %e %e\n", i+idn1, Gtm->gnodes[i].x, Gtm->gnodes[i].y, Gtm->gnodes[i].z);
  fprintf(out, "end Coordinates\n");
}


/**
  The function exports element 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 i    - element id
  @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 write_gid_element(FILE *out, long i, long idn1, long ide1)
{
  long j;
  fprintf(out, "%ld ", i+ide1);
  for (j=0; j<Gtm->gelements[i].give_nne(); j++)
    fprintf(out, "%ld ", Gtm->gelements[i].nodes[j]+idn1);
  fprintf(out, "%ld\n", Mt->elements[i].idm[0]+1);
/*  for (j=0; j < Mt->elements[i].nm; j++)
    fprintf(out, "%ld%ld", Mt->elements[i].tm[j],Mt->elements[i].idm[j]+1);
  fprintf(out, "\n");*/
}

/**
  The function exports 3D brick element as 2D plane element to the file given by parameter out in GiD format. 
  Purpose of the function is to perform cuts of the domain which consists of the regular mesh
  of brick elements.

  Parameters :
  @param out  - pointer to the opened text file where the output will be produced
  @param i    - brick element id
  @param id1  - brick element index of the first node for plane element
  @param nne  - brick element index of the last node for plane element
  @param icut - index of the cut - will be used in the element property
  @param di   - the start index of the plane elements
  @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 write_gid_2delement(FILE *out, long i, long id1, long nne, long icut, long di, long idn1, long ide1)
{
  long j;
  fprintf(out, "%ld ", i+di+ide1);
  for (j=id1; j<nne; j++)
    fprintf(out, "%ld ", Gtm->gelements[i].nodes[j]+idn1);
  fprintf(out, "%ld%ld\n", icut+1, Mt->elements[i].idm[0]+1);
/*  for (j=0; j < Mt->elements[i].nm; j++)
    fprintf(out, "%ld%ld", Mt->elements[i].tm[j],Mt->elements[i].idm[j]+1);
  fprintf(out, "\n");*/
}

/**
  The function writes displacement vector for all nodes to the file given by parameter out in GiD format. 
  The results are printed for all nodes with no dependency on the outdriver selection.

  Parameters :
  @param out - pointer to the opened text file where the output will be produced
  @param lcid - load case id
  @param desclcid - string with description of loadcase

  Return :
  The function does not return anything
*/
void write_gid_displ(FILE *out, long lcid, char *desclcid)
{
  long i, j;
  vector l,g;
  matrix tm(3,3);

  fprintf(out, "\nResult \"Displacements\" \"%ld\" %s Vector OnNodes\n", lcid, desclcid);
  fprintf(out, "ComponentNames \"r_1\" \"r_2\" \"r_3\"\nValues\n");
  for (i=0; i < Mt->nn; i++)
  {
    allocv (Gtm->gnodes[i].ndofn,l);
    noddispl (lcid, l.a, i);

    if (Mt->nodes[i].transf>0){

      //fprintf (Outm,"\n uzel %ld   %ld",i,Mt->nodes[i].transf);

      allocv (Gtm->gnodes[i].ndofn,g);
      
      tm[0][0]=Mt->nodes[i].e1[0];
      tm[1][0]=Mt->nodes[i].e1[1];
      tm[2][0]=Mt->nodes[i].e1[2];
      
      tm[0][1]=Mt->nodes[i].e2[0];
      tm[1][1]=Mt->nodes[i].e2[1];
      tm[2][1]=Mt->nodes[i].e2[2];
      
      tm[0][2]=Mt->nodes[i].e3[0];
      tm[1][2]=Mt->nodes[i].e3[1];
      tm[2][2]=Mt->nodes[i].e3[2];
      
      mxv(tm,l,g);
      copyv (g,l);
      destrv (g);
    }
    
    if (Mp->tprob == growing_mech_structure){
      fprintf(out, "%ld", i+Outdm->idn1);
      for (j = 0; j < Gtm->gnodes[i].ndofn; j++)
	//for (j = 0; j < Gtm->gnodes[i].ndofn-1; j++)
	fprintf(out, " % e", l[j]+Mt->nodes[i].nvgraph[j]);
      fprintf(out, "\n");
      destrv (l);
    }
    else{
      fprintf(out, "%ld", i+Outdm->idn1);
      for (j = 0; j < Gtm->gnodes[i].ndofn; j++)
      //for (j = 0; j < 3; j++)
	fprintf(out, " % e", l[j]);
      fprintf(out, "\n");
      destrv (l);
    }
  }
  fprintf(out, "End Values\n");
}

/**
  The function writes a scalar value given by parameter scal on all nodes to the file given by parameter out in GiD format. 
  The results are printed at all nodes with no dependency of the nodal selection specified in the
  outdriver. 

  Parameters :
  @param out - pointer to the opened text file where the output will be produced
  @param scal - specifies type of required scalar quantity (strain/stres/other)
  @param lcid - load case id
  @param dir - specifies which component of the quantity array will be printed
  @param desclcid - string with description of loadcase

  Return :
  The function does not return anything
*/
void write_gid_nodscalar(FILE *out, strastre scal, long lcid, long dir, char *desclcid)
{
  char *sig = "";
  long i, ncompstr, ncompother;

  if (scal == stress)
  {
    if (Mt->nodes[0].ncompstr < 6)
    {
      switch (dir)
      {
        case 0:
          sig = "sig_n_1";
          break;
        case 1:
          sig = "sig_n_2";
          break;
        case 2:
          sig = "sig_n_3";
          break;
        case 3:
          sig = "sig_n_4";
          break;
        case 4:
          sig = "sig_n_5";
          break;
        default :
          fprintf(stderr, "\n\n Error - unknown direction of stress");
          fprintf(stderr, "\n in function wr_nodscalar, file %s, line %d\n", __FILE__, __LINE__);
      }
    }
    if (Mt->nodes[0].ncompstr == 6)
    {
      switch (dir)
      {
        case 0:
          sig = "sig_n_x";
          break;
        case 1:
          sig = "sig_n_y";
          break;
        case 2:
          sig = "sig_n_z";
          break;
        case 3:
          sig = "tau_n_yz";
          break;
        case 4:
          sig = "tau_n_xz";
          break;
        case 5:
          sig = "tau_n_xy";
          break;
        default :
          fprintf(stderr, "\n\n Error - unknown direction of stress");
          fprintf(stderr, "\n in function write_gid_nodscalar, file %s, line %d\n", __FILE__, __LINE__);
      }
    }
  }
  if (scal == strain)
  {
    if (Mt->nodes[0].ncompstr < 6)
    {
      switch (dir)
      {
        case 0:
          sig = "eps_n_1";
          break;
        case 1:
          sig = "eps_n_2";
          break;
        case 2:
          sig = "eps_n_3";
          break;
        case 3:
          sig = "eps_n_4";
          break;
        case 4:
          sig = "eps_n_5";