From 08015b45cda3414977cb44b07188e1b657523608 Mon Sep 17 00:00:00 2001
From: Corentin Prigent <corentin.prigent@u-bordeaux.fr>
Date: Wed, 31 Jan 2024 16:11:51 +0100
Subject: [PATCH 1/2] Added output messages in npy reading

---
 src/hashhexa.c | 2 +-
 src/inout.c    | 8 ++++++--
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/hashhexa.c b/src/hashhexa.c
index befdc76..5777e5b 100644
--- a/src/hashhexa.c
+++ b/src/hashhexa.c
@@ -249,7 +249,7 @@ int H2T_hashQuad(MMG5_pMesh mesh,MMG5_Hash *hash) {
  */
 int H2T_hashGetRef(MMG5_pMesh mesh,MMG5_Hash *hash) {
   MMG5_pTetra pt;
-  int         *adja;
+  MMG5_int         *adja;
   int         ie,i,nt,ia,ib,ic,k,ref;
 
   /* Create tetrahedra adjacency */
diff --git a/src/inout.c b/src/inout.c
index cd75f26..bc8d0a7 100644
--- a/src/inout.c
+++ b/src/inout.c
@@ -25,8 +25,8 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   FILE* inm;
   unsigned char buffer = 0x00;
   char* str = NULL;
-  int pos1, pos2, dim = 0, t[3], nhex;
-  MMG5_int np, ne, i, j, k, ref, pos;
+  int pos1, pos2, dim = 0, t[3];
+  MMG5_int nhex, np, ne, i, j, k, ref, pos;
 
   /* Input data and creation of the hexa array */
   if( !(inm = fopen(mmgMesh->namein,"rb")) ) {
@@ -84,6 +84,7 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   /* Point coordinates from grid indices */
   ref = 0;
   pos = 0;
+  fprintf(stdout,"  READING %" MMG5_PRId " VERTICES\n",np);
   for (i=0;i<t[0];i++) {
     for (j=0;j<t[1];j++) {
       for (k=0;k<t[2];k++) {
@@ -96,6 +97,7 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
 
   /* Hexahedra */
   pos = 1;
+  fprintf(stdout,"  READING %" MMG5_PRId " HEXA\n",nhex);
   for ( i=0; i<t[0]-1; ++i ) {
     for ( j=0; j<t[1]-1; ++j ) {
       for ( k=0; k<t[2]-1; ++k ) {
@@ -121,6 +123,7 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   /* Edges */
   ref = 0;
   pos = 0;
+  fprintf(stdout,"  READING %" MMG5_PRId " EDGES\n",ne);
   for (i=0;i<t[0]-1; ++i) {
     int np0, np1;
     np0 = H2T_npy_point_index(i  ,0     ,0  ,t);
@@ -179,6 +182,7 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   }
 
   /* Ridges */
+  fprintf(stdout,"  READING %" MMG5_PRId " RIDGES\n",ne);
   for (i=1;i<=ne;i++) {
     MMG3D_Set_ridge(mmgMesh,i);
   }

From e5f6edf9d5299f2af93acfe101ac92ed66c1ce0e Mon Sep 17 00:00:00 2001
From: Corentin Prigent <corentin.prigent@u-bordeaux.fr>
Date: Thu, 1 Feb 2024 15:45:58 +0100
Subject: [PATCH 2/2] Proper input of initial data from .npy file : number of
 entities in every direction increased by one

---
 src/inout.c | 78 +++++++++++++++++++++++++++--------------------------
 1 file changed, 40 insertions(+), 38 deletions(-)

diff --git a/src/inout.c b/src/inout.c
index bc8d0a7..059f068 100644
--- a/src/inout.c
+++ b/src/inout.c
@@ -16,8 +16,8 @@
 #include <math.h>
 
 static inline
-int H2T_npy_point_index(int i,int j,int k, int *n) {
-  return n[1]*n[2]*i + n[2]*j + k + 1;
+MMG5_int H2T_npy_point_index(int i,int j,int k, int *n) {
+  return (n[1]+1)*(n[2]+1)*i + (n[2]+1)*j + k + 1;
 }
 
 int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
@@ -71,9 +71,9 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
     fread(&buffer,sizeof buffer,1,inm);
   }
 
-  np   = t[0]*t[1]*t[2];
-  nhex = (t[0]-1)*(t[1]-1)*(t[2]-1);
-  ne   = 4 * ((t[0]-1)+(t[1]-1)+(t[2]-1));
+  np   = (t[0]+1)*(t[1]+1)*(t[2]+1);
+  nhex = t[0]*t[1]*t[2];
+  ne   = 4 * (t[0]+t[1]+t[2]);
 
   if ( H2T_Set_meshSize(mmgMesh,np,nhex,0,ne) != 1 ) {
     return -1;
@@ -85,9 +85,9 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   ref = 0;
   pos = 0;
   fprintf(stdout,"  READING %" MMG5_PRId " VERTICES\n",np);
-  for (i=0;i<t[0];i++) {
-    for (j=0;j<t[1];j++) {
-      for (k=0;k<t[2];k++) {
+  for (i=0;i<=t[0];i++) {
+    for (j=0;j<=t[1];j++) {
+      for (k=0;k<=t[2];k++) {
         if ( H2T_Set_vertex(mmgMesh,(double)i ,(double)j ,(double)k,ref,++pos) != 1 ) {
           return -1;
         }
@@ -98,10 +98,10 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   /* Hexahedra */
   pos = 1;
   fprintf(stdout,"  READING %" MMG5_PRId " HEXA\n",nhex);
-  for ( i=0; i<t[0]-1; ++i ) {
-    for ( j=0; j<t[1]-1; ++j ) {
-      for ( k=0; k<t[2]-1; ++k ) {
-        int iadr = 9*pos;
+  for ( i=0; i<t[0]; ++i ) {
+    for ( j=0; j<t[1]; ++j ) {
+      for ( k=0; k<t[2]; ++k ) {
+        MMG5_int iadr = 9*pos;
 
         /* Hexa vertices */
         (*tabhex)[iadr+0] = H2T_npy_point_index(i  ,j  ,k  ,t);
@@ -114,7 +114,7 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
         (*tabhex)[iadr+7] = H2T_npy_point_index(i  ,j+1,k+1,t);
 
         /* Hexa references */
-        fread(&(*tabhex)[iadr+8],sizeof(int16_t),1,inm);
+        fread(&(*tabhex)[iadr+8],sizeof(uint32_t),1,inm);
         ++pos;
       }
     }
@@ -124,60 +124,62 @@ int H2T_loadNpy(MMG5_pMesh mmgMesh, int** tabhex, char* filename) {
   ref = 0;
   pos = 0;
   fprintf(stdout,"  READING %" MMG5_PRId " EDGES\n",ne);
-  for (i=0;i<t[0]-1; ++i) {
-    int np0, np1;
+  for (i=0;i<t[0]; ++i) {
+    MMG5_int np0, np1;
     np0 = H2T_npy_point_index(i  ,0     ,0  ,t);
     np1 = H2T_npy_point_index(i+1,0     ,0  ,t);
-    H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
+    if (H2T_Set_edge(mmgMesh,np0,np1,ref,++pos) != 1) {
+      return -1;
+    }
 
-    np0 = H2T_npy_point_index(i  ,t[1]-1,0  ,t);
-    np1 = H2T_npy_point_index(i+1,t[1]-1,0  ,t);
+    np0 = H2T_npy_point_index(i  ,t[1],0  ,t);
+    np1 = H2T_npy_point_index(i+1,t[1],0  ,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(i  ,0     ,t[2]-1,t);
-    np1 = H2T_npy_point_index(i+1,0     ,t[2]-1,t);
+    np0 = H2T_npy_point_index(i  ,0     ,t[2],t);
+    np1 = H2T_npy_point_index(i+1,0     ,t[2],t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(i  ,t[1]-1,t[2]-1,t);
-    np1 = H2T_npy_point_index(i+1,t[1]-1,t[2]-1,t);
+    np0 = H2T_npy_point_index(i  ,t[1],t[2],t);
+    np1 = H2T_npy_point_index(i+1,t[1],t[2],t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
   }
 
-  for (j=0;j<t[1]-1; ++j) {
-    int np0, np1;
+  for (j=0;j<t[1]; ++j) {
+    MMG5_int np0, np1;
     np0 = H2T_npy_point_index(0  ,j  ,0  ,t);
     np1 = H2T_npy_point_index(0  ,j+1,0  ,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(t[0]-1,j  ,0  ,t);
-    np1 = H2T_npy_point_index(t[0]-1,j+1,0  ,t);
+    np0 = H2T_npy_point_index(t[0],j  ,0  ,t);
+    np1 = H2T_npy_point_index(t[0],j+1,0  ,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(0  ,j  ,t[2]-1,t);
-    np1 = H2T_npy_point_index(0  ,j+1,t[2]-1,t);
+    np0 = H2T_npy_point_index(0  ,j  ,t[2],t);
+    np1 = H2T_npy_point_index(0  ,j+1,t[2],t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(t[0]-1,j  ,t[2]-1,t);
-    np1 = H2T_npy_point_index(t[0]-1,j+1,t[2]-1,t);
+    np0 = H2T_npy_point_index(t[0],j  ,t[2],t);
+    np1 = H2T_npy_point_index(t[0],j+1,t[2],t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
   }
 
-  for (k=0;k<t[2]-1; ++k) {
-    int np0, np1;
+  for (k=0;k<t[2]; ++k) {
+    MMG5_int np0, np1;
     np0 = H2T_npy_point_index(0  ,0  ,k  ,t);
     np1 = H2T_npy_point_index(0  ,0  ,k+1,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(t[0]-1,0  ,k  ,t);
-    np1 = H2T_npy_point_index(t[0]-1,0  ,k+1,t);
+    np0 = H2T_npy_point_index(t[0],0  ,k  ,t);
+    np1 = H2T_npy_point_index(t[0],0  ,k+1,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(0  ,t[1]-1,k  ,t);
-    np1 = H2T_npy_point_index(0  ,t[1]-1,k+1,t);
+    np0 = H2T_npy_point_index(0  ,t[1],k  ,t);
+    np1 = H2T_npy_point_index(0  ,t[1],k+1,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
 
-    np0 = H2T_npy_point_index(t[0]-1,t[1]-1,k  ,t);
-    np1 = H2T_npy_point_index(t[0]-1,t[1]-1,k+1,t);
+    np0 = H2T_npy_point_index(t[0],t[1],k  ,t);
+    np1 = H2T_npy_point_index(t[0],t[1],k+1,t);
     H2T_Set_edge(mmgMesh,np0,np1,ref,++pos);
   }