Skip to content
Snippets Groups Projects
Select Git revision
  • d8b14b3a0c8b112933980fa4a65b53835279bcee
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

GModelIO.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    GModelIO.cpp 31.05 KiB
    #include <map>
    #include <string>
    
    #include "Message.h"
    #include "GmshDefines.h"
    #include "gmshRegion.h"
    #include "gmshFace.h"
    #include "gmshEdge.h"
    #include "MElement.h"
    #include "SBoundingBox3d.h"
    
    static int getNumVerticesForElementTypeMSH(int type)
    {
      switch (type) {
      case PNT : return 1;
      case LGN1: return 2;
      case LGN2: return 2 + 1;
      case TRI1: return 3;
      case TRI2: return 3 + 3;
      case QUA1: return 4;
      case QUA2: return 4 + 4 + 1;
      case TET1: return 4;
      case TET2: return 4 + 6;
      case HEX1: return 8;
      case HEX2: return 8 + 12 + 6 + 1;
      case PRI1: return 6;
      case PRI2: return 6 + 9 + 3;
      case PYR1: return 5;
      case PYR2: return 5 + 8 + 1;
      default: return 0;
      }
    }
    
    template<class T>
    void copyElements(std::vector<T*> &dst, const std::vector<MElement*> &src)
    {
      dst.resize(src.size());
      for(unsigned int i = 0; i < src.size(); i++) dst[i] = (T*)src[i];
    }
    
    static void storeElementsInEntities(GModel *m, int type, 
    				    std::map<int, std::vector<MElement*> > &map)
    {
      std::map<int, std::vector<MElement*> >::const_iterator it = map.begin();
      std::map<int, std::vector<MElement*> >::const_iterator ite = map.end();
      for(; it != ite; ++it){
        switch(type){
        case LGN1:     
          {
    	GEdge *e = m->edgeByTag(it->first);
    	if(!e){
    	  e = new gmshEdge(m, it->first);
    	  m->add(e);
    	}
    	copyElements(e->lines, it->second);
          }
          break;
        case TRI1: case QUA1: 
          {
    	GFace *f = m->faceByTag(it->first);
    	if(!f){
    	  f = new gmshFace(m, it->first);
    	  m->add(f);
    	}
    	if(type == TRI1) copyElements(f->triangles, it->second);
    	else if(type == QUA1) copyElements(f->quadrangles, it->second);
          }
          break;
        case TET1: case HEX1: case PRI1: case PYR1:
          {
    	GRegion *r = m->regionByTag(it->first);
    	if(!r){
    	  r = new gmshRegion(m, it->first);
    	  m->add(r);
    	}
    	if(type == TET1) copyElements(r->tetrahedra, it->second);
    	else if(type == HEX1) copyElements(r->hexahedra, it->second);
    	else if(type == PRI1) copyElements(r->prisms, it->second);
    	else if(type == PYR1) copyElements(r->pyramids, it->second);
          }
          break;
        }
      }
    }
    
    template<class T>
    static void associateEntityWithVertices(GEntity *ge, std::vector<T*> &elements)
    {
      for(unsigned int i = 0; i < elements.size(); i++)
        for(int j = 0; j < elements[i]->getNumVertices(); j++)
          elements[i]->getVertex(j)->setEntity(ge);
    }
    
    static void associateEntityWithVertices(GModel *m)
    {
      // loop on regions, then on faces, edges and vertices and store the
      // entity pointer in the the elements' vertices (this way we
      // associate the entity of lowest geometrical degree with each
      // vertex)
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
        associateEntityWithVertices(*it, (*it)->tetrahedra);
        associateEntityWithVertices(*it, (*it)->hexahedra);
        associateEntityWithVertices(*it, (*it)->prisms);
        associateEntityWithVertices(*it, (*it)->pyramids);
      }
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
        associateEntityWithVertices(*it, (*it)->triangles);
        associateEntityWithVertices(*it, (*it)->quadrangles);
      }
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
        associateEntityWithVertices(*it, (*it)->lines);
      }
      for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){
        (*it)->mesh_vertices[0]->setEntity(*it);
      }
    }
    
    static void storeVerticesInEntities(std::map<int, MVertex*> &vertices)
    {
      std::map<int, MVertex*>::const_iterator it = vertices.begin();
      std::map<int, MVertex*>::const_iterator ite = vertices.end();
      for(; it != ite; ++it){
        MVertex *v = it->second;
        GEntity *ge = v->onWhat();
        if(ge) 
          ge->mesh_vertices.push_back(v);
        else
          delete v; // we delete all unused vertices
      }
    }
    
    static void storePhysicalTagsInEntities(GModel *m, int dim,
    					std::map<int, std::map<int, std::string> > &map)
    {
      std::map<int, std::map<int, std::string> >::const_iterator it = map.begin();
      std::map<int, std::map<int, std::string> >::const_iterator ite = map.end();
      for(; it != ite; ++it){
        GEntity *ge = 0;
        switch(dim){
        case 0: ge = m->vertexByTag(it->first); break;
        case 1: ge = m->edgeByTag(it->first); break;
        case 2: ge = m->faceByTag(it->first); break;
        case 3: ge = m->regionByTag(it->first); break;
        }
        if(ge){
          std::map<int, std::string>::const_iterator it2 = it->second.begin();
          std::map<int, std::string>::const_iterator ite2 = it->second.end();
          for(; it2 != ite2; ++it2)
    	ge->physicals.push_back(it2->first);
        }
      }
    }
    
    int GModel::readMSH(const std::string &name)
    {
      FILE *fp = fopen(name.c_str(), "r");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      int elementTypes[7] = {LGN1, TRI1, QUA1, TET1, HEX1, PRI1, PYR1};
      double version = 1.0;
      char str[256];
      SBoundingBox3d bbox;
      std::map<int, MVertex*> vertices;
      std::map<int, std::vector<MVertex*> > points;
      std::map<int, std::vector<MElement*> > elements[7];
      std::map<int, std::map<int, std::string> > physicals[4];
    
      while(1) {
    
        do {
          if(!fgets(str, sizeof(str), fp) || feof(fp))
            break;
        } while(str[0] != '$');
    
        if(feof(fp))
          break;
    
        if(!strncmp(&str[1], "MeshFormat", 10)) {
    
          int format, size;
          fscanf(fp, "%lf %d %d\n", &version, &format, &size);
    
        }
        else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5)) {
    
          int numVertices;
          fscanf(fp, "%d", &numVertices);
          Msg(INFO, "%d vertices", numVertices);
    
          int progress = (numVertices > 100000) ? numVertices / 50 : numVertices / 10;
          for(int i = 0; i < numVertices; i++) {
    	int num;
    	double x, y, z;
            fscanf(fp, "%d %lf %lf %lf", &num, &x, &y, &z);
    	bbox += SPoint3(x, y, z);
    	if(vertices.count(num))
    	  Msg(WARNING, "Skipping duplicate vertex %d", num);
    	else
    	  vertices[num] = new MVertex(x, y, z);
    	if(progress && (i % progress == progress - 1))
    	  Msg(PROGRESS, "Read %d vertices", i + 1);
          }
          Msg(PROGRESS, "");
    
        }
        else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) {
    
          int numElements;
          fscanf(fp, "%d", &numElements);
          Msg(INFO, "%d elements", numElements);
    
          int progress = (numElements > 100000) ? numElements / 50 : numElements / 10;
          for(int i = 0; i < numElements; i++) {
    	int num, type, physical = 1, elementary = 1, partition = 1, numVertices;
    	if(version <= 1.0){
    	  fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices);
    	  int check = getNumVerticesForElementTypeMSH(type);
    	  if(!check){
    	    Msg(GERROR, "Unknown type for element %d", num); 
    	    continue;
    	  }
    	  if(numVertices != check){
    	    Msg(GERROR, "Wrong number of vertices (%d) for element %d", numVertices, num);
    	    continue;
    	  }
    	}
    	else{
    	  int numTags;
    	  fscanf(fp, "%d %d %d", &num, &type, &numTags);
    	  for(int j = 0; j < numTags; j++){
    	    int tag;
    	    fscanf(fp, "%d", &tag);	    
    	    if(j == 0)      physical = tag;
    	    else if(j == 1) elementary = tag;
    	    else if(j == 2) partition = tag;
    	    // ignore any other tags for now
    	  }
    	  numVertices = getNumVerticesForElementTypeMSH(type);
    	  if(!numVertices){
    	    Msg(GERROR, "Unknown type (%d) for element %d", type, num); 
    	    continue;
    	  }
    	}
    	int n[30];
            for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &n[j]);
    	int dim = 0;
            switch (type) {
            case PNT:
    	  points[elementary].push_back(vertices[n[0]]);
    	  dim = 0;
              break;
            case LGN1:
    	  elements[0][elementary].push_back
    	    (new MLine(vertices[n[0]], vertices[n[1]], num, partition));
    	  dim = 1;
              break;
            case TRI1:
    	  elements[1][elementary].push_back
    	    (new MTriangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			   num, partition));
    	  dim = 2;
              break;
            case QUA1:
    	  elements[2][elementary].push_back
    	    (new MQuadrangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			     vertices[n[3]], num, partition));
    	  dim = 2;
              break;
    	case TET1:
    	  elements[3][elementary].push_back
    	    (new MTetrahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			      vertices[n[3]], num, partition));
    	  dim = 3; 
    	  break;
    	case HEX1:
    	  elements[4][elementary].push_back
    	    (new MHexahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			     vertices[n[3]], vertices[n[4]], vertices[n[5]], 
    			     vertices[n[6]], vertices[n[7]], num, partition));
    	  dim = 3; 
    	  break;
    	case PRI1: 
    	  elements[5][elementary].push_back
    	    (new MPrism(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			vertices[n[3]], vertices[n[4]], vertices[n[5]], 
    			num, partition));
    	  dim = 3; 
    	  break;
    	case PYR1: 
    	  elements[6][elementary].push_back
    	    (new MPyramid(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
    			  vertices[n[3]], vertices[n[4]], num, partition));
    	  dim = 3; 
    	  break;
            default:
    	  Msg(GERROR, "Unknown type (%d) for element %d", type, num); 
              break;
            }
    
    	if(!physicals[dim].count(elementary) || !physicals[dim][elementary].count(physical))
    	  physicals[dim][elementary][physical] = "unnamed";
    	
    	if(progress && (i % progress == progress - 1))
    	  Msg(PROGRESS, "Read %d elements", i + 1);
          }
          Msg(PROGRESS, "");
    
        }
    
        do {
          if(!fgets(str, sizeof(str), fp) || feof(fp))
    	Msg(GERROR, "Prematured end of mesh file");
        } while(str[0] != '$');
    
      }
    
      // store the elements in their associated elementary entity. If the
      // entity does not exist, create a new one.
      for(int i = 0; i < 7; i++)
        storeElementsInEntities(this, elementTypes[i], elements[i]);
    
      // treat points separately
      {
        std::map<int, std::vector<MVertex*> >::const_iterator it = points.begin();
        std::map<int, std::vector<MVertex*> >::const_iterator ite = points.end();
        for(; it != ite; ++it){
          GVertex *v = vertexByTag(it->first);
          if(!v){
    	v = new gmshVertex(this, it->first);
    	add(v);
          }
          v->mesh_vertices.push_back(it->second[0]);
        }
      }
    
      // associate the correct geometrical entity with each mesh vertex
      associateEntityWithVertices(this);
    
      // special case for geometry vertices: now that the correct
      // geometrical entity has been associated with the vertices, we
      // reset mesh_vertices so that it can be filled again below
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        (*it)->mesh_vertices.clear();
    
      // store the vertices in their associated geometrical entity
      storeVerticesInEntities(vertices);
    
      // store the physical tags
      for(int i = 0; i < 4; i++)  
        storePhysicalTagsInEntities(this, i, physicals[i]);
    
      boundingBox += bbox;
    
      fclose(fp);
      return 1;
    }
    
    template<class T>
    static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, int saveAll, 
    			     double version, int &num, int elementary, 
    			     std::vector<int> &physicals)
    {
      for(unsigned int i = 0; i < ele.size(); i++)
        if(saveAll)
          ele[i]->writeMSH(fp, version, ++num, elementary, elementary);
        else
          for(unsigned int j = 0; j < physicals.size(); j++)
    	ele[i]->writeMSH(fp, version, ++num, elementary, physicals[j]);
    }
    
    int GModel::writeMSH(const std::string &name, double version, bool saveAll, 
    		     double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      // if there are no physicals we save all the elements
      if(noPhysicals()) saveAll = true;
    
      // get the number of vertices and renumber the vertices in a
      // continuous sequence
      int numVertices = renumberMeshVertices();
      
      // get the number of elements
      int numElements = 0;
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
          (*it)->mesh_vertices.size();
      for(eiter it = firstEdge(); it != lastEdge(); ++it)
        numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
          (*it)->lines.size();
      for(fiter it = firstFace(); it != lastFace(); ++it)
        numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
          ((*it)->triangles.size() + (*it)->quadrangles.size());
      for(riter it = firstRegion(); it != lastRegion(); ++it)
        numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
          ((*it)->tetrahedra.size() + (*it)->hexahedra.size() +
           (*it)->prisms.size() + (*it)->pyramids.size());
    
      if(version > 2.0){
        fprintf(fp, "$MeshFormat\n");
        fprintf(fp, "%g %d %d\n", version, 0, (int)sizeof(double));
        fprintf(fp, "$EndMeshFormat\n");
        fprintf(fp, "$Nodes\n");
      }
      else
        fprintf(fp, "$NOD\n");
     
      fprintf(fp, "%d\n", numVertices);
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
      for(eiter it = firstEdge(); it != lastEdge(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
          (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
      for(fiter it = firstFace(); it != lastFace(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
      for(riter it = firstRegion(); it != lastRegion(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
    
      if(version > 2.0){
        fprintf(fp, "$EndNodes\n");
        fprintf(fp, "$Elements\n");
      }
      else{
        fprintf(fp, "$ENDNOD\n");
        fprintf(fp, "$ELM\n");
      }
    
      fprintf(fp, "%d\n", numElements);
      int num = 0;
    
      for(viter it = firstVertex(); it != lastVertex(); ++it){
        writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
      }
      for(eiter it = firstEdge(); it != lastEdge(); ++it){
        writeElementsMSH(fp, (*it)->lines, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
      }
      for(fiter it = firstFace(); it != lastFace(); ++it){
        writeElementsMSH(fp, (*it)->triangles, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
        writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
      }
      for(riter it = firstRegion(); it != lastRegion(); ++it){
        writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
        writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
        writeElementsMSH(fp, (*it)->prisms, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
        writeElementsMSH(fp, (*it)->pyramids, saveAll, version, num,
    		     (*it)->tag(), (*it)->physicals);
      }
      
      if(version > 2.0){
        fprintf(fp, "$EndElements\n");
      }
      else{
        fprintf(fp, "$ENDELM\n");
      }
    
      return 1;
    }
    
    int GModel::writePOS(const std::string &name, double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      if(numRegion()){
        fprintf(fp, "View \"Volumes\" {\n");
        fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", "
    	    "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8));
        for(riter it = firstRegion(); it != lastRegion(); ++it) {
          for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
    	(*it)->tetrahedra[i]->writePOS(fp, scalingFactor, (*it)->tag());
          for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
    	(*it)->hexahedra[i]->writePOS(fp, scalingFactor, (*it)->tag());
          for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
    	(*it)->prisms[i]->writePOS(fp, scalingFactor, (*it)->tag());
          for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
    	(*it)->pyramids[i]->writePOS(fp, scalingFactor, (*it)->tag());
        }
        fprintf(fp, "};\n");
      }
      
      if(numFace()){
        fprintf(fp, "View \"Surfaces\" {\n");
        fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", "
    	    "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8));
        for(fiter it = firstFace(); it != lastFace(); ++it) {
          for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
    	(*it)->triangles[i]->writePOS(fp, scalingFactor, (*it)->tag());
          for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
    	(*it)->quadrangles[i]->writePOS(fp, scalingFactor, (*it)->tag());
        }
        fprintf(fp, "};\n");
      }
    
      if(numEdge()){
        fprintf(fp, "View \"Lines\" {\n");
        fprintf(fp, "T2(1.e5,30,%d){\"Elementary Entity\", \"Element Number\", "
    	    "\"Gamma\", \"Eta\", \"Rho\"};\n", (1<<16)|(4<<8));
        for(eiter it = firstEdge(); it != lastEdge(); ++it) {
          for(unsigned int i = 0; i < (*it)->lines.size(); i++)
     	(*it)->lines[i]->writePOS(fp, scalingFactor, (*it)->tag());
        }
        fprintf(fp, "};\n");
      }
    
      fclose(fp);
      return 1;
    }
    
    static void swapBytes(char *array, int size, int n)
    {
      char *x = new char[size];
      for(int i = 0; i < n; i++) {
        char *a = &array[i * size];
        memcpy(x, a, size);
        for(int c = 0; c < size; c++)
          a[size - 1 - c] = x[c];
      }
      delete [] x;
    }
    
    int GModel::readSTL(const std::string &name, double tolerance)
    {
      FILE *fp = fopen(name.c_str(), "rb");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      // read all facets and store triplets of points
      std::vector<SPoint3> points;
      SBoundingBox3d bbox;
    
      // "solid" or binary data header
      char buffer[256];
      fgets(buffer, sizeof(buffer), fp);
    
      if(!strncmp(buffer, "solid", 5)) { 
        // ASCII STL
        while(!feof(fp)) {
          // "facet normal x y z" or "endsolid"
          if(!fgets(buffer, sizeof(buffer), fp)) break;
          if(!strncmp(buffer, "endsolid", 8)) break;
          char s1[256], s2[256];
          float x, y, z;
          sscanf(buffer, "%s %s %f %f %f", s1, s2, &x, &y, &z);
          // "outer loop"
          if(!fgets(buffer, sizeof(buffer), fp)) break;
          // "vertex x y z"
          for(int i = 0; i < 3; i++){
    	if(!fgets(buffer, sizeof(buffer), fp)) break;
    	sscanf(buffer, "%s %f %f %f", s1, &x, &y, &z);
    	SPoint3 p(x, y, z);
    	points.push_back(p);
    	bbox += p;
          }
          // "endloop"
          if(!fgets(buffer, sizeof(buffer), fp)) break;
          // "endfacet"
          if(!fgets(buffer, sizeof(buffer), fp)) break;
        }
      }
      else{
        // Binary STL
        rewind(fp);
        char header[80];
        if(fread(header, sizeof(char), 80, fp)){
          unsigned int nfacets = 0;
          size_t ret = fread(&nfacets, sizeof(unsigned int), 1, fp);
          bool swap = false;
          if(nfacets > 10000000){
    	Msg(INFO, "Swapping bytes from binary file");
    	swap = true;
    	swapBytes((char*)&nfacets, sizeof(unsigned int), 1);
          }
          if(ret && nfacets){
    	char *data = new char[nfacets * 50 * sizeof(char)];
    	ret = fread(data, sizeof(char), nfacets * 50, fp);
    	if(ret == nfacets * 50){
    	  for(unsigned int i = 0; i < nfacets; i++) {
    	    float *xyz = (float *)&data[i * 50 * sizeof(char)];
    	    if(swap) swapBytes((char*)xyz, sizeof(float), 12);
    	    for(int j = 0; j < 3; j++){
    	      SPoint3 p(xyz[3 + 3 * j], xyz[3 + 3 * j + 1], xyz[3 + 3 * j + 2]);
    	      points.push_back(p);
    	      bbox += p;
    	    }
    	  }
    	}
    	delete data;
          }
        }
      }
    
      if(!points.size()){
        Msg(GERROR, "No facets found in STL file");
        return 0;
      }
      
      if(points.size() % 3){
        Msg(GERROR, "Wrong number of points in STL file");
        return 0;
      }
    
      Msg(INFO, "%d facets", points.size() / 3);
    
      // create face
      GFace *face = new gmshFace(this, numFace() + 1);
      add(face);
    
      // create (unique) vertices and triangles
      SPoint3 min = bbox.min();
      SPoint3 max = bbox.max();
      double dx = max.x() - min.x(), dy = max.y() - min.y(), dz = max.z() - min.z();
      double lc = sqrt(dx * dx + dy * dy + dz * dz);
      MVertexLessThanLexicographic::tolerance = lc * tolerance;
      std::set<MVertex*, MVertexLessThanLexicographic> vertices;
      for(unsigned int i = 0; i < points.size(); i += 3){
        MVertex *v[3];
        for(int j = 0; j < 3; j++){
          double x = points[i + j].x(), y = points[i + j].y(), z = points[i + j].z();
          MVertex w(x, y, z);
          std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = vertices.find(&w);
          if(it != vertices.end()) {
            v[j] = *it;
          }
          else {
            v[j] = new MVertex(x, y, z);
    	vertices.insert(v[j]);
    	face->mesh_vertices.push_back(v[j]);
          }
        }
        face->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
      }
    
      boundingBox += bbox;
    
      fclose(fp);
      return 1;
    }
    
    int GModel::writeSTL(const std::string &name, double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      fprintf(fp, "solid Created by Gmsh\n");
      for(fiter it = firstFace(); it != lastFace(); ++it) {
        for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
          (*it)->triangles[i]->writeSTL(fp, scalingFactor);
        for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
          (*it)->quadrangles[i]->writeSTL(fp, scalingFactor);
      }
      fprintf(fp, "endsolid Created by Gmsh\n");
    
      fclose(fp);
      return 1;
    }
    
    int GModel::writeVRML(const std::string &name, double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      renumberMeshVertices();
    
      fprintf(fp, "#VRML V1.0 ascii\n");
      fprintf(fp, "#created by Gmsh\n");
      fprintf(fp, "Coordinate3 {\n");
      fprintf(fp, "  point [\n");
    
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
      for(eiter it = firstEdge(); it != lastEdge(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
          (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
      for(fiter it = firstFace(); it != lastFace(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
    
      fprintf(fp, "  ]\n");
      fprintf(fp, "}\n");
    
      for(eiter it = firstEdge(); it != lastEdge(); ++it){
        fprintf(fp, "DEF Curve%d IndexedLineSet {\n", (*it)->tag());
        fprintf(fp, "  coordIndex [\n");
        for(unsigned int i = 0; i < (*it)->lines.size(); i++)
          (*it)->lines[i]->writeVRML(fp);
        fprintf(fp, "  ]\n");
        fprintf(fp, "}\n");
      }
    
      for(fiter it = firstFace(); it != lastFace(); ++it){
        fprintf(fp, "DEF Surface%d IndexedFaceSet {\n", (*it)->tag());
        fprintf(fp, "  coordIndex [\n");
        for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
          (*it)->triangles[i]->writeVRML(fp);
        for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
          (*it)->quadrangles[i]->writeVRML(fp);
        fprintf(fp, "  ]\n");
        fprintf(fp, "}\n");
      }
      
      fclose(fp);
      return 1;
    }
    
    static void addInGroupOfNodesUNV(FILE *fp, std::set<int> &nodes, int num)
    {
      std::set<int>::iterator it = nodes.find(num);
      if(it == nodes.end()){
        nodes.insert(num);
        fprintf(fp, "%10d%10d%2d%2d%2d%2d%2d%2d\n", num, 1, 0, 1, 0, 0, 0, 0);
        fprintf(fp, "   0.0000000000000000D+00   1.0000000000000000D+00"
    	    "   0.0000000000000000D+00\n");
        fprintf(fp, "   0.0000000000000000D+00   0.0000000000000000D+00"
    	    "   0.0000000000000000D+00\n");
        fprintf(fp, "%10d%10d%10d%10d%10d%10d\n", 0, 0, 0, 0, 0, 0);
      }
    }
    
    template<class T>
    static void addInGroupOfNodesUNV(FILE *fp, std::set<int> &nodes, 
    				 std::vector<T*> &elements)
    {
      for(unsigned int i = 0; i < elements.size(); i++)
        for(int j = 0; j < elements[i]->getNumVertices(); j++)
          addInGroupOfNodesUNV(fp, nodes, elements[i]->getVertex(j)->getNum());
    }
    
    int GModel::writeUNV(const std::string &name, double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      // IDEAS records
      const int NODES=2411, ELEMENTS=2412, GROUPOFNODES=790;
    
      // IDEAS elements
      const int BEAM=21, THINSHLL=91, QUAD=94, SOLIDFEM=111, WEDGE=112, BRICK=115;
      //const int BEAM2=24, THINSHLL2=92, QUAD2=95/*?*/, SOLIDFEM2=118;
    
      renumberMeshVertices();
      
      fprintf(fp, "%6d\n", -1);
      fprintf(fp, "%6d\n", NODES);
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
      for(eiter it = firstEdge(); it != lastEdge(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
          (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
      for(fiter it = firstFace(); it != lastFace(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
      for(riter it = firstRegion(); it != lastRegion(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeUNV(fp, scalingFactor);
      fprintf(fp, "%6d\n", -1);  
    
      fprintf(fp, "%6d\n", -1);
      fprintf(fp, "%6d\n", ELEMENTS);
      for(eiter it = firstEdge(); it != lastEdge(); ++it){
        for(unsigned int i = 0; i < (*it)->lines.size(); i++)
          (*it)->lines[i]->writeUNV(fp, BEAM, (*it)->tag());
      }
      for(fiter it = firstFace(); it != lastFace(); ++it){
        for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
          (*it)->triangles[i]->writeUNV(fp, THINSHLL, (*it)->tag());
        for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
          (*it)->quadrangles[i]->writeUNV(fp, QUAD, (*it)->tag());
      }
      for(riter it = firstRegion(); it != lastRegion(); ++it){
        for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
          (*it)->tetrahedra[i]->writeUNV(fp, SOLIDFEM, (*it)->tag());
        for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
          (*it)->hexahedra[i]->writeUNV(fp, BRICK, (*it)->tag());
        for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
          (*it)->prisms[i]->writeUNV(fp, WEDGE, (*it)->tag());
      }
      fprintf(fp, "%6d\n", -1);
    
      std::map<int, std::vector<GEntity*> > physicals[4];
      getPhysicalGroups(physicals);
    
      for(int dim = 0; dim < 4; dim++){
        std::map<int, std::vector<GEntity*> >::const_iterator it = physicals[dim].begin();
        std::map<int, std::vector<GEntity*> >::const_iterator ite = physicals[dim].end();
        for(; it != ite; ++it){
          fprintf(fp, "%6d\n", -1);
          fprintf(fp, "%6d\n", GROUPOFNODES);
          fprintf(fp, "%10d%10d\n", it->first, 1);
          fprintf(fp, "LOAD SET %2d\n", 1);
          std::set<int> nodes;
          for(unsigned int i = 0; i < it->second.size(); i++){
    	// we could also do this using the mesh_vertices of the entity
    	// and all the entities in the closure
    	GVertex *v;
    	switch(dim){
    	case 0: 
    	  v = (GVertex*)it->second[i];
    	  for(unsigned int j = 0; j < v->mesh_vertices.size(); j++)
    	    addInGroupOfNodesUNV(fp, nodes, v->mesh_vertices[j]->getNum());
    	  break;
    	case 1: 
    	  addInGroupOfNodesUNV(fp, nodes, ((GEdge*)it->second[i])->lines);
    	  break;
    	case 2: 
    	  addInGroupOfNodesUNV(fp, nodes, ((GFace*)it->second[i])->triangles);
    	  addInGroupOfNodesUNV(fp, nodes, ((GFace*)it->second[i])->quadrangles);
    	  break;
    	case 3: 
    	  addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->tetrahedra);
    	  addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->hexahedra);
    	  addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->prisms);
    	  addInGroupOfNodesUNV(fp, nodes, ((GRegion*)it->second[i])->pyramids);
    	  break;
    	}
          }
          fprintf(fp, "%6d\n", -1);
        }
      }
    
      fclose(fp);
      return 1;
    }
    
    int GModel::readMESH(const std::string &name)
    {
      FILE *fp = fopen(name.c_str(), "r");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      char buffer[256];
      fgets(buffer, sizeof(buffer), fp);
    
      char str[256];
      int format;
      sscanf(buffer, "%s %d", str, &format);
    
      if(format != 1){
        Msg(GERROR, "Non-ASCII INRIA mesh import is not yet available");
        return 0;
      }
    
      std::map<int, MVertex*> vertices;
      int elementTypes[2] = {TRI1, QUA1};
      std::map<int, std::vector<MElement*> > elements[2];
      SBoundingBox3d bbox;
    
      while(!feof(fp)) {
        fgets(buffer, sizeof(buffer), fp);
        if(buffer[0] != '#'){ // skip comments
          sscanf(buffer, "%s", str);
          if(!strcmp(str, "Dimension")){
    	fgets(buffer, sizeof(buffer), fp);
          }
          else if(!strcmp(str, "Vertices")){
    	fgets(buffer, sizeof(buffer), fp);
    	int nbv;
    	sscanf(buffer, "%d", &nbv);
    	Msg(INFO, "%d vertices", nbv);
    	for(int i = 0; i < nbv; i++) {
    	  fgets(buffer, sizeof(buffer), fp);
    	  int cl;
    	  double x, y, z;
    	  sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &cl);
    	  bbox += SPoint3(x, y, z);
    	  vertices[i + 1] = new MVertex(x, y, z);
    	}
          }
          else if(!strcmp(str, "Triangles")){
    	fgets(buffer, sizeof(buffer), fp);
    	int nbt;
    	sscanf(buffer, "%d", &nbt);
    	Msg(INFO, "%d triangles", nbt);
    	for(int i = 0; i < nbt; i++) {
    	  fgets(buffer, sizeof(buffer), fp);
    	  int n1, n2, n3, cl;
    	  sscanf(buffer, "%d %d %d %d", &n1, &n2, &n3, &cl);
    	  elements[0][cl].push_back
    	    (new MTriangle(vertices[n1], vertices[n2], vertices[n3]));
    	}
          }
          else if(!strcmp(str, "Quadrilaterals")) {
    	fgets(buffer, sizeof(buffer), fp);
    	int nbq;
    	sscanf(buffer, "%d", &nbq);
    	Msg(INFO, "%d quadrangles", nbq);
    	for(int i = 0; i < nbq; i++) {
    	  fgets(buffer, sizeof(buffer), fp);
    	  int n1, n2, n3, n4, cl;
    	  sscanf(buffer, "%d %d %d %d %d", &n1, &n2, &n3, &n4, &cl);
    	  elements[1][cl].push_back
    	    (new MQuadrangle(vertices[n1], vertices[n2], vertices[n3], vertices[n4]));
    	}
          }
        }
      }
    
      // store the elements in their associated elementary entity. If the
      // entity does not exist, create a new one.
      for(int i = 0; i < 2; i++)
        storeElementsInEntities(this, elementTypes[i], elements[i]);
    
      // associate the correct geometrical entity with each mesh vertex
      associateEntityWithVertices(this);
    
      // store the vertices in their associated geometrical entity
      storeVerticesInEntities(vertices);
    
      boundingBox += bbox;
    
      fclose(fp);
      return 1;
    }
    
    int GModel::writeMESH(const std::string &name, double scalingFactor)
    {
      FILE *fp = fopen(name.c_str(), "w");
      if(!fp){
        Msg(GERROR, "Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      fprintf(fp, " MeshVersionFormatted 1\n");
      fprintf(fp, " Dimension\n");
      fprintf(fp, " 3\n");
    
      int numVertices = renumberMeshVertices();
      fprintf(fp, " Vertices\n");
      fprintf(fp, " %d\n", numVertices);
      for(viter it = firstVertex(); it != lastVertex(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
      for(eiter it = firstEdge(); it != lastEdge(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
          (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
      for(fiter it = firstFace(); it != lastFace(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
      for(riter it = firstRegion(); it != lastRegion(); ++it)
        for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
          (*it)->mesh_vertices[i]->writeMESH(fp, scalingFactor);
      
      int numTriangles = 0, numQuadrangles = 0;
      for(fiter it = firstFace(); it != lastFace(); ++it){
        numTriangles += (*it)->triangles.size();
        numQuadrangles += (*it)->quadrangles.size();
      }
      if(numTriangles){
        fprintf(fp, " Triangles\n");
        fprintf(fp, " %d\n", numTriangles);
        for(fiter it = firstFace(); it != lastFace(); ++it)
          for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
    	(*it)->triangles[i]->writeMESH(fp, (*it)->tag());
      }
      if(numQuadrangles){
        fprintf(fp, " Quadrilaterals\n");
        fprintf(fp, " %d\n", numQuadrangles);
        for(fiter it = firstFace(); it != lastFace(); ++it)
          for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
    	(*it)->quadrangles[i]->writeMESH(fp, (*it)->tag());
      }
    
      fprintf(fp, " End\n");
      
      fclose(fp);
      return 1;
    }