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