From efa4f073a657c22ce66639e51a186709a55248e2 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 3 Oct 2009 11:30:43 +0000
Subject: [PATCH] fix VisualFEA output

---
 Common/CommandLine.cpp |   6 +-
 Common/CreateFile.cpp  |  14 ++--
 Common/GmshDefines.h   |   2 +-
 Geo/GModel.h           |   6 +-
 Geo/GModelIO_Mesh.cpp  | 151 ++++++++++++++++-------------------------
 Geo/MElement.cpp       |  12 ++++
 Geo/MElement.h         |   2 +
 7 files changed, 87 insertions(+), 106 deletions(-)

diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 33511f2280..d26c848a6e 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -451,9 +451,9 @@ void GetOptions(int argc, char *argv[])
             CTX::instance()->mesh.format = FORMAT_DIFF;
           else if(!strcmp(argv[i], "med"))
             CTX::instance()->mesh.format = FORMAT_MED;
-           else if(!strcmp(argv[i], "fea"))
-            CTX::instance()->mesh.format = FORMAT_VISUALFEA;
-         else
+          else if(!strcmp(argv[i], "fea"))
+            CTX::instance()->mesh.format = FORMAT_FEA;
+          else
             Msg::Fatal("Unknown mesh format");
           i++;
         }
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 7d369001c5..8abd69e6b4 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -38,7 +38,7 @@ int GuessFileFormatFromFileName(std::string fileName)
   else if(ext == ".stl")  return FORMAT_STL;
   else if(ext == ".cgns") return FORMAT_CGNS;
   else if(ext == ".med")  return FORMAT_MED;
-  else if(ext == ".fea")  return FORMAT_VISUALFEA;
+  else if(ext == ".fea")  return FORMAT_FEA;
   else if(ext == ".mesh") return FORMAT_MESH;
   else if(ext == ".bdf")  return FORMAT_BDF;
   else if(ext == ".diff") return FORMAT_DIFF;
@@ -79,7 +79,7 @@ std::string GetDefaultFileName(int format)
   case FORMAT_STL:  name += ".stl"; break;
   case FORMAT_CGNS: name += ".cgns"; break;
   case FORMAT_MED:  name += ".med"; break;
-  case FORMAT_VISUALFEA: name += ".fea"; break;
+  case FORMAT_FEA:  name += ".fea"; break;
   case FORMAT_MESH: name += ".mesh"; break;
   case FORMAT_BDF:  name += ".bdf"; break;
   case FORMAT_DIFF: name += ".diff"; break;
@@ -95,9 +95,9 @@ std::string GetDefaultFileName(int format)
   case FORMAT_SVG:  name += ".svg"; break;
   case FORMAT_PPM:  name += ".ppm"; break;
   case FORMAT_YUV:  name += ".yuv"; break;
-  case FORMAT_BREP:  name += ".brep"; break;
-  case FORMAT_IGES:  name += ".iges"; break;
-  case FORMAT_STEP:  name += ".step"; break;
+  case FORMAT_BREP: name += ".brep"; break;
+  case FORMAT_IGES: name += ".iges"; break;
+  case FORMAT_STEP: name += ".step"; break;
   default: break;
   }
   return name;
@@ -207,8 +207,8 @@ void CreateOutputFile(std::string fileName, int format)
        CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
     break;
 
-  case FORMAT_VISUALFEA:
-    GModel::current()->writeVisualFEA
+  case FORMAT_FEA:
+    GModel::current()->writeFEA
       (fileName, CTX::instance()->mesh.saveElementTagType, 
        CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
     break;
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index bbd732ae8b..7b04811ea6 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -40,7 +40,7 @@
 #define FORMAT_BREP          35
 #define FORMAT_STEP          36
 #define FORMAT_IGES          37
-#define FORMAT_VISUALFEA     38
+#define FORMAT_FEA           38
 
 // Element types
 #define TYPE_PNT 1
diff --git a/Geo/GModel.h b/Geo/GModel.h
index c7fffad427..96a913fd61 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -344,9 +344,9 @@ class GModel
   int writeMSH(const std::string &name, double version=1.0, bool binary=false,
                bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0);
 
-  // VisualFEA file format
-  int writeVisualFEA(const std::string &name, int elementTagType, 
-		     bool saveAll, double scalingFactor);
+  // Visual FEA file format
+  int writeFEA(const std::string &name, int elementTagType, 
+               bool saveAll, double scalingFactor);
 
   // mesh statistics (saved as a Gmsh post-processing view)
   int writePOS(const std::string &name, bool printElementary,
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 78969fc777..ee25d5e1fe 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1610,98 +1610,6 @@ int GModel::readMESH(const std::string &name)
   return 1;
 }
 
-
-int GModel::writeVisualFEA(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);
-  int numHexahedra = 0, numTriangles = 0, numQuadrangles = 0, numTetrahedra = 0;
-  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();
-    }
-  }
-
-  fprintf(fp,"33\n");
-  fprintf(fp,"%d %d %d\n", numVertices, numTriangles+numQuadrangles, numTetrahedra+numHexahedra);
-  
-  std::vector<GEntity*> entities;
-  getEntities(entities);
-  for(unsigned int i = 0; i < entities.size(); i++)
-    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
-      if (entities[i]->mesh_vertices[j]->getIndex() >= 0)
-	fprintf(fp,"%d %g %g %g\n",entities[i]->mesh_vertices[j]->getIndex(),
-		entities[i]->mesh_vertices[j]->x()*scalingFactor,
-		entities[i]->mesh_vertices[j]->y()*scalingFactor,
-		entities[i]->mesh_vertices[j]->z()*scalingFactor);  
-  int iElement = 1;
-  for(fiter it = firstFace(); it != lastFace(); ++it){
-    int numPhys = (*it)->physicals.size();
-    if(saveAll || numPhys){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
-	MTriangle *t = (*it)->triangles[i];
-	fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
-		t->getVertex(0)->getNum(),
-		t->getVertex(1)->getNum(),
-		t->getVertex(2)->getNum());
-      }
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
-	MQuadrangle *t = (*it)->quadrangles[i];
-	fprintf(fp,"%d %d 4 %d %d %d\n",iElement++,(*it)->tag(),
-		t->getVertex(0)->getNum(),
-		t->getVertex(1)->getNum(),
-		t->getVertex(2)->getNum(),
-		t->getVertex(3)->getNum());
-      }
-    }
-  }
-
-  for(riter it = firstRegion(); it != lastRegion(); ++it){
-    int numPhys = (*it)->physicals.size();
-    if(saveAll || numPhys){
-      for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++){
-	MTetrahedron *t = (*it)->tetrahedra[i];
-	fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
-		t->getVertex(0)->getNum(),
-		t->getVertex(1)->getNum(),
-		t->getVertex(2)->getNum(),
-		t->getVertex(3)->getNum());
-      }
-      for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++){
-	MHexahedron *t = (*it)->hexahedra[i];
-	fprintf(fp,"%d %d 3 %d %d %d\n",iElement++,(*it)->tag(),
-		t->getVertex(0)->getNum(),
-		t->getVertex(1)->getNum(),
-		t->getVertex(2)->getNum(),
-		t->getVertex(3)->getNum(),
-		t->getVertex(4)->getNum(),
-		t->getVertex(5)->getNum(),
-		t->getVertex(6)->getNum(),
-		t->getVertex(7)->getNum());
-      }
-    }
-  }
-  
-  fclose(fp);
-  return 1;
-}
-
-
 int GModel::writeMESH(const std::string &name, int elementTagType, 
                       bool saveAll, double scalingFactor)
 {
@@ -1800,6 +1708,65 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
   return 1;
 }
 
