Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_10_5
  • master default protected
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
40 results

t15.geo

Blame
  • Print_Mesh.cpp 52.07 KiB
    // $Id: Print_Mesh.cpp,v 1.75 2006-05-14 00:48:20 geuzaine Exp $
    //
    // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
    //
    // This program is free software; you can redistribute it and/or modify
    // it under the terms of the GNU General Public License as published by
    // the Free Software Foundation; either version 2 of the License, or
    // (at your option) any later version.
    //
    // This program is distributed in the hope that it will be useful,
    // but WITHOUT ANY WARRANTY; without even the implied warranty of
    // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    // GNU General Public License for more details.
    //
    // You should have received a copy of the GNU General Public License
    // along with this program; if not, write to the Free Software
    // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
    // USA.
    // 
    // Please report all bugs and problems to <gmsh@geuz.org>.
    
    #include "Gmsh.h"
    #include "Numeric.h"
    #include "Geo.h"
    #include "CAD.h"
    #include "Mesh.h"
    #include "Create.h"
    #include "Context.h"
    #include <list>
    #include <map>
    
    extern Context_T CTX;
    extern Mesh *THEM;
    
    // Write mesh in native MSH format
    
    #define LINE            1
    #define TRIANGLE        2
    #define QUADRANGLE      3
    #define TETRAHEDRON     4
    #define HEXAHEDRON      5
    #define PRISM           6
    #define PYRAMID         7
    #define LINE_2          8
    #define TRIANGLE_2      9
    #define QUADRANGLE_2   10
    #define TETRAHEDRON_2  11
    #define HEXAHEDRON_2   12
    #define PRISM_2        13
    #define PYRAMID_2      14
    #define POINT          15
    
    static int MSH_VOL_NUM, MSH_SUR_NUM, MSH_LIN_NUM;
    static int MSH_NODE_NUM, MSH_ELEMENT_NUM, MSH_ADD;
    static int MSH_PHYSICAL_NUM, MSH_PHYSICAL_ORI;
    static FILE *MSHFILE;
    
    static void _msh_print_node(void *a, void *b)
    {
      Vertex *V = *(Vertex **) a;
    
      MSH_NODE_NUM++;
      if(CTX.mesh.renumber_nodes_continuous)
        V->Num = MSH_NODE_NUM;
    
      fprintf(MSHFILE, "%d %.16g %.16g %.16g\n",
              V->Num,
              V->Pos.X * CTX.mesh.scaling_factor,
              V->Pos.Y * CTX.mesh.scaling_factor,
              V->Pos.Z * CTX.mesh.scaling_factor);
    }
    
    static void _msh_process_nodes(Mesh *M)
    {
      int i, j, Num;
      PhysicalGroup *p;
      Vertex *pv, **ppv, v;
    
      for(i = 0; i < List_Nbr(M->PhysicalGroups); i++) {
        List_Read(M->PhysicalGroups, i, &p);
        if(p->Typ == MSH_PHYSICAL_POINT) {
          for(j = 0; j < List_Nbr(p->Entities); j++) {
            List_Read(p->Entities, j, &Num);
            pv = &v;
            pv->Num = abs(Num);
            if(!Tree_Search(M->Vertices, &pv)) {
              if((ppv = (Vertex **) Tree_PQuery(M->Points, &pv)))
                Tree_Add(M->Vertices, ppv);
            }
          }
        }
      }
    
      if(CTX.mesh.msh_file_version == 2.0)
        fprintf(MSHFILE, "$Nodes\n");
      else
        fprintf(MSHFILE, "$NOD\n");
      fprintf(MSHFILE, "%d\n", Tree_Nbr(M->Vertices));
    
      MSH_NODE_NUM = 0;
      Tree_Action(M->Vertices, _msh_print_node);
      if(CTX.mesh.msh_file_version == 2.0)
        fprintf(MSHFILE, "$EndNodes\n");
      else
        fprintf(MSHFILE, "$ENDNOD\n");
    }
    
    int _get_partition(int index)
    {
      MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, index);
      if(!part)
        return -1; // OK, no partitions available
      else
        return (*part)->Num; // partition number
    }
    
    static void _msh_print_simplex(void *a, void *b)
    {
      int i, type, nbn, nbs = 0;
    
      SimplexBase *s = *(SimplexBase **) a;
    
      if(MSH_VOL_NUM && (MSH_VOL_NUM != s->iEnt))
        return;
    
      if(MSH_SUR_NUM && (MSH_SUR_NUM != s->iEnt))
        return;
    
      if(MSH_LIN_NUM && (MSH_LIN_NUM != s->iEnt))
        return;
    
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      if(!s->V[2]) {
        nbn = 2;
        if(s->VSUP) {
          type = LINE_2;
          nbs = 1;
        }
        else
          type = LINE;
      }
      else if(!s->V[3]) {
        nbn = 3;
        if(s->VSUP) {
          type = TRIANGLE_2;
          nbs = 3;
        }
        else
          type = TRIANGLE;
      }
      else {
        nbn = 4;
        if(s->VSUP) {
          type = TETRAHEDRON_2;
          nbs = 6;
          if(s->Volume_Simplexe() < 0) {
    	Vertex *temp;
    	temp = s->V[0];	s->V[0] = s->V[1]; s->V[1] = temp;
    	temp = s->VSUP[1]; s->VSUP[1] = s->VSUP[2]; s->VSUP[2] = temp;
    	temp = s->VSUP[5]; s->VSUP[5] = s->VSUP[3]; s->VSUP[3] = temp;
          }
        }
        else{
          type = TETRAHEDRON;
          if(s->Volume_Simplexe() < 0) {
    	Vertex *temp;
    	temp = s->V[0]; s->V[0] = s->V[1]; s->V[1] = temp;
          }
        }
      }
    
      if(CTX.mesh.msh_file_version == 2.0){
        int part = _get_partition(s->iPart);
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : s->iEnt, s->iEnt);
        if(part >= 0) fprintf(MSHFILE, " %d", part);
      }
      else
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : s->iEnt, s->iEnt,
    	    nbn + nbs);
    
      if(MSH_PHYSICAL_ORI > 0) {
        for(i = 0; i < nbn; i++)
          fprintf(MSHFILE, " %d", s->V[i]->Num);
        for(i = 0; i < nbs; i++)
          fprintf(MSHFILE, " %d", s->VSUP[i]->Num);
      }
      else {
        for(i = 0; i < nbn; i++)
          fprintf(MSHFILE, " %d", s->V[nbn - i - 1]->Num);
        for(i = 0; i < nbs; i++)
          fprintf(MSHFILE, " %d", s->VSUP[nbs - i - 1]->Num);
      }
    
      fprintf(MSHFILE, "\n");
    }
    
    static void _msh_print_quadrangle(void *a, void *b)
    {
      int i, type, nbn, nbs = 0;
    
      Quadrangle *q = *(Quadrangle **) a;
    
      if(MSH_SUR_NUM && (MSH_SUR_NUM != q->iEnt))
        return;
    
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      nbn = 4;
      if(q->VSUP) {
        type = QUADRANGLE_2;
        nbs = 4 + 1;
      }
      else
        type = QUADRANGLE;
    
      if(CTX.mesh.msh_file_version == 2.0){
        int part = _get_partition(q->iPart);
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : q->iEnt, q->iEnt);
        if(part >= 0) fprintf(MSHFILE, " %d", part);
      }
      else
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : q->iEnt, q->iEnt,
    	    nbn + nbs);
    
      if(MSH_PHYSICAL_ORI > 0) {
        for(i = 0; i < nbn; i++)
          fprintf(MSHFILE, " %d", q->V[i]->Num);
        for(i = 0; i < nbs; i++)
          fprintf(MSHFILE, " %d", q->VSUP[i]->Num);
      }
      else {
        for(i = 0; i < nbn; i++)
          fprintf(MSHFILE, " %d", q->V[nbn - i - 1]->Num);
        if(nbs){
          for(i = 0; i < nbs - 1; i++)
    	fprintf(MSHFILE, " %d", q->VSUP[nbs - i - 2]->Num);
          fprintf(MSHFILE, " %d", q->VSUP[nbs - 1]->Num);
        }
      }
    
      fprintf(MSHFILE, "\n");
    }
    
    static void _msh_print_hexahedron(void *a, void *b)
    {
      int i, type, nbn, nbs = 0;
    
      Hexahedron *h = *(Hexahedron **) a;
    
      if(MSH_VOL_NUM && (MSH_VOL_NUM != h->iEnt))
        return;
    
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      nbn = 8;
      if(h->VSUP) {
        type = HEXAHEDRON_2;
        nbs = 12 + 6 + 1;
        if(h->Orientation() < 0) {
          Vertex *temp;
          temp = h->V[0]; h->V[0] = h->V[2]; h->V[2] = temp;
          temp = h->V[4]; h->V[4] = h->V[6]; h->V[6] = temp;
          Vertex *old[12];
          for(int i = 0; i < 12; i++) old[i] = h->VSUP[i];
          h->VSUP[0] = old[3]; h->VSUP[1] = old[5]; h->VSUP[2] = old[6];
          h->VSUP[3] = old[0]; h->VSUP[4] = old[4]; h->VSUP[5] = old[1];
          h->VSUP[6] = old[2]; h->VSUP[7] = old[7]; h->VSUP[8] = old[10];
          h->VSUP[9] = old[11]; h->VSUP[10] = old[8]; h->VSUP[11] = old[9];
        }
      }
      else {
        type = HEXAHEDRON;
        if(h->Orientation() < 0) {
          Vertex *temp;
          temp = h->V[0]; h->V[0] = h->V[2]; h->V[2] = temp;
          temp = h->V[4]; h->V[4] = h->V[6]; h->V[6] = temp;
        }
      }
    
      if(CTX.mesh.msh_file_version == 2.0){
        int part = _get_partition(h->iPart);
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : h->iEnt, h->iEnt);
        if(part >= 0) fprintf(MSHFILE, " %d", part);
      }
      else
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : h->iEnt, h->iEnt,
    	    nbn + nbs);
    
      for(i = 0; i < nbn; i++)
        fprintf(MSHFILE, " %d", h->V[i]->Num);
      for(i = 0; i < nbs; i++)
        fprintf(MSHFILE, " %d", h->VSUP[i]->Num);
    
      fprintf(MSHFILE, "\n");
    }
    
    static void _msh_print_prism(void *a, void *b)
    {
      int i, type, nbn, nbs = 0;
    
      Prism *p = *(Prism **) a;
    
      if(MSH_VOL_NUM && (MSH_VOL_NUM != p->iEnt))
        return;
    
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      nbn = 6;
      if(p->VSUP) {
        type = PRISM_2;
        nbs = 9 + 3;
        if(p->Orientation() < 0) {
          Vertex *temp;
          temp = p->V[0]; p->V[0] = p->V[1]; p->V[1] = temp;
          temp = p->V[3]; p->V[3] = p->V[4]; p->V[4] = temp;
          temp = p->VSUP[1]; p->VSUP[1] = p->VSUP[3]; p->VSUP[3] = temp;
          temp = p->VSUP[2]; p->VSUP[2] = p->VSUP[4]; p->VSUP[4] = temp;
          temp = p->VSUP[7]; p->VSUP[7] = p->VSUP[8]; p->VSUP[8] = temp;
        }
      }
      else {
        type = PRISM;
        if(p->Orientation() < 0) {
          Vertex *temp;
          temp = p->V[0]; p->V[0] = p->V[1]; p->V[1] = temp;
          temp = p->V[3]; p->V[3] = p->V[4]; p->V[4] = temp;
        }
      }
    
      if(CTX.mesh.msh_file_version == 2.0){
        int part = _get_partition(p->iPart);
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt);
        if(part >= 0) fprintf(MSHFILE, " %d", part);
      }
      else
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt,
    	    nbn + nbs);
    
      for(i = 0; i < nbn; i++)
        fprintf(MSHFILE, " %d", p->V[i]->Num);
      for(i = 0; i < nbs; i++)
        fprintf(MSHFILE, " %d", p->VSUP[i]->Num);
    
      fprintf(MSHFILE, "\n");
    }
    
    static void _msh_print_pyramid(void *a, void *b)
    {
      int i, type, nbn, nbs = 0;
    
      Pyramid *p = *(Pyramid **) a;
    
      if(MSH_VOL_NUM && (MSH_VOL_NUM != p->iEnt))
        return;
    
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      nbn = 5;
      if(p->VSUP) {
        type = PYRAMID_2;
        nbs = 8 + 1;
        if(p->Orientation() < 0) {
          Vertex *temp;
          temp = p->V[0]; p->V[0] = p->V[2]; p->V[2] = temp;
          temp = p->VSUP[0]; p->VSUP[0] = p->VSUP[3]; p->VSUP[3] = temp;
          temp = p->VSUP[1]; p->VSUP[1] = p->VSUP[5]; p->VSUP[5] = temp;
          temp = p->VSUP[2]; p->VSUP[2] = p->VSUP[6]; p->VSUP[6] = temp;
        }
      }
      else {
        type = PYRAMID;
        if(p->Orientation() < 0) {
          Vertex *temp;
          temp = p->V[0]; p->V[0] = p->V[2]; p->V[2] = temp;
        }
      }
    
      if(CTX.mesh.msh_file_version == 2.0){
        int part = _get_partition(p->iPart);
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt);
        if(part >= 0) fprintf(MSHFILE, " %d", part);
      }
      else
        fprintf(MSHFILE, "%d %d %d %d %d",
    	    MSH_ELEMENT_NUM++, type,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt,
    	    nbn + nbs);
    
      for(i = 0; i < nbn; i++)
        fprintf(MSHFILE, " %d", p->V[i]->Num);
      for(i = 0; i < nbs; i++)
        fprintf(MSHFILE, " %d", p->VSUP[i]->Num);
    
      fprintf(MSHFILE, "\n");
    }
    
    static void _msh_print_point(Vertex *V)
    {
      if(!MSH_ADD) {
        MSH_ELEMENT_NUM++;
        return;
      }
    
      if(CTX.mesh.msh_file_version == 2.0)
        fprintf(MSHFILE, "%d %d 2 %d %d %d\n",
    	    MSH_ELEMENT_NUM++, POINT, MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : V->Num, V->Num,
    	    V->Num);
      else
        fprintf(MSHFILE, "%d %d %d %d 1 %d\n",
    	    MSH_ELEMENT_NUM++, POINT,
    	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : V->Num, V->Num, V->Num);
    }
    
    static void _msh_print_elements(Mesh *M)
    {
      int i, j, k, Num;
    
      PhysicalGroup *p;
      Volume *pV;
      Surface *ps;
      Curve *pc;
      Vertex *pv, v;
    
      List_T *ListCurves = Tree2List(M->Curves);
      List_T *ListSurfaces = Tree2List(M->Surfaces);
      List_T *ListVolumes = Tree2List(M->Volumes);
    
      for(i = 0; i < List_Nbr(M->PhysicalGroups); i++) {
        List_Read(M->PhysicalGroups, i, &p);
        MSH_PHYSICAL_NUM = p->Num;
        MSH_VOL_NUM = MSH_SUR_NUM = MSH_LIN_NUM = 0;
    
        switch (p->Typ) {
    
        case MSH_PHYSICAL_POINT:
          for(j = 0; j < List_Nbr(p->Entities); j++) {
            pv = &v;
            List_Read(p->Entities, j, &Num);
            pv->Num = abs(Num);
            MSH_PHYSICAL_ORI = sign(Num);
            if(Tree_Query(M->Vertices, &pv))
              _msh_print_point(pv);
          }
          break;
    
        case MSH_PHYSICAL_LINE:
          if(CTX.mesh.oldxtrude) {  //for old extrusion mesh generator
            for(k = 0; k < List_Nbr(ListVolumes); k++) {
              List_Read(ListVolumes, k, &pV);
              for(j = 0; j < List_Nbr(p->Entities); j++) {
                List_Read(p->Entities, j, &Num);
                MSH_LIN_NUM = abs(Num);
                MSH_PHYSICAL_ORI = sign(Num);
                Tree_Action(pV->Lin_Surf, _msh_print_simplex);
              }
            }
          }
          else{
    	for(k = 0; k < List_Nbr(ListCurves); k++) {
    	  List_Read(ListCurves, k, &pc);
    	  for(j = 0; j < List_Nbr(p->Entities); j++) {
    	    List_Read(p->Entities, j, &Num);
    	    MSH_LIN_NUM = abs(Num);
    	    MSH_PHYSICAL_ORI = sign(Num);
    	    Tree_Action(pc->Simplexes, _msh_print_simplex);
    	    Tree_Action(pc->SimplexesBase, _msh_print_simplex);
    	  }
    	}
          }
          break;
    
        case MSH_PHYSICAL_SURFACE:
          if(CTX.mesh.oldxtrude) {  //for old extrusion mesh generator
            for(k = 0; k < List_Nbr(ListVolumes); k++) {
              List_Read(ListVolumes, k, &pV);
              for(j = 0; j < List_Nbr(p->Entities); j++) {
                List_Read(p->Entities, j, &Num);
                MSH_SUR_NUM = abs(Num);
                MSH_PHYSICAL_ORI = sign(Num);
                Tree_Action(pV->Simp_Surf, _msh_print_simplex);
                Tree_Action(pV->Quad_Surf, _msh_print_quadrangle);
              }
            }
          }
          else{
    	for(k = 0; k < List_Nbr(ListSurfaces); k++) {
    	  List_Read(ListSurfaces, k, &ps);
    	  for(j = 0; j < List_Nbr(p->Entities); j++) {
    	    List_Read(p->Entities, j, &Num);
    	    MSH_SUR_NUM = abs(Num);
    	    MSH_PHYSICAL_ORI = sign(Num);
    	    Tree_Action(ps->Simplexes, _msh_print_simplex);
    	    Tree_Action(ps->SimplexesBase, _msh_print_simplex);
    	    Tree_Action(ps->Quadrangles, _msh_print_quadrangle);
    	  }
    	}
          }
          break;
    
        case MSH_PHYSICAL_VOLUME:
          for(k = 0; k < List_Nbr(ListVolumes); k++) {
            List_Read(ListVolumes, k, &pV);
            for(j = 0; j < List_Nbr(p->Entities); j++) {
              List_Read(p->Entities, j, &Num);
              MSH_VOL_NUM = abs(Num);
              MSH_PHYSICAL_ORI = sign(Num);
              Tree_Action(pV->Simplexes, _msh_print_simplex);
              Tree_Action(pV->SimplexesBase, _msh_print_simplex);
              Tree_Action(pV->Hexahedra, _msh_print_hexahedron);
              Tree_Action(pV->Prisms, _msh_print_prism);
              Tree_Action(pV->Pyramids, _msh_print_pyramid);
            }
          }
          break;
    
        default:
          Msg(GERROR, "Unknown type of physical group");
          break;
        }
    
      }
    
      List_Delete(ListCurves);
      List_Delete(ListSurfaces);
      List_Delete(ListVolumes);
    }
    
    
    static void _get_all_model_points ( std::list<Vertex*> &mp )
    {
      List_T *curves = Tree2List(THEM->Curves);
      std::set<Vertex*> points;
      for(int i = 0; i < List_Nbr(curves); i++){
        Curve *c;
        List_Read(curves, i, &c);
        if (c->Num >=0)
          {
    	if (c->beg && points.find(c->beg) == points.end())
    	  {
    	    points.insert(c->beg);
    	    mp.push_back(c->beg); 
    	  }
    	if (c->end && points.find(c->end) == points.end())
    	  {
    	    points.insert(c->end);
    	    mp.push_back(c->end); 
    	  }
          }
      }
    }
    
    static void _msh_print_all_modelpoints()
    {
      std::list<Vertex*> mp;
      _get_all_model_points ( mp );
      std::list<Vertex*>::iterator it = mp.begin();
      while (it != mp.end())
        {
          Vertex *v = (*it);
          _msh_print_point(v);
          it++;
        }
    }
    
    static void _msh_print_all_curves(void *a, void *b)
    {
      Curve *c = *(Curve **) a; 
    
      Tree_Action(c->Simplexes, _msh_print_simplex);
      Tree_Action(c->SimplexesBase, _msh_print_simplex);
    }
    
    static void _msh_print_all_surfaces(void *a, void *b)
    {
      Surface *s = *(Surface **) a;
      Tree_Action(s->Simplexes, _msh_print_simplex);
      Tree_Action(s->SimplexesBase, _msh_print_simplex);
      Tree_Action(s->Quadrangles, _msh_print_quadrangle);
    }
    
    static void _msh_print_all_simpsurf(void *a, void *b)
    {
      Volume *v = *(Volume **) a;
      Tree_Action(v->Simp_Surf, _msh_print_simplex);
      Tree_Action(v->Quad_Surf, _msh_print_quadrangle);
    }
    
    static void _msh_print_all_linsurf(void *a, void *b)
    {
      Volume *v = *(Volume **) a;
      Tree_Action(v->Lin_Surf, _msh_print_simplex);
    }
    
    static void _msh_print_all_volumes(void *a, void *b)
    {
      Volume *v = *(Volume **) a;
      Tree_Action(v->Simplexes, _msh_print_simplex);
      Tree_Action(v->SimplexesBase, _msh_print_simplex);
      Tree_Action(v->Hexahedra, _msh_print_hexahedron);
      Tree_Action(v->Prisms, _msh_print_prism);
      Tree_Action(v->Pyramids, _msh_print_pyramid);
    }
    
    static void _msh_print_all_elements(Mesh *M)
    {
      MSH_PHYSICAL_NUM = 0;
      MSH_PHYSICAL_ORI = 1;
      MSH_LIN_NUM = MSH_SUR_NUM = MSH_VOL_NUM = 0;
    
    
      _msh_print_all_modelpoints();
    
      if(CTX.mesh.oldxtrude) {
        Tree_Action(M->Volumes, _msh_print_all_simpsurf);
        Tree_Action(M->Volumes, _msh_print_all_linsurf);
      }
      else {
        Tree_Action(M->Curves, _msh_print_all_curves);
        Tree_Action(M->Surfaces, _msh_print_all_surfaces);
      }
    
      Tree_Action(M->Volumes, _msh_print_all_volumes);
    }
    
    static void _msh_process_elements(Mesh *M)
    {
      MSH_ADD = 0;
      MSH_ELEMENT_NUM = 1;
    
      if(!List_Nbr(M->PhysicalGroups) || CTX.mesh.save_all) {
        Msg(INFO, "Saving all elements (discarding physical groups)");
        _msh_print_all_elements(M);
      }
      else
        _msh_print_elements(M);
    
      if(CTX.mesh.msh_file_version == 2.0)
        fprintf(MSHFILE, "$Elements\n");
      else
        fprintf(MSHFILE, "$ELM\n");
    
      fprintf(MSHFILE, "%d\n", MSH_ELEMENT_NUM - 1);
    
      if(MSH_ELEMENT_NUM == 1)
        Msg(WARNING, "No elements to save");
    
      MSH_ADD = 1;
      MSH_ELEMENT_NUM = 1;
      if(!List_Nbr(M->PhysicalGroups) || CTX.mesh.save_all)
        _msh_print_all_elements(M);
      else
        _msh_print_elements(M);
    
      if(CTX.mesh.msh_file_version == 2.0)
        fprintf(MSHFILE, "$EndElements\n");
      else
        fprintf(MSHFILE, "$ENDELM\n");
    }
    
    void Print_Mesh_MSH(Mesh *M, FILE *fp)
    {
      MSHFILE = fp;
      if(CTX.mesh.msh_file_version == 1.0){
        // OK, no header
      }
      else if(CTX.mesh.msh_file_version == 2.0){
        fprintf(MSHFILE, "$MeshFormat\n");
        fprintf(MSHFILE, "%g %d %d\n", CTX.mesh.msh_file_version,
    	    LIST_FORMAT_ASCII, (int)sizeof(double));
        fprintf(MSHFILE, "$EndMeshFormat\n");
      }
      else{
        Msg(GERROR, "Unknown MSH file version to generate (%g)", 
    	CTX.mesh.msh_file_version);
        return;
      }
      _msh_process_nodes(M);
      _msh_process_elements(M);
      Msg(INFO, "%d nodes", MSH_NODE_NUM);
      Msg(INFO, "%d elements", MSH_ELEMENT_NUM - 1);
    }
    
    // Write mesh in UNV format
    
    // IDEAS records
    #define HEADER       151
    #define UNITS        164
    #define NODES        2411
    #define ELEMENTS     2412
    #define RESNODE      55
    #define RESELEM      56
    #define RESVECT      57
    #define GROUPOFNODES 790
    
    // IDEAS elements
    #define BEAM         21
    #define BEAM2        24
    #define THINSHLL     91
    #define THINSHLL2    92
    #define QUAD         94
    #define QUAD2        95 // Ca c'est une impro !!!
    #define SOLIDFEM     111
    #define WEDGE        112
    #define BRICK        115
    #define SOLIDFEM2    118
    
    static int ELEMENT_ID;
    static Tree_T *tree;
    static int UNV_VOL_NUM;
    static FILE *UNVFILE;
    
    static void _unv_process_nodes(Mesh *M)
    {
      Vertex *v;
      List_T *Nodes = Tree2List(M->Vertices);
    
      fprintf(UNVFILE, "%6d\n", -1);
      fprintf(UNVFILE, "%6d\n", NODES);
      int nbnod = List_Nbr(Nodes);
    
      for(int i = 0; i < nbnod; i++) {
        List_Read(Nodes, i, &v);
        int idnod = v->Num;
        double x = v->Pos.X * CTX.mesh.scaling_factor;
        double y = v->Pos.Y * CTX.mesh.scaling_factor;
        double z = v->Pos.Z * CTX.mesh.scaling_factor;
        fprintf(UNVFILE, "%10d%10d%10d%10d\n", idnod, 1, 1, 11);
        char tmp[128];
        // ugly hack to print number with 'D+XX' exponents
        sprintf(tmp, "%25.16E%25.16E%25.16E\n", x, y, z);
        tmp[21] = 'D';
        tmp[46] = 'D';
        tmp[71] = 'D';
        fprintf(UNVFILE, tmp);    
      }
    
      List_Delete(Nodes);
      fprintf(UNVFILE, "%6d\n", -1);
    }
    
    static void _unv_print_record(int num, int fetyp, int geo, int n, int nsup, 
    			      Vertex **v, Vertex **vsup)
    {
      fprintf(UNVFILE, "%10d%10d%10d%10d%10d%10d\n",
    	  num, fetyp, geo, geo, 7, n + nsup);
      if(fetyp == BEAM || fetyp == BEAM2)
        fprintf(UNVFILE, "%10d%10d%10d\n", 0, 0, 0);
      int ntot = 0;
      for(int k = 0; k < n; k++) {
        fprintf(UNVFILE, "%10d", v[k]->Num);
        if(ntot % 8 == 7)
          fprintf(UNVFILE, "\n");
        ntot++;
      }
      for(int k = 0; k < nsup; k++) {
        fprintf(UNVFILE, "%10d", vsup[k]->Num);
        if(ntot % 8 == 7)
          fprintf(UNVFILE, "\n");
        ntot++;
      }
      if(ntot - 1 % 8 != 7)
        fprintf(UNVFILE, "\n");
    }
    
    static void _unv_process_1D_elements(Mesh *m)
    {
      List_T *ListCurves = Tree2List(m->Curves);
      List_T *Elements;
      SimplexBase *sx;
      Curve *c;
    
      for(int i = 0; i < List_Nbr(ListCurves); i++) {
        List_Read(ListCurves, i, &c);
        for(int simtype = 0; simtype < 2; simtype ++){
          Elements = (!simtype) ? Tree2List(c->Simplexes) : Tree2List(c->SimplexesBase);
          for(int j = 0; j < List_Nbr(Elements); j++) {
    	List_Read(Elements, j, &sx);
    	if(sx->VSUP)
    	  _unv_print_record(sx->Num, BEAM2, c->Num, 2, 1, &sx->V[0], sx->VSUP);
    	else 
    	  _unv_print_record(sx->Num, BEAM, c->Num, 2, 0, &sx->V[0], NULL);
          }
          List_Delete(Elements);
        }
      }
    
      List_Delete(ListCurves);
    }
    
    static void _unv_process_2D_elements(Mesh *m)
    {
      List_T *ListSurfaces = Tree2List(m->Surfaces);
      List_T *Elements;
      Surface *s;
      SimplexBase *sx;
      Quadrangle *qx;
    
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
          
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          Elements = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Elements); j++) {
    	List_Read(Elements, j, &sx);
    	if(sx->VSUP)
    	  _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 3, &sx->V[0], sx->VSUP);
    	else
    	  _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 0, &sx->V[0], NULL);
          }
          List_Delete(Elements);
        }
        
        // quadrangles
        Elements = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Elements); j++) {
          List_Read(Elements, j, &qx);
          if(qx->VSUP)
    	_unv_print_record(abs(qx->Num), QUAD, s->Num, 4, 4+1, &qx->V[0], qx->VSUP);
          else
    	_unv_print_record(abs(qx->Num), QUAD2, s->Num, 4, 0, &qx->V[0], NULL);
        }
        List_Delete(Elements);
      }
      List_Delete(ListSurfaces);
    }
    
    static void _unv_process_3D_elements(Mesh *m)
    {
      List_T *ListVolumes = Tree2List(m->Volumes);
      List_T *Elements;
      SimplexBase *sx;
      Hexahedron *hx;
      Prism *px;
      Volume *v;
    
      for(int i = 0; i < List_Nbr(ListVolumes); i++) {
        List_Read(ListVolumes, i, &v);
    
        // Tets
        for(int simtype = 0; simtype < 2; simtype++){
          Elements = (!simtype) ? Tree2List(v->Simplexes) : Tree2List(v->SimplexesBase);
          for(int j = 0; j < List_Nbr(Elements); j++) {
    	List_Read(Elements, j, &sx);
    	if(sx->Volume_Simplexe() < 0) {
    	  Vertex *temp;
    	  temp = sx->V[0];
    	  sx->V[0] = sx->V[1];
    	  sx->V[1] = temp;
    	  if(sx->VSUP){
    	    temp = sx->VSUP[1];
    	    sx->VSUP[1] = sx->VSUP[2];
    	    sx->VSUP[2] = temp;
    	    temp = sx->VSUP[5];
    	    sx->VSUP[5] = sx->VSUP[3];
    	    sx->VSUP[3] = temp;
    	  }
    	}
    	if(sx->VSUP)
    	  _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 6, &sx->V[0], sx->VSUP);
    	else
    	  _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 0, &sx->V[0], NULL);
          }
          List_Delete(Elements);
        }
    
        // Prisms
        Elements = Tree2List(v->Prisms);
        for(int j = 0; j < List_Nbr(Elements); j++) {
          List_Read(Elements, j, &px);
          if(px->VSUP)
    	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 9+3, &px->V[0], px->VSUP);
          else
    	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 0, &px->V[0], NULL);
        }
        List_Delete(Elements);
    
        // Hexas
        Elements = Tree2List(v->Hexahedra);
        for(int j = 0; j < List_Nbr(Elements); j++) {
          List_Read(Elements, j, &hx);
          if(hx->VSUP)
    	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 12+6+1, &hx->V[0], hx->VSUP);
          else
    	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 0, &hx->V[0], NULL);
        }
        List_Delete(Elements);
      }
      List_Delete(ListVolumes);
    }
    
    static void _unv_add_vertex(void *a, void *b)
    {
      Vertex *v = *(Vertex **) a;
      if(Tree_Search(tree, &v->Num))
        return;
      Tree_Add(tree, &v->Num);
      fprintf(UNVFILE, "%10d%10d%2d%2d%2d%2d%2d%2d\n", v->Num, 1, 0, 1, 0, 0, 0, 0);
      fprintf(UNVFILE, "   0.0000000000000000D+00   1.0000000000000000D+00   0.0000000000000000D+00\n");
      fprintf(UNVFILE, "   0.0000000000000000D+00   0.0000000000000000D+00   0.0000000000000000D+00\n");
      fprintf(UNVFILE, "%10d%10d%10d%10d%10d%10d\n", 0, 0, 0, 0, 0, 0);
    }
    
    static void _unv_add_simplex_vertices(void *a, void *b)
    {
      SimplexBase *s = *(SimplexBase **) a;
      if(s->iEnt != UNV_VOL_NUM)
        return;
      for(int i = 0; i < 4; i++)
        _unv_add_vertex(&s->V[i], NULL);
    }
    
    static void _unv_add_hexahedron_vertices(void *a, void *b)
    {
      Hexahedron *h = *(Hexahedron **) a;
      if(h->iEnt != UNV_VOL_NUM)
        return;
      for(int i = 0; i < 8; i++)
        _unv_add_vertex(&h->V[i], NULL);
    }
    
    static void _unv_add_prism_vertices(void *a, void *b)
    {
      Prism *p = *(Prism **) a;
      if(p->iEnt != UNV_VOL_NUM)
        return;
      for(int i = 0; i < 6; i++)
        _unv_add_vertex(&p->V[i], NULL);
    }
    
    static void _unv_add_pyramid_vertices(void *a, void *b)
    {
      Pyramid *p = *(Pyramid **) a;
      if(p->iEnt != UNV_VOL_NUM)
        return;
      for(int i = 0; i < 5; i++)
        _unv_add_vertex(&p->V[i], NULL);
    }
    
    static void _unv_process_groups(Mesh *m)
    {
      Volume *pV;
      Surface *ps, s;
      Curve *pc, c;
      Vertex *pv, v;
      PhysicalGroup *p;
      List_T *ListVolumes;
    
      for(int i = 0; i < List_Nbr(m->PhysicalGroups); i++) {
    
        List_Read(m->PhysicalGroups, i, &p);
    
        fprintf(UNVFILE, "%6d\n", -1);
        fprintf(UNVFILE, "%6d\n", GROUPOFNODES);
        fprintf(UNVFILE, "%10d%10d\n", p->Num, 1);
        fprintf(UNVFILE, "LOAD SET %2d\n", 1);
    
        tree = Tree_Create(sizeof(int), fcmp_absint);
    
        switch (p->Typ) {
        case MSH_PHYSICAL_VOLUME:
          ListVolumes = Tree2List(m->Volumes);
          for(int k = 0; k < List_Nbr(ListVolumes); k++) {
            List_Read(ListVolumes, k, &pV);
            for(int j = 0; j < List_Nbr(p->Entities); j++) {
              List_Read(p->Entities, j, &UNV_VOL_NUM);
              Tree_Action(pV->Simplexes, _unv_add_simplex_vertices);
              Tree_Action(pV->SimplexesBase, _unv_add_simplex_vertices);
              Tree_Action(pV->Hexahedra, _unv_add_hexahedron_vertices);
              Tree_Action(pV->Prisms, _unv_add_prism_vertices);
              Tree_Action(pV->Pyramids, _unv_add_pyramid_vertices);
            }
          }
          List_Delete(ListVolumes);
          break;
        case MSH_PHYSICAL_SURFACE:
          for(int j = 0; j < List_Nbr(p->Entities); j++) {
            ps = &s;
            List_Read(p->Entities, j, &ps->Num);
            if(Tree_Query(m->Surfaces, &ps))
              Tree_Action(ps->Vertices, _unv_add_vertex);
          }
          break;
        case MSH_PHYSICAL_LINE:
          for(int j = 0; j < List_Nbr(p->Entities); j++) {
            pc = &c;
            List_Read(p->Entities, j, &pc->Num);
            if(Tree_Query(m->Curves, &pc))
              for(int k = 0; k < List_Nbr(pc->Vertices); k++)
                _unv_add_vertex(List_Pointer(pc->Vertices, k), NULL);
          }
          break;
        case MSH_PHYSICAL_POINT:
          for(int j = 0; j < List_Nbr(p->Entities); j++) {
            pv = &v;
            List_Read(p->Entities, j, &pv->Num);
            if(Tree_Query(m->Vertices, &pv))
              _unv_add_vertex(&pv, NULL);
          }
          break;
        }
    
        Tree_Delete(tree);
    
        fprintf(UNVFILE, "%6d\n", -1);
      }
    }
    
    void Print_Mesh_UNV(Mesh *M, FILE *fp)
    {
      UNVFILE = fp;
      _unv_process_nodes(M);
      fprintf(UNVFILE, "%6d\n", -1);
      fprintf(UNVFILE, "%6d\n", ELEMENTS);
      ELEMENT_ID = 1;
      _unv_process_3D_elements(M);
      _unv_process_2D_elements(M);
      _unv_process_1D_elements(M);
      fprintf(UNVFILE, "%6d\n", -1);
      _unv_process_groups(M);
    }
    
    // Write mesh in Gref format
    
    static FILE *GREFFILE;
    
    static void _gref_consecutive_nodes(Mesh *M, Tree_T *ConsecutiveNTree,
    				    Tree_T *ConsecutiveETree)
    {
      SimplexBase *sx;
      Quadrangle *qx;
      Surface *s;
      int nbnod, nbedges, nbdof;
    
      int newnum = 0;
    
      List_T *ListSurfaces = Tree2List(M->Surfaces);
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Triangles); j++) {
    	List_Read(Triangles, j, &sx);
    	for(int k = 0; k < 3; k++) {
    	  if(sx->V[k]->Frozen >= 0) {
    	    sx->V[k]->Frozen = --newnum;
    	    Tree_Insert(ConsecutiveNTree, &(sx->V[k]));
    	  }
    	}
          }
          List_Delete(Triangles);
        }
        // quads
        List_T *Quadrangles = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Quadrangles); j++) {
          List_Read(Quadrangles, j, &qx);
          for(int k = 0; k < 4; k++) {
            if(qx->V[k]->Frozen >= 0) {
              qx->V[k]->Frozen = --newnum;
              Tree_Insert(ConsecutiveNTree, &(qx->V[k]));
            }
          }
        }
        List_Delete(Quadrangles);
      }
    
      nbnod = -newnum;
    
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Triangles); j++) {
    	List_Read(Triangles, j, &sx);
    	for(int k = 0; k < 3; k++) {
    	  if(sx->VSUP[k]->Frozen >= 0) {
    	    sx->VSUP[k]->Frozen = --newnum;
    	    Tree_Insert(ConsecutiveETree, &(sx->VSUP[k]));
    	  }
    	}
          }
          List_Delete(Triangles);
        }
        // quads
        List_T *Quadrangles = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Quadrangles); j++) {
          List_Read(Quadrangles, j, &qx);
          for(int k = 0; k < 4; k++) {
    	if(qx->VSUP[k]->Frozen >= 0) {
    	  qx->VSUP[k]->Frozen = --newnum;
    	  Tree_Insert(ConsecutiveETree, &(qx->VSUP[k]));
    	}
          }
        }
        List_Delete(Quadrangles);
      }
    
      List_Delete(ListSurfaces);
    
      nbedges = -newnum - nbnod;
      nbdof = nbnod + nbedges;
      Msg(INFO, "%d Dofs", nbdof);
    }
    
    static void _gref_end_consecutive_nodes(Mesh *M)
    {
      SimplexBase *sx;
      Quadrangle *qx;
      Surface *s;
    
      List_T *ListSurfaces = Tree2List(M->Surfaces);
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Triangles); j++) {
    	List_Read(Triangles, j, &sx);
    	for(int k = 0; k < 3; k++)
    	  sx->V[k]->Frozen = 0;
    	for(int k = 0; k < 3; k++)
    	  sx->VSUP[k]->Frozen = 0;
          }
          List_Delete(Triangles);
        }
        // quads
        List_T *Quadrangles = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Quadrangles); j++) {
          List_Read(Quadrangles, j, &qx);
          for(int k = 0; k < 4; k++)
            qx->V[k]->Frozen = 0;
          for(int k = 0; k < 4; k++)
    	qx->VSUP[k]->Frozen = 0;
        }
        List_Delete(Quadrangles);
      }
      List_Delete(ListSurfaces);
    }
    
    static int _gref_compare_frozen(const void *a, const void *b)
    {
      Vertex *q, *w;
      q = *(Vertex **) a;
      w = *(Vertex **) b;
      return w->Frozen - q->Frozen;
    }
    
    static int _gref_process_nodes(Mesh *M, Tree_T *ConsecutiveNTree, 
    			       Tree_T *ConsecutiveETree)
    {
      int i, nbtriqua;
      Vertex *v;
      Surface *s;
      List_T *Nodes;
    
      Tree_Action(M->Curves, Degre2_Curve);
      Tree_Action(M->Surfaces, Degre2_Surface);
     
      List_T *ListSurfaces = Tree2List(M->Surfaces);
      nbtriqua = 0;
      for(i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        nbtriqua += Tree_Nbr(s->Simplexes) + Tree_Nbr(s->SimplexesBase);
        nbtriqua += Tree_Nbr(s->Quadrangles);
      }
      List_Delete(ListSurfaces);
    
      _gref_consecutive_nodes(M, ConsecutiveNTree, ConsecutiveETree);
    
      fprintf(GREFFILE, "%d %d %d\n", nbtriqua, Tree_Nbr(ConsecutiveNTree),
              Tree_Nbr(ConsecutiveNTree) + Tree_Nbr(ConsecutiveETree));
    
      Nodes = Tree2List(ConsecutiveNTree);
      for(i = 0; i < List_Nbr(Nodes); i++) {
        List_Read(Nodes, i, &v);
        fprintf(GREFFILE, "%21.16e ", v->Pos.X * CTX.mesh.scaling_factor);
        if(i % 3 == 2)
          fprintf(GREFFILE, "\n");
      }
      if((List_Nbr(Nodes) - 1) % 3 != 2)
        fprintf(GREFFILE, "\n");
      for(i = 0; i < List_Nbr(Nodes); i++) {
        List_Read(Nodes, i, &v);
        fprintf(GREFFILE, "%21.16e ", v->Pos.Y * CTX.mesh.scaling_factor);
        if(i % 3 == 2)
          fprintf(GREFFILE, "\n");
      }
      if((List_Nbr(Nodes) - 1) % 3 != 2)
        fprintf(GREFFILE, "\n");
      i = Tree_Nbr(ConsecutiveNTree);
      List_Delete(Nodes);
      return i;
    }
    
    static int _gref_find_physical(Vertex *v, Mesh *m)
    {
      PhysicalGroup *p;
      Curve *c;
      int i, j;
      for(i = 0; i < List_Nbr(m->PhysicalGroups); i++) {
        List_Read(m->PhysicalGroups, i, &p);
        if(p->Typ == MSH_PHYSICAL_POINT) {
          if(List_Search(p->Entities, &v->Num, fcmp_absint)) {
            return p->Num;
          }
        }
      }
    
      if(v->ListCurves) {
        for(i = 0; i < List_Nbr(m->PhysicalGroups); i++) {
          List_Read(m->PhysicalGroups, i, &p);
          if(p->Typ == MSH_PHYSICAL_LINE) {
            for(j = 0; j < List_Nbr(v->ListCurves); j++) {
              List_Read(v->ListCurves, j, &c);
              if(List_Search(p->Entities, &c->Num, fcmp_absint)) {
                return p->Num;
              }
            }
          }
        }
      }
      return 0;
    }
    
    static void _gref_process_boundary_conditions(Mesh *M, Tree_T *TRN, Tree_T *TRE)
    {
      int i, ent;
      Vertex *v;
    
      List_T *Nodes = Tree2List(TRN);
      for(i = 0; i < List_Nbr(Nodes); i++) {
        List_Read(Nodes, i, &v);
        ent = _gref_find_physical(v, M);
        fprintf(GREFFILE, "%d %d ", ent, ent);
        if(i % 3 == 2)
          fprintf(GREFFILE, "\n");
      }
      if((List_Nbr(Nodes) - 1) % 3 != 2)
        fprintf(GREFFILE, "\n");
      List_Delete(Nodes);
    
      Nodes = Tree2List(TRE);
      for(i = 0; i < List_Nbr(Nodes); i++) {
        List_Read(Nodes, i, &v);
        ent = _gref_find_physical(v, M);
        fprintf(GREFFILE, "%d %d ", ent, ent);
        if(i % 3 == 2)
          fprintf(GREFFILE, "\n");
      }
      if((List_Nbr(Nodes) - 1) % 3 != 2)
        fprintf(GREFFILE, "\n");
      List_Delete(Nodes);
    }
    
    static void _gref_process_elements(Mesh *M, int nn)
    {
      SimplexBase *sx;
      Quadrangle *qx;
      Surface *s;
    
      List_T *ListSurfaces = Tree2List(M->Surfaces);
    
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Triangles); j++) {
    	List_Read(Triangles, j, &sx);
    	fprintf(GREFFILE, "%d %d %d\n", -sx->V[0]->Frozen,
    		-sx->V[1]->Frozen, -sx->V[2]->Frozen);
          }
          List_Delete(Triangles);
        }
        // quads
        List_T *Quadrangles = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Quadrangles); j++) {
          List_Read(Quadrangles, j, &qx);
          fprintf(GREFFILE, "%d %d %d %d\n", -qx->V[0]->Frozen,
    	      -qx->V[1]->Frozen, -qx->V[2]->Frozen, -qx->V[3]->Frozen);
        }
        List_Delete(Quadrangles);
      }
    
      // Degres de Liberte
      for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &s);
        // triangles
        for(int simtype = 0; simtype < 2; simtype++){
          List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
          for(int j = 0; j < List_Nbr(Triangles); j++) {
    	List_Read(Triangles, j, &sx);
    	fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d\n",
    		-2 * sx->V[0]->Frozen - 1,
    		-2 * sx->V[0]->Frozen,
    		-2 * sx->VSUP[0]->Frozen - 1,
    		-2 * sx->VSUP[0]->Frozen,
    		-2 * sx->V[1]->Frozen - 1,
    		-2 * sx->V[1]->Frozen,
    		-2 * sx->VSUP[1]->Frozen - 1,
    		-2 * sx->VSUP[1]->Frozen,
    		-2 * sx->V[2]->Frozen - 1,
    		-2 * sx->V[2]->Frozen,
    		-2 * sx->VSUP[2]->Frozen - 1, -2 * sx->VSUP[2]->Frozen);
          }
          List_Delete(Triangles);
        }
        // quads
        List_T *Quadrangles = Tree2List(s->Quadrangles);
        for(int j = 0; j < List_Nbr(Quadrangles); j++) {
          List_Read(Quadrangles, j, &qx);
          fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
    	      -2 * qx->V[0]->Frozen - 1,
    	      -2 * qx->V[0]->Frozen,
    	      -2 * qx->VSUP[0]->Frozen - 1,
    	      -2 * qx->VSUP[0]->Frozen,
    	      -2 * qx->V[1]->Frozen - 1,
    	      -2 * qx->V[1]->Frozen,
    	      -2 * qx->VSUP[1]->Frozen - 1,
    	      -2 * qx->VSUP[1]->Frozen,
    	      -2 * qx->V[2]->Frozen - 1,
    	      -2 * qx->V[2]->Frozen,
    	      -2 * qx->VSUP[2]->Frozen - 1,
    	      -2 * qx->VSUP[2]->Frozen,
    	      -2 * qx->V[3]->Frozen - 1,
    	      -2 * qx->V[3]->Frozen,
    	      -2 * qx->VSUP[3]->Frozen - 1, -2 * qx->VSUP[3]->Frozen);
        }
        List_Delete(Quadrangles);
      }
      List_Delete(ListSurfaces);
    }
    
    void Print_Mesh_GREF(Mesh *M, FILE *fp)
    {
      GREFFILE = fp;
      Tree_T *TRN = Tree_Create(sizeof(Vertex *), _gref_compare_frozen);
      Tree_T *TRE = Tree_Create(sizeof(Vertex *), _gref_compare_frozen);
      _gref_process_nodes(M, TRN, TRE);
      _gref_process_elements(M, Tree_Nbr(TRN));
      _gref_process_boundary_conditions(M, TRN, TRE);
      Tree_Delete(TRN);
      Tree_Delete(TRE);
      _gref_end_consecutive_nodes(M);
    }
    
    // Write mesh in VRML 1 format
    
    static List_T *wrlnodes = NULL;
    static FILE *WRLFILE;
    
    static void _wrl_print_node(void *a, void *b)
    {
      Vertex *V = *(Vertex **) a;
      fprintf(WRLFILE, "%.16g %.16g %.16g,\n",
              V->Pos.X * CTX.mesh.scaling_factor,
              V->Pos.Y * CTX.mesh.scaling_factor,
              V->Pos.Z * CTX.mesh.scaling_factor);
      List_Add(wrlnodes, &V->Num);
    }
    
    static void _wrl_process_nodes(Mesh *M)
    {
      if(!wrlnodes)
        wrlnodes = List_Create(Tree_Size(M->Vertices), 100, sizeof(int));
      else
        List_Reset(wrlnodes);
      fprintf(WRLFILE, "#VRML V1.0 ascii\n");
      fprintf(WRLFILE, "#created by Gmsh\n");
      fprintf(WRLFILE, "Coordinate3 {\n");
      fprintf(WRLFILE, "  point [\n");
      Tree_Action(M->Vertices, _wrl_print_node);
      fprintf(WRLFILE, "  ]\n");
      fprintf(WRLFILE, "}\n");
    }
    
    static void _wrl_print_line(void *a, void *b)
    {
      SimplexBase *s = *(SimplexBase **) a;
      for(int i = 0; i < 2; i++){
        if(s->V[i]){
          int j = List_ISearch(wrlnodes, &s->V[i]->Num, fcmp_int);
          if(j < 0)
    	Msg(GERROR, "Unknown node %d in line %d", s->V[i]->Num, s->Num);
          else
    	fprintf(WRLFILE, "%d,", j);
        }
      }
      fprintf(WRLFILE, "-1,\n");
    }
    
    static void _wrl_print_triangle(void *a, void *b)
    {
      SimplexBase *s = *(SimplexBase **) a;
      for(int i = 0; i < 3; i++){
        if(s->V[i]){
          int j = List_ISearch(wrlnodes, &s->V[i]->Num, fcmp_int);
          if(j < 0)
    	Msg(GERROR, "Unknown node %d in triangle %d", s->V[i]->Num, s->Num);
          else
    	fprintf(WRLFILE, "%d,", j);
        }
      }
      fprintf(WRLFILE, "-1,\n");
    }
    
    static void _wrl_print_quadrangle(void *a, void *b)
    {
      Quadrangle *q = *(Quadrangle **) a;
      for(int i = 0; i < 4; i++){
        if(q->V[i]){
          int j = List_ISearch(wrlnodes, &q->V[i]->Num, fcmp_int);
          if(j < 0)
    	Msg(GERROR, "Unknown node %d in quadrangle %d", q->V[i]->Num, q->Num);
          else
    	fprintf(WRLFILE, "%d,", j);
        }
      }
      fprintf(WRLFILE, "-1,\n");
    }
    
    static void _wrl_print_all_curves(void *a, void *b)
    {
      Curve *c = *(Curve **) a;
      if(c->Num < 0)
        return;
      fprintf(WRLFILE, "DEF Curve%d IndexedLineSet {\n", c->Num);
      fprintf(WRLFILE, "  coordIndex [\n");
      Tree_Action(c->Simplexes, _wrl_print_line);
      Tree_Action(c->SimplexesBase, _wrl_print_line);
      fprintf(WRLFILE, "  ]\n");
      fprintf(WRLFILE, "}\n");
    }
    
    static void _wrl_print_all_surfaces(void *a, void *b)
    {
      Surface *s = *(Surface **) a;
      fprintf(WRLFILE, "DEF Surface%d IndexedFaceSet {\n", s->Num);
      fprintf(WRLFILE, "  coordIndex [\n");
      Tree_Action(s->Simplexes, _wrl_print_triangle);
      Tree_Action(s->SimplexesBase, _wrl_print_triangle);
      Tree_Action(s->Quadrangles, _wrl_print_quadrangle);
      fprintf(WRLFILE, "  ]\n");
      fprintf(WRLFILE, "}\n");
    }
    
    static void _wrl_process_elements(Mesh *M)
    {
      if(!wrlnodes)
        Msg(GERROR, "VRML node list does not exist");
      else {
        Tree_Action(M->Curves, _wrl_print_all_curves);
        Tree_Action(M->Surfaces, _wrl_print_all_surfaces);
      }
    }
    
    void Print_Mesh_WRL(Mesh *M, FILE *fp)
    {
      WRLFILE = fp;
      _wrl_process_nodes(M);
      _wrl_process_elements(M);
    }
    
    // Write surface mesh in STL format
    
    static FILE *STLFILE;
    
    static void _stl_print_triangle(void *a, void *b)
    {
      SimplexBase *s = *(SimplexBase **) a;
      
      if(!s->V[2]) return;
      
      double n[3];
      normal3points(s->V[0]->Pos.X, s->V[0]->Pos.Y, s->V[0]->Pos.Z, 
    		s->V[1]->Pos.X, s->V[1]->Pos.Y, s->V[1]->Pos.Z, 
    		s->V[2]->Pos.X, s->V[2]->Pos.Y, s->V[2]->Pos.Z, n);
    
      fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
      fprintf(STLFILE, "  outer loop\n");
      fprintf(STLFILE, "    vertex %g %g %g\n", s->V[0]->Pos.X, s->V[0]->Pos.Y, s->V[0]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", s->V[1]->Pos.X, s->V[1]->Pos.Y, s->V[1]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", s->V[2]->Pos.X, s->V[2]->Pos.Y, s->V[2]->Pos.Z);
      fprintf(STLFILE, "  endloop\n");
      fprintf(STLFILE, "endfacet\n");
    }
    
    static void _stl_print_quadrangle(void *a, void *b)
    {
      Quadrangle *q = *(Quadrangle **) a;
      
      double n[3];
      normal3points(q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z, 
    		q->V[1]->Pos.X, q->V[1]->Pos.Y, q->V[1]->Pos.Z, 
    		q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z, n);
    
      fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
      fprintf(STLFILE, "  outer loop\n");
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[1]->Pos.X, q->V[1]->Pos.Y, q->V[1]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z);
      fprintf(STLFILE, "  endloop\n");
      fprintf(STLFILE, "endfacet\n");
      fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
      fprintf(STLFILE, "  outer loop\n");
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z);
      fprintf(STLFILE, "    vertex %g %g %g\n", q->V[3]->Pos.X, q->V[3]->Pos.Y, q->V[3]->Pos.Z);
      fprintf(STLFILE, "  endloop\n");
      fprintf(STLFILE, "endfacet\n");
    }
    
    static void _stl_print_all_surfaces(void *a, void *b)
    {
      Surface *s = *(Surface **) a;
      Tree_Action(s->Simplexes, _stl_print_triangle);
      Tree_Action(s->SimplexesBase, _stl_print_triangle);
      Tree_Action(s->Quadrangles, _stl_print_quadrangle);
    }
    
    void Print_Mesh_STL(Mesh *M, FILE *fp)
    {
      STLFILE = fp;
      fprintf(STLFILE, "solid Created by Gmsh\n");
      Tree_Action(M->Surfaces, _stl_print_all_surfaces);
      fprintf(STLFILE, "endsolid Created by Gmsh\n");
    }
    
    // Write mesh in DMG format
    
    static int _dmg_is_topologic(Vertex *v, List_T *curves)
    {
      Curve *c;
      for(int i = 0; i < List_Nbr(curves); i++) {
        List_Read(curves, i, &c);
        if(!compareVertex(&v, &c->beg))
          return 1;
      }
      return 0;
    }
    
    void Print_Mesh_DMG(Mesh *m, FILE *fp)
    {
      int i, j;
      List_T *ll, *l;
      Vertex *v;
      Curve *c;
      Surface *s;
      int k;
    
      l = Tree2List(m->Points);
      ll = Tree2List(m->Curves);
    
      k = 0;
      for(i = 0; i < List_Nbr(l); i++) {
        List_Read(l, i, &v);
        if(_dmg_is_topologic(v, ll)) {
          k++;
        }
      }
    
      // write first the global infos 
    
      fprintf(fp, "%d %d %d %d \n", Tree_Nbr(m->Volumes),
              Tree_Nbr(m->Surfaces),
              Tree_Nbr(m->Curves) / 2,     // the 2 is for the reverse curves
              k);
    
      // then write the bounding box
    
      fprintf(fp, "%12.5E %12.5E %12.5E \n", 
    	  CTX.min[0], CTX.min[1], CTX.min[2]);
      fprintf(fp, "%12.5E %12.5E %12.5E \n", 
    	  CTX.max[0], CTX.max[1], CTX.max[2]);
    
      // write the points
      k = 0;
      for(i = 0; i < List_Nbr(l); i++) {
        List_Read(l, i, &v);
        if(_dmg_is_topologic(v, ll)) {
          v->Frozen = k++;
          fprintf(fp, "%d %12.5E %12.5E %12.5E \n", 
    	      v->Frozen, v->Pos.X, v->Pos.Y, v->Pos.Z);
        }
      }
      List_Delete(l);
    
      // write the curves
      l = ll;
      k = 0;
      for(i = 0; i < List_Nbr(l); i++) {
        List_Read(l, i, &c);
        if(c->Num > 0) {
          c->ipar[3] = k;
          Curve *cinv = FindCurve(-c->Num, m);
          cinv->ipar[3] = k++;
          fprintf(fp, "%d %d %d \n", 
    	      c->ipar[3], c->beg->Frozen, c->end->Frozen);
        }
      }
    
      List_Delete(l);
    
      // write the surfaces
      l = Tree2List(m->Surfaces);
    
      for(i = 0; i < List_Nbr(l); i++) {
        List_Read(l, i, &s);
    
        int numEdgeLoop[2000], iLoop = 0;
        Vertex *beg = NULL;
        numEdgeLoop[iLoop] = 0;
        int deb = 1;
        for(j = 0; j < List_Nbr(s->Generatrices); j++) {
          List_Read(s->Generatrices, j, &c);
          if(deb) {
            beg = c->beg;
            deb = 0;
          }
          Msg(INFO, "beg->%d end->%d", c->beg->Num, c->end->Num);
          (numEdgeLoop[iLoop])++;
          if(c->end == beg) {
            iLoop++;
            numEdgeLoop[iLoop] = 0;
            deb = 1;
          }
        }
        s->ipar[3] = i;
        fprintf(fp, "%d %d\n", i, iLoop);
        int iEdge = 0;
        for(k = 0; k < iLoop; k++) {
          fprintf(fp, "%d ", numEdgeLoop[k]);
          for(j = 0; j < numEdgeLoop[k]; j++) {
            List_Read(s->Generatrices, iEdge++, &c);
            fprintf(fp, "%d %d ", abs(c->ipar[3]), (c->Num > 0) ? 1 : -1);
          }
          fprintf(fp, "\n");
        }
      }
      List_Delete(l);
    
      // write the volumes (2 b continued...)
    }
    
    // Write mesh in Plot3D format, ASCII structured, multi-zone
    // FIXME: still need to implement this for extruded grids
    
    static FILE *P3DFILE;
    
    int _p3d_cmp_entities(const void *a, const void *b)
    {
      Element **e1 = (Element **) a;
      Element **e2 = (Element **) b;
    
      return (*e1)->iEnt - (*e2)->iEnt;
    }
    
    int _p3d_cmp_surf_num(const void *a, const void *b)
    {
      Surface **e1 = (Surface **) a;
      Surface **e2 = (Surface **) b;
    
      return (*e1)->Num - (*e2)->Num;
    }
    
    int _p3d_cmp_vol_num(const void *a, const void *b)
    {
      Volume **e1 = (Volume **) a;
      Volume **e2 = (Volume **) b;
    
      return (*e1)->Num - (*e2)->Num;
    }
    
    void _p3d_print_quads(List_T *ListQuads, int Nu, int Nv)
    {
      int i = 0, j = 0, c = 1, curQuad = 0;
      double coord;
      Quadrangle *pQ;
    
      for (c = 0; c < 3; c++) {
        for (j = 0; j < (Nu - 1); j++) {
          for (i = 0; i < (Nv - 1); i++) {
            curQuad = i + (j * (Nv - 1));
            List_Read(ListQuads, curQuad, &pQ);
            coord = (c==0) ? pQ->V[0]->Pos.X : ((c==1) ? pQ->V[0]->Pos.Y : pQ->V[0]->Pos.Z);
            fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
          }
          coord = (c==0) ? pQ->V[3]->Pos.X : ((c==1) ? pQ->V[3]->Pos.Y : pQ->V[3]->Pos.Z);
          fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
        } j--;
        for (i = 0; i < (Nv - 1); i++) {
          curQuad = i + (j * (Nv - 1));
          List_Read(ListQuads, curQuad, &pQ);
          coord = (c==0) ? pQ->V[1]->Pos.X : ((c==1) ? pQ->V[1]->Pos.Y : pQ->V[1]->Pos.Z);
          fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
        }
        coord = (c==0) ? pQ->V[2]->Pos.X : ((c==1) ? pQ->V[2]->Pos.Y : pQ->V[2]->Pos.Z);
        fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
      }
    }
    
    void _p3d_print_hex(List_T *ListHex, int Nu, int Nv, int Nw)
    {
      int i = 0, j = 0, k = 0, c = 0, curHex = 0;
      double coord;
      Hexahedron *pH;
    
      for (c = 0; c < 3; c++) {
        for (k = 0; k < (Nu - 1); k++) {
          for (j = 0; j < (Nv - 1); j++) {
            for (i = 0; i < (Nw - 1); i++) {
              curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
              List_Read(ListHex, curHex, &pH);
              coord = (c==0) ? pH->V[0]->Pos.X : ((c==1) ? pH->V[0]->Pos.Y : pH->V[0]->Pos.Z);
              fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
            }
           coord = (c==0) ? pH->V[4]->Pos.X : ((c==1) ? pH->V[4]->Pos.Y : pH->V[4]->Pos.Z);
           fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
          } j--;
          for (i = 0; i < (Nw - 1); i++) {
            curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
            List_Read(ListHex, curHex, &pH);
            coord = (c==0) ? pH->V[3]->Pos.X : ((c==1) ? pH->V[3]->Pos.Y : pH->V[3]->Pos.Z);
            fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
          }
          coord = (c==0) ? pH->V[7]->Pos.X : ((c==1) ? pH->V[7]->Pos.Y : pH->V[7]->Pos.Z);
          fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
        } k--;
        for (j = 0; j < (Nv - 1); j++) {
          for (i = 0; i < (Nw - 1); i++) {
            curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
            List_Read(ListHex, curHex, &pH);
            coord = (c==0) ? pH->V[1]->Pos.X : ((c==1) ? pH->V[1]->Pos.Y : pH->V[1]->Pos.Z);
            fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
          }
          coord = (c==0) ? pH->V[5]->Pos.X : ((c==1) ? pH->V[5]->Pos.Y : pH->V[5]->Pos.Z);
          fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
        } j--;
        for (i = 0; i < (Nw - 1); i++) {
          curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
          List_Read(ListHex, curHex, &pH);
          coord = (c==0) ? pH->V[2]->Pos.X : ((c==1) ? pH->V[2]->Pos.Y : pH->V[2]->Pos.Z);
          fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
        }
        coord = (c==0) ? pH->V[6]->Pos.X : ((c==1) ? pH->V[6]->Pos.Y : pH->V[6]->Pos.Z);
        fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
      }
    }
    
    void _p3d_print_all_elements(Mesh *M)
    {
      int i, Nv, Nu, ElemCnt = 0;
      Volume *pV;
      Surface *pS;
      List_T *ListQuads;
      List_T *ListHex;
    
      List_T *ListSurfaces = Tree2List(M->Surfaces);
      List_T *ListVolumes = Tree2List(M->Volumes);
    
      List_T *ListStructSurf = List_Create(1,1,sizeof(Surface *));
      List_T *ListStructVol = List_Create(1,1,sizeof(Volume *));
    
      List_Sort(ListSurfaces, _p3d_cmp_surf_num);
      List_Sort(ListVolumes, _p3d_cmp_vol_num);
    
      for(i = 0; i < List_Nbr(ListSurfaces); i++) {
        List_Read(ListSurfaces, i, &pS);
        if (pS->Nu &&
            pS->Nv &&
            (Tree_Nbr(pS->Quadrangles) == (pS->Nu - 1) * (pS->Nv - 1))){
            ElemCnt++;
            List_Add(ListStructSurf, &pS);
        }
      }
    
      for(i = 0; i < List_Nbr(ListVolumes); i++) {
        List_Read(ListVolumes, i, &pV);
        List_Read(pV->Surfaces, 0, &pS);
        Nu = pS->Nu;
        Nv = pS->Nv;
        List_Read(pV->Surfaces, 1, &pS);
        if (Nu &&
            Nv &&
            pS->Nv &&
            (Tree_Nbr(pV->Hexahedra) == (Nu - 1) * (Nv - 1) * (pS->Nv - 1))){
            ElemCnt++;
            List_Add(ListStructVol, &pV);
        }
      }
    
      fprintf(P3DFILE, "%d\n", ElemCnt);
    
      for(i = 0; i < List_Nbr(ListStructSurf); i++) {
        List_Read(ListStructSurf, i, &pS);
        fprintf(P3DFILE, "%d %d 1\n", pS->Nv, pS->Nu);
      }
    
      for(i = 0; i < List_Nbr(ListStructVol); i++) {
        List_Read(ListStructVol, i, &pV);
        List_Read(pV->Surfaces, 0, &pS);
        Nu = pS->Nu;
        Nv = pS->Nv;
        List_Read(pV->Surfaces, 1, &pS);
        fprintf(P3DFILE, "%d %d %d\n", pS->Nv, Nv, Nu);
      }
    
      for(i = 0; i < List_Nbr(ListStructSurf); i++) {
        List_Read(ListStructSurf, i, &pS);
        ListQuads = Tree2List(pS->Quadrangles);
        List_Sort(ListQuads, _p3d_cmp_entities);
        _p3d_print_quads(ListQuads, pS->Nu, pS->Nv);
        List_Delete(ListQuads);
      }
    
      for(i = 0; i < List_Nbr(ListStructVol); i++) {
        List_Read(ListStructVol, i, &pV);
        List_Read(pV->Surfaces, 0, &pS);
        Nu = pS->Nu;
        Nv = pS->Nv;
        List_Read(pV->Surfaces, 1, &pS);
        ListHex = Tree2List(pV->Hexahedra);
        List_Sort(ListHex, _p3d_cmp_entities);
        _p3d_print_hex(ListHex, Nu, Nv, pS->Nv);
        List_Delete(ListHex);
      }
    
      List_Delete(ListSurfaces);
      List_Delete(ListVolumes);
      List_Delete(ListStructSurf);
      List_Delete(ListStructVol);
    }
    
    void _p3d_print_elements(Mesh *M)
    {
      // FIXME: to do
    }
    
    void Print_Mesh_P3D(Mesh *M, FILE *fp)
    {
      P3DFILE = fp;
    
      if(!List_Nbr(M->PhysicalGroups) || CTX.mesh.save_all) {
        Msg(INFO, "Saving all elements (discarding physical groups)");
        _p3d_print_all_elements(M);
      }
      else{
        Msg(WARNING, "No Plot3d format of physicals yet, saving all elements!");
        //_p3d_print_elements(M);
        _p3d_print_all_elements(M);
      }
    }
    
    // Public Print_Mesh routine
    
    void GetDefaultMeshFileName(Mesh *M, int Type, char *name)
    {
      char ext[10] = "";
      strcpy(name, M->name);
      switch(Type){
      case FORMAT_MSH:  strcpy(ext, ".msh"); break;
      case FORMAT_VRML: strcpy(ext, ".wrl"); break;
      case FORMAT_UNV:  strcpy(ext, ".unv"); break;
      case FORMAT_GREF: strcpy(ext, ".Gref"); break;
      case FORMAT_DMG:  strcpy(ext, ".dmg"); break;
      case FORMAT_STL:  strcpy(ext, ".stl"); break;
      case FORMAT_P3D:  strcpy(ext, ".p3d"); break;
      default: Msg(GERROR, "Unknown mesh file format %d", Type); break;
      }
      strcat(name, ext);
    }
    
    void Print_Mesh(Mesh *M, char *filename, int Type)
    {
      char name[256];
    
      if(CTX.threads_lock) {
        Msg(INFO, "I'm busy! Ask me that later...");
        return;
      }
    
      CTX.threads_lock = 1;
    
      if(!filename)
        GetDefaultMeshFileName(M, Type, name);
      else
        strcpy(name, filename);
    
      Msg(INFO, "Writing mesh file '%s'", name);
    
      FILE *fp = fopen(name, "w");
      if(!fp) {
        Msg(GERROR, "Unable to open file '%s'", name);
        CTX.threads_lock = 0;
        return;
      }
    
      switch(Type){
      case FORMAT_MSH:  Print_Mesh_MSH(M, fp); break;
      case FORMAT_VRML: Print_Mesh_WRL(M, fp); break;
      case FORMAT_UNV:  Print_Mesh_UNV(M, fp); break;
      case FORMAT_GREF: Print_Mesh_GREF(M, fp); break;
      case FORMAT_DMG:  Print_Mesh_DMG(M, fp); break;
      case FORMAT_STL:  Print_Mesh_STL(M, fp); break;
      case FORMAT_P3D:  Print_Mesh_P3D(M, fp); break;
      default:
        break;
      }
    
      fclose(fp);
      Msg(INFO, "Wrote mesh file '%s'", name);
      Msg(STATUS2N, "Wrote '%s'", name);
    
      CTX.threads_lock = 0;
    }