From 2f0dec3154472e650708fd518d30c508707c0394 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 15 May 2011 15:55:38 +0000
Subject: [PATCH] merge DIFFPACK patch from Jacques Lechelle

---
 Common/OpenFile.cpp   |   3 ++
 Geo/GModelIO_Mesh.cpp | 102 +++++++++++++++++++++---------------------
 2 files changed, 55 insertions(+), 50 deletions(-)

diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index fd4ef3a651..1a3c451364 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -299,6 +299,9 @@ int MergeFile(std::string fileName, bool warnIfMissing)
   else if(ext == ".mesh" || ext == ".MESH"){
     status = GModel::current()->readMESH(fileName);
   }
+  else if(ext == ".diff" || ext == ".DIFF"){
+    status = GModel::current()->readDIFF(fileName);
+  }
   else if(ext == ".med" || ext == ".MED" || ext == ".mmed" || ext == ".MMED" ||
           ext == ".rmed" || ext == ".RMED"){
     status = GModel::readMED(fileName);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index da5b4aadcb..3953916624 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -3058,14 +3058,13 @@ int GModel::readDIFF(const std::string &name)
   std::map<int, std::map<int, std::string> > physicals[4];
   std::map<int, MVertex*> vertexMap;
   std::vector<MVertex*> vertexVector;
-
-  while(1) {
-
+ 
+  {
     while(strstr(str, "Number of space dim. =") == NULL){
       if(!fgets(str, sizeof(str), fp) || feof(fp))
         break;
     }
-
+    
     int dim;
     if(sscanf(str, "%*s %*s %*s %*s %*s %d", &dim) != 1) return 0;
     Msg::Info("dimension %d", dim);
@@ -3096,24 +3095,27 @@ int GModel::readDIFF(const std::string &name)
     }
     if(sscanf(str, "%*s %*s %*s %*s %*s %*s %*s %d", &numVerticesPerElement) != 1)
       return 0;
-    Msg::Info("numVerticesPerElement %d",numVerticesPerElement);
+    Msg::Info("numVerticesPerElement %d", numVerticesPerElement);
 
     bool several_subdomains;
     if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
-    while(strstr(str, "Only one subdomain el              :") == NULL){
-      if(!fgets(str, sizeof(str), fp) || feof(fp))
-        break;
+    if(!strncmp(&str[2], "Only one material", 17) || 
+       !strncmp(&str[2], "Only one subdomain", 18)){
+      if(!strncmp(&str[37], "dpTRUE", 6) || !strncmp(&str[37], "true", 4) || 
+         !strncmp(&str[36], "dpTRUE", 6) || !strncmp(&str[36], "true", 4)){
+        several_subdomains = false;
+      }
+      else{
+        several_subdomains = true;
+      }
+      Msg::Info("several_subdomains %x %s", several_subdomains, str);
     }
-    if(!strncmp(&str[39], "dpFALSE", 6))
-      several_subdomains = true;
-    else
-      several_subdomains = false;
-    Msg::Info("several_subdomains %x", several_subdomains);
-
+    
     int nbi;
     std::vector<int> bi;
     if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
-    while(strstr(str, "Boundary indicators:") == NULL){
+    while(strstr(str, "Boundary indicators:") == NULL &&
+          strstr(str, "boundary indicators:") == NULL){
       if(!fgets(str, sizeof(str), fp) || feof(fp))
         break;
     }
@@ -3132,7 +3134,7 @@ int GModel::readDIFF(const std::string &name)
       if(sscanf(str, format_read_bi.c_str(), &bi[i]) != 1) return 0;
       Msg::Info("bi[%d]=%d", i, bi[i]);
     }
-
+    
     while(str[0] != '#'){
       if(!fgets(str, sizeof(str), fp) || feof(fp))
         break;
@@ -3148,8 +3150,8 @@ int GModel::readDIFF(const std::string &name)
       if(!fgets(str, sizeof(str), fp)) return 0;
       double xyz[3];
       int tmp;
-      if(sscanf(str, "%d ( %lf , %lf , %lf ) [%d]", &num, &xyz[0], &xyz[1], &xyz[2],
-                &tmp) != 5) return 0;
+      if(sscanf(str, "%d ( %lf , %lf , %lf ) [%d]", &num, 
+                &xyz[0], &xyz[1], &xyz[2], &tmp) != 5) return 0;
       elementary[i].resize(tmp + 1);
       elementary[i][0] = tmp;
       minVertex = std::min(minVertex, num);
@@ -3158,16 +3160,16 @@ int GModel::readDIFF(const std::string &name)
         Msg::Warning("Skipping duplicate vertex %d", num);
       else
         vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
-      if(numVertices > 100000)
+      if(numVertices > 100000) 
         Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
-      // If the vertex numbering is dense, transfer the map into a
+      // If the vertex numbering is dense, tranfer the map into a
       // vector to speed up element creation
-      if((int)vertexMap.size() == numVertices &&
+      if((int)vertexMap.size() == numVertices && 
          ((minVertex == 1 && maxVertex == numVertices) ||
           (minVertex == 0 && maxVertex == numVertices - 1))){
         Msg::Info("Vertex numbering is dense");
         vertexVector.resize(vertexMap.size() + 1);
-        if(minVertex == 1)
+        if(minVertex == 1) 
           vertexVector[0] = 0;
         else
           vertexVector[numVertices] = 0;
@@ -3176,7 +3178,7 @@ int GModel::readDIFF(const std::string &name)
           vertexVector[it->first] = it->second;
         vertexMap.clear();
       }
-      Msg::Info("%d ( %lf , %lf , %lf ) [%d]",i, xyz[0], xyz[1], xyz[2],
+      Msg::Info("%d ( %lf , %lf , %lf ) [%d]",i, xyz[0], xyz[1], xyz[2], 
                 elementary[i][0]);
       std::string format_read_bi = "%*d ( %*lf , %*lf , %*lf ) [%*d]";
       for(int j = 0; j < elementary[i][0]; j++){
@@ -3197,7 +3199,7 @@ int GModel::readDIFF(const std::string &name)
     }
     std::vector<int> material(numElements);
     std::vector<std::vector<int> > ElementsNodes(numElements);
-    for(int i = 0; i < numVertices; i++){
+    for(int i = 0; i < numElements; i++){
       ElementsNodes[i].resize(numVerticesPerElement);
     }
     char eleTypec[20]="";
@@ -3212,31 +3214,31 @@ int GModel::readDIFF(const std::string &name)
       eleType = std::string(eleTypec);
       int k2; // local number for the element
       int NoVertices; // number of vertices per element
-      if(eleType == "ElmT3n2D"){
+      if(eleType == "ElmT3n2D"){ 
         NoVertices = 3;
         static int map[3] = {0, 1, 2}; // identical to gmsh
         mapping=std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_TRI_3;
       }
-      else if(eleType == "ElmT6n2D"){
+      else if(eleType == "ElmT6n2D"){ 
         NoVertices = 6;
         static int map[6] = {0, 1, 2, 3, 4, 5}; // identical to gmsh
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_TRI_6;
       }
-      else if(eleType == "ElmB4n2D"){
+      else if(eleType == "ElmB4n2D"){ 
         NoVertices = 4;
         static int map[4] = {0, 1, 3, 2}; // local numbering
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_QUA_4;
       }
-      else if(eleType == "ElmB8n2D"){
+      else if(eleType == "ElmB8n2D"){ 
         NoVertices = 8;
         static int map[8] = {0, 1, 3, 2, 4, 6, 7, 5}; // local numbering
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_QUA_8;
       }
-      else if(eleType == "ElmB9n2D"){
+      else if(eleType == "ElmB9n2D"){ 
         NoVertices = 9;
         static int map[9] = {0, 4, 1, 7, 8, 5, 3, 6, 2}; // local numbering
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
@@ -3247,8 +3249,8 @@ int GModel::readDIFF(const std::string &name)
         static int map[4] = {0, 1, 2, 3}; // identical to gmsh
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_TET_4;
-      }
-      else if(eleType == "ElmT10n3D"){
+      } 
+      else if(eleType == "ElmT10n3D"){ 
         NoVertices = 10;
         static int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; // local numbering
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
@@ -3259,21 +3261,21 @@ int GModel::readDIFF(const std::string &name)
         static int map[8] = {4, 5, 0, 1, 7, 6, 3, 2};
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_HEX_8;
-      }
+      } 
       else if(eleType == "ElmB20n3D"){
         NoVertices = 20;
-        static int map[20] = {4, 5, 0, 1, 7, 6, 3, 2, 16, 8, 19, 13, 15, 12,
+        static int map[20] = {4, 5, 0, 1, 7, 6, 3, 2, 16, 8, 19, 13, 15, 12, 
                               14, 17, 18, 9, 11};
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_HEX_20;
-      }
+      } 
       else if(eleType == "ElmB27n3D"){
         NoVertices = 27;
-        static int map[27] = {4, 16, 5, 10, 21, 12, 0, 8, 1, 17, 25, 18, 22,
+        static int map[27] = {4, 16, 5, 10, 21, 12, 0, 8, 1, 17, 25, 18, 22, 
                               26, 23, 9, 20, 11, 7, 19, 6, 15, 24, 14, 3, 13, 2};
         mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
         type = MSH_HEX_27;
-      }
+      } 
       else{
         return 0;
       }
@@ -3282,7 +3284,7 @@ int GModel::readDIFF(const std::string &name)
         if(format_read_vertices[format_read_vertices.size()-2] != '*') {
           format_read_vertices[format_read_vertices.size()-1] = '*';
           format_read_vertices += "d %d";
-        }
+        } 
         else
           format_read_vertices += " %d";
         k2 = mapping[k];
@@ -3294,21 +3296,21 @@ int GModel::readDIFF(const std::string &name)
         indices[j] = ElementsNodes[i - 1][j];
       std::vector<MVertex*> vertices;
       if(vertexVector.size()){
-        if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices))
+        if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices)) 
           return 0;
       }
       else{
-        if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices))
+        if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices)) 
           return 0;
       }
-      createElementMSH(this, num, type, physical, elementary[i-1][1], partition,
-                       vertices, elements, physicals);
+      createElementMSH(this, num, type, physical, elementary[i-1][1], partition, 
+                       vertices, elements, physicals); 
       // trouble if elementary[i-1][0]>1 nodal post-processing needed ?
-      if(numElements > 100000)
+      if(numElements > 100000) 
         Msg::ProgressMeter(i + 1, numElements, "Reading elements");
     }
   }