+int GModel::writeFEA(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), num2D = 0, num3D = 0;
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+    if(saveAll || (*it)->physicals.size()) 
+      num2D += (*it)->getNumMeshElements();
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    if(saveAll || (*it)->physicals.size())
+      num3D += (*it)->getNumMeshElements();
+
+  fprintf(fp,"33\n");
+  if(num2D && num3D)
+    fprintf(fp,"%d %d %d\n", numVertices, num2D, num3D);
+  else
+    fprintf(fp,"%d %d\n", numVertices, num2D ? num2D : num3D);
+  
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
+      if(entities[i]->mesh_vertices[j]->getIndex() >= 0)
+	fprintf(fp,"%d %g %g %g\n", entities[i]->mesh_vertices[j]->getIndex(),
+		entities[i]->mesh_vertices[j]->x() * scalingFactor,
+		entities[i]->mesh_vertices[j]->y() * scalingFactor,
+		entities[i]->mesh_vertices[j]->z() * scalingFactor);  
+
+  int iElement = 1;
+  for(fiter it = firstFace(); it != lastFace(); ++it){
+    int numPhys = (*it)->physicals.size();
+    if(saveAll || numPhys)
+      for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++)
+        (*it)->getMeshElement(i)->writeFEA
+          (fp, elementTagType, iElement++, (*it)->tag(), 
+           numPhys ? (*it)->physicals[0] : 0);
+  }
+
+  for(riter it = firstRegion(); it != lastRegion(); ++it){
+    int numPhys = (*it)->physicals.size();
+    if(saveAll || numPhys)
+      for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++)
+        (*it)->getMeshElement(i)->writeFEA
+          (fp, elementTagType, iElement++, (*it)->tag(), 
+           numPhys ? (*it)->physicals[0] : 0);
+  }
+  
+  fclose(fp);
+  return 1;
+}
+
 static int getFormatBDF(char *buffer, int &keySize)
 {
   if(buffer[keySize] == '*'){ keySize++; return 2; } // long fields
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 14c715ed05..01b7398403 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -604,6 +604,18 @@ void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
           (elementTagType == 2) ? physical : elementary);
 }
 
+void MElement::writeFEA(FILE *fp, int elementTagType, int num, int elementary, 
+                        int physical)
+{
+  int numVert = getNumVertices();
+  setVolumePositive();
+  fprintf(fp, "%d %d %d", num, (elementTagType == 3) ? _partition : 
+          (elementTagType == 2) ? physical : elementary, numVert);
+  for(int i = 0; i < numVert; i++)
+    fprintf(fp, " %d", getVertex(i)->getIndex());
+  fprintf(fp, "\n");
+}
+
 void MElement::writeBDF(FILE *fp, int format, int elementTagType, int elementary,
                         int physical)
 {
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 8407e8b1f5..6f98c42918 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -234,6 +234,8 @@ class MElement
   virtual void writeVTK(FILE *fp, bool binary=false, bool bigEndian=false);
   virtual void writeMESH(FILE *fp, int elementTagType=1, int elementary=1, 
                          int physical=0);
+  virtual void writeFEA(FILE *fp, int elementTagType, int num, int elementary, 
+                        int physical);
   virtual void writeBDF(FILE *fp, int format=0, int elementTagType=1, 
                         int elementary=1, int physical=0);
   virtual void writeDIFF(FILE *fp, int num, bool binary=false, 
-- 
GitLab