1 #include "hxt_mesh.h"
2 #include "hxt_tools.h"
3 #include <stdio.h>
4 #include <float.h>
5 
6 /* Element types, same as Gmsh */
7 #define POINTID 15
8 #define LINEID 1
9 #define TRIID 2
10 #define TETID 4
11 
hxtMeshCreate(HXTMesh ** mesh)12 HXTStatus  hxtMeshCreate ( HXTMesh** mesh) {
13   HXT_CHECK( hxtMalloc (mesh, sizeof(HXTMesh)) );
14   if (*mesh == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
15 
16   // This line, using a compound literal, will initialize everything
17   // to zero (and pointers to NULL)
18   **mesh = (HXTMesh) {0};
19 
20   return HXT_STATUS_OK;
21 }
22 
hxtMeshDelete(HXTMesh ** mesh)23 HXTStatus hxtMeshDelete ( HXTMesh** mesh) {
24   // free on NULL pointer is well defined (it does nothing)
25 
26   // vertices
27   HXT_CHECK( hxtAlignedFree(&(*mesh)->vertices.coord) );
28 
29   // tetrahedra
30   HXT_CHECK( hxtAlignedFree(&(*mesh)->tetrahedra.color) );
31   HXT_CHECK( hxtAlignedFree(&(*mesh)->tetrahedra.flag) );
32   HXT_CHECK( hxtAlignedFree(&(*mesh)->tetrahedra.node) );
33   HXT_CHECK( hxtAlignedFree(&(*mesh)->tetrahedra.neigh) );
34   // HXT_CHECK( hxtAlignedFree(&(*mesh)->tetrahedra.subdet) );
35 
36   // triangles
37   HXT_CHECK( hxtAlignedFree(&(*mesh)->triangles.node) );
38   HXT_CHECK( hxtAlignedFree(&(*mesh)->triangles.color) );
39 
40   // lines
41   HXT_CHECK( hxtAlignedFree(&(*mesh)->lines.node) );
42   HXT_CHECK( hxtAlignedFree(&(*mesh)->lines.color) );
43 
44   // points
45   HXT_CHECK( hxtAlignedFree(&(*mesh)->points.node) );
46   HXT_CHECK( hxtAlignedFree(&(*mesh)->points.color) );
47 
48   // boundary representation
49   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.numSurfacesPerVolume) );
50   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.surfacesPerVolume) );
51   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.numCurvesPerSurface) );
52   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.curvesPerSurface) );
53   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.endPointsOfCurves) );
54   HXT_CHECK( hxtAlignedFree(&(*mesh)->brep.points) );
55 
56   // reset every fields to zero
57   **mesh = (HXTMesh) {0};
58 
59   HXT_CHECK( hxtFree(mesh) );
60   return HXT_STATUS_OK;
61 }
62 
63 
64 // TODO: more checking of fgets
ReadNodesFromGmsh(FILE * fp,HXTMesh * m)65 static HXTStatus ReadNodesFromGmsh(FILE *fp,  HXTMesh* m)
66 {
67   rewind (fp);
68   size_t n;
69   char buf[BUFSIZ];
70   // scan for Nodes
71   m->vertices.num = 0;
72   while( fgets(buf, BUFSIZ, fp )!=NULL){
73 
74     if(strstr(buf, "$Nodes")){
75       if(fgets(buf, BUFSIZ, fp )==NULL)
76         return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
77 
78       m->vertices.num = atoi(buf);
79       HXT_CHECK( hxtAlignedMalloc(&m->vertices.coord, sizeof( double )*4*m->vertices.num) );
80       if (m->vertices.coord == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
81       m->vertices.size = m->vertices.num;
82       for(n=0;n<m->vertices.num;++n){
83         if(fgets(buf, BUFSIZ, fp )==NULL)
84           return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
85         sscanf(buf, "%*d %lf %lf %lf", &m->vertices.coord[(size_t) 4*n+0], &m->vertices.coord[(size_t) 4*n+1], &m->vertices.coord[(size_t) 4*n+2]);
86       }
87       break;
88     }
89   }
90   return HXT_STATUS_OK;
91 }
92 
93 
ReadMeshFormatFromGmsh(FILE * fp)94 static HXTStatus ReadMeshFormatFromGmsh(FILE *fp)
95 {
96   char buf[BUFSIZ];
97   int found = 0;
98   while(fgets(buf, BUFSIZ, fp )){
99     if(strstr(buf, "$MeshFormat")){
100       if(fgets(buf, BUFSIZ, fp )==NULL)
101         return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to format");
102       float version;
103       int match = sscanf(buf, " %f 0 8 ", &version);
104       if(match!=1 || version!=2.2f) {
105         return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "HXT only support Gmsh's MSH2 ASCII file format %f %d", version, match);
106       }
107       found = 1;
108       break;
109     }
110   }
111 
112   if(!found)
113     return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Unrecognized format. HXT only support Gmsh's MSH2 ASCII file format");
114 
115   return HXT_STATUS_OK;
116 }
117 
118 
119 // TODO: add possibility to read the elements from an array
120 
121 // TODO: more checking of fgets
122 // TODO: why set to zero AND use malloc ? => either set to zero and use realloc, either trust it is zero and malloc
ReadElementsFromGmsh(FILE * fp,HXTMesh * m)123 static HXTStatus ReadElementsFromGmsh(FILE *fp, HXTMesh* m)
124 {
125   int k;
126   char buf[BUFSIZ];
127 
128   rewind (fp);
129 
130   m->lines.num = 0;
131   m->points.num = 0;
132   m->triangles.num = 0;
133   m->tetrahedra.num = 0;
134 
135   while(fgets(buf, BUFSIZ, fp )){
136     if(strstr(buf, "$Elements")){
137       if(fgets(buf, BUFSIZ, fp )==NULL)
138         return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
139       int tmpK = atoi(buf);
140       for(k=0;k<tmpK;++k){
141         int etype = 0;
142         if(fgets(buf, BUFSIZ, fp )==NULL)
143           return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
144         sscanf(buf, "%*d %d", &etype);
145         if(etype==TETID){ // tets
146           ++(m->tetrahedra.num);
147         }
148         else if(etype==TRIID){ // triangles
149           ++(m->triangles.num);
150         }
151         else if(etype==LINEID){ // lines
152           ++(m->lines.num);
153         }
154         else if(etype==POINTID){ // points
155           //          ++(m->lines.num);
156           ++(m->points.num);
157         }
158         else {
159           return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "unsupported element. HXT only support points, lines, triangles and tetrahedra.");
160         }
161       }
162       break;
163     }
164   }
165 
166   rewind (fp);
167 
168 
169   // prescan for elements
170   if (m->tetrahedra.num){
171     HXT_CHECK( hxtAlignedMalloc(&m->tetrahedra.node, (m->tetrahedra.num)*4*sizeof(uint32_t)) );
172     if (m->tetrahedra.node == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
173     HXT_CHECK( hxtAlignedMalloc(&m->tetrahedra.color, (m->tetrahedra.num)*sizeof(uint32_t)) );
174     if (m->tetrahedra.color == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
175     m->tetrahedra.size = m->tetrahedra.num;
176   }
177   if (m->triangles.num){
178     HXT_CHECK( hxtAlignedMalloc(&m->triangles.node, (m->triangles.num)*3*sizeof(uint32_t)) );
179     if (m->triangles.node == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
180     HXT_CHECK( hxtAlignedMalloc(&m->triangles.color, (m->triangles.num)*sizeof(uint32_t)) );
181     if (m->triangles.color == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
182     m->triangles.size = m->triangles.num;
183   }
184   if (m->lines.num){
185     HXT_CHECK( hxtAlignedMalloc(&m->lines.node, (m->lines.num)*2*sizeof(uint32_t)) );
186     if (m->lines.node == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
187     HXT_CHECK( hxtAlignedMalloc(&m->lines.color, (m->lines.num)*sizeof(uint32_t)) );
188     if (m->lines.color == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
189     m->lines.size = m->lines.num;
190   }
191   if (m->points.num){
192     HXT_CHECK( hxtAlignedMalloc(&m->points.node, (m->points.num)*sizeof(uint32_t)) );
193     if (m->points.node == NULL)return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
194     m->points.size = m->points.num;
195   }
196 
197   while( fgets(buf, BUFSIZ, fp )){
198 
199     if(strstr(buf, "$Elements")){
200       if(fgets(buf, BUFSIZ, fp )==NULL)
201         return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
202 
203       int tmpK = atoi(buf);
204       m->tetrahedra.num=0;
205       m->triangles.num=0;
206       m->lines.num=0;
207       m->points.num=0;
208 
209       for(k=0;k<tmpK;++k){
210         int etype = 0, ntags;
211         if(fgets(buf, BUFSIZ, fp )==NULL)
212           return HXT_ERROR_MSG(HXT_STATUS_READ_ERROR, "Failed to read line");
213         sscanf(buf, "%*d %d %d", &etype, &ntags);
214         if(etype==TETID){ // tets
215           if(ntags==2){ //
216             int a, b, c, d, color;
217             sscanf(buf, "%*d %*d %*d %*d %d  %d %d %d %d",
218                    &color,&a, &b, &c, &d);
219             m->tetrahedra.node[4*m->tetrahedra.num+0] = a-1;
220             m->tetrahedra.node[4*m->tetrahedra.num+1] = b-1;
221             m->tetrahedra.node[4*m->tetrahedra.num+2] = c-1;
222             m->tetrahedra.node[4*m->tetrahedra.num+3] = d-1;
223             m->tetrahedra.color[m->tetrahedra.num] = color;
224           }
225           ++m->tetrahedra.num;
226         }
227         else if(etype==TRIID){ // triangles
228           if(ntags==2){ //
229             int a, b, c, color;
230             sscanf(buf, "%*d %*d %*d %*d %d  %d %d %d",
231                    &color,&a, &b, &c);
232             m->triangles.node[3*m->triangles.num+0] = a-1;
233             m->triangles.node[3*m->triangles.num+1] = b-1;
234             m->triangles.node[3*m->triangles.num+2] = c-1;
235             m->triangles.color[m->triangles.num] = color;
236           }
237           ++m->triangles.num;
238         }
239         else if(etype==LINEID){ // lines
240           if(ntags==2){ //
241             int a, b, color;
242             sscanf(buf, "%*d %*d %*d %*d %d  %d %d",
243             &color, &a, &b);
244             m->lines.node[2*m->lines.num+0] = a-1;
245             m->lines.node[2*m->lines.num+1] = b-1;
246             m->lines.color[m->lines.num] = color;
247           }
248           ++m->lines.num;
249         }
250         else if(etype==POINTID){ // points
251           if(ntags==2){ //
252             int a;
253             sscanf(buf, "%*d %*d %*d %*d %*d  %d",&a);
254             m->points.node[m->points.num] = a-1;
255             //            m->lines.node[2*m->lines.num+0] = a-1;
256             //            m->lines.node[2*m->lines.num+1] = b-1;
257           }
258           ++m->points.num;
259         }
260         else return HXT_STATUS_READ_ERROR;
261       }
262       break;
263     }
264   }
265 
266   return HXT_STATUS_OK;
267 }
268 
hxtMeshReadGmsh(HXTMesh * m,const char * filename)269 HXTStatus  hxtMeshReadGmsh  ( HXTMesh* m , const char *filename) {
270   FILE *f = fopen (filename, "r");
271   if (f==NULL)
272     return HXT_ERROR_MSG(HXT_STATUS_FILE_CANNOT_BE_OPENED,
273       "Cannot open mesh file \"%s\"",(filename==NULL)?"(null)":filename);
274 
275   HXT_CHECK(
276     ReadMeshFormatFromGmsh(f));
277 
278   HXT_CHECK(
279     ReadNodesFromGmsh(f,m));
280 
281   HXT_CHECK(
282     ReadElementsFromGmsh(f, m));
283 
284   fclose (f);
285   return HXT_STATUS_OK;
286 }
287 
288 
289 
290 // writes the gmsh stuff
hxtMeshWriteGmsh(HXTMesh * mesh,const char * filename)291 HXTStatus  hxtMeshWriteGmsh  ( HXTMesh* mesh , const char *filename) {
292   FILE* file = fopen(filename,"w");
293   if(file==NULL)
294     return HXT_ERROR_MSG(HXT_STATUS_FILE_CANNOT_BE_OPENED,
295       "Cannot open mesh file \"%s\"",(filename==NULL)?"(null)":filename);
296 
297   /* format for gmsh */
298   fprintf(file,"$MeshFormat\n"
299                "2.2 0 %u\n"
300                "$EndMeshFormat\n"
301                "$Nodes\n"
302                "%u\n",(unsigned) sizeof(double), mesh->vertices.num);
303 
304   { /* print the nodes */
305     uint32_t i;
306     for (i=0; i<mesh->vertices.num; i++)
307       fprintf(file,"%u %.10E %.10E %.10E\n",i+1, mesh->vertices.coord[(size_t) 4*i+0], mesh->vertices.coord[(size_t) 4*i+1], mesh->vertices.coord[(size_t) 4*i+2]);
308   }
309 
310   // count non-ghost tetrahedra:
311   uint64_t numRealTet = 0;
312   uint64_t i;
313   #pragma omp parallel for reduction(+:numRealTet)
314   for (i=0; i<mesh->tetrahedra.num; i++)
315   {
316     if(mesh->tetrahedra.node[i*4 + 3]!=UINT32_MAX){
317       uint32_t myColor = mesh->tetrahedra.color!=NULL ? mesh->tetrahedra.color[i] : 0;
318       // color = HXT_COLOR_OUT --> outside the domain
319       if (myColor != HXT_COLOR_OUT)
320         ++numRealTet;
321     }
322   }
323 
324   fprintf(file,"$EndNodes\n"
325           "$Elements\n"
326           "%" HXTu64 "\n",
327           numRealTet
328           + mesh->points.num
329           + mesh->lines.num
330           + mesh->triangles.num
331           );
332 
333   { /* print the elements */
334 
335     uint64_t index = 0;
336     for (i=0; i<mesh->points.num; i++){
337       fprintf(file,"%lu %u 2 0 %u %u\n", ++index,POINTID,
338               mesh->points.node[i]+1,
339               mesh->points.node[i]+1);
340     }
341     for (i=0; i<mesh->lines.num; i++){
342       uint32_t myColor = mesh->lines.color ? mesh->lines.color[i] : 0;
343       if (myColor != HXT_COLOR_OUT) {
344         fprintf(file, "%" HXTu64 " %u 2 0 %u %u %u\n", ++index,LINEID,
345                 myColor,
346                 mesh->lines.node[i*2]+1,
347                 mesh->lines.node[i*2 + 1]+1);
348       }
349     }
350     for (i=0; i<mesh->triangles.num; i++){
351       uint32_t myColor = mesh->triangles.color ? mesh->triangles.color[i] : 0;
352       if (myColor != HXT_COLOR_OUT) {
353         fprintf(file, "%" HXTu64 " %u 2 0 %u %u %u %u\n", ++index,TRIID,
354                 myColor,
355                 mesh->triangles.node[i*3]+1,
356                 mesh->triangles.node[i*3 + 1]+1,
357                 mesh->triangles.node[i*3 + 2]+1);
358       }
359     }
360     for (i=0; i<mesh->tetrahedra.num; i++){
361       if(mesh->tetrahedra.node[i*4 + 3]!=UINT32_MAX){
362         uint32_t myColor = mesh->tetrahedra.color ? mesh->tetrahedra.color[i] : 0;
363         if (myColor != HXT_COLOR_OUT) {
364           fprintf(file, "%" HXTu64 " %u 2 0 %u %u %u %u %u\n", ++index,TETID,
365               myColor,
366               mesh->tetrahedra.node[i*4]+1,
367               mesh->tetrahedra.node[i*4 + 1]+1,
368               mesh->tetrahedra.node[i*4 + 2]+1,
369               mesh->tetrahedra.node[i*4 + 3]+1);
370         }
371       }
372     }
373   }
374 
375   fputs("$EndElements\n",file);
376   fclose(file);
377   return HXT_STATUS_OK;
378 }
379 
380