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