Skip to content
Snippets Groups Projects
Select Git revision
  • 74600b626df8c41ff216ffda666f083054469f52
  • 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

TriNedelecBasis.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    GModelIO_MSH2.cpp 38.68 KiB
    // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    
    #include <stdlib.h>
    #include <string.h>
    #include <map>
    #include <string>
    #include <sstream>
    #include <cassert>
    #include <iomanip>
    #include "GModel.h"
    #include "OS.h"
    #include "GmshDefines.h"
    #include "MPoint.h"
    #include "MLine.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MTetrahedron.h"
    #include "MHexahedron.h"
    #include "MPrism.h"
    #include "MPyramid.h"
    #include "MElementCut.h"
    #include "StringUtils.h"
    #include "GmshMessage.h"
    #include "Context.h"
    #include "OS.h"
    
    // periodic nodes and entities backported from MSH3 format
    extern void writeMSHPeriodicNodes(FILE *fp, std::vector<GEntity *> &entities,
                                      bool renumber, bool saveAll);
    extern void readMSHPeriodicNodes(FILE *fp, GModel *gm);
    extern void writeMSHEntities(FILE *fp, GModel *gm);
    
    static bool getMeshVertices(int num, int *indices,
                                std::map<int, MVertex *> &map,
                                std::vector<MVertex *> &vertices)
    {
      for(int i = 0; i < num; i++) {
        if(!map.count(indices[i])) {
          Msg::Error("Wrong vertex index %d", indices[i]);
          return false;
        }
        else
          vertices.push_back(map[indices[i]]);
      }
      return true;
    }
    
    static bool getMeshVertices(int num, int *indices, std::vector<MVertex *> &vec,
                                std::vector<MVertex *> &vertices, int minVertex = 0)
    {
      for(int i = 0; i < num; i++) {
        if(indices[i] < minVertex ||
           indices[i] > (int)(vec.size() - 1 + minVertex)) {
          Msg::Error("Wrong vertex index %d", indices[i]);
          return false;
        }
        else
          vertices.push_back(vec[indices[i]]);
      }
      return true;
    }
    
    static MElement *
    createElementMSH2(GModel *m, int num, int typeMSH, int physical, int reg,
                      unsigned int part, std::vector<MVertex *> &v,
                      std::map<int, std::vector<MElement *> > elements[10],
                      std::map<int, std::map<int, std::string> > physicals[4],
                      bool owner = false, MElement *parent = 0, MElement *d1 = 0,
                      MElement *d2 = 0)
    {
      if(CTX::instance()->mesh.switchElementTags) {
        int tmp = reg;
        reg = physical;
        physical = tmp;
      }
    
      MElementFactory factory;
      MElement *e = factory.create(typeMSH, v, num, part, owner, 0, parent, d1, d2);
    
      if(!e) {
        Msg::Error("Unknown type of element %d", typeMSH);
        return NULL;
      }
    
      int dim = e->getDim();
      if(physical &&
         (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
        physicals[dim][reg][physical] = "unnamed";
    
      if(part > m->getNumPartitions()) m->setNumPartitions(part);
      return e;
    }
    
    int GModel::_readMSH2(const std::string &name)
    {
      FILE *fp = Fopen(name.c_str(), "rb");
      if(!fp) {
        Msg::Error("Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      char str[256] = "XXX";
      double version = 1.0;
      bool binary = false, swap = false, postpro = false;
      std::map<int, std::vector<MElement *> > elements[10];
      std::map<int, std::map<int, std::string> > physicals[4];
      std::map<int, MVertex *> vertexMap;
      std::vector<MVertex *> vertexVector;
      int minVertex = 0;
    
      while(1) {
        while(str[0] != '$') {
          if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
        }
    
        if(feof(fp)) break;
    
        if(!strncmp(&str[1], "MeshFormat", 10)) {
          if(!fgets(str, sizeof(str), fp)) {
            fclose(fp);
            return 0;
          }
          int format, size;
          if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) {
            fclose(fp);
            return 0;
          }
          if(format) {
            binary = true;
            Msg::Debug("Mesh is in binary format");
            int one;
            if(fread(&one, sizeof(int), 1, fp) != 1) {
              fclose(fp);
              return 0;
            }
            if(one != 1) {
              swap = true;
              Msg::Debug("Swapping bytes from binary file");
            }
          }
        }
        else if(!strncmp(&str[1], "PhysicalNames", 13)) {
          if(!fgets(str, sizeof(str), fp)) {
            fclose(fp);
            return 0;
          }
          int numNames;
          if(sscanf(str, "%d", &numNames) != 1) {
            fclose(fp);
            return 0;
          }
          for(int i = 0; i < numNames; i++) {
            int dim = -1, num;
            if(version > 2.0) {
              if(fscanf(fp, "%d", &dim) != 1) {
                fclose(fp);
                return 0;
              }
            }
            if(fscanf(fp, "%d", &num) != 1) {
              fclose(fp);
              return 0;
            }
            if(!fgets(str, sizeof(str), fp)) {
              fclose(fp);
              return 0;
            }
            std::string name = ExtractDoubleQuotedString(str, 256);
            if(name.size()) setPhysicalName(name, dim, num);
          }
        }
        else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5) ||
                !strncmp(&str[1], "ParametricNodes", 15)) {
          const bool parametric = !strncmp(&str[1], "ParametricNodes", 15);
          if(!fgets(str, sizeof(str), fp)) {
            fclose(fp);
            return 0;
          }
          int numVertices = -1;
          if(sscanf(str, "%d", &numVertices) != 1) {
            fclose(fp);
            return 0;
          }
          Msg::Info("%d vertices", numVertices);
          Msg::ResetProgressMeter();
          vertexVector.clear();
          vertexMap.clear();
          minVertex = numVertices + 1;
          int maxVertex = -1;
          for(int i = 0; i < numVertices; i++) {
            int num;
            double xyz[3], uv[2];
            MVertex *newVertex = 0;
            if(!parametric) {
              if(!binary) {
                if(fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) !=
                   4) {
                  fclose(fp);
                  return 0;
                }
              }
              else {
                if(fread(&num, sizeof(int), 1, fp) != 1) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)&num, sizeof(int), 1);
                if(fread(xyz, sizeof(double), 3, fp) != 3) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
              }
              newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
            }
            else {
              int iClasDim, iClasTag;
              if(!binary) {
                if(fscanf(fp, "%d %lf %lf %lf %d %d", &num, &xyz[0], &xyz[1],
                          &xyz[2], &iClasDim, &iClasTag) != 6) {
                  fclose(fp);
                  return 0;
                }
              }
              else {
                if(fread(&num, sizeof(int), 1, fp) != 1) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)&num, sizeof(int), 1);
                if(fread(xyz, sizeof(double), 3, fp) != 3) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
                if(fread(&iClasDim, sizeof(int), 1, fp) != 1) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)&iClasDim, sizeof(int), 1);
                if(fread(&iClasTag, sizeof(int), 1, fp) != 1) {
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)&iClasTag, sizeof(int), 1);
              }
              if(iClasDim == 0) {
                GVertex *gv = getVertexByTag(iClasTag);
                if(gv) gv->mesh_vertices.clear();
                newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
              }
              else if(iClasDim == 1) {
                GEdge *ge = getEdgeByTag(iClasTag);
                if(!binary) {
                  if(fscanf(fp, "%lf", &uv[0]) != 1) {
                    fclose(fp);
                    return 0;
                  }
                }
                else {
                  if(fread(uv, sizeof(double), 1, fp) != 1) {
                    fclose(fp);
                    return 0;
                  }
                  if(swap) SwapBytes((char *)uv, sizeof(double), 1);
                }
                newVertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, uv[0], num);
              }
              else if(iClasDim == 2) {
                GFace *gf = getFaceByTag(iClasTag);
                if(!binary) {
                  if(fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) {
                    fclose(fp);
                    return 0;
                  }
                }
                else {
                  if(fread(uv, sizeof(double), 2, fp) != 2) {
                    fclose(fp);
                    return 0;
                  }
                  if(swap) SwapBytes((char *)uv, sizeof(double), 2);
                }
                newVertex =
                  new MFaceVertex(xyz[0], xyz[1], xyz[2], gf, uv[0], uv[1], num);
              }
              else if(iClasDim == 3) {
                GRegion *gr = getRegionByTag(iClasTag);
                newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
              }
            }
            minVertex = std::min(minVertex, num);
            maxVertex = std::max(maxVertex, num);
            if(vertexMap.count(num))
              Msg::Warning("Skipping duplicate vertex %d", num);
            vertexMap[num] = newVertex;
            if(numVertices > 100000)
              Msg::ProgressMeter(i + 1, numVertices, true, "Reading nodes");
          }
          // If the vertex numbering is dense, transfer the map into a
          // vector to speed up element creation
          if((int)vertexMap.size() == numVertices &&
             ((minVertex == 1 && maxVertex == numVertices) ||
              (minVertex == 0 && maxVertex == numVertices - 1))) {
            Msg::Debug("Vertex numbering is dense");
            vertexVector.resize(vertexMap.size() + 1);
            if(minVertex == 1)
              vertexVector[0] = 0;
            else
              vertexVector[numVertices] = 0;
            std::map<int, MVertex *>::const_iterator it = vertexMap.begin();
            for(; it != vertexMap.end(); ++it) vertexVector[it->first] = it->second;
            vertexMap.clear();
          }
        }
        else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) {
          if(!fgets(str, sizeof(str), fp)) {
            fclose(fp);
            return 0;
          }
          int numElements;
    
          std::map<int, MElement *> elems;
          std::map<int, MElement *> parents;
          std::map<int, MElement *> domains;
          std::map<int, int> elemreg;
          std::map<int, int> elemphy;
    
          std::set<MElement *> parentsOwned;
          sscanf(str, "%d", &numElements);
          Msg::Info("%d elements", numElements);
          Msg::ResetProgressMeter();
          if(!binary) {
            for(int i = 0; i < numElements; i++) {
              int num, type, physical = 0, elementary = 0, partition = 0,
                             parent = 0;
              int dom1 = 0, dom2 = 0, numVertices;
              std::vector<short> ghosts;
              if(version <= 1.0) {
                if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary,
                          &numVertices) != 5) {
                  fclose(fp);
                  return 0;
                }
                if(numVertices != (int)MElement::getInfoMSH(type)) {
                  fclose(fp);
                  return 0;
                }
              }
              else {
                int numTags;
                if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) {
                  fclose(fp);
                  return 0;
                }
                int numPartitions = 0;
                for(int j = 0; j < numTags; j++) {
                  int tag;
                  if(fscanf(fp, "%d", &tag) != 1) {
                    fclose(fp);
                    return 0;
                  }
                  if(j == 0)
                    physical = tag;
                  else if(j == 1)
                    elementary = tag;
                  else if(version < 2.2 && j == 2)
                    partition = tag;
                  else if(version >= 2.2 && j == 2 && numTags > 3)
                    numPartitions = tag;
                  else if(version >= 2.2 && j == 3)
                    partition = tag;
                  else if(j >= 4 && j < 4 + numPartitions - 1)
                    ghosts.push_back(-tag);
                  else if(j == 3 + numPartitions && (numTags == 4 + numPartitions))
                    parent = tag;
                  else if(j == 3 + numPartitions &&
                          (numTags == 5 + numPartitions)) {
                    dom1 = tag;
                    j++;
                    if(fscanf(fp, "%d", &dom2) != 1) {
                      fclose(fp);
                      return 0;
                    }
                  }
                }
                if(!(numVertices = MElement::getInfoMSH(type))) {
                  if(type != MSH_POLYG_ && type != MSH_POLYH_ &&
                     type != MSH_POLYG_B) {
                    fclose(fp);
                    return 0;
                  }
                  if(fscanf(fp, "%d", &numVertices) != 1) {
                    fclose(fp);
                    return 0;
                  }
                }
              }
              int *indices = new int[numVertices];
              for(int j = 0; j < numVertices; j++) {
                if(fscanf(fp, "%d", &indices[j]) != 1) {
                  delete[] indices;
                  fclose(fp);
                  return 0;
                }
              }
              std::vector<MVertex *> vertices;
              if(vertexVector.size()) {
                if(!getMeshVertices(numVertices, indices, vertexVector, vertices,
                                    minVertex)) {
                  delete[] indices;
                  fclose(fp);
                  return 0;
                }
              }
              else {
                if(!getMeshVertices(numVertices, indices, vertexMap, vertices)) {
                  delete[] indices;
                  fclose(fp);
                  return 0;
                }
              }
              MElement *p = NULL;
              bool own = false;
    
              // search parent element
              if(parent != 0) {
                std::map<int, MElement *>::iterator ite = elems.find(parent);
                if(ite == elems.end())
                  Msg::Error(
                    "Parent element (ascii) %d not found for element %d of type %d",
                    parent, num, type);
                else {
                  p = ite->second;
                  parents[parent] = p;
                }
                std::set<MElement *>::iterator itpo = parentsOwned.find(p);
                if(itpo == parentsOwned.end()) {
                  own = true;
                  parentsOwned.insert(p);
                }
                assert(p != NULL);
              }
    
              // search domains
              MElement *doms[2] = {NULL, NULL};
              if(dom1) {
                std::map<int, MElement *>::iterator ite = elems.find(dom1);
                if(ite == elems.end())
                  Msg::Error("Domain element %d not found for element %d", dom1,
                             num);
                else
                  doms[0] = ite->second;
    
                ite = elems.find(dom2);
                if(ite != elems.end()) doms[1] = ite->second;
    
                if(!doms[0])
                  Msg::Error("Domain element %d not found for element %d", dom1,
                             num);
                if(dom2 && !doms[1])
                  Msg::Error("Domain element %d not found for element %d", dom2,
                             num);
              }
              delete[] indices;
    
              if(elementary < 0) continue;
              MElement *e = createElementMSH2(this, num, type, physical, elementary,
                                              partition, vertices, elements,
                                              physicals, own, p, doms[0], doms[1]);
              elems[num] = e;
              elemreg[num] = elementary;
              elemphy[num] = physical;
              for(std::size_t j = 0; j < ghosts.size(); j++)
                _ghostCells.insert(std::pair<MElement *, short>(e, ghosts[j]));
              if(numElements > 100000)
                Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
            }
          }
          else {
            int numElementsPartial = 0;
            while(numElementsPartial < numElements) {
              int header[3];
              if(fread(header, sizeof(int), 3, fp) != 3) {
                fclose(fp);
                return 0;
              }
              if(swap) SwapBytes((char *)header, sizeof(int), 3);
              int type = header[0];
              int numElms = header[1];
              int numTags = header[2];
              int numVertices = MElement::getInfoMSH(type);
              unsigned int n = 1 + numTags + numVertices;
              int *data = new int[n];
              for(int i = 0; i < numElms; i++) {
                if(fread(data, sizeof(int), n, fp) != n) {
                  delete[] data;
                  fclose(fp);
                  return 0;
                }
                if(swap) SwapBytes((char *)data, sizeof(int), n);
                int num = data[0];
                int physical = (numTags > 0) ? data[1] : 0;
                int elementary = (numTags > 1) ? data[2] : 0;
                int numPartitions = (version >= 2.2 && numTags > 3) ? data[3] : 0;
                int partition = (version < 2.2 && numTags > 2) ?
                                  data[3] :
                                  (version >= 2.2 && numTags > 3) ? data[4] : 0;
                int parent = (version < 2.2 && numTags > 3) ||
                                 (version >= 2.2 && numPartitions &&
                                  numTags > 3 + numPartitions) ||
                                 (version >= 2.2 && !numPartitions && numTags > 2) ?
                               data[numTags] :
                               0;
                int *indices = &data[numTags + 1];
                std::vector<MVertex *> vertices;
                if(vertexVector.size()) {
                  if(!getMeshVertices(numVertices, indices, vertexVector, vertices,
                                      minVertex)) {
                    delete[] data;
                    fclose(fp);
                    return 0;
                  }
                }
                else {
                  if(!getMeshVertices(numVertices, indices, vertexMap, vertices)) {
                    delete[] data;
                    fclose(fp);
                    return 0;
                  }
                }
                MElement *p = NULL;
                bool own = false;
                if(parent) {
                  std::map<int, MElement *>::iterator ite = elems.find(parent);
                  if(ite == elems.end())
                    Msg::Error(
                      "Parent (binary) element %d not found for element %d", parent,
                      num);
                  else {
                    p = ite->second;
                    parents[parent] = p;
                  }
                  std::set<MElement *>::iterator itpo = parentsOwned.find(p);
                  if(itpo == parentsOwned.end()) {
                    own = true;
                    parentsOwned.insert(p);
                  }
                  assert(p != NULL);
                }
                MElement *e = createElementMSH2(this, num, type, physical,
                                                elementary, partition, vertices,
                                                elements, physicals, own, p);
                elems[num] = e;
                elemreg[num] = elementary;
                elemphy[num] = physical;
                if(numPartitions > 1)
                  for(int j = 0; j < numPartitions - 1; j++)
                    _ghostCells.insert(
                      std::pair<MElement *, short>(e, -data[5 + j]));
                if(numElements > 100000)
                  Msg::ProgressMeter(numElementsPartial + i + 1, numElements, true,
                                     "Reading elements");
              }
              delete[] data;
              numElementsPartial += numElms;
            }
          }
          for(int i = 0; i < 10; i++) elements[i].clear();
    
          std::map<int, MElement *>::iterator ite;
          for(ite = elems.begin(); ite != elems.end(); ite++) {
            int num = ite->first;
            MElement *e = ite->second;
            if(parents.find(num) == parents.end()) {
              int reg;
              if(CTX::instance()->mesh.switchElementTags) {
                reg = elemphy[num];
              }
              else {
                reg = elemreg[num];
              }
              switch(e->getType()) {
              case TYPE_PNT: elements[0][reg].push_back(e); break;
              case TYPE_LIN: elements[1][reg].push_back(e); break;
              case TYPE_TRI: elements[2][reg].push_back(e); break;
              case TYPE_QUA: elements[3][reg].push_back(e); break;
              case TYPE_TET: elements[4][reg].push_back(e); break;
              case TYPE_HEX: elements[5][reg].push_back(e); break;
              case TYPE_PRI: elements[6][reg].push_back(e); break;
              case TYPE_PYR: elements[7][reg].push_back(e); break;
              case TYPE_POLYG: elements[8][reg].push_back(e); break;
              case TYPE_POLYH: elements[9][reg].push_back(e); break;
              default:
                Msg::Error("Wrong type of element");
                fclose(fp);
                return 0;
              }
            }
          }
        }
        else if(!strncmp(&str[1], "NodeData", 8)) {
          // there's some nodal post-processing data to read later on, so
          // cache the vertex indexing data
          if(vertexVector.size())
            _vertexVectorCache = vertexVector;
          else
            _vertexMapCache = vertexMap;
          postpro = true;
          break;
        }
        else if(!strncmp(&str[1], "ElementData", 11) ||
                !strncmp(&str[1], "ElementNodeData", 15)) {
          // there's some element post-processing data to read later on
          postpro = true;
          break;
        }
    
        do {
          if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
        } while(str[0] != '$');
      }
    
      // store the elements in their associated elementary entity. If the
      // entity does not exist, create a new (discrete) one.
      for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
        _storeElementsInEntities(elements[i]);
    
      // associate the correct geometrical entity with each mesh vertex
      _associateEntityWithMeshVertices();
    
      // store the vertices in their associated geometrical entity
      if(vertexVector.size())
        _storeVerticesInEntities(vertexVector);
      else
        _storeVerticesInEntities(vertexMap);
    
      // store the physical tags
      for(int i = 0; i < 4; i++) _storePhysicalTagsInEntities(i, physicals[i]);
    
      // copying periodic information from the mesh
    
      if(!CTX::instance()->mesh.ignorePeriodicity) {
        rewind(fp);
    
        while(1) {
          while(str[0] != '$') {
            if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
          }
          if(feof(fp)) break;
    
          if(!strncmp(&str[1], "Periodic", 8) &&
             strncmp(&str[1], "PeriodicNodes", 13)) {
            readMSHPeriodicNodes(fp, this);
            break;
          }
          do {
            if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
          } while(str[0] != '$');
        }
      }
    
      fclose(fp);
    
      if(!CTX::instance()->mesh.ignorePeriodicity) alignPeriodicBoundaries();
    
      return postpro ? 2 : 1;
    }
    
    template <class T>
    static void writeElementMSH(FILE *fp, GModel *model, GEntity *ge, T *ele,
                                bool saveAll, double version, bool binary, int &num,
                                int elementary, std::vector<int> &physicals,
                                int parentNum = 0, int dom1Num = 0, int dom2Num = 0)
    {
      if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
         && ge->getParentEntity()->dim() > ge->dim())
        return; // ignore partition boundaries
    
      std::vector<short> ghosts;
      if(model->getGhostCells().size()) {
        std::pair<std::multimap<MElement *, short>::iterator,
                  std::multimap<MElement *, short>::iterator>
          itp = model->getGhostCells().equal_range(ele);
        for(std::multimap<MElement *, short>::iterator it = itp.first;
            it != itp.second; it++)
          ghosts.push_back(it->second);
      }
    
      if(saveAll)
        ele->writeMSH2(fp, version, binary, ++num, elementary, 0, parentNum,
                       dom1Num, dom2Num, &ghosts);
      else {
        if(parentNum) parentNum = parentNum - physicals.size() + 1;
        for(std::size_t j = 0; j < physicals.size(); j++) {
          ele->writeMSH2(fp, version, binary, ++num, elementary, physicals[j],
                         parentNum, dom1Num, dom2Num, &ghosts);
          if(parentNum) parentNum++;
        }
      }
    
      model->setMeshElementIndex(ele, num); // should really be a multimap...
    
      if(CTX::instance()->mesh.saveTri && ele->getNumChildren())
        num += ele->getNumChildren() - 1;
    }
    
    template <class T>
    static void writeElementsMSH(FILE *fp, GModel *model, GEntity *ge,
                                 std::vector<T *> &ele,
                                 bool saveAll, int saveSinglePartition,
                                 double version, bool binary, int &num,
                                 int elementary, std::vector<int> &physicals)
    {
      if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
         && ge->getParentEntity()->dim() > ge->dim())
        return; // ignore partition boundaries
    
      // Hack to save each partition as a separate physical entity
      if(saveSinglePartition < 0) {
        if(ele.empty()) return;
        int offset = -saveSinglePartition;
        int maxPhysical = model->getMaxPhysicalNumber(ele[0]->getDim());
        for(std::size_t i = 0; i < ele.size(); i++) {
          std::vector<int> newPhysicals;
          int newElementary = elementary;
          if(elementary > 0) { // classical entity
            newElementary = elementary * offset + ele[i]->getPartition();
            for(std::size_t j = 0; j < physicals.size(); j++)
              newPhysicals.push_back(physicals[j] * offset +
                                     ele[i]->getPartition());
          }
          else if(elementary < 0) { // partition boundary
            newPhysicals.push_back((maxPhysical - elementary) * offset);
          }
          ele[i]->setPartition(0);
          writeElementMSH(fp, model, ge, ele[i], saveAll, version, binary, num,
                          newElementary, newPhysicals);
        }
        return;
      }
    
      for(std::size_t i = 0; i < ele.size(); i++) {
        if(saveSinglePartition && ele[i]->getPartition() != saveSinglePartition)
          continue;
        if(ele[i]->getDomain(0)) continue;
        int parentNum = 0;
        MElement *parent = ele[i]->getParent();
        if(parent) parentNum = model->getMeshElementIndex(parent);
        writeElementMSH(fp, model, ge, ele[i], saveAll, version, binary, num,
                        elementary, physicals, parentNum);
      }
    }
    
    static int getNumElementsMSH(GEntity *ge, bool saveAll, int saveSinglePartition)
    {
      if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
         && ge->getParentEntity()->dim() > ge->dim())
        return 0; // ignore partition boundaries
    
      int n = 0, p = saveAll ? 1 : ge->physicals.size();
    
      if(saveSinglePartition < 0 && ge->tag() < 0) p = 1; // partition boundary
    
      if(saveSinglePartition <= 0)
        n = p * ge->getNumMeshElements();
      else
        for(std::size_t i = 0; i < ge->getNumMeshElements(); i++)
          if(ge->getMeshElement(i)->getPartition() == saveSinglePartition) n += p;
      return n;
    }
    
    static int getNumElementsMSH(GModel *m, bool saveAll, int saveSinglePartition)
    {
      int n = 0;
      for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(!CTX::instance()->mesh.saveTri) {
          for(std::size_t i = 0; i < (*it)->points.size(); i++)
            if((*it)->points[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
        }
      }
      for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(!CTX::instance()->mesh.saveTri) {
          for(std::size_t i = 0; i < (*it)->lines.size(); i++)
            if((*it)->lines[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
        }
      }
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(CTX::instance()->mesh.saveTri) {
          for(std::size_t i = 0; i < (*it)->polygons.size(); i++) {
            int nbC = (*it)->polygons[i]->getNumChildren() - 1;
            n += (saveAll ? nbC : nbC * (*it)->physicals.size());
          }
        }
        else {
          for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
            if((*it)->triangles[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
          for(std::size_t i = 0; i < (*it)->polygons.size(); i++)
            if((*it)->polygons[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
        }
      }
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
        n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
        if(CTX::instance()->mesh.saveTri) {
          for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++) {
            int nbC = (*it)->polyhedra[i]->getNumChildren() - 1;
            n += (saveAll ? nbC : nbC * (*it)->physicals.size());
          }
        }
        else {
          for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++)
            if((*it)->tetrahedra[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
          for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++)
            if((*it)->polyhedra[i]->ownsParent())
              n += (saveAll ? 1 : (*it)->physicals.size());
        }
        n -= (*it)->trihedra.size();
      }
      return n;
    }
    
    static int _getElementary(GEntity *ge)
    {
      if(CTX::instance()->mesh.partitionOldStyleMsh2 && ge->getParentEntity()
         && ge->getParentEntity()->dim() == ge->dim()) {
        // hack for backward compatibility of partitioned meshes in MSH2 format: use
        // elementary tag of parent entity if they are of the same dimension
        // (i.e. if they are not partition boundaries)
        return ge->getParentEntity()->tag();
      }
      else {
        return ge->tag();
      }
    }
    
    int GModel::_writeMSH2(const std::string &name, double version, bool binary,
                           bool saveAll, bool saveParametric, double scalingFactor,
                           int elementStartNum, int saveSinglePartition,
                           bool append, bool renumberVertices)
    {
      FILE *fp;
      if(append)
        fp = Fopen(name.c_str(), binary ? "ab" : "a");
      else
        fp = Fopen(name.c_str(), binary ? "wb" : "w");
      if(!fp) {
        Msg::Error("Unable to open file '%s'", name.c_str());
        return 0;
      }
    
      // binary format exists only in version 2
      if(version > 1 || binary)
        version = 2.2;
      else
        version = 1.0;
    
      // if there are no physicals we save all the elements
      if(noPhysicalGroups()) saveAll = true;
    
      // get the number of vertices and index the vertices in a continuous
      // sequence
    
      int numVertices =
        indexMeshVertices(saveAll, saveSinglePartition, renumberVertices);
    
      // get the number of elements we need to save
      int numElements = getNumElementsMSH(this, saveAll, saveSinglePartition);
    
      if(version >= 2.0) {
        fprintf(fp, "$MeshFormat\n");
        fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
        if(binary) {
          int one = 1;
          fwrite(&one, sizeof(int), 1, fp);
          fprintf(fp, "\n");
        }
        fprintf(fp, "$EndMeshFormat\n");
    
        if(numPhysicalNames()) {
          fprintf(fp, "$PhysicalNames\n");
          fprintf(fp, "%d\n", numPhysicalNames());
          for(piter it = firstPhysicalName(); it != lastPhysicalName(); it++) {
            std::string name = it->second;
            if(name.size() > 128) name.resize(128);
            fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
                    name.c_str());
          }
          fprintf(fp, "$EndPhysicalNames\n");
        }
    
        if(CTX::instance()->mesh.saveTopology) writeMSHEntities(fp, this);
    
        if(saveParametric)
          fprintf(fp, "$ParametricNodes\n");
        else
          fprintf(fp, "$Nodes\n");
      }
      else
        fprintf(fp, "$NOD\n");
    
      fprintf(fp, "%d\n", numVertices);
    
      std::vector<GEntity *> entities;
      getEntities(entities);
      for(std::size_t i = 0; i < entities.size(); i++)
        for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
          entities[i]->mesh_vertices[j]->writeMSH2(fp, binary, saveParametric,
                                                   scalingFactor);
    
      if(binary) fprintf(fp, "\n");
    
      if(version >= 2.0) {
        if(saveParametric)
          fprintf(fp, "$EndParametricNodes\n");
        else
          fprintf(fp, "$EndNodes\n");
        fprintf(fp, "$Elements\n");
      }
      else {
        fprintf(fp, "$ENDNOD\n");
        fprintf(fp, "$ELM\n");
      }
    
      fprintf(fp, "%d\n", numElements);
      int num = elementStartNum;
    
      _elementIndexCache.clear();
    
      // parents
      if(!CTX::instance()->mesh.saveTri) {
        for(viter it = firstVertex(); it != lastVertex(); ++it) {
          for(std::size_t i = 0; i < (*it)->points.size(); i++)
            if((*it)->points[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->points[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
        for(eiter it = firstEdge(); it != lastEdge(); ++it) {
          for(std::size_t i = 0; i < (*it)->lines.size(); i++)
            if((*it)->lines[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->lines[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
        for(fiter it = firstFace(); it != lastFace(); ++it) {
          for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
            if((*it)->triangles[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->triangles[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
        for(riter it = firstRegion(); it != lastRegion(); ++it) {
          for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++)
            if((*it)->tetrahedra[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->tetrahedra[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
        for(fiter it = firstFace(); it != lastFace(); ++it) {
          for(std::size_t i = 0; i < (*it)->polygons.size(); i++)
            if((*it)->polygons[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->polygons[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
        for(riter it = firstRegion(); it != lastRegion(); ++it) {
          for(std::size_t i = 0; i < (*it)->polyhedra.size(); i++)
            if((*it)->polyhedra[i]->ownsParent())
              writeElementMSH(fp, this, *it, (*it)->polyhedra[i]->getParent(), saveAll,
                              version, binary, num, _getElementary(*it),
                              (*it)->physicals);
        }
      }
      // points
      for(viter it = firstVertex(); it != lastVertex(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->points, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // lines
      for(eiter it = firstEdge(); it != lastEdge(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->lines, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // triangles
      for(fiter it = firstFace(); it != lastFace(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->triangles, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // quads
      for(fiter it = firstFace(); it != lastFace(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->quadrangles, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // polygons
      for(fiter it = firstFace(); it != lastFace(); it++) {
        writeElementsMSH(fp, this, *it, (*it)->polygons, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // tets
      for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->tetrahedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // hexas
      for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->hexahedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // prisms
      for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->prisms, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // pyramids
      for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->pyramids, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // polyhedra
      for(riter it = firstRegion(); it != lastRegion(); ++it) {
        writeElementsMSH(fp, this, *it, (*it)->polyhedra, saveAll, saveSinglePartition,
                         version, binary, num, _getElementary(*it),
                         (*it)->physicals);
      }
      // level set faces
      for(fiter it = firstFace(); it != lastFace(); ++it) {
        for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
          MTriangle *t = (*it)->triangles[i];
          if(t->getDomain(0))
            writeElementMSH(fp, this, *it, t, saveAll, version, binary, num,
                            _getElementary(*it), (*it)->physicals, 0,
                            getMeshElementIndex(t->getDomain(0)),
                            getMeshElementIndex(t->getDomain(1)));
        }
        for(std::size_t i = 0; i < (*it)->polygons.size(); i++) {
          MPolygon *p = (*it)->polygons[i];
          if(p->getDomain(0))
            writeElementMSH(fp, this, *it, p, saveAll, version, binary, num,
                            _getElementary(*it), (*it)->physicals, 0,
                            getMeshElementIndex(p->getDomain(0)),
                            getMeshElementIndex(p->getDomain(1)));
        }
      }
      // level set lines
      for(eiter it = firstEdge(); it != lastEdge(); ++it) {
        for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
          MLine *l = (*it)->lines[i];
          if(l->getDomain(0))
            writeElementMSH(fp, this, *it, l, saveAll, version, binary, num,
                            _getElementary(*it), (*it)->physicals, 0,
                            getMeshElementIndex(l->getDomain(0)),
                            getMeshElementIndex(l->getDomain(1)));
        }
      }
    
      if(binary) fprintf(fp, "\n");
    
      if(version >= 2.0) {
        fprintf(fp, "$EndElements\n");
      }
      else {
        fprintf(fp, "$ENDELM\n");
      }
    
      writeMSHPeriodicNodes(fp, entities, renumberVertices, saveAll);
    
      fclose(fp);
    
      return 1;
    }
    
    int GModel::_writePartitionedMSH2(const std::string &baseName, bool binary,
                                      bool saveAll, bool saveParametric,
                                      double scalingFactor)
    {
      int numElements;
      int startNum = 0;
      for(std::size_t partition = 1; partition <= getNumPartitions();
          partition++) {
        std::ostringstream sstream;
        sstream << baseName << "_" << partition << ".msh";
    
        numElements = getNumElementsMSH(this, saveAll, partition);
        Msg::Info("Writing partition %d in file '%s'", partition,
                  sstream.str().c_str());
        _writeMSH2(sstream.str(), 2.2, binary, saveAll, saveParametric,
                   scalingFactor, startNum, partition, false,
                   false); // NO RENUMBERING!
        startNum += numElements; // update for next iteration in the loop
      }
    
    #if 0
      if(_ghostCells.size()){
        Msg::Info("Writing ghost cells in debug file 'ghosts.pos'");
        FILE *fp = Fopen("ghosts.pos", "w");
        fprintf(fp, "View \"ghosts\"{\n");
        for(std::multimap<MElement*, short>::iterator it = _ghostCells.begin();
            it != _ghostCells.end(); it++)
          it->first->writePOS(fp, false, true, false, false, false, false);
        fprintf(fp, "};\n");
        fclose(fp);
      }
    #endif
    
      return 1;
    }