// 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. // // Contributor(s): // Anthony Royer #include <cstdio> #include <fstream> #include <vector> #include <map> #include <algorithm> #include <sstream> #include <string> #include <cstdlib> #include <limits> #include "GmshDefines.h" #include "OS.h" #include "Context.h" #include "GModel.h" #include "GEntity.h" #include "partitionRegion.h" #include "partitionFace.h" #include "partitionEdge.h" #include "partitionVertex.h" #include "discreteEdge.h" #include "ghostFace.h" #include "ghostEdge.h" #include "ghostRegion.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 "MTrihedron.h" #include "StringUtils.h" static bool readMSH4Physicals(GModel *const model, FILE *fp, GEntity *const entity, bool binary, char *str, bool swap) { std::size_t numPhy = 0; if(binary) { if(fread(&numPhy, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numPhy, sizeof(std::size_t), 1); std::vector<int> phyTags(numPhy); if(fread(&phyTags[0], sizeof(int), numPhy, fp) != numPhy) { return false; } if(swap) SwapBytes((char *)&phyTags[0], sizeof(int), numPhy); for(std::size_t i = 0; i < numPhy; i++) { entity->addPhysicalEntity(phyTags[i]); } } else { sscanf(str, "%lu %[0-9- ]", &numPhy, str); for(std::size_t i = 0; i < numPhy; i++) { int phyTag = 0; if(i == numPhy - 1 && entity->dim() == 0) { if(sscanf(str, "%d", &phyTag) != 1) { return false; } } else { if(sscanf(str, "%d %[0-9- ]", &phyTag, str) != 2) { return false; } } entity->addPhysicalEntity(phyTag); } } return true; } static bool readMSH4BoundingEntities(GModel *const model, FILE *fp, GEntity *const entity, bool binary, char *str, bool swap) { std::size_t numBrep = 0; std::vector<GEntity *> boundingEntities; std::vector<int> boundingSign; if(binary) { if(fread(&numBrep, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numBrep, sizeof(std::size_t), 1); std::vector<int> brepTags(numBrep); if(fread(&brepTags[0], sizeof(int), numBrep, fp) != numBrep) { return false; } if(swap) SwapBytes((char *)&brepTags[0], sizeof(int), numBrep); for(std::size_t i = 0; i < numBrep; i++) { GEntity *brep = model->getEntityByTag(entity->dim() - 1, std::abs(brepTags[i])); if(!brep) { Msg::Warning("Entity %d not found in the Brep of entity %d", brepTags[i], entity->tag()); } else { boundingEntities.push_back(brep); boundingSign.push_back((std::abs(brepTags[i]) == brepTags[i] ? 1 : -1)); } } } else { sscanf(str, "%lu %[0-9- ]", &numBrep, str); for(std::size_t i = 0; i < numBrep; i++) { int entityTag = 0; if(i != numBrep - 1) { if(sscanf(str, "%d %[0-9- ]", &entityTag, str) != 2) { return false; } } else { if(sscanf(str, "%d", &entityTag) != 1) { return false; } } GEntity *brep = model->getEntityByTag(entity->dim() - 1, std::abs(entityTag)); if(!brep) { Msg::Warning("Entity %d not found in the Brep of entity %d", entityTag, entity->tag()); } else { boundingEntities.push_back(brep); boundingSign.push_back((std::abs(entityTag) == entityTag ? 1 : -1)); } } } switch(entity->dim()) { case 1: if(boundingEntities.size() == 2) { reinterpret_cast<discreteEdge *>(entity)->setBeginVertex( reinterpret_cast<GVertex *>(boundingEntities[0])); reinterpret_cast<discreteEdge *>(entity)->setEndVertex( reinterpret_cast<GVertex *>(boundingEntities[1])); } else if(boundingEntities.size() == 1) { if(boundingSign[0] == 1) { reinterpret_cast<discreteEdge *>(entity)->setBeginVertex( reinterpret_cast<GVertex *>(boundingEntities[0])); } else { reinterpret_cast<discreteEdge *>(entity)->setEndVertex( reinterpret_cast<GVertex *>(boundingEntities[0])); } } break; case 2: { std::vector<int> tags(boundingEntities.size()); for(std::size_t i = 0; i < boundingEntities.size(); i++) tags[i] = std::abs(boundingEntities[i]->tag()); reinterpret_cast<discreteFace *>(entity)->setBoundEdges(tags, boundingSign); } break; case 3: { std::vector<int> tags(boundingEntities.size()); for(std::size_t i = 0; i < boundingEntities.size(); i++) tags[i] = std::abs(boundingEntities[i]->tag()); reinterpret_cast<discreteRegion *>(entity)->setBoundFaces(tags, boundingSign); } break; default: break; } return true; } static bool readMSH4EntityInfo(FILE *fp, bool binary, char *str, int sizeofstr, bool swap, double version, bool partition, int dim, int &tag, int &parentDim, int &parentTag, std::vector<int> &partitions, double &minX, double &minY, double &minZ, double &maxX, double &maxY, double &maxZ) { if(partition) { if(binary) { int dataInt[3]; if(fread(dataInt, sizeof(int), 3, fp) != 3) { return false; } if(swap) SwapBytes((char *)dataInt, sizeof(int), 3); tag = dataInt[0]; parentDim = dataInt[1]; parentTag = dataInt[2]; std::size_t numPart; if(fread(&numPart, sizeof(std::size_t), 1, fp) != 1) { return false; } partitions.resize(numPart, 0); if(fread(&partitions[0], sizeof(int), numPart, fp) != numPart) { return false; } if(swap) SwapBytes((char *)&partitions[0], sizeof(int), numPart); double dataDouble[6]; const std::size_t nbb = (dim > 0) ? 6 : 3; if(fread(dataDouble, sizeof(double), nbb, fp) != nbb) { return false; } if(swap) SwapBytes((char *)dataDouble, sizeof(double), nbb); minX = dataDouble[0]; minY = dataDouble[1]; minZ = dataDouble[2]; maxX = dataDouble[(nbb == 6) ? 3 : 0]; maxY = dataDouble[(nbb == 6) ? 4 : 1]; maxZ = dataDouble[(nbb == 6) ? 5 : 2]; } else { std::size_t numPart = 0; if(fscanf(fp, "%d %d %d %lu", &tag, &parentDim, &parentTag, &numPart) != 4) { return false; } partitions.resize(numPart, 0); for(std::size_t i = 0; i < numPart; i++) { if(fscanf(fp, "%d", &partitions[i]) != 1) { return false; } } if(version < 4.1 || dim > 0) { if(fscanf(fp, "%lf %lf %lf %lf %lf %lf", &minX, &minY, &minZ, &maxX, &maxY, &maxZ) != 6) { return false; } } else { if(fscanf(fp, "%lf %lf %lf", &minX, &minY, &minZ) != 3) { return false; } maxX = minX; maxY = minY; maxZ = minZ; } if(!fgets(str, sizeofstr, fp)) { return false; } } } else { if(binary) { int dataInt; if(fread(&dataInt, sizeof(int), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&dataInt, sizeof(int), 1); double dataDouble[6]; const std::size_t nbb = (dim > 0) ? 6 : 3; if(fread(dataDouble, sizeof(double), nbb, fp) != nbb) { return false; } if(swap) SwapBytes((char *)dataDouble, sizeof(double), nbb); tag = dataInt; minX = dataDouble[0]; minY = dataDouble[1]; minZ = dataDouble[2]; maxX = dataDouble[(nbb == 6) ? 3 : 0]; maxY = dataDouble[(nbb == 6) ? 4 : 1]; maxZ = dataDouble[(nbb == 6) ? 5 : 2]; } else { if(version < 4.1 || dim > 0) { if(fscanf(fp, "%d %lf %lf %lf %lf %lf %lf", &tag, &minX, &minY, &minZ, &maxX, &maxY, &maxZ) != 7) { return false; } } else { if(fscanf(fp, "%d %lf %lf %lf", &tag, &minX, &minY, &minZ) != 4) { return false; } maxX = minX; maxY = minY; maxZ = minZ; } if(!fgets(str, sizeofstr, fp)) { return false; } } } return true; } static bool readMSH4Entities(GModel *const model, FILE *fp, bool partition, bool binary, bool swap, double version) { if(partition) { std::size_t numPartitions = 0; std::size_t ghostSize = 0; std::vector<int> ghostTags; if(binary) { if(fread(&numPartitions, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numPartitions, sizeof(std::size_t), 1); if(fread(&ghostSize, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&ghostSize, sizeof(std::size_t), 1); if(ghostSize) { ghostTags.resize(2 * ghostSize); if(fread(&ghostTags[0], sizeof(int), 2 * ghostSize, fp) != 2 * ghostSize) { return false; } if(swap) SwapBytes((char *)&ghostTags[0], sizeof(int), 2 * ghostSize); } } else { if(fscanf(fp, "%lu", &numPartitions) != 1) { return false; } if(fscanf(fp, "%lu", &ghostSize) != 1) { return false; } if(ghostSize) { ghostTags.resize(2 * ghostSize); for(std::size_t i = 0; i < 2 * ghostSize; i += 2) { if(fscanf(fp, "%d %d", &ghostTags[i], &ghostTags[i + 1]) != 2) { return false; } } } } model->setNumPartitions(numPartitions); Msg::Info("%lu partitions", model->getNumPartitions()); for(std::size_t i = 0; i < 2 * ghostSize; i += 2) { switch(model->getDim()) { case 1: { ghostEdge *ghostEntities = new ghostEdge(model, ghostTags[i], ghostTags[i + 1]); model->add(ghostEntities); } break; case 2: { ghostFace *ghostEntities = new ghostFace(model, ghostTags[i], ghostTags[i + 1]); model->add(ghostEntities); } break; case 3: { ghostRegion *ghostEntities = new ghostRegion(model, ghostTags[i], ghostTags[i + 1]); model->add(ghostEntities); } break; default: break; } } } std::size_t numEntities[4] = {0, 0, 0, 0}; if(binary) { if(fread(numEntities, sizeof(std::size_t), 4, fp) != 4) { return false; } if(swap) SwapBytes((char *)numEntities, sizeof(std::size_t), 4); } else { if(fscanf(fp, "%lu %lu %lu %lu", &numEntities[0], &numEntities[1], &numEntities[2], &numEntities[3]) != 4) { return false; } } // max length of line for ascii input file (should be large enough to handle // entities with many entities on their boundary) int nume = numEntities[0] + numEntities[1] + numEntities[2] + numEntities[3]; std::size_t strl = std::max(4096, 128 * nume); char *str = new char[strl]; for(int dim = 0; dim < 4; dim++) { for(std::size_t i = 0; i < numEntities[dim]; i++) { int tag = 0, parentDim = 0, parentTag = 0; std::vector<int> parts; double minX = 0., minY = 0., minZ = 0., maxX = 0., maxY = 0., maxZ = 0.; if(!readMSH4EntityInfo(fp, binary, str, strl, swap, version, partition, dim, tag, parentDim, parentTag, parts, minX, minY, minZ, maxX, maxY, maxZ)) { delete[] str; return false; } // FIXME: rewrite partition classes in terms of int std::vector<unsigned int> partitions(parts.begin(), parts.end()); switch(dim) { case 0: { GVertex *gv = model->getVertexByTag(tag); if(!gv) { if(partition) { gv = new partitionVertex(model, tag, partitions); if(parentTag) static_cast<partitionVertex *>(gv)->setParentEntity( model->getEntityByTag(parentDim, parentTag)); } else { gv = new discreteVertex(model, tag, minX, minY, minZ); } model->add(gv); } if(!readMSH4Physicals(model, fp, gv, binary, str, swap)) { delete[] str; return false; } } break; case 1: { GEdge *ge = model->getEdgeByTag(tag); if(!ge) { if(partition) { ge = new partitionEdge(model, tag, 0, 0, partitions); if(parentTag) static_cast<partitionEdge *>(ge)->setParentEntity( model->getEntityByTag(parentDim, parentTag)); } else { ge = new discreteEdge(model, tag, 0, 0); } model->add(ge); } if(!readMSH4Physicals(model, fp, ge, binary, str, swap)) { delete[] str; return false; } if(!readMSH4BoundingEntities(model, fp, ge, binary, str, swap)) { delete[] str; return false; } } break; case 2: { GFace *gf = model->getFaceByTag(tag); if(!gf) { if(partition) { gf = new partitionFace(model, tag, partitions); if(parentTag) static_cast<partitionFace *>(gf)->setParentEntity( model->getEntityByTag(parentDim, parentTag)); } else { gf = new discreteFace(model, tag); } model->add(gf); } if(!readMSH4Physicals(model, fp, gf, binary, str, swap)) { delete[] str; return false; } if(!readMSH4BoundingEntities(model, fp, gf, binary, str, swap)) { delete[] str; return false; } } break; case 3: { GRegion *gr = model->getRegionByTag(tag); if(!gr) { if(partition) { gr = new partitionRegion(model, tag, partitions); if(parentTag) static_cast<partitionRegion *>(gr)->setParentEntity( model->getEntityByTag(parentDim, parentTag)); } else { gr = new discreteRegion(model, tag); } model->add(gr); } if(!readMSH4Physicals(model, fp, gr, binary, str, swap)) { delete[] str; return false; } if(!readMSH4BoundingEntities(model, fp, gr, binary, str, swap)) { delete[] str; return false; } } break; } } } delete[] str; return true; } static std::pair<std::size_t, MVertex *> * readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense, std::size_t &totalNumNodes, std::size_t &maxNodeNum, bool swap, double version) { std::size_t numBlock = 0, minTag = 0, maxTag = 0; totalNumNodes = 0; maxNodeNum = 0; if(binary) { std::size_t data[4]; if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return 0; } if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4); numBlock = data[0]; totalNumNodes = data[1]; minTag = data[2]; maxTag = data[3]; } else { if(version >= 4.1) { if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumNodes, &minTag, &maxTag) != 4) { return 0; } } else { if(fscanf(fp, "%lu %lu", &numBlock, &totalNumNodes) != 2) { return 0; } } } std::size_t nodeRead = 0; std::size_t minNodeNum = std::numeric_limits<std::size_t>::max(); std::pair<std::size_t, MVertex *> *vertexCache = new std::pair<std::size_t, MVertex *>[totalNumNodes]; Msg::Info("%lu nodes", totalNumNodes); for(std::size_t i = 0; i < numBlock; i++) { int parametric = 0; int entityTag = 0, entityDim = 0; std::size_t numNodes = 0; if(binary) { int data[3]; if(fread(data, sizeof(int), 3, fp) != 3) { delete[] vertexCache; return 0; } if(swap) SwapBytes((char *)data, sizeof(int), 3); entityDim = data[0]; entityTag = data[1]; parametric = data[2]; if(fread(&numNodes, sizeof(std::size_t), 1, fp) != 1) { delete[] vertexCache; return 0; } if(swap) SwapBytes((char *)&numNodes, sizeof(std::size_t), 1); } else { if(version >= 4.1) { if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, ¶metric, &numNodes) != 4) { delete[] vertexCache; return 0; } } else { if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, ¶metric, &numNodes) != 4) { delete[] vertexCache; return 0; } } } GEntity *entity = model->getEntityByTag(entityDim, entityTag); if(!entity) { switch(entityDim) { case 0: { Msg::Info("Creating discrete point %d", entityTag); GVertex *gv = new discreteVertex(model, entityTag); GModel::current()->add(gv); entity = gv; break; } case 1: { Msg::Info("Creating discrete curve %d", entityTag); GEdge *ge = new discreteEdge(model, entityTag, 0, 0); GModel::current()->add(ge); entity = ge; break; } case 2: { Msg::Info("Creating discrete surface %d", entityTag); GFace *gf = new discreteFace(model, entityTag); GModel::current()->add(gf); entity = gf; break; } case 3: { Msg::Info("Creating discrete volume %d", entityTag); GRegion *gr = new discreteRegion(model, entityTag); GModel::current()->add(gr); entity = gr; break; } default: { Msg::Error("Invalid dimension %d to create discrete entity", entityDim); delete[] vertexCache; return 0; } } } std::size_t n = 3; if(parametric) n += entityDim; std::vector<std::size_t> tags(numNodes); if(binary) { if(fread(&tags[0], sizeof(std::size_t), numNodes, fp) != numNodes) { delete[] vertexCache; return 0; } if(swap) SwapBytes((char *)&tags[0], sizeof(std::size_t), numNodes); std::vector<double> coord(n * numNodes); if(fread(&coord[0], sizeof(double), n * numNodes, fp) != n * numNodes) { delete[] vertexCache; return 0; } if(swap) SwapBytes((char *)&coord[0], sizeof(double), n * numNodes); std::size_t k = 0; for(std::size_t j = 0; j < numNodes; j++) { MVertex *mv = 0; std::size_t tagNode = tags[j]; if(n == 5) { mv = new MFaceVertex(coord[k], coord[k + 1], coord[k + 2], entity, coord[k + 3], coord[k + 4], tagNode); } else if(n == 4) { mv = new MEdgeVertex(coord[k], coord[k + 1], coord[k + 2], entity, coord[k + 3], tagNode); } else { mv = new MVertex(coord[k], coord[k + 1], coord[k + 2], entity, tagNode); } k += n; entity->addMeshVertex(mv); mv->setEntity(entity); minNodeNum = std::min(minNodeNum, tagNode); maxNodeNum = std::max(maxNodeNum, tagNode); vertexCache[nodeRead] = std::pair<std::size_t, MVertex *>(tagNode, mv); nodeRead++; if(totalNumNodes > 100000) Msg::ProgressMeter(nodeRead, totalNumNodes, true, "Reading nodes"); } } else { if(version >= 4.1) { for(std::size_t j = 0; j < numNodes; j++) { if(fscanf(fp, "%lu", &tags[j]) != 1) { delete[] vertexCache; return 0; } } } for(std::size_t j = 0; j < numNodes; j++) { std::size_t tagNode = 0; if(version >= 4.1) { tagNode = tags[j]; } else { if(fscanf(fp, "%lu", &tagNode) != 1) { delete[] vertexCache; return 0; } } MVertex *mv = 0; if(n == 5) { double x, y, z, u, v; if(fscanf(fp, "%lf %lf %lf %lf %lf", &x, &y, &z, &u, &v) != 5) { delete[] vertexCache; return 0; } mv = new MFaceVertex(x, y, z, entity, u, v, tagNode); } else if(n == 4) { double x, y, z, u; if(fscanf(fp, "%lf %lf %lf %lf", &x, &y, &z, &u) != 4) { delete[] vertexCache; return 0; } mv = new MEdgeVertex(x, y, z, entity, u, tagNode); } else { double x, y, z; if(fscanf(fp, "%lf %lf %lf", &x, &y, &z) != 3) { delete[] vertexCache; return 0; } // discard extra parametric coordinates, as Gmsh does not use them for(std::size_t k = 3; k < n; k++) { double dummy; if(fscanf(fp, "%lf", &dummy) != 1) { delete[] vertexCache; return 0; } } mv = new MVertex(x, y, z, entity, tagNode); } entity->addMeshVertex(mv); mv->setEntity(entity); minNodeNum = std::min(minNodeNum, tagNode); maxNodeNum = std::max(maxNodeNum, tagNode); vertexCache[nodeRead] = std::pair<std::size_t, MVertex *>(tagNode, mv); nodeRead++; if(totalNumNodes > 100000) Msg::ProgressMeter(nodeRead, totalNumNodes, true, "Reading nodes"); } } } if(version >= 4.1) { // consistency check if(minTag != minNodeNum || maxTag != maxNodeNum) Msg::Warning("Min/Max node tags reported in section header are wrong: " "(%d/%d) != (%d/%d)", minTag, maxTag, minNodeNum, maxNodeNum); } // if the vertex numbering is (fairly) dense, we fill the vector cache, // otherwise we fill the map cache if(minNodeNum == 1 && maxNodeNum == totalNumNodes) { Msg::Debug("Vertex numbering is dense"); dense = true; } else if(maxNodeNum < 10 * totalNumNodes) { Msg::Debug( "Vertex numbering is fairly dense - still caching with a vector"); dense = true; } else { Msg::Debug("Vertex numbering is not dense"); dense = false; } return vertexCache; } static std::pair<std::size_t, MElement *> * readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense, std::size_t &totalNumElements, std::size_t &maxElementNum, bool swap, double version) { char str[1024]; std::size_t numBlock = 0, minTag = 0, maxTag = 0; totalNumElements = 0; maxElementNum = 0; if(binary) { std::size_t data[4]; if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return 0; } if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4); numBlock = data[0]; totalNumElements = data[1]; minTag = data[2]; maxTag = data[3]; } else { if(version >= 4.1) { if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumElements, &minTag, &maxTag) != 4) { return 0; } } else { if(fscanf(fp, "%lu %lu", &numBlock, &totalNumElements) != 2) { return 0; } } } std::size_t elementRead = 0; std::size_t minElementNum = std::numeric_limits<std::size_t>::max(); std::pair<std::size_t, MElement *> *elementCache = new std::pair<std::size_t, MElement *>[totalNumElements]; Msg::Info("%lu elements", totalNumElements); for(std::size_t i = 0; i < numBlock; i++) { int entityTag = 0, entityDim = 0, elmType = 0; std::size_t numElements = 0; if(binary) { int data[3]; if(fread(data, sizeof(int), 3, fp) != 3) { delete[] elementCache; return 0; } if(swap) SwapBytes((char *)data, sizeof(int), 3); entityDim = data[0]; entityTag = data[1]; elmType = data[2]; if(fread(&numElements, sizeof(std::size_t), 1, fp) != 1) { delete[] elementCache; return 0; } if(swap) SwapBytes((char *)&numElements, sizeof(std::size_t), 1); } else { if(version >= 4.1) { if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, &elmType, &numElements) != 4) { delete[] elementCache; return 0; } } else { if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, &elmType, &numElements) != 4) { delete[] elementCache; return 0; } } } GEntity *entity = model->getEntityByTag(entityDim, entityTag); if(!entity) { Msg::Error("Unknown entity %d of dimension %d", entityTag, entityDim); delete[] elementCache; return 0; } if(entity->geomType() == GEntity::GhostCurve) { static_cast<ghostEdge *>(entity)->haveMesh(true); } else if(entity->geomType() == GEntity::GhostSurface) { static_cast<ghostFace *>(entity)->haveMesh(true); } else if(entity->geomType() == GEntity::GhostVolume) { static_cast<ghostRegion *>(entity)->haveMesh(true); } const int numVertPerElm = MElement::getInfoMSH(elmType); if(binary) { std::size_t n = 1 + numVertPerElm; std::vector<std::size_t> data(numElements * n); if(fread(&data[0], sizeof(std::size_t), numElements * n, fp) != numElements * n) { delete[] elementCache; return 0; } if(swap) SwapBytes((char *)&data[0], sizeof(std::size_t), numElements * n); std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)0); for(std::size_t j = 0; j < numElements * n; j += n) { for(int k = 0; k < numVertPerElm; k++) { vertices[k] = model->getMeshVertexByTag(data[j + k + 1]); if(!vertices[k]) { Msg::Error("Unknown node %lu in element %lu", data[j + k + 1], data[j]); delete[] elementCache; return 0; } } MElementFactory elementFactory; MElement *element = elementFactory.create(elmType, vertices, data[j], 0, false, 0, 0, 0, 0); if(!element) { Msg::Error("Could not create element %lu of type %d", data[j], elmType); delete[] elementCache; return 0; } if(entity->geomType() != GEntity::GhostCurve && entity->geomType() != GEntity::GhostSurface && entity->geomType() != GEntity::GhostVolume) { entity->addElement(element->getType(), element); } minElementNum = std::min(minElementNum, data[j]); maxElementNum = std::max(maxElementNum, data[j]); elementCache[elementRead] = std::pair<std::size_t, MElement *>(data[j], element); elementRead++; if(totalNumElements > 100000) Msg::ProgressMeter(elementRead, totalNumElements, true, "Reading elements"); } } else { for(std::size_t j = 0; j < numElements; j++) { std::size_t elmTag = 0; if(fscanf(fp, "%lu", &elmTag) != 1) { delete[] elementCache; return 0; } if(!fgets(str, sizeof(str), fp)) { delete[] elementCache; return 0; } std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)0); for(int k = 0; k < numVertPerElm; k++) { std::size_t vertexTag = 0; if(k != numVertPerElm - 1) { if(sscanf(str, "%lu %[0-9- ]", &vertexTag, str) != 2) { delete[] elementCache; return 0; } } else { if(sscanf(str, "%lu", &vertexTag) != 1) { delete[] elementCache; return 0; } } vertices[k] = model->getMeshVertexByTag(vertexTag); if(!vertices[k]) { Msg::Error("Unknown node %lu in element %lu", vertexTag, elmTag); delete[] elementCache; return 0; } } MElementFactory elementFactory; MElement *element = elementFactory.create(elmType, vertices, elmTag, 0, false, 0, 0, 0, 0); if(!element) { Msg::Error("Could not create element %lu of type %d", elmTag, elmType); delete[] elementCache; return 0; } if(entity->geomType() != GEntity::GhostCurve && entity->geomType() != GEntity::GhostSurface && entity->geomType() != GEntity::GhostVolume) { entity->addElement(element->getType(), element); } minElementNum = std::min(minElementNum, elmTag); maxElementNum = std::max(maxElementNum, elmTag); elementCache[elementRead] = std::pair<std::size_t, MElement *>(elmTag, element); elementRead++; if(totalNumElements > 100000) Msg::ProgressMeter(elementRead, totalNumElements, true, "Reading elements"); } } } // if the vertex numbering is dense, we fill the vector cache, otherwise we // fill the map cache if(minElementNum == 1 && maxElementNum == totalNumElements) { Msg::Debug("Element numbering is dense"); dense = true; } else if(maxElementNum < 10 * totalNumElements) { Msg::Debug( "Element numbering is fairly dense - still caching with a vector"); dense = true; } else { Msg::Debug("Element numbering is not dense"); dense = false; } return elementCache; } static bool readMSH4PeriodicNodes(GModel *const model, FILE *fp, bool binary, bool swap, double version) { std::size_t numPeriodicLinks = 0; if(binary) { if(fread(&numPeriodicLinks, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numPeriodicLinks, sizeof(std::size_t), 1); } else { if(fscanf(fp, "%lu", &numPeriodicLinks) != 1) { return false; } } for(std::size_t i = 0; i < numPeriodicLinks; i++) { int slaveDim = 0, slaveTag = 0, masterTag = 0; if(binary) { int data[3]; if(fread(data, sizeof(int), 3, fp) != 3) { return false; } if(swap) SwapBytes((char *)data, sizeof(int), 3); slaveDim = data[0]; slaveTag = data[1]; masterTag = data[2]; } else { if(fscanf(fp, "%d %d %d", &slaveDim, &slaveTag, &masterTag) != 3) { return false; } } GEntity *slave = 0, *master = 0; switch(slaveDim) { case 0: slave = model->getVertexByTag(slaveTag); master = model->getVertexByTag(masterTag); break; case 1: slave = model->getEdgeByTag(slaveTag); master = model->getEdgeByTag(masterTag); break; case 2: slave = model->getFaceByTag(masterTag); master = model->getFaceByTag(masterTag); break; } if(!slave) { Msg::Info("Could not find periodic entity %d of dimension %d", slaveTag, slaveDim); continue; } if(!master) { Msg::Info("Could not find periodic source entity %d of dimension %d", masterTag, slaveDim); continue; } std::size_t correspondingVertexSize = 0; if(binary) { std::size_t numAffine; if(fread(&numAffine, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numAffine, sizeof(std::size_t), 1); if(numAffine) { std::vector<double> tfo(numAffine); if(fread(&tfo[0], sizeof(double), numAffine, fp) != numAffine) { return false; } if(swap) SwapBytes((char *)&tfo[0], sizeof(double), numAffine); slave->setMeshMaster(master, tfo); } else { slave->setMeshMaster(master); } if(fread(&correspondingVertexSize, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&correspondingVertexSize, sizeof(std::size_t), 1); } else { if(version >= 4.1) { std::size_t numAffine; if(!fscanf(fp, "%lu", &numAffine)) { return false; } if(numAffine) { std::vector<double> tfo(numAffine); for(std::size_t i = 0; i < numAffine; i++) { if(fscanf(fp, "%lf", &tfo[i]) != 1) { return false; } } slave->setMeshMaster(master, tfo); } else { slave->setMeshMaster(master); } if(fscanf(fp, "%lu", &correspondingVertexSize) != 1) { return false; } } else { char affine[256]; if(!fscanf(fp, "%s", affine)) { return false; } if(!strncmp(affine, "Affine", 6)) { if(!fgets(affine, sizeof(affine), fp)) { return false; } std::vector<double> tfo(16); if(sscanf(affine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " "%lf %lf %lf %lf", &tfo[0], &tfo[1], &tfo[2], &tfo[3], &tfo[4], &tfo[5], &tfo[6], &tfo[7], &tfo[8], &tfo[9], &tfo[10], &tfo[11], &tfo[12], &tfo[13], &tfo[14], &tfo[15]) != 16) { return false; } slave->setMeshMaster(master, tfo); if(fscanf(fp, "%lu", &correspondingVertexSize) != 1) { return false; } } else { slave->setMeshMaster(master); if(sscanf(affine, "%lu", &correspondingVertexSize) != 1) { return false; } } } } for(std::size_t j = 0; j < correspondingVertexSize; j++) { std::size_t v1 = 0, v2 = 0; if(binary) { std::size_t data[2]; if(fread(data, sizeof(std::size_t), 2, fp) != 2) { return false; } if(swap) SwapBytes((char *)data, sizeof(std::size_t), 2); v1 = data[0]; v2 = data[1]; } else { if(fscanf(fp, "%lu %lu", &v1, &v2) != 2) { return false; } } MVertex *mv1 = model->getMeshVertexByTag(v1); MVertex *mv2 = model->getMeshVertexByTag(v2); if(mv1 && mv2) { slave->correspondingVertices[mv1] = mv2; } else { if(!mv1) { Msg::Info("Could not find periodic node %d", v1); } else if(!mv2) { Msg::Info("Could not find periodic node %d", v2); } } } } return true; } static bool readMSH4GhostElements(GModel *const model, FILE *fp, bool binary, bool swap) { std::size_t numGhostCells = 0; if(binary) { if(fread(&numGhostCells, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numGhostCells, sizeof(std::size_t), 1); } else { if(fscanf(fp, "%lu", &numGhostCells) != 1) { return false; } } std::multimap<std::pair<MElement *, int>, int> ghostCells; for(std::size_t i = 0; i < numGhostCells; i++) { std::size_t elmTag = 0; int partNum = 0; std::size_t numGhostPartitions = 0; char str[1024]; if(binary) { if(fread(&elmTag, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&elmTag, sizeof(std::size_t), 1); if(fread(&partNum, sizeof(int), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&partNum, sizeof(int), 1); if(fread(&numGhostPartitions, sizeof(std::size_t), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&numGhostPartitions, sizeof(std::size_t), 1); } else { if(fscanf(fp, "%lu %d %lu", &elmTag, &partNum, &numGhostPartitions) != 3) { return false; } if(!fgets(str, sizeof(str), fp)) { return false; } } MElement *elm = model->getMeshElementByTag(elmTag); if(!elm) { Msg::Error("No element with tag %lu", elmTag); continue; } for(std::size_t j = 0; j < numGhostPartitions; j++) { int ghostPartition = 0; if(binary) { if(fread(&ghostPartition, sizeof(int), 1, fp) != 1) { return false; } if(swap) SwapBytes((char *)&ghostPartition, sizeof(int), 1); } else { if(j == numGhostPartitions - 1) { if(sscanf(str, "%d", &ghostPartition) != 1) { return false; } } else { if(sscanf(str, "%d %[0-9- ]", &ghostPartition, str) != 2) { return false; } } } ghostCells.insert(std::pair<std::pair<MElement *, int>, int>( std::pair<MElement *, int>(elm, partNum), ghostPartition)); } } std::vector<GEntity *> ghostEntities(model->getNumPartitions() + 1, 0); std::vector<GEntity *> entities; model->getEntities(entities); for(std::size_t i = 0; i < entities.size(); i++) { GEntity *ge = entities[i]; int partNum = -1; if(ge->geomType() == GEntity::GhostCurve) partNum = static_cast<ghostEdge *>(ge)->getPartition(); else if(ge->geomType() == GEntity::GhostSurface) partNum = static_cast<ghostFace *>(ge)->getPartition(); else if(ge->geomType() == GEntity::GhostVolume) partNum = static_cast<ghostRegion *>(ge)->getPartition(); if(partNum >= 0 && partNum < (int)ghostEntities.size()) ghostEntities[partNum] = ge; } for(std::multimap<std::pair<MElement *, int>, int>::iterator it = ghostCells.begin(); it != ghostCells.end(); ++it) { if(it->second >= (int)ghostEntities.size()) { Msg::Error("Invalid partition %d in ghost elements", it->second); return false; } GEntity *ge = ghostEntities[it->second]; if(!ge) { Msg::Warning("Missing ghost entity on partition %d", it->second); } else if(ge->geomType() == GEntity::GhostCurve) { static_cast<ghostEdge *>(ge)->addElement( it->first.first->getType(), it->first.first, it->first.second); } else if(ge->geomType() == GEntity::GhostSurface) { static_cast<ghostFace *>(ge)->addElement( it->first.first->getType(), it->first.first, it->first.second); } else if(ge->geomType() == GEntity::GhostVolume) { static_cast<ghostRegion *>(ge)->addElement( it->first.first->getType(), it->first.first, it->first.second); } } return true; } static bool readMSH4Parametrizations(GModel *const model, FILE *fp, bool binary) { std::size_t nParamE, nParamF; if(binary){ if(fread(&nParamE, sizeof(std::size_t), 1, fp) != 1) { return false; } if(fread(&nParamF, sizeof(std::size_t), 1, fp) != 1) { return false; } } else{ if(fscanf(fp, "%lu %lu", &nParamE, &nParamF) != 2) { return false; } } for(std::size_t edge = 0; edge < nParamE; edge++) { int tag; if(binary){ if(fread(&tag, sizeof(int), 1, fp) != 1) { return false; } } else{ if(fscanf(fp, "%d", &tag) != 1) { return false; } } GEdge *ge = model->getEdgeByTag(tag); if(ge) { discreteEdge *de = dynamic_cast<discreteEdge *>(ge); if(de) { if(!de->readParametrization(fp, binary)) return false; } } } for(std::size_t face = 0; face < nParamF; face++) { int tag; if(binary){ if(fread(&tag, sizeof(int), 1, fp) != 1) { return false; } } else{ if(fscanf(fp, "%d", &tag) != 1) { return false; } } GFace *gf = model->getFaceByTag(tag); if(gf) { discreteFace *df = dynamic_cast<discreteFace *>(gf); if(df) { if(!df->readParametrization(fp, binary)) return false; } } } return true; } int GModel::_readMSH4(const std::string &name) { bool partitioned = false; FILE *fp = Fopen(name.c_str(), "rb"); if(!fp) { Msg::Error("Unable to open file '%s'", name.c_str()); return 0; } char str[1024] = "x"; double version = 1.0; bool binary = false, swap = false, postpro = false; while(1) { while(str[0] != '$') { if(!fgets(str, sizeof(str), fp) || feof(fp)) break; } std::string sectionName(&str[1]); std::string endSectionName = "End" + sectionName; if(feof(fp)) break; if(!strncmp(&str[1], "MeshFormat", 10)) { if(!fgets(str, sizeof(str), fp) || feof(fp)) { fclose(fp); return 0; } int format; std::size_t size; if(sscanf(str, "%lf %d %lu", &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"); } } if(binary && size != sizeof(std::size_t)) { Msg::Error("Binary file has sizeof(size_t) = %d, not matching " "machine sizeof(size_t) = %d", size, sizeof(std::size_t)); return false; } if(binary && version < 4.1) { Msg::Error("Can only read MSH 4.0 format in ASCII mode"); return false; } } else if(!strncmp(&str[1], "PhysicalNames", 13)) { if(!fgets(str, sizeof(str), fp) || feof(fp)) { fclose(fp); return 0; } int numPhysicalNames = 0; if(sscanf(str, "%d", &numPhysicalNames) != 1) { fclose(fp); return 0; } std::vector<GModel::piter> iterators; getInnerPhysicalNamesIterators(iterators); for(int i = 0; i < numPhysicalNames; i++) { int dim = 0, tag = 0; if(fscanf(fp, "%d %d", &dim, &tag) != 2) { fclose(fp); return 0; } char name[256]; if(!fgets(name, sizeof(name), fp)) { fclose(fp); return 0; } std::string physicalName = ExtractDoubleQuotedString(name, 256); if(physicalName.size()) iterators[dim] = setPhysicalName(iterators[dim], physicalName, dim, tag); } } else if(!strncmp(&str[1], "Entities", 8)) { if(!readMSH4Entities(this, fp, false, binary, swap, version)) { Msg::Error("Could not read entities"); fclose(fp); return 0; } } else if(!strncmp(&str[1], "PartitionedEntities", 19)) { if(!readMSH4Entities(this, fp, true, binary, swap, version)) { Msg::Error("Could not read partitioned entities"); fclose(fp); return 0; } partitioned = true; } else if(!strncmp(&str[1], "Nodes", 5)) { _vertexVectorCache.clear(); _vertexMapCache.clear(); Msg::ResetProgressMeter(); bool dense = false; std::size_t totalNumNodes = 0, maxNodeNum; std::pair<std::size_t, MVertex *> *vertexCache = readMSH4Nodes( this, fp, binary, dense, totalNumNodes, maxNodeNum, swap, version); if(!vertexCache) { Msg::Error("Could not read nodes"); fclose(fp); return false; } if(dense) { _vertexVectorCache.resize(maxNodeNum + 1, 0); for(std::size_t i = 0; i < totalNumNodes; i++) { if(!_vertexVectorCache[vertexCache[i].first]) { _vertexVectorCache[vertexCache[i].first] = vertexCache[i].second; } else { Msg::Info("Skipping duplicate node %d", vertexCache[i].first); } } } else { for(std::size_t i = 0; i < totalNumNodes; i++) { if(_vertexMapCache.count(vertexCache[i].first) == 0) { _vertexMapCache[vertexCache[i].first] = vertexCache[i].second; } else { Msg::Info("Skipping duplicate node %d", vertexCache[i].first); } } } delete[] vertexCache; } else if(!strncmp(&str[1], "Elements", 8)) { Msg::ResetProgressMeter(); bool dense = false; std::size_t totalNumElements = 0, maxElementNum = 0; std::pair<std::size_t, MElement *> *elementCache = readMSH4Elements(this, fp, binary, dense, totalNumElements, maxElementNum, swap, version); if(!elementCache) { Msg::Error("Could not read elements"); fclose(fp); return 0; } if(dense) { _elementVectorCache.resize(maxElementNum + 1, (MElement *)0); for(std::size_t i = 0; i < totalNumElements; i++) { if(!_elementVectorCache[elementCache[i].first]) { _elementVectorCache[elementCache[i].first] = elementCache[i].second; } else { Msg::Info("Skipping duplicate element %d", elementCache[i].first); } } } else { for(std::size_t i = 0; i < totalNumElements; i++) { if(_elementMapCache.count(elementCache[i].first) == 0) { _elementMapCache[elementCache[i].first] = elementCache[i].second; } else { Msg::Info("Skipping duplicate element %d", elementCache[i].first); } } } delete[] elementCache; } else if(!strncmp(&str[1], "Periodic", 8)) { if(!readMSH4PeriodicNodes(this, fp, binary, swap, version)) { Msg::Error("Could not read periodic section"); fclose(fp); return 0; } } else if(!strncmp(&str[1], "GhostElements", 13)) { if(!readMSH4GhostElements(this, fp, binary, swap)) { Msg::Error("Could not read ghost elements"); fclose(fp); return 0; } } else if(!strncmp(&str[1], "Parametrizations", 16)) { if(!readMSH4Parametrizations(this, fp, binary)) { Msg::Error("Could not read parametrizations"); fclose(fp); return 0; } } else if(!strncmp(&str[1], "NodeData", 8) || !strncmp(&str[1], "ElementData", 11) || !strncmp(&str[1], "ElementNodeData", 15)) { postpro = true; break; } while(strncmp(&str[1], endSectionName.c_str(), endSectionName.size())) { if(!fgets(str, sizeof(str), fp) || feof(fp)) { break; } } str[0] = 'a'; } fclose(fp); if(partitioned) { // This part is added to ensure the compatibility between the new // partitioning and the old one. std::vector<GEntity *> entities; getEntities(entities); for(std::size_t i = 0; i < entities.size(); i++) { if(entities[i]->geomType() == GEntity::PartitionPoint) { partitionVertex *pv = static_cast<partitionVertex *>(entities[i]); if(pv->numPartitions() == 1) { const int part = pv->getPartition(0); for(std::size_t j = 0; j < pv->getNumMeshElements(); j++) { pv->getMeshElement(j)->setPartition(part); } } } else if(entities[i]->geomType() == GEntity::PartitionCurve) { partitionEdge *pe = static_cast<partitionEdge *>(entities[i]); if(pe->numPartitions() == 1) { const int part = pe->getPartition(0); for(std::size_t j = 0; j < pe->getNumMeshElements(); j++) { pe->getMeshElement(j)->setPartition(part); } } } else if(entities[i]->geomType() == GEntity::PartitionSurface) { partitionFace *pf = static_cast<partitionFace *>(entities[i]); if(pf->numPartitions() == 1) { const int part = pf->getPartition(0); for(std::size_t j = 0; j < pf->getNumMeshElements(); j++) { pf->getMeshElement(j)->setPartition(part); } } } else if(entities[i]->geomType() == GEntity::PartitionVolume) { partitionRegion *pr = static_cast<partitionRegion *>(entities[i]); if(pr->numPartitions() == 1) { const int part = pr->getPartition(0); for(std::size_t j = 0; j < pr->getNumMeshElements(); j++) { pr->getMeshElement(j)->setPartition(part); } } } } } return postpro ? 2 : 1; } static void writeMSH4Physicals(FILE *fp, GEntity *const entity, bool binary) { if(binary) { std::vector<int> phys = entity->getPhysicalEntities(); std::size_t phySize = phys.size(); fwrite(&phySize, sizeof(std::size_t), 1, fp); for(std::size_t i = 0; i < phys.size(); i++) { int phy = phys[i]; fwrite(&phy, sizeof(int), 1, fp); } } else { std::vector<int> phys = entity->getPhysicalEntities(); fprintf(fp, "%lu", phys.size()); for(std::size_t i = 0; i < phys.size(); i++) { fprintf(fp, " %d", phys[i]); } fprintf(fp, " "); } } static void writeMSH4BoundingBox(SBoundingBox3d boundBox, FILE *fp, double scalingFactor, bool binary, int dim, double version) { double bb[6] = {0., 0., 0., 0., 0., 0.}; if(!boundBox.empty()) { boundBox *= scalingFactor; bb[0] = boundBox.min().x(); bb[1] = boundBox.min().y(); bb[2] = boundBox.min().z(); bb[3] = boundBox.max().x(); bb[4] = boundBox.max().y(); bb[5] = boundBox.max().z(); } if(binary) { fwrite(bb, sizeof(double), (dim > 0) ? 6 : 3, fp); } else { std::size_t n = (version < 4.1 || dim > 0) ? 6 : 3; for(std::size_t i = 0; i < n; i++) fprintf(fp, "%.16g ", bb[i]); } } static void writeMSH4Entities(GModel *const model, FILE *fp, bool partition, bool binary, double scalingFactor, double version) { std::set<GEntity *, GEntityLessThan> ghost; std::set<GRegion *, GEntityLessThan> regions; std::set<GFace *, GEntityLessThan> faces; std::set<GEdge *, GEntityLessThan> edges; std::set<GVertex *, GEntityLessThan> vertices; if(partition) { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) { if((*it)->geomType() == GEntity::PartitionPoint) vertices.insert(*it); } for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { if((*it)->geomType() == GEntity::PartitionCurve) edges.insert(*it); if((*it)->geomType() == GEntity::GhostCurve) ghost.insert(*it); } for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) { if((*it)->geomType() == GEntity::PartitionSurface) faces.insert(*it); if((*it)->geomType() == GEntity::GhostSurface) ghost.insert(*it); } for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) { if((*it)->geomType() == GEntity::PartitionVolume) regions.insert(*it); if((*it)->geomType() == GEntity::GhostVolume) ghost.insert(*it); } } else { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) if((*it)->geomType() != GEntity::PartitionPoint) vertices.insert(*it); for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) if((*it)->geomType() != GEntity::PartitionCurve && (*it)->geomType() != GEntity::GhostCurve) edges.insert(*it); for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) if((*it)->geomType() != GEntity::PartitionSurface && (*it)->geomType() != GEntity::GhostSurface) faces.insert(*it); for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) if((*it)->geomType() != GEntity::PartitionVolume && (*it)->geomType() != GEntity::GhostVolume) regions.insert(*it); } if(binary) { if(partition) { std::size_t numPartitions = model->getNumPartitions(); fwrite(&numPartitions, sizeof(std::size_t), 1, fp); // write the ghostentities' tag std::size_t ghostSize = ghost.size(); std::vector<int> tags; if(ghostSize) { tags.resize(2 * ghostSize); int index = 0; for(std::set<GEntity *, GEntityLessThan>::iterator it = ghost.begin(); it != ghost.end(); ++it) { if((*it)->geomType() == GEntity::GhostCurve) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostEdge *>(*it)->getPartition(); } else if((*it)->geomType() == GEntity::GhostSurface) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostFace *>(*it)->getPartition(); } else if((*it)->geomType() == GEntity::GhostVolume) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostRegion *>(*it)->getPartition(); } index++; } } fwrite(&ghostSize, sizeof(std::size_t), 1, fp); if(ghostSize) { fwrite(&tags[0], sizeof(int), 2 * ghostSize, fp); } } std::size_t verticesSize = vertices.size(); std::size_t edgesSize = edges.size(); std::size_t facesSize = faces.size(); std::size_t regionsSize = regions.size(); fwrite(&verticesSize, sizeof(std::size_t), 1, fp); fwrite(&edgesSize, sizeof(std::size_t), 1, fp); fwrite(&facesSize, sizeof(std::size_t), 1, fp); fwrite(®ionsSize, sizeof(std::size_t), 1, fp); for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it) { int entityTag = (*it)->tag(); fwrite(&entityTag, sizeof(int), 1, fp); if(partition) { partitionVertex *pv = static_cast<partitionVertex *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pv->getParentEntity()) { parentEntityDim = pv->getParentEntity()->dim(); parentEntityTag = pv->getParentEntity()->tag(); } fwrite(&parentEntityDim, sizeof(int), 1, fp); fwrite(&parentEntityTag, sizeof(int), 1, fp); std::vector<int> partitions(pv->getPartitions().begin(), pv->getPartitions().end()); // FIXME std::size_t numPart = partitions.size(); fwrite(&numPart, sizeof(std::size_t), 1, fp); fwrite(&partitions[0], sizeof(int), partitions.size(), fp); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 0, version); writeMSH4Physicals(fp, *it, binary); } for(GModel::eiter it = edges.begin(); it != edges.end(); ++it) { std::vector<GVertex *> vertices; std::vector<int> ori; if((*it)->getBeginVertex()) { vertices.push_back((*it)->getBeginVertex()); ori.push_back(1); } if((*it)->getEndVertex()) { // convention that the end vertex is negative vertices.push_back((*it)->getEndVertex()); ori.push_back(-1); } std::size_t verticesSize = vertices.size(); int entityTag = (*it)->tag(); fwrite(&entityTag, sizeof(int), 1, fp); if(partition) { partitionEdge *pe = static_cast<partitionEdge *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pe->getParentEntity()) { parentEntityDim = pe->getParentEntity()->dim(); parentEntityTag = pe->getParentEntity()->tag(); } fwrite(&parentEntityDim, sizeof(int), 1, fp); fwrite(&parentEntityTag, sizeof(int), 1, fp); std::vector<int> partitions(pe->getPartitions().begin(), pe->getPartitions().end()); // FIXME std::size_t numPart = partitions.size(); fwrite(&numPart, sizeof(std::size_t), 1, fp); fwrite(&partitions[0], sizeof(int), partitions.size(), fp); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 1, version); writeMSH4Physicals(fp, *it, binary); fwrite(&verticesSize, sizeof(std::size_t), 1, fp); int oriI = 0; for(std::vector<GVertex *>::const_iterator itv = vertices.begin(); itv != vertices.end(); itv++) { int brepTag = ori[oriI] * (*itv)->tag(); fwrite(&brepTag, sizeof(int), 1, fp); oriI++; } } for(GModel::fiter it = faces.begin(); it != faces.end(); ++it) { std::vector<GEdge *> const &edges = (*it)->edges(); std::vector<int> const &ori = (*it)->edgeOrientations(); std::size_t edgesSize = edges.size(); int entityTag = (*it)->tag(); fwrite(&entityTag, sizeof(int), 1, fp); if(partition) { partitionFace *pf = static_cast<partitionFace *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pf->getParentEntity()) { parentEntityDim = pf->getParentEntity()->dim(); parentEntityTag = pf->getParentEntity()->tag(); } fwrite(&parentEntityDim, sizeof(int), 1, fp); fwrite(&parentEntityTag, sizeof(int), 1, fp); std::vector<int> partitions(pf->getPartitions().begin(), pf->getPartitions().end()); // FIXME std::size_t numPart = partitions.size(); fwrite(&numPart, sizeof(std::size_t), 1, fp); fwrite(&partitions[0], sizeof(int), partitions.size(), fp); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 2, version); writeMSH4Physicals(fp, *it, binary); fwrite(&edgesSize, sizeof(std::size_t), 1, fp); std::vector<int> tags, signs; for(std::vector<GEdge *>::const_iterator ite = edges.begin(); ite != edges.end(); ite++) tags.push_back((*ite)->tag()); signs.insert(signs.end(), ori.begin(), ori.end()); if(tags.size() == signs.size()) { for(std::size_t i = 0; i < tags.size(); i++) tags[i] *= (signs[i] > 0 ? 1 : -1); } for(std::size_t i = 0; i < tags.size(); i++) { int brepTag = tags[i]; fwrite(&brepTag, sizeof(int), 1, fp); } } for(GModel::riter it = regions.begin(); it != regions.end(); ++it) { std::vector<GFace *> faces = (*it)->faces(); std::vector<int> const &ori = (*it)->faceOrientations(); std::size_t facesSize = faces.size(); int entityTag = (*it)->tag(); fwrite(&entityTag, sizeof(int), 1, fp); if(partition) { partitionRegion *pr = static_cast<partitionRegion *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pr->getParentEntity()) { parentEntityDim = pr->getParentEntity()->dim(); parentEntityTag = pr->getParentEntity()->tag(); } fwrite(&parentEntityDim, sizeof(int), 1, fp); fwrite(&parentEntityTag, sizeof(int), 1, fp); std::vector<int> partitions(pr->getPartitions().begin(), pr->getPartitions().end()); // FIXME std::size_t numPart = partitions.size(); fwrite(&numPart, sizeof(std::size_t), 1, fp); fwrite(&partitions[0], sizeof(int), partitions.size(), fp); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 3, version); writeMSH4Physicals(fp, *it, binary); fwrite(&facesSize, sizeof(std::size_t), 1, fp); std::vector<int> tags, signs; for(std::vector<GFace *>::iterator itf = faces.begin(); itf != faces.end(); itf++) tags.push_back((*itf)->tag()); for(std::vector<int>::const_iterator itf = ori.begin(); itf != ori.end(); itf++) signs.push_back(*itf); if(tags.size() == signs.size()) { for(std::size_t i = 0; i < tags.size(); i++) tags[i] *= (signs[i] > 0 ? 1 : -1); } for(std::size_t i = 0; i < tags.size(); i++) { int brepTag = tags[i]; fwrite(&brepTag, sizeof(int), 1, fp); } } fprintf(fp, "\n"); } else { if(partition) { fprintf(fp, "%lu\n", model->getNumPartitions()); // write the ghostentities' tag std::size_t ghostSize = ghost.size(); std::vector<int> tags; if(ghostSize) { tags.resize(2 * ghostSize); int index = 0; for(std::set<GEntity *, GEntityLessThan>::iterator it = ghost.begin(); it != ghost.end(); ++it) { if((*it)->geomType() == GEntity::GhostCurve) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostEdge *>(*it)->getPartition(); } else if((*it)->geomType() == GEntity::GhostSurface) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostFace *>(*it)->getPartition(); } else if((*it)->geomType() == GEntity::GhostVolume) { tags[index] = (*it)->tag(); tags[++index] = static_cast<ghostRegion *>(*it)->getPartition(); } index++; } } fprintf(fp, "%lu\n", ghostSize); if(ghostSize) { for(std::size_t i = 0; i < 2 * ghostSize; i += 2) { fprintf(fp, "%d %d\n", tags[i], tags[i + 1]); } } } fprintf(fp, "%lu %lu %lu %lu\n", vertices.size(), edges.size(), faces.size(), regions.size()); for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it) { fprintf(fp, "%d ", (*it)->tag()); if(partition) { partitionVertex *pv = static_cast<partitionVertex *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pv->getParentEntity()) { parentEntityDim = pv->getParentEntity()->dim(); parentEntityTag = pv->getParentEntity()->tag(); } fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag); std::vector<int> partitions(pv->getPartitions().begin(), pv->getPartitions().end()); // FIXME fprintf(fp, "%lu ", partitions.size()); for(std::size_t i = 0; i < partitions.size(); i++) fprintf(fp, "%d ", partitions[i]); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 0, version); writeMSH4Physicals(fp, *it, binary); fprintf(fp, "\n"); } for(GModel::eiter it = edges.begin(); it != edges.end(); ++it) { std::vector<GVertex *> vertices; std::vector<int> ori; if((*it)->getBeginVertex()) { vertices.push_back((*it)->getBeginVertex()); ori.push_back(1); } if((*it)->getEndVertex()) { // I use the convention that the end vertex is // negative vertices.push_back((*it)->getEndVertex()); ori.push_back(-1); } fprintf(fp, "%d ", (*it)->tag()); if(partition) { partitionEdge *pe = static_cast<partitionEdge *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pe->getParentEntity()) { parentEntityDim = pe->getParentEntity()->dim(); parentEntityTag = pe->getParentEntity()->tag(); } fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag); std::vector<int> partitions(pe->getPartitions().begin(), pe->getPartitions().end()); // FIXME fprintf(fp, "%lu ", partitions.size()); for(std::size_t i = 0; i < partitions.size(); i++) fprintf(fp, "%d ", partitions[i]); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 1, version); writeMSH4Physicals(fp, *it, binary); fprintf(fp, "%lu ", vertices.size()); int oriI = 0; for(std::vector<GVertex *>::iterator itv = vertices.begin(); itv != vertices.end(); itv++) { fprintf(fp, "%d ", ori[oriI] * (*itv)->tag()); oriI++; } fprintf(fp, "\n"); } for(GModel::fiter it = faces.begin(); it != faces.end(); ++it) { std::vector<GEdge *> const &edges = (*it)->edges(); std::vector<int> const &ori = (*it)->edgeOrientations(); fprintf(fp, "%d ", (*it)->tag()); if(partition) { partitionFace *pf = static_cast<partitionFace *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pf->getParentEntity()) { parentEntityDim = pf->getParentEntity()->dim(); parentEntityTag = pf->getParentEntity()->tag(); } fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag); std::vector<int> partitions(pf->getPartitions().begin(), pf->getPartitions().end()); // FIXME fprintf(fp, "%lu ", partitions.size()); for(std::size_t i = 0; i < partitions.size(); i++) fprintf(fp, "%d ", partitions[i]); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 2, version); writeMSH4Physicals(fp, *it, binary); fprintf(fp, "%lu ", edges.size()); std::vector<int> tags, signs; for(std::vector<GEdge *>::const_iterator ite = edges.begin(); ite != edges.end(); ite++) tags.push_back((*ite)->tag()); for(std::vector<int>::const_iterator ite = ori.begin(); ite != ori.end(); ite++) signs.push_back(*ite); if(tags.size() == signs.size()) { for(std::size_t i = 0; i < tags.size(); i++) tags[i] *= (signs[i] > 0 ? 1 : -1); } for(std::size_t i = 0; i < tags.size(); i++) fprintf(fp, "%d ", tags[i]); fprintf(fp, "\n"); } for(GModel::riter it = regions.begin(); it != regions.end(); ++it) { std::vector<GFace *> const &faces = (*it)->faces(); std::vector<int> const &ori = (*it)->faceOrientations(); fprintf(fp, "%d ", (*it)->tag()); if(partition) { partitionRegion *pr = static_cast<partitionRegion *>(*it); int parentEntityDim = 0, parentEntityTag = 0; if(pr->getParentEntity()) { parentEntityDim = pr->getParentEntity()->dim(); parentEntityTag = pr->getParentEntity()->tag(); } fprintf(fp, "%d %d ", parentEntityDim, parentEntityTag); std::vector<int> partitions(pr->getPartitions().begin(), pr->getPartitions().end()); // FIXME fprintf(fp, "%lu ", partitions.size()); for(std::size_t i = 0; i < partitions.size(); i++) fprintf(fp, "%d ", partitions[i]); } writeMSH4BoundingBox((*it)->bounds(), fp, scalingFactor, binary, 3, version); writeMSH4Physicals(fp, *it, binary); fprintf(fp, "%lu ", faces.size()); // TODO C++11 std::transform or similiar std::vector<int> tags; tags.reserve(faces.size()); for(std::vector<GFace *>::const_iterator itf = faces.begin(); itf != faces.end(); itf++) tags.push_back((*itf)->tag()); std::vector<int> signs(ori.begin(), ori.end()); if(tags.size() == signs.size()) { for(std::size_t i = 0; i < tags.size(); i++) tags[i] *= (signs[i] > 0 ? 1 : -1); } for(std::size_t i = 0; i < tags.size(); i++) fprintf(fp, "%d ", tags[i]); fprintf(fp, "\n"); } } } static void writeMSH4EntityNodes(GEntity *ge, FILE *fp, bool binary, int saveParametric, double scalingFactor, double version) { int parametric = saveParametric; if(ge->dim() != 1 && ge->dim() != 2) parametric = 0; // Gmsh only stores parametric coordinates for dim 1 and 2 if(binary) { int entityDim = ge->dim(); int entityTag = ge->tag(); std::size_t numVerts = ge->getNumMeshVertices(); fwrite(&entityDim, sizeof(int), 1, fp); fwrite(&entityTag, sizeof(int), 1, fp); fwrite(¶metric, sizeof(int), 1, fp); fwrite(&numVerts, sizeof(std::size_t), 1, fp); } else { fprintf(fp, "%d %d %d %lu\n", (version >= 4.1) ? ge->dim() : ge->tag(), (version >= 4.1) ? ge->tag() : ge->dim(), parametric, ge->getNumMeshVertices()); } std::size_t N = ge->getNumMeshVertices(); std::size_t n = 3; if(parametric) n += ge->dim(); if(binary) { std::vector<size_t> tags(N); for(std::size_t i = 0; i < N; i++) tags[i] = ge->getMeshVertex(i)->getNum(); fwrite(&tags[0], sizeof(std::size_t), N, fp); std::vector<double> coord(n * N); std::size_t j = 0; for(std::size_t i = 0; i < N; i++) { MVertex *mv = ge->getMeshVertex(i); coord[j++] = mv->x() * scalingFactor; coord[j++] = mv->y() * scalingFactor; coord[j++] = mv->z() * scalingFactor; if(n >= 4) mv->getParameter(0, coord[j++]); if(n == 5) mv->getParameter(1, coord[j++]); } fwrite(&coord[0], sizeof(double), n * N, fp); } else { if(version >= 4.1) { for(std::size_t i = 0; i < N; i++) fprintf(fp, "%lu\n", ge->getMeshVertex(i)->getNum()); } for(std::size_t i = 0; i < N; i++) { MVertex *mv = ge->getMeshVertex(i); double x = mv->x() * scalingFactor; double y = mv->y() * scalingFactor; double z = mv->z() * scalingFactor; if(version < 4.1) fprintf(fp, "%lu ", mv->getNum()); if(n == 5) { double u, v; mv->getParameter(0, u); mv->getParameter(1, v); fprintf(fp, "%.16g %.16g %.16g %.16g %.16g\n", x, y, z, u, v); } else if(n == 4) { double u; mv->getParameter(0, u); fprintf(fp, "%.16g %.16g %.16g %.16g\n", x, y, z, u); } else { fprintf(fp, "%.16g %.16g %.16g\n", x, y, z); } } } } static std::size_t getAdditionalEntities(std::set<GRegion *, GEntityLessThan> ®ions, std::set<GFace *, GEntityLessThan> &faces, std::set<GEdge *, GEntityLessThan> &edges, std::set<GVertex *, GEntityLessThan> &vertices) { std::size_t numVertices = 0; for(std::set<GVertex *, GEntityLessThan>::iterator it = vertices.begin(); it != vertices.end(); ++it) { numVertices += (*it)->getNumMeshVertices(); } for(std::set<GEdge *, GEntityLessThan>::iterator it = edges.begin(); it != edges.end(); ++it) { numVertices += (*it)->getNumMeshVertices(); for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) { for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++) { if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) { GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat(); switch(entity->dim()) { case 0: if(vertices.find(static_cast<GVertex *>(entity)) == vertices.end()) { vertices.insert(static_cast<GVertex *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 1: if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) { edges.insert(static_cast<GEdge *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 2: if(faces.find(static_cast<GFace *>(entity)) == faces.end()) { faces.insert(static_cast<GFace *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 3: if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) { regions.insert(static_cast<GRegion *>(entity)); numVertices += entity->getNumMeshVertices(); } break; default: break; } } } } } for(std::set<GFace *, GEntityLessThan>::iterator it = faces.begin(); it != faces.end(); ++it) { numVertices += (*it)->getNumMeshVertices(); for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) { for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++) { if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) { GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat(); switch(entity->dim()) { case 0: if(vertices.find(static_cast<GVertex *>(entity)) == vertices.end()) { vertices.insert(static_cast<GVertex *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 1: if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) { edges.insert(static_cast<GEdge *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 2: if(faces.find(static_cast<GFace *>(entity)) == faces.end()) { faces.insert(static_cast<GFace *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 3: if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) { regions.insert(static_cast<GRegion *>(entity)); numVertices += entity->getNumMeshVertices(); } break; default: break; } } } } } for(std::set<GRegion *, GEntityLessThan>::iterator it = regions.begin(); it != regions.end(); ++it) { numVertices += (*it)->getNumMeshVertices(); for(std::size_t i = 0; i < (*it)->getNumMeshElements(); i++) { for(std::size_t j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++) { if((*it)->getMeshElement(i)->getVertex(j)->onWhat() != (*it)) { GEntity *entity = (*it)->getMeshElement(i)->getVertex(j)->onWhat(); switch(entity->dim()) { case 0: if(vertices.find(static_cast<GVertex *>(entity)) == vertices.end()) { vertices.insert(static_cast<GVertex *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 1: if(edges.find(static_cast<GEdge *>(entity)) == edges.end()) { edges.insert(static_cast<GEdge *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 2: if(faces.find(static_cast<GFace *>(entity)) == faces.end()) { faces.insert(static_cast<GFace *>(entity)); numVertices += entity->getNumMeshVertices(); } break; case 3: if(regions.find(static_cast<GRegion *>(entity)) == regions.end()) { regions.insert(static_cast<GRegion *>(entity)); numVertices += entity->getNumMeshVertices(); } break; default: break; } } } } } return numVertices; } static std::size_t getEntitiesForNodes(GModel *const model, bool partitioned, bool saveAll, std::set<GRegion *, GEntityLessThan> ®ions, std::set<GFace *, GEntityLessThan> &faces, std::set<GEdge *, GEntityLessThan> &edges, std::set<GVertex *, GEntityLessThan> &vertices) { if(partitioned) { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) if((*it)->geomType() == GEntity::PartitionPoint) vertices.insert(*it); for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { if((*it)->geomType() == GEntity::PartitionCurve) edges.insert(*it); else if((*it)->geomType() == GEntity::GhostCurve) if(static_cast<ghostEdge *>(*it)->saveMesh()) edges.insert(*it); } for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) { if((*it)->geomType() == GEntity::PartitionSurface) faces.insert(*it); else if((*it)->geomType() == GEntity::GhostSurface) if(static_cast<ghostFace *>(*it)->saveMesh()) faces.insert(*it); } for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) { if((*it)->geomType() == GEntity::PartitionVolume) regions.insert(*it); else if((*it)->geomType() == GEntity::GhostVolume) if(static_cast<ghostRegion *>(*it)->saveMesh()) regions.insert(*it); } } else { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) if((*it)->geomType() != GEntity::PartitionPoint && (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0))) vertices.insert(*it); for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) if((*it)->geomType() != GEntity::PartitionCurve && (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) || (*it)->geomType() == GEntity::GhostCurve)) edges.insert(*it); for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) if((*it)->geomType() != GEntity::PartitionSurface && (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) || (*it)->geomType() == GEntity::GhostSurface)) faces.insert(*it); for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) if((*it)->geomType() != GEntity::PartitionVolume && (saveAll || (!saveAll && (*it)->getPhysicalEntities().size() != 0) || (*it)->geomType() == GEntity::GhostVolume)) regions.insert(*it); } std::size_t numVertices = model->getNumMeshVertices(); if(!saveAll && !partitioned) { numVertices = getAdditionalEntities(regions, faces, edges, vertices); } return numVertices; } static void writeMSH4Nodes(GModel *const model, FILE *fp, bool partitioned, bool binary, int saveParametric, double scalingFactor, bool saveAll, double version) { std::set<GRegion *, GEntityLessThan> regions; std::set<GFace *, GEntityLessThan> faces; std::set<GEdge *, GEntityLessThan> edges; std::set<GVertex *, GEntityLessThan> vertices; std::size_t numNodes = getEntitiesForNodes(model, partitioned, saveAll, regions, faces, edges, vertices); std::size_t minTag = std::numeric_limits<std::size_t>::max(), maxTag = 0; for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it) { for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) { minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum()); maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum()); } } for(GModel::eiter it = edges.begin(); it != edges.end(); ++it) { for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) { minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum()); maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum()); } } for(GModel::fiter it = faces.begin(); it != faces.end(); ++it) { for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) { minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum()); maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum()); } } for(GModel::riter it = regions.begin(); it != regions.end(); ++it) { for(std::size_t i = 0; i < (*it)->getNumMeshVertices(); i++) { minTag = std::min(minTag, (*it)->getMeshVertex(i)->getNum()); maxTag = std::max(maxTag, (*it)->getMeshVertex(i)->getNum()); } } if(binary) { std::size_t numSection = vertices.size() + edges.size() + faces.size() + regions.size(); fwrite(&numSection, sizeof(std::size_t), 1, fp); fwrite(&numNodes, sizeof(std::size_t), 1, fp); fwrite(&minTag, sizeof(std::size_t), 1, fp); fwrite(&maxTag, sizeof(std::size_t), 1, fp); } else { if(version >= 4.1) { fprintf(fp, "%lu %lu %lu %lu\n", vertices.size() + edges.size() + faces.size() + regions.size(), numNodes, minTag, maxTag); } else { fprintf(fp, "%lu %lu\n", vertices.size() + edges.size() + faces.size() + regions.size(), numNodes); } } for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it) { writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor, version); } for(GModel::eiter it = edges.begin(); it != edges.end(); ++it) { writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor, version); } for(GModel::fiter it = faces.begin(); it != faces.end(); ++it) { writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor, version); } for(GModel::riter it = regions.begin(); it != regions.end(); ++it) { writeMSH4EntityNodes(*it, fp, binary, saveParametric, scalingFactor, version); } if(binary) fprintf(fp, "\n"); } static void writeMSH4Elements(GModel *const model, FILE *fp, bool partitioned, bool binary, bool saveAll, double version) { std::set<GRegion *, GEntityLessThan> regions; std::set<GFace *, GEntityLessThan> faces; std::set<GEdge *, GEntityLessThan> edges; std::set<GVertex *, GEntityLessThan> vertices; if(partitioned) { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) if((*it)->geomType() == GEntity::PartitionPoint) vertices.insert(*it); for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { if((*it)->geomType() == GEntity::PartitionCurve) edges.insert(*it); else if((*it)->geomType() == GEntity::GhostCurve) if(static_cast<ghostEdge *>(*it)->saveMesh()) edges.insert(*it); } for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) { if((*it)->geomType() == GEntity::PartitionSurface) faces.insert(*it); else if((*it)->geomType() == GEntity::GhostSurface) if(static_cast<ghostFace *>(*it)->saveMesh()) faces.insert(*it); } for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) { if((*it)->geomType() == GEntity::PartitionVolume) regions.insert(*it); else if((*it)->geomType() == GEntity::GhostVolume) if(static_cast<ghostRegion *>(*it)->saveMesh()) regions.insert(*it); } } else { for(GModel::viter it = model->firstVertex(); it != model->lastVertex(); ++it) if((*it)->geomType() != GEntity::PartitionPoint) vertices.insert(*it); for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) if((*it)->geomType() != GEntity::PartitionCurve && (*it)->geomType() != GEntity::GhostCurve) edges.insert(*it); for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) if((*it)->geomType() != GEntity::PartitionSurface && (*it)->geomType() != GEntity::GhostSurface) faces.insert(*it); for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) if((*it)->geomType() != GEntity::PartitionVolume && (*it)->geomType() != GEntity::GhostVolume) regions.insert(*it); } std::map<std::pair<int, int>, std::vector<MElement *> > elementsByType[4]; std::size_t numElements = 0; for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it) { if(!saveAll && (*it)->physicals.size() == 0) continue; numElements += (*it)->points.size(); for(std::size_t i = 0; i < (*it)->points.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->points[i]->getTypeForMSH()); elementsByType[0][p].push_back((*it)->points[i]); } } for(GModel::eiter it = edges.begin(); it != edges.end(); ++it) { if(!saveAll && (*it)->physicals.size() == 0 && (*it)->geomType() != GEntity::GhostCurve) continue; numElements += (*it)->lines.size(); for(std::size_t i = 0; i < (*it)->lines.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->lines[i]->getTypeForMSH()); elementsByType[1][p].push_back((*it)->lines[i]); } } for(GModel::fiter it = faces.begin(); it != faces.end(); ++it) { if(!saveAll && (*it)->physicals.size() == 0 && (*it)->geomType() != GEntity::GhostSurface) continue; numElements += (*it)->triangles.size(); for(std::size_t i = 0; i < (*it)->triangles.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->triangles[i]->getTypeForMSH()); elementsByType[2][p].push_back((*it)->triangles[i]); } numElements += (*it)->quadrangles.size(); for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->quadrangles[i]->getTypeForMSH()); elementsByType[2][p].push_back((*it)->quadrangles[i]); } } for(GModel::riter it = regions.begin(); it != regions.end(); ++it) { if(!saveAll && (*it)->physicals.size() == 0 && (*it)->geomType() != GEntity::GhostVolume) continue; numElements += (*it)->tetrahedra.size(); for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->tetrahedra[i]->getTypeForMSH()); elementsByType[3][p].push_back((*it)->tetrahedra[i]); } numElements += (*it)->hexahedra.size(); for(std::size_t i = 0; i < (*it)->hexahedra.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->hexahedra[i]->getTypeForMSH()); elementsByType[3][p].push_back((*it)->hexahedra[i]); } numElements += (*it)->prisms.size(); for(std::size_t i = 0; i < (*it)->prisms.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->prisms[i]->getTypeForMSH()); elementsByType[3][p].push_back((*it)->prisms[i]); } numElements += (*it)->pyramids.size(); for(std::size_t i = 0; i < (*it)->pyramids.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->pyramids[i]->getTypeForMSH()); elementsByType[3][p].push_back((*it)->pyramids[i]); } numElements += (*it)->trihedra.size(); for(std::size_t i = 0; i < (*it)->trihedra.size(); i++) { std::pair<int, int> p((*it)->tag(), (*it)->trihedra[i]->getTypeForMSH()); elementsByType[3][p].push_back((*it)->trihedra[i]); } } std::size_t numSection = 0; for(int dim = 0; dim <= 3; dim++) numSection += elementsByType[dim].size(); std::size_t minTag = std::numeric_limits<std::size_t>::max(), maxTag = 0; for(int dim = 0; dim <= 3; dim++) { for(std::map<std::pair<int, int>, std::vector<MElement *> >::iterator it = elementsByType[dim].begin(); it != elementsByType[dim].end(); ++it) { for(std::size_t i = 0; i < it->second.size(); i++) { minTag = std::min(minTag, it->second[i]->getNum()); maxTag = std::max(maxTag, it->second[i]->getNum()); } } } if(binary) { fwrite(&numSection, sizeof(std::size_t), 1, fp); fwrite(&numElements, sizeof(std::size_t), 1, fp); fwrite(&minTag, sizeof(std::size_t), 1, fp); fwrite(&maxTag, sizeof(std::size_t), 1, fp); } else { if(version >= 4.1) fprintf(fp, "%lu %lu %lu %lu\n", numSection, numElements, minTag, maxTag); else fprintf(fp, "%lu %lu\n", numSection, numElements); } for(int dim = 0; dim <= 3; dim++) { for(std::map<std::pair<int, int>, std::vector<MElement *> >::iterator it = elementsByType[dim].begin(); it != elementsByType[dim].end(); ++it) { int entityTag = it->first.first; int elmType = it->first.second; std::size_t numElm = it->second.size(); if(binary) { fwrite(&dim, sizeof(int), 1, fp); fwrite(&entityTag, sizeof(int), 1, fp); fwrite(&elmType, sizeof(int), 1, fp); fwrite(&numElm, sizeof(std::size_t), 1, fp); } else { fprintf(fp, "%d %d %d %lu\n", (version >= 4.1) ? dim : entityTag, (version >= 4.1) ? entityTag : dim, elmType, numElm); } std::size_t N = it->second.size(); if(binary) { const int numVertPerElm = MElement::getInfoMSH(elmType); std::size_t n = 1 + numVertPerElm; std::vector<std::size_t> tags(N * n); std::size_t k = 0; for(std::size_t i = 0; i < N; i++) { MElement *e = it->second[i]; tags[k] = e->getNum(); for(int j = 0; j < numVertPerElm; j++) { tags[k + 1 + j] = e->getVertex(j)->getNum(); } k += n; } fwrite(&tags[0], sizeof(std::size_t), N * n, fp); } else { for(std::size_t i = 0; i < N; i++) { MElement *e = it->second[i]; fprintf(fp, "%lu ", e->getNum()); for(std::size_t i = 0; i < e->getNumVertices(); i++) { fprintf(fp, "%lu ", e->getVertex(i)->getNum()); } fprintf(fp, "\n"); } } } } if(binary) fprintf(fp, "\n"); } static void writeMSH4PeriodicNodes(GModel *const model, FILE *fp, bool partitioned, bool binary, double version) { // To avoid saving correspondances bwteen nodes that are not saved (either in // the same file or not at all, e.g. in the partitioned case, or simply if // some physical entities are not defined), we could only apply the code below // to the entities returned by getEntitiesForNodes(). // The current choice is to save everything, and not complain if a node is not // found when reading the file. This should be reevaluated at some point. std::size_t count = 0; std::vector<GEntity *> entities; model->getEntities(entities); for(std::size_t i = 0; i < entities.size(); i++) if(entities[i]->getMeshMaster() != entities[i]) count++; if(!count) return; fprintf(fp, "$Periodic\n"); if(binary) { fwrite(&count, sizeof(std::size_t), 1, fp); } else { fprintf(fp, "%lu\n", count); } for(std::size_t i = 0; i < entities.size(); i++) { GEntity *g_slave = entities[i]; GEntity *g_master = g_slave->getMeshMaster(); if(g_slave != g_master) { if(binary) { int gSlaveDim = g_slave->dim(); int gSlaveTag = g_slave->tag(); int gMasterTag = g_master->tag(); fwrite(&gSlaveDim, sizeof(int), 1, fp); fwrite(&gSlaveTag, sizeof(int), 1, fp); fwrite(&gMasterTag, sizeof(int), 1, fp); if(g_slave->affineTransform.size() == 16) { std::size_t numAffine = 16; fwrite(&numAffine, sizeof(std::size_t), 1, fp); for(int j = 0; j < 16; j++) { double value = g_slave->affineTransform[j]; fwrite(&value, sizeof(double), 1, fp); } } else { std::size_t numAffine = 0; fwrite(&numAffine, sizeof(std::size_t), 1, fp); } std::size_t corrVertSize = g_slave->correspondingVertices.size(); fwrite(&corrVertSize, sizeof(std::size_t), 1, fp); for(std::map<MVertex *, MVertex *>::iterator it = g_slave->correspondingVertices.begin(); it != g_slave->correspondingVertices.end(); ++it) { std::size_t numFirst = it->first->getNum(); std::size_t numSecond = it->second->getNum(); fwrite(&numFirst, sizeof(std::size_t), 1, fp); fwrite(&numSecond, sizeof(std::size_t), 1, fp); } } else { fprintf(fp, "%d %d %d\n", g_slave->dim(), g_slave->tag(), g_master->tag()); if(version >= 4.1) { if(g_slave->affineTransform.size() == 16) { fprintf(fp, "16"); for(int j = 0; j < 16; j++) fprintf(fp, " %.16g", g_slave->affineTransform[j]); fprintf(fp, "\n"); } else { fprintf(fp, "0\n"); } } else { if(g_slave->affineTransform.size() == 16) { fprintf(fp, "Affine"); for(int j = 0; j < 16; j++) fprintf(fp, " %.16g", g_slave->affineTransform[j]); fprintf(fp, "\n"); } } fprintf(fp, "%lu\n", g_slave->correspondingVertices.size()); for(std::map<MVertex *, MVertex *>::iterator it = g_slave->correspondingVertices.begin(); it != g_slave->correspondingVertices.end(); ++it) { fprintf(fp, "%lu %lu\n", it->first->getNum(), it->second->getNum()); } } } } if(binary) fprintf(fp, "\n"); fprintf(fp, "$EndPeriodic\n"); } static void writeMSH4GhostCells(GModel *const model, FILE *fp, bool binary) { std::vector<GEntity *> entities; model->getEntities(entities); std::map<MElement *, std::vector<int> > ghostCells; for(std::size_t i = 0; i < entities.size(); i++) { std::map<MElement *, unsigned int> ghostElements; // FIXME int partition; if(entities[i]->geomType() == GEntity::GhostCurve) { ghostElements = static_cast<ghostEdge *>(entities[i])->getGhostCells(); partition = static_cast<ghostEdge *>(entities[i])->getPartition(); } else if(entities[i]->geomType() == GEntity::GhostSurface) { ghostElements = static_cast<ghostFace *>(entities[i])->getGhostCells(); partition = static_cast<ghostFace *>(entities[i])->getPartition(); } else if(entities[i]->geomType() == GEntity::GhostVolume) { ghostElements = static_cast<ghostRegion *>(entities[i])->getGhostCells(); partition = static_cast<ghostRegion *>(entities[i])->getPartition(); } for(std::map<MElement *, unsigned int>::iterator it = ghostElements.begin(); it != ghostElements.end(); ++it) { if(ghostCells[it->first].size() == 0) ghostCells[it->first].push_back(it->second); ghostCells[it->first].push_back(partition); } } if(ghostCells.size() != 0) { fprintf(fp, "$GhostElements\n"); if(binary) { std::size_t ghostCellsSize = ghostCells.size(); fwrite(&ghostCellsSize, sizeof(std::size_t), 1, fp); for(std::map<MElement *, std::vector<int> >::iterator it = ghostCells.begin(); it != ghostCells.end(); ++it) { std::size_t elmTag = it->first->getNum(); int partNum = it->second[0]; std::size_t numGhostPartitions = it->second.size() - 1; fwrite(&elmTag, sizeof(std::size_t), 1, fp); fwrite(&partNum, sizeof(int), 1, fp); fwrite(&numGhostPartitions, sizeof(std::size_t), 1, fp); for(std::size_t i = 1; i < it->second.size(); i++) { fwrite(&it->second[i], sizeof(int), 1, fp); } } fprintf(fp, "\n"); } else { fprintf(fp, "%ld\n", ghostCells.size()); for(std::map<MElement *, std::vector<int> >::iterator it = ghostCells.begin(); it != ghostCells.end(); ++it) { fprintf(fp, "%lu %d %ld", it->first->getNum(), it->second[0], it->second.size() - 1); for(std::size_t i = 1; i < it->second.size(); i++) { fprintf(fp, " %d", it->second[i]); } fprintf(fp, "\n"); } } fprintf(fp, "$EndGhostElements\n"); } } static void writeMSH4Parametrizations(GModel *const model, FILE *fp, bool binary) { fprintf(fp, "$Parametrizations\n"); std::size_t nParamE = 0, nParamF = 0; for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { discreteEdge *de = dynamic_cast<discreteEdge *>(*it); if(de && de->haveParametrization()) { nParamE++; } } for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) { discreteFace *df = dynamic_cast<discreteFace *>(*it); if(df && df->haveParametrization()) { nParamF++; } } if(binary){ fwrite(&nParamE, sizeof(std::size_t), 1, fp); fwrite(&nParamF, sizeof(std::size_t), 1, fp); } else{ fprintf(fp, "%lu %lu\n", nParamE, nParamF); } for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { discreteEdge *de = dynamic_cast<discreteEdge *>(*it); if(de && de->haveParametrization()) { int t = de->tag(); if(binary) fwrite(&t, sizeof(int), 1, fp); else fprintf(fp, "%d\n", t); de->writeParametrization(fp, binary); } } for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) { discreteFace *df = dynamic_cast<discreteFace *>(*it); if(df && df->haveParametrization()) { int t = df->tag(); if(binary) fwrite(&t, sizeof(int), 1, fp); else fprintf(fp, "%d\n", t); df->writeParametrization(fp, binary); } } if(binary) fprintf(fp, "\n"); fprintf(fp, "$EndParametrizations\n"); } int GModel::_writeMSH4(const std::string &name, double version, bool binary, bool saveAll, bool saveParametric, double scalingFactor, bool append) { FILE *fp = 0; 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; } if(version < 4.1 && binary) { Msg::Error("Can only write MSH 4.0 format in ASCII mode"); return 0; } // if there are no physicals we save all the elements if(noPhysicalGroups()) saveAll = true; // header fprintf(fp, "$MeshFormat\n"); fprintf(fp, "%g %d %lu\n", version, (binary ? 1 : 0), sizeof(std::size_t)); if(binary) { int one = 1; fwrite(&one, sizeof(int), 1, fp); // swapping byte fprintf(fp, "\n"); } fprintf(fp, "$EndMeshFormat\n"); // physicals if(numPhysicalNames() > 0) { 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"); } // entities fprintf(fp, "$Entities\n"); writeMSH4Entities(this, fp, false, binary, scalingFactor, version); fprintf(fp, "$EndEntities\n"); // partitioned entities if(getNumPartitions() > 0) { fprintf(fp, "$PartitionedEntities\n"); writeMSH4Entities(this, fp, true, binary, scalingFactor, version); fprintf(fp, "$EndPartitionedEntities\n"); } // nodes fprintf(fp, "$Nodes\n"); writeMSH4Nodes(this, fp, getNumPartitions() == 0 ? false : true, binary, saveParametric ? 1 : 0, scalingFactor, saveAll, version); fprintf(fp, "$EndNodes\n"); // elements fprintf(fp, "$Elements\n"); writeMSH4Elements(this, fp, getNumPartitions() == 0 ? false : true, binary, saveAll, version); fprintf(fp, "$EndElements\n"); // periodic writeMSH4PeriodicNodes(this, fp, getNumPartitions() == 0 ? false : true, binary, version); // ghostCells writeMSH4GhostCells(this, fp, binary); // parametrizations writeMSH4Parametrizations(this, fp, binary); fclose(fp); return 1; } static void associateVertices(GModel *model) { for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(0); } } (*it)->mesh_vertices.clear(); } for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(0); } } (*it)->mesh_vertices.clear(); } for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(0); } } (*it)->mesh_vertices.clear(); } for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(0); } } (*it)->mesh_vertices.clear(); } for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it); (*it)->mesh_vertices.push_back( (*it)->getMeshElement(j)->getVertex(k)); } } } } for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it); (*it)->mesh_vertices.push_back( (*it)->getMeshElement(j)->getVertex(k)); } } } } for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it); (*it)->mesh_vertices.push_back( (*it)->getMeshElement(j)->getVertex(k)); } } } } for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it) { for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) { for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices(); k++) { if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) { (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it); (*it)->mesh_vertices.push_back( (*it)->getMeshElement(j)->getVertex(k)); } } } } } int GModel::_writePartitionedMSH4(const std::string &baseName, double version, bool binary, bool saveAll, bool saveParametric, double scalingFactor) { std::vector<GEntity *> ghostEntities; std::vector<GEntity *> entities; getEntities(entities); for(std::size_t i = 0; i < entities.size(); i++) { if(entities[i]->geomType() == GEntity::GhostCurve || entities[i]->geomType() == GEntity::GhostSurface || entities[i]->geomType() == GEntity::GhostVolume) { ghostEntities.push_back(entities[i]); } } // Create a temporary model GModel *tmp = new GModel(); tmp->setPhysicalNames(getPhysicalNames()); tmp->setNumPartitions(getNumPartitions()); for(std::size_t i = 1; i <= getNumPartitions(); i++) { std::set<GEntity *> entitiesSet; GEntity *ghostEntity = 0; for(std::size_t j = 0; j < entities.size(); j++) { GEntity *ge = entities[j]; switch(ge->geomType()) { case GEntity::PartitionVolume: { partitionRegion *pr = static_cast<partitionRegion *>(ge); if(std::find(pr->getPartitions().begin(), pr->getPartitions().end(), i) != pr->getPartitions().end()) { tmp->add(pr); if(ghostEntities.size()) entitiesSet.insert(pr); } } break; case GEntity::PartitionSurface: { partitionFace *pf = static_cast<partitionFace *>(ge); if(std::find(pf->getPartitions().begin(), pf->getPartitions().end(), i) != pf->getPartitions().end()) { tmp->add(pf); if(ghostEntities.size()) entitiesSet.insert(pf); } } break; case GEntity::PartitionCurve: { partitionEdge *pe = static_cast<partitionEdge *>(ge); if(std::find(pe->getPartitions().begin(), pe->getPartitions().end(), i) != pe->getPartitions().end()) { tmp->add(pe); if(ghostEntities.size()) entitiesSet.insert(pe); } } break; case GEntity::PartitionPoint: { partitionVertex *pv = static_cast<partitionVertex *>(ge); if(std::find(pv->getPartitions().begin(), pv->getPartitions().end(), i) != pv->getPartitions().end()) { tmp->add(pv); if(ghostEntities.size()) entitiesSet.insert(pv); } } break; case GEntity::GhostCurve: if(i == static_cast<ghostEdge *>(ge)->getPartition()) { static_cast<ghostEdge *>(ge)->saveMesh(true); tmp->add(static_cast<ghostEdge *>(ge)); if(ghostEntities.size()) entitiesSet.insert(ge); ghostEntity = ge; } break; case GEntity::GhostSurface: if(i == static_cast<ghostFace *>(ge)->getPartition()) { static_cast<ghostFace *>(ge)->saveMesh(true); tmp->add(static_cast<ghostFace *>(ge)); if(ghostEntities.size()) entitiesSet.insert(ge); ghostEntity = ge; } break; case GEntity::GhostVolume: if(i == static_cast<ghostRegion *>(ge)->getPartition()) { static_cast<ghostRegion *>(ge)->saveMesh(true); tmp->add(static_cast<ghostRegion *>(ge)); if(ghostEntities.size()) entitiesSet.insert(ge); ghostEntity = ge; } break; default: switch(ge->dim()) { case 0: tmp->add(static_cast<GVertex *>(ge)); break; case 1: tmp->add(static_cast<GEdge *>(ge)); break; case 2: tmp->add(static_cast<GFace *>(ge)); break; case 3: tmp->add(static_cast<GRegion *>(ge)); break; } if(ghostEntities.size()) entitiesSet.insert(ge); break; } } associateVertices(tmp); if(ghostEntity) { for(std::size_t j = 0; j < ghostEntity->getNumMeshElements(); j++) { MElement *e = ghostEntity->getMeshElement(j); for(std::size_t k = 0; k < e->getNumVertices(); k++) { MVertex *v = e->getVertex(k); if(entitiesSet.find(v->onWhat()) == entitiesSet.end()) { v->setEntity(ghostEntity); ghostEntity->addMeshVertex(v); } } } } std::ostringstream sstream; sstream << baseName << "_" << i << ".msh"; if(getNumPartitions() > 100) { if(i % 100 == 1) { Msg::Info("Writing partition %d/%d in file '%s'", i, getNumPartitions(), sstream.str().c_str()); } } else { Msg::Info("Writing partition %d in file '%s'", i, sstream.str().c_str()); } tmp->_writeMSH4(sstream.str(), version, binary, saveAll, saveParametric, scalingFactor, false); tmp->remove(); } delete tmp; associateVertices(this); return 1; } static bool getPhyscialNameInfo(const std::string &name, int &parentPhysicalTag, std::vector<int> &partitions) { if(name[0] != '_') return false; const std::string part = "_part{"; const std::string physical = "_physical{"; size_t firstPart = name.find(part) + part.size(); size_t lastPart = name.find_first_of('}', firstPart); const std::string partString = name.substr(firstPart, lastPart - firstPart); size_t firstPhysical = name.find(physical) + physical.size(); size_t lastPhysical = name.find_first_of('}', firstPhysical); const std::string physicalString = name.substr(firstPhysical, lastPhysical - firstPhysical); std::string number; for(size_t i = 0; i < partString.size(); ++i) { if(partString[i] == ',') { partitions.push_back(atoi(number.c_str())); number.clear(); } else { number += partString[i]; } } partitions.push_back(atoi(number.c_str())); parentPhysicalTag = atoi(physicalString.c_str()); return true; } int GModel::writePartitionedTopology(std::string &name) { Msg::Info("Writing '%s'", name.c_str()); std::vector<std::map<int, std::pair<int, std::vector<int> > > > allParts(4); std::vector<GEntity *> entities; getEntities(entities); for(size_t i = 0; i < entities.size(); i++) { std::vector<int> physicals = entities[i]->getPhysicalEntities(); for(size_t j = 0; j < physicals.size(); ++j) { const std::string phyName = this->getPhysicalName(entities[i]->dim(), physicals[j]); int parentPhysicalTag; std::vector<int> partitions; if(getPhyscialNameInfo(phyName, parentPhysicalTag, partitions)) { allParts[entities[i]->dim()].insert( std::pair<int, std::pair<int, std::vector<int> > >( physicals[j], std::pair<int, std::vector<int> >(parentPhysicalTag, partitions))); } } } FILE *fp = Fopen(name.c_str(), "w"); if(!fp) { Msg::Error("Could not open file '%s'", name.c_str()); return 0; } #if __cplusplus < 201103L char intToChar[20]; #endif fprintf(fp, "Group{\n"); fprintf(fp, " // Part~{dim}~{parentPhysicalTag}~{part1}~{part2}~...\n\n"); std::vector<std::map<int, std::string> > tagToString(4); for(size_t i = 4; i > 0; --i) { fprintf(fp, " // Dim %lu\n", i - 1); for(std::multimap<int, std::pair<int, std::vector<int> > >::iterator it = allParts[i - 1].begin(); it != allParts[i - 1].end(); ++it) { #if __cplusplus >= 201103L std::string partName = "Part~{" + std::to_string(i - 1) + "}~{" + std::to_string(it->second.first) + "}"; #else std::string partName = "Part~{"; sprintf(intToChar, "%lu", i - 1); partName += intToChar; partName += "}~{"; sprintf(intToChar, "%d", it->second.first); partName += intToChar; partName += "}"; #endif fprintf(fp, " Part~{%lu}~{%d}", i - 1, it->second.first); for(size_t j = 0; j < it->second.second.size(); ++j) { #if __cplusplus >= 201103L partName += "~{" + std::to_string(it->second.second[j]) + "}"; #else partName += "~{"; sprintf(intToChar, "%d", it->second.second[j]); partName += intToChar; partName += "}"; #endif fprintf(fp, "~{%d}", it->second.second[j]); } tagToString[i - 1].insert( std::pair<int, std::string>(it->first, partName)); fprintf(fp, " = Region[{%d}];\n", it->first); } fprintf(fp, "\n"); } fprintf(fp, " // Global names\n\n"); std::map<int, std::vector<int> > omegas; std::map<std::pair<int, int>, std::vector<int> > sigmasij; std::map<int, std::vector<int> > sigmas; std::map<int, std::set<int> > neighbors; std::size_t omegaDim = 0; for(size_t i = 4; i > 0; --i) { if(allParts[i - 1].size() != 0) { omegaDim = i - 1; break; } } // omega for(std::multimap<int, std::pair<int, std::vector<int> > >::iterator it = allParts[omegaDim].begin(); it != allParts[omegaDim].end(); ++it) { if(it->second.second.size() == 1) { omegas[it->second.second[0]].push_back(it->first); } } fprintf(fp, " // Omega\n"); for(std::map<int, std::vector<int> >::iterator it = omegas.begin(); it != omegas.end(); ++it) { fprintf(fp, " Omega~{%d} = Region[{", it->first); for(size_t j = 0; j < it->second.size(); ++j) { if(j == 0) fprintf(fp, "%s", tagToString[omegaDim][it->second[j]].c_str()); else fprintf(fp, ", %s", tagToString[omegaDim][it->second[j]].c_str()); } fprintf(fp, "}];\n"); } fprintf(fp, "\n"); if(omegaDim > 0) { // sigma for(std::multimap<int, std::pair<int, std::vector<int> > >::iterator it = allParts[omegaDim - 1].begin(); it != allParts[omegaDim - 1].end(); ++it) { if(it->second.second.size() == 2) { sigmasij[std::pair<int, int>(it->second.second[0], it->second.second[1])] .push_back(it->first); sigmasij[std::pair<int, int>(it->second.second[1], it->second.second[0])] .push_back(it->first); sigmas[it->second.second[0]].push_back(it->first); sigmas[it->second.second[1]].push_back(it->first); } } fprintf(fp, " // Sigma\n"); for(std::map<std::pair<int, int>, std::vector<int> >::iterator it = sigmasij.begin(); it != sigmasij.end(); ++it) { fprintf(fp, " Sigma~{%d}~{%d} = Region[{", it->first.first, it->first.second); for(size_t j = 0; j < it->second.size(); ++j) { if(j == 0) fprintf(fp, "%s", tagToString[omegaDim - 1][it->second[j]].c_str()); else fprintf(fp, ", %s", tagToString[omegaDim - 1][it->second[j]].c_str()); } fprintf(fp, "}];\n"); } fprintf(fp, "\n"); for(std::map<int, std::vector<int> >::iterator it = sigmas.begin(); it != sigmas.end(); ++it) { fprintf(fp, " Sigma~{%d} = Region[{", it->first); for(size_t j = 0; j < it->second.size(); ++j) { if(j == 0) fprintf(fp, "%s", tagToString[omegaDim - 1][it->second[j]].c_str()); else fprintf(fp, ", %s", tagToString[omegaDim - 1][it->second[j]].c_str()); } fprintf(fp, "}];\n"); } fprintf(fp, "\n"); } // D fprintf(fp, " D() = {"); for(size_t i = 1; i <= getNumPartitions(); ++i) { if(i != 1) fprintf(fp, ", "); fprintf(fp, "%lu", i); } fprintf(fp, "};\n"); if(omegaDim > 0) { // D~{i} for(std::multimap<int, std::pair<int, std::vector<int> > >::iterator it = allParts[omegaDim - 1].begin(); it != allParts[omegaDim - 1].end(); ++it) { if(it->second.second.size() == 2) { neighbors[it->second.second[0]].insert(it->second.second[1]); neighbors[it->second.second[1]].insert(it->second.second[0]); } } for(size_t i = 1; i <= getNumPartitions(); ++i) { fprintf(fp, " D~{%lu}() = {", i); for(std::set<int>::iterator it = neighbors[i].begin(); it != neighbors[i].end(); ++it) { if(it != neighbors[i].begin()) fprintf(fp, ", "); fprintf(fp, "%d", *it); } fprintf(fp, "};\n"); } } fprintf(fp, "}\n\n"); fclose(fp); Msg::Info("Done writing '%s'", name.c_str()); return 1; }