Select Git revision
GModelIO_MESH.cpp
-
Christophe Geuzaine authoredChristophe Geuzaine authored
GModelIO_MESH.cpp 11.94 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#include <stdlib.h>
#include <string.h>
#include "GModel.h"
#include "OS.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "Context.h"
static bool getMeshVertices(int num, int *indices, std::vector<MVertex *> &vec,
std::vector<MVertex *> &vertices)
{
for(int i = 0; i < num; i++) {
if(indices[i] < 0 || indices[i] > (int)(vec.size() - 1)) {
Msg::Error("Wrong vertex index %d", indices[i]);
return false;
}
else
vertices.push_back(vec[indices[i]]);
}
return true;
}
int GModel::readMESH(const std::string &name)
{
FILE *fp = Fopen(name.c_str(), "r");
if(!fp) {
Msg::Error("Unable to open file '%s'", name.c_str());
return 0;
}
char buffer[256];
if(!fgets(buffer, sizeof(buffer), fp)) {
fclose(fp);
return 0;
}
char str[256];
int format;
sscanf(buffer, "%s %d", str, &format);
if(format == 3) {
Msg::Error("Medit mesh import only available for ASCII files");
fclose(fp);
return 0;
}
std::vector<MVertex *> vertexVector;
std::map<int, std::vector<MElement *> > elements[5];
while(!feof(fp)) {
if(!fgets(buffer, 256, fp)) break;
if(buffer[0] != '#') { // skip comments and empty lines
str[0] = '\0';
sscanf(buffer, "%s", str);
if(!strncmp(buffer, "Dimension 3", 11)) {
// alternative single-line 'Dimension' field used by CGAL
}
else if(!strcmp(str, "Dimension")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
}
else if(!strcmp(str, "Vertices")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbv;
sscanf(buffer, "%d", &nbv);
Msg::Info("%d vertices", nbv);
vertexVector.resize(nbv);
for(int i = 0; i < nbv; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int dum;
double x, y, z;
sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &dum);
vertexVector[i] = new MVertex(x, y, z);
}
}
else if(!strcmp(str, "Edges")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d edges", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[2], cl;
sscanf(buffer, "%d %d %d", &n[0], &n[1], &cl);
for(int j = 0; j < 2; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(2, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[0][cl].push_back(new MLine(vertices));
}
}
else if(!strcmp(str, "EdgesP2")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d edges", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[3], cl;
sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &cl);
for(int j = 0; j < 3; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(3, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[0][cl].push_back(new MLine3(vertices));
}
}
else if(!strcmp(str, "Triangles")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d triangles", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[3], cl;
sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &cl);
for(int j = 0; j < 3; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(3, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[1][cl].push_back(new MTriangle(vertices));
}
}
else if(!strcmp(str, "TrianglesP2")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d triangles", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[6], cl;
sscanf(buffer, "%d %d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3],
&n[4], &n[5], &cl);
for(int j = 0; j < 6; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(6, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[1][cl].push_back(new MTriangle6(vertices));
}
}
else if(!strcmp(str, "Quadrilaterals")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d quadrangles", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[4], cl;
sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(4, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[2][cl].push_back(new MQuadrangle(vertices));
}
}
else if(!strcmp(str, "Tetrahedra")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d tetrahedra", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[4], cl;
sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(4, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[3][cl].push_back(new MTetrahedron(vertices));
}
}
else if(!strcmp(str, "TetrahedraP2")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d tetrahedra", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[10], cl;
sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d %d", &n[0], &n[1],
&n[2], &n[3], &n[4], &n[5], &n[6], &n[7], &n[9], &n[8], &cl);
for(int j = 0; j < 10; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(10, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[3][cl].push_back(new MTetrahedron10(vertices));
}
}
else if(!strcmp(str, "Hexahedra")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
sscanf(buffer, "%d", &nbe);
Msg::Info("%d hexahedra", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[8], cl;
sscanf(buffer, "%d %d %d %d %d %d %d %d %d", &n[0], &n[1], &n[2],
&n[3], &n[4], &n[5], &n[6], &n[7], &cl);
for(int j = 0; j < 8; j++) n[j]--;
std::vector<MVertex *> vertices;
if(!getMeshVertices(8, n, vertexVector, vertices)) {
fclose(fp);
return 0;
}
elements[4][cl].push_back(new MHexahedron(vertices));
}
}
}
}
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertexVector);
fclose(fp);
return 1;
}
int GModel::writeMESH(const std::string &name, int elementTagType, bool saveAll,
double scalingFactor)
{
FILE *fp = Fopen(name.c_str(), "w");
if(!fp) {
Msg::Error("Unable to open file '%s'", name.c_str());
return 0;
}
if(noPhysicalGroups()) saveAll = true;
int numVertices = indexMeshVertices(saveAll);
fprintf(fp, " MeshVersionFormatted 2\n");
fprintf(fp, " Dimension\n");
fprintf(fp, " 3\n");
fprintf(fp, " Vertices\n");
fprintf(fp, " %d\n", numVertices);
std::vector<GEntity *> entities;
getEntities(entities);
for(std::size_t i = 0; i < entities.size(); i++)
for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeMESH(fp, scalingFactor);
int numEdges = 0, numTriangles = 0, numQuadrangles = 0;
int numTetrahedra = 0, numHexahedra = 0;
for(eiter it = firstEdge(); it != lastEdge(); ++it) {
if(saveAll || (*it)->physicals.size()) {
numEdges += (*it)->lines.size();
}
}
for(fiter it = firstFace(); it != lastFace(); ++it) {
if(saveAll || (*it)->physicals.size()) {
numTriangles += (*it)->triangles.size();
numQuadrangles += (*it)->quadrangles.size();
}
}
for(riter it = firstRegion(); it != lastRegion(); ++it) {
if(saveAll || (*it)->physicals.size()) {
numTetrahedra += (*it)->tetrahedra.size();
numHexahedra += (*it)->hexahedra.size();
}
}
if(numEdges) {
if(CTX::instance()->mesh.order == 2) // FIXME (check getPolynomialOrder())
fprintf(fp, " EdgesP2\n");
else
fprintf(fp, " Edges\n");
fprintf(fp, " %d\n", numEdges);
for(eiter it = firstEdge(); it != lastEdge(); ++it) {
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys) {
for(std::size_t i = 0; i < (*it)->lines.size(); i++)
(*it)->lines[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
}
if(numTriangles) {
if(CTX::instance()->mesh.order == 2) // FIXME (check getPolynomialOrder())
fprintf(fp, " TrianglesP2\n");
else
fprintf(fp, " Triangles\n");
fprintf(fp, " %d\n", numTriangles);
for(fiter it = firstFace(); it != lastFace(); ++it) {
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys) {
for(std::size_t i = 0; i < (*it)->triangles.size(); i++)
(*it)->triangles[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
}
if(numQuadrangles) {
fprintf(fp, " Quadrilaterals\n");
fprintf(fp, " %d\n", numQuadrangles);
for(fiter it = firstFace(); it != lastFace(); ++it) {
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys) {
for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++)
(*it)->quadrangles[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
}
if(numTetrahedra) {
if(CTX::instance()->mesh.order == 2)
fprintf(fp, " TetrahedraP2\n"); // FIXME (check getPolynomialOrder())
else
fprintf(fp, " Tetrahedra\n");
fprintf(fp, " %d\n", numTetrahedra);
for(riter it = firstRegion(); it != lastRegion(); ++it) {
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys) {
for(std::size_t i = 0; i < (*it)->tetrahedra.size(); i++)
(*it)->tetrahedra[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
}
if(numHexahedra) {
fprintf(fp, " Hexahedra\n");
fprintf(fp, " %d\n", numHexahedra);
for(riter it = firstRegion(); it != lastRegion(); ++it) {
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys) {
for(std::size_t i = 0; i < (*it)->hexahedra.size(); i++)
(*it)->hexahedra[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
}
fprintf(fp, " End\n");
fclose(fp);
return 1;
}