unsorted nodes in vtk output

Author: gtheler

ChangeLog.md

@@ -1,5 +1,6 @@
 # Current (under development)
 
+ * VTK/VTU writer now supports sparse or unsorted node tags
  * basic example for a large-deformation mechanical case (NAFEMS GNL-5 benchmar problem)
  * outputs from `WRITE_RESULTS` or `WRITE_MESH` in `vtu` or `vtk` for transient problems create a `.pvd` file
  * instruction `PROBLEM_SOLVE` is not mandatory anymore, now FeenoX can guess where it should be called

src/Makefile-base.am

@@ -83,6 +83,7 @@ feenox_SOURCES = \
 ./mesh/calculix.c \
 ./mesh/vtk.c \
 ./mesh/vtu.c \
+./mesh/tag_index_map.c \
 ./mesh/elements/line2.c \
 ./mesh/elements/line3.c \
 ./mesh/elements/triang3.c \

src/feenox.h

@@ -1277,6 +1277,21 @@ struct bc_data_t {
 };
 
 
+// hash map entry for sparse node ids
+typedef struct {
+  size_t gmsh_tag;     // Gmsh node tag (key)
+  size_t vtk_idx;      // VTK point index (value)
+  UT_hash_handle hh;
+} tag_index_hash_t;
+
+// struct to encapsulate either mapping
+typedef struct {
+  int use_array;    // 1 for array, 0 for hashmap
+  size_t min_id, max_id;
+  size_t*array;       // if use_array
+  tag_index_hash_t *hash; // if !use_array
+} tag_index_map_t;
+
 // unstructured mesh
 struct mesh_t {
   file_t *file;
@@ -1312,6 +1327,9 @@ struct mesh_t {
   size_t *tag2index;     // array to map tags to indexes (we need a -1)
   // pointer to tag_min entries before so we can use the tag as an offset
   size_t *tag2index_from_tag_min; 
+  size_t node_tag_min;
+  size_t node_tag_max;
+  tag_index_map_t *index2tag;
 
   enum {
     data_type_element,
@@ -2284,6 +2302,7 @@ extern double feenox_gsl_function(double x, void *params);
 // mesh.c
 extern int feenox_instruction_mesh_read(void *arg);
 extern int feenox_mesh_create_nodes_argument(mesh_t *);
+extern int feenox_mesh_create_index2tag(mesh_t *this);
 extern int feenox_mesh_free(mesh_t *);
 
 extern int feenox_mesh_read_vtk(mesh_t *);
@@ -2392,6 +2411,12 @@ extern int feenox_mesh_init_nodal_indexes(mesh_t *, int dofs);
 // cell.c
 extern int feenox_mesh_element2cell(mesh_t *);
 
+// tag_index_map.c
+extern int tag_index_map_init(tag_index_map_t *map, size_t min_id, size_t max_id, size_t n_nodes, double threshold);
+extern int tag_index_map_insert(tag_index_map_t *map, size_t gmsh_id, size_t vtk_idx);
+extern size_t tag_index_map_lookup(tag_index_map_t *map, size_t gmsh_id);
+extern int tag_index_map_free(tag_index_map_t *map);
+
 // vtk.c
 extern int feenox_mesh_read_vtk_field_node(mesh_t *mesh, FILE *fp, const char *name, unsigned int size);
 extern int feenox_mesh_write_unstructured_mesh_vtk(mesh_t *mesh, FILE *file);

src/mesh/gmsh.c

@@ -1,7 +1,7 @@
 /*------------ -------------- -------- --- ----- ---   --       -            -
  *  feenox's mesh-related gmsh routines
  *
- *  Copyright (C) 2014--2024 Jeremy Theler
+ *  Copyright (C) 2014--2025 Jeremy Theler
  *
  *  This file is part of feenox.
  *
@@ -399,9 +399,7 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
           }
 
           // en msh2 si no es sparse, los tags son indices pero si es sparse es otro cantar
-          if (j+1 != tag) {
-            this->sparse = 1;
-          }
+          this->sparse |= (j+1 != tag);
 
           this->node[j].tag = tag;
           if (tag < tag_min) {
@@ -416,6 +414,7 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
 
         // finished reading the node data, check if they are sparse
         // if they are, build tag2index
+        // TODO: use the vtk hashmap
         if (this->sparse) {
           if (tag_max < this->n_nodes) {
             feenox_push_error_message("maximum node tag %d is less that number of nodes %d", tag_max, this->n_nodes);
@@ -433,9 +432,19 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
 
       } else if (version_maj == 4) {
         // number of blocks and nodes
-        size_t num_blocks;
+        size_t num_blocks = 0;
         if (version_min == 0) {
-          // in 4.0 there is no min/ max
+          // in 4.0 there's no min/max, I asked for this in
+          // https://onelab.info/pipermail/gmsh/2018/012432.html
+          //
+          // > I would need to know the biggest tag of the nodes or elements before
+          // > readuing them, not just the number of non-zero elements. Otherwise I
+          // > need to either double-parse the file or reallocate the array each time.
+          //  
+          // we could indeed add min/max vertex/element tags in the section header
+          // in a future revision of the format
+          //
+          // see also https://gitlab.onelab.info/gmsh/gmsh/-/issues/444
           size_t data[2] = {0,0};
           feenox_call(feenox_gmsh_read_data_size_t(this, 2, data, binary));
           num_blocks = data[0];
@@ -457,7 +466,8 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
         feenox_check_alloc(this->node = calloc(this->n_nodes, sizeof(node_t)));
 
         if (tag_max != 0) {
-          // we can do this as a one single pass because we have tag_max (I asked for this in v4.1)
+          // we can do this as a one single pass because we have tag_max
+          // (I asked for this in v4.1, see links above)
           if (tag_max < this->n_nodes) {
             feenox_push_error_message("maximum node tag %d is less that number of nodes %d", tag_max, this->n_nodes);
           }
@@ -496,12 +506,15 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
                 return FEENOX_ERROR;
               }
 
+              if (this->node[node_index].tag < tag_min) {
+                tag_min = this->node[node_index].tag;
+              }
               if (this->node[node_index].tag > tag_max) {
                 tag_max = this->node[node_index].tag;
               }
-              
-              // in msh4 the tags are the indices of the global mesh
               this->node[node_index].index_mesh = node_index;
+              // we don't have tag2index yet, but we can check if it is sparse
+              this->sparse |= (this->node[node_index].tag != node_index+1);
               node_index++;
             }
           } else {
@@ -530,14 +543,18 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
                 return FEENOX_ERROR;
               }
               this->tag2index_from_tag_min[this->node[node_index].tag] = node_index;
+              this->sparse |= (this->node[node_index].tag != node_index+1);
               node_index++;
             }
             feenox_free(node_coords);
           }
         }
 
+        this->node_tag_min = tag_min;
+        this->node_tag_max = tag_max;
         if (version_min == 0) {
-          // for v4.0 we need an extra loop because we did not have the actual tag_max
+          // for v4.0 we need an extra loop because we did not have
+          // the actual tag_min/tag_max before reading
           feenox_mesh_tag2index_alloc(this, tag_min, tag_max);
 
           for (size_t j = 0; j < this->n_nodes; j++) {
@@ -668,7 +685,6 @@ int feenox_mesh_read_gmsh(mesh_t *this) {
       } else if (version_maj == 4) {
         size_t num_blocks;
         if (version_min == 0) {
-          // in 4.0 there's no min/max (I asked for this)
           tag_min = 0;
           tag_max = 0;
           if (fscanf(fp, "%lu %lu", &num_blocks, &this->n_elements) < 2) {

src/mesh/mesh.c

@@ -671,8 +671,12 @@ int feenox_mesh_free(mesh_t *mesh) {
     // feenox_free(physical_group->element);
     // feenox_free(physical_group);
   }
-*/    
+*/ 
 
+  if (mesh->index2tag != NULL) {
+    feenox_call(tag_index_map_free(mesh->index2tag));
+    feenox_free(mesh->index2tag);
+  }  
   mesh->initialized = 0;
 
   return FEENOX_OK;
@@ -700,3 +704,22 @@ int feenox_mesh_create_nodes_argument(mesh_t *this) {
 
   return FEENOX_OK;
 }
+
+int feenox_mesh_create_index2tag(mesh_t *this) {
+  if (this->index2tag != NULL) {
+    feenox_call(tag_index_map_free(this->index2tag));
+    feenox_free(this->index2tag);
+  }
+  this->index2tag = calloc(1, sizeof(tag_index_map_t));
+  tag_index_map_init(this->index2tag, this->node_tag_min, this->node_tag_max, this->n_nodes, 4);
+  size_t index = 0;
+  size_t n_nodes = 0;
+  for (size_t tag = this->node_tag_min; tag <= this->node_tag_max; tag++) {
+    if ((index = this->tag2index_from_tag_min[tag]) != SIZE_MAX) {
+      tag_index_map_insert(this->index2tag, tag, index);
+      n_nodes++;
+    }
+  }
+  
+  return FEENOX_OK;
+}

src/mesh/tag_index_map.c

@@ -0,0 +1,94 @@
+/*------------ -------------- -------- --- ----- ---   --       -            -
+ *  feenox's routines to perform array or hash-based lookups
+ *
+ *  Copyright (C) 2025 Jeremy Theler
+ *
+ *  This file is part of feenox.
+ *
+ *  feenox is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  FeenoX is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with feenox.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * 
+ * Portions of this file were developed with the assistance of GitHub Copilot.
+ *------------------- ------------  ----    --------  --     -       -         -
+ */
+
+#include "feenox.h"
+
+int tag_index_map_init(tag_index_map_t *map, size_t min_id, size_t max_id, size_t n_nodes, double threshold) {
+  map->min_id = min_id;
+  map->max_id = max_id;
+  size_t range = max_id - min_id + 1;
+  if (range <= threshold * n_nodes) {
+    map->use_array = 1;
+    feenox_check_alloc(map->array = (size_t*)malloc(range * sizeof(size_t)));
+    for (size_t i = 0; i < range; ++i) {
+      map->array[i] = SIZE_MAX;
+    }
+    map->hash = NULL;
+  } else {
+    map->use_array = 0;
+    map->array = NULL;
+    map->hash = NULL;
+  }
+  return FEENOX_OK;
+}
+
+int tag_index_map_insert(tag_index_map_t *map, size_t gmsh_tag, size_t vtk_idx) {
+  if (map->use_array) {
+    map->array[gmsh_tag - map->min_id] = vtk_idx;
+  } else {
+    tag_index_hash_t *entry = NULL;
+    HASH_FIND(hh, map->hash, &gmsh_tag, sizeof(size_t), entry);
+    if (!entry) {
+      entry = (tag_index_hash_t*)malloc(sizeof(tag_index_hash_t));
+      entry->gmsh_tag = gmsh_tag;
+      entry->vtk_idx = vtk_idx;
+      HASH_ADD(hh, map->hash, gmsh_tag, sizeof(size_t), entry);
+    } else {
+      entry->vtk_idx = vtk_idx; // Overwrite if needed
+    }
+  }
+  return FEENOX_OK;
+}
+
+size_t tag_index_map_lookup(tag_index_map_t *map, size_t gmsh_tag) {
+  if (map) {
+    if (map->use_array) {
+      return map->array[gmsh_tag - map->min_id];
+    } else {
+      tag_index_hash_t *entry = NULL;
+      HASH_FIND_INT(map->hash, &gmsh_tag, entry);
+      if (entry) {
+        return entry->vtk_idx;
+      } else {
+        return SIZE_MAX; // Not found
+      }
+    }
+  } else {
+    return gmsh_tag-1;
+  }
+}
+
+int tag_index_map_free(tag_index_map_t *map) {
+  if (map->use_array) {
+    feenox_free(map->array);
+  } else {
+    tag_index_hash_t *cur, *tmp;
+    HASH_ITER(hh, map->hash, cur, tmp) {
+      HASH_DEL(map->hash, cur);
+      feenox_free(cur);
+    }
+  }
+  return FEENOX_OK;
+}