-
+  
   // store the elements in their associated elementary entity. If the
   // entity does not exist, create a new (discrete) one.
   for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
@@ -3344,7 +3346,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
     Msg::Error("Unable to open file '%s'", name.c_str());
     return 0;
   }
-
+  
   if(noPhysicalGroups()) saveAll = true;
 
   // get the number of vertices and index the vertices in a continuous
@@ -3390,7 +3392,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
     if(entities[i]->physicals.size() || saveAll)
       for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
         dim = std::max(dim, entities[i]->getMeshElement(j)->getDim());
-
+  
   // loop over all elements we need to save
   int numElements = 0, maxNumNodesPerElement = 0;
   for(unsigned int i = 0; i < entities.size(); i++){
@@ -3412,13 +3414,13 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   fprintf(fp, " Number of nodes      =  %d\n\n", numVertices);
   fprintf(fp, " All elements are of the same type : dpTRUE\n");
   fprintf(fp, " Max number of nodes in an element: %d \n", maxNumNodesPerElement);
-  fprintf(fp, " Only one subdomain el              : dpFALSE\n");
+  fprintf(fp, " Only one subdomain               : dpFALSE\n");
   fprintf(fp, " Lattice data                     ? 0\n\n\n\n");
   fprintf(fp, " %d Boundary indicators:  ", (int)boundaryIndicators.size());
   for(std::list<int>::iterator it = boundaryIndicators.begin();
       it != boundaryIndicators.end(); it++)
     fprintf(fp, " %d", *it);
-
+  
   fprintf(fp, "\n\n\n");
   fprintf(fp,"  Nodal coordinates and nodal boundary indicators,\n");
   fprintf(fp,"  the columns contain:\n");
@@ -3427,7 +3429,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
   fprintf(fp,"   - no of boundary indicators that are set (ON)\n");
   fprintf(fp,"   - the boundary indicators that are set (ON) if any.\n");
   fprintf(fp,"#\n");
-
+  
   // write mesh vertices
   for(unsigned int i = 0; i < entities.size(); i++){
     for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
@@ -3442,7 +3444,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
       }
     }
   }
-
+  
   fprintf(fp, "\n");
   fprintf(fp, "\n");
   fprintf(fp,     "  Element types and connectivity\n");
@@ -3465,7 +3467,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
     }
   }
   fprintf(fp, "\n");
-
+  
   fclose(fp);
   return 1;
 }
-- 
GitLab