mesh.c 21 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
/*
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
2
 * MigFlow - Copyright (C) <2010-2018>
3
4
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
Michel Henry's avatar
Michel Henry committed
5
 *
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
 * List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
 * Description and complete License: see LICENSE file.
Michel Henry's avatar
Michel Henry committed
8
9
10
 *
 * This program (MigFlow) is free software:
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
11
12
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
Michel Henry's avatar
Michel Henry committed
13
 *
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
14
15
16
 * This program 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
17
 * GNU Lesser General Public License for more details.
Michel Henry's avatar
Michel Henry committed
18
 *
19
 * You should have received a copy of the GNU Lesser General Public License
Michel Henry's avatar
Michel Henry committed
20
 * along with this program (see COPYING and COPYING.LESSER files).  If not,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
21
22
23
 * see <http://www.gnu.org/licenses/>.
 */

24
#include "mesh.h"
Michel Henry's avatar
Michel Henry committed
25
#include "gmsh_mesh.h"
26
27
28
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
29
#include "tools.h"
30

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
32
33
34
35
36
37
38
39
typedef struct {
  int n[DIMENSION];
  int element;
  int pos; // or tag if element == -1
}HalfEdge;

int hedgecmp(const void *pe0, const void *pe1) {
  int *e0 = ((HalfEdge*)pe0)->n;
  int *e1 = ((HalfEdge*)pe1)->n;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
40
  if (e0[0] != e1[0]) return e0[0]-e1[0];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
#if DIMENSION==2
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
42
  return e0[1]-e1[1];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
43
44
45
46
#else
  if (e0[1] != e1[1]) return e0[1]-e1[1];
  return e0[2]-e1[2];
#endif
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
47
48
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49
50
51
52
53
54
static void swapint(int *a, int *b) {
  int c = *a;
  *a = *b;
  *b = c;
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
55
56
57
58
59
60
61
62
63
64
65
static void sort_edge_nodes(int *n) {
  if(n[0]>n[1])
    swapint(n,n+1);
#if DIMENSION == 3
  if(n[0]>n[2])
    swapint(n,n+2);
  if(n[1]>n[2])
    swapint(n+1,n+2);
#endif
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66
void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
67
68
69
  int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements;
  HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges);
  HalfEdge *he = halfedges; // current
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
70
  int id = 0;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
71
  for (int i=0; i<m->n_elements; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
72
73
74
    for (int j = 0; j < DIMENSION+1; ++j) {
      for (int k = 0; k < DIMENSION; ++k) {
        he->n[k] = m->elements[i*(DIMENSION+1)+elbnd[j][k]];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
75
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76
77
78
79
      sort_edge_nodes(he->n);
      he->element = i;
      he->pos = j;
      he++;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
    }
  }
  for (int i=0; i<n_bnd_elements; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
83
    for (int k = 0; k < DIMENSION; ++k){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
84
      he->n[k] = bnd_elements[i*DIMENSION+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
85
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
86
87
    sort_edge_nodes(he->n);
    he->element = -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
    he->pos = bnd_tags[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89
    he ++;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
90
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
91
  qsort(halfedges,n_half_edges,sizeof(HalfEdge),hedgecmp);
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
92
93
  int n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
94
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)) n_edges++;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
95
96
97
98
  }
  int *edges = malloc(n_edges*4*sizeof(int));
  n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
99
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)){
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
100
      n_edges += 1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
102
      edges[(n_edges-1)*4+0] = halfedges[i].element;
      edges[(n_edges-1)*4+1] = halfedges[i].pos;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
103
104
105
106
      edges[(n_edges-1)*4+2] = -1;
      edges[(n_edges-1)*4+3] = -1;
    }
    else {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
107
108
      edges[(n_edges-1)*4+2] = halfedges[i].element;
      edges[(n_edges-1)*4+3] = halfedges[i].pos;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
109
110
111
112
113
    }
  }
  m->n_interfaces = n_edges;
  m->interfaces = edges;
  free(halfedges);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
114
  for (int i = 0; i < n_edges; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
115
    int *e = edges + i*4;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
    // ensure boundary tag is on the second part of the closure
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
117
118
119
    if (e[0] < 0) {
      swapint(&e[0],&e[2]);
      swapint(&e[1],&e[3]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
120
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
122
123
124
    // compute the closure id of the second element
    if (e[2] >= 0) {
#if DIMENSION==2
      e[3] += 3;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
125
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
126
127
128
129
130
      int k = 0;
      const int *e0 = &m->elements[e[0]*4];
      const int *e1 = &m->elements[e[2]*4];
      e[3] += 12;
      while (e0[elbnd[e[1]][0]] != e1[elbnd[e[3]][k]]) e[3] += 4;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
131
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
132
133
    }
  }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
134
135
}

136
137
138
static int check_word_in_file(FILE *f, const char *w) {
  char word[256];
  if(fscanf(f,"%255s", word) != 1){
139
    printf("Invalid mesh file (\"%s\" expected).\n", w);
140
141
142
    return -1;
  }
  if(strcmp(word, w) != 0){
143
    printf("Invalid mesh file (\"%s\" expected).\n", w);
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
    return -1;
  }
  return 0;
}

static int read_quoted_string(FILE *f, char *w, size_t maxl) {
  char c;
  size_t pos = 0;
  c = fgetc(f);
  while (c != '\"' && !feof(f)) {
    c = fgetc(f);
  }
  do {
    c = fgetc(f);
    if (c != '\"')
      w[pos++] = c;
  }while (c != '\"' && !feof(f) && pos < maxl);
  if (pos == maxl)
    return printf("string too long in mesh file\n");
  if(feof(f))
    return printf("string not terminated in mesh file\n");
  w[pos] = '\0';
  return 0;
}

static int int_cmp(const void *p0, const void *p1)
{
  int i0 = *(int*)p0;
  int i1 = *(int*)p1;
  if (i0 == i1) return 0;
  if (i0 < i1) return -1;
  return 1;
}

178
179
180
181
182
183
184
185
186
187
static void sortDi(int i[DIMENSION]) {
#if DIMENSION == 2
   if (i[0]>i[1]) {
     int n = i[0];
     i[0] = i[1];
     i[1] = n;
   }
#else
   int n[3] = {i[0],i[1],i[2]};
   if (n[0]>n[1]) {
Michel Henry's avatar
Michel Henry committed
188
189
      i[1]=n[0];
      i[0]=n[1];
190
   } else {
Michel Henry's avatar
Michel Henry committed
191
192
193
194
195
196
197
      i[1]=n[1];
      i[0]=n[0];
   }
   if (i[1]>n[2]) {
      i[2]=i[1];
      if(i[0]>n[2]){
         i[1]=i[0];
198
         i[0]=n[2];
Michel Henry's avatar
Michel Henry committed
199
200
      } else i[1]=n[2];
   } else i[2]=n[2];
201
202
#endif
}
Michel Henry's avatar
Michel Henry committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
Mesh *mesh_load_msh(const char *filename)
{
  GmshMesh* gmsh_m = gmsh_mesh_new();
  printf("Try to read !\n");
  gmsh_mesh_read(gmsh_m, filename);
  printf("Mesh read !\n");
  Mesh *m = malloc(sizeof(Mesh));
  m->interfaces = NULL;
  // Boundary
  printf("Add Boundary ?\t");
  m->n_boundaries = gmsh_m->n_physical_names;
  m->boundary_names = malloc(sizeof(char*)*m->n_boundaries);
  for(int i = 0; i < m->n_boundaries; ++i){
    m->boundary_names[i] = strdup(gmsh_m->physical_names[i]);
  }
  printf("Boundary stored\n");
  // Nodes
  printf("Add Nodes ?\t");
  m->n_nodes = 0;
  for (int d = 0; d < 4; d ++){
    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
      const Entity *e = gmsh_m->entities[d][ie];
      m->n_nodes += e->n_nodes;
    }
  }
  printf("n_nodes : %d\n", m->n_nodes);
  m->x = malloc(sizeof(double)*3*m->n_nodes);
  int i_nodes = 0;
  for (int d = 0; d < 4; ++d){
    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
      const Entity *e = gmsh_m->entities[d][ie];
      for(size_t i = 0; i < e->n_nodes; ++i){
        // e->x[i*3+0] e->x[i*3+1] e->x[i*3+2]
        //printf("Entity[%d][%d]:\t Nodes %d :\t (%f %f %f)\n",d,ie,i,e->x[i*3+0],e->x[i*3+1],e->x[i*3+2]);
        m->x[i_nodes+0] = e->x[i*3+0];
        m->x[i_nodes+1] = e->x[i*3+1];
        m->x[i_nodes+2] = e->x[i*3+2];
        i_nodes += 3;
      }
    }
  }
  for(int i = 0; i < m->n_nodes; ++i){
    printf("Nodes %d :\t (%f %f %f)\n",i,m->x[i*3+0],m->x[i*3+1],m->x[i*3+2]);
  }
  printf("Nodes stored\n");
  // elements
  printf("Add Elements ?\t");
  m->n_elements = 0;
  int n_nodes_by_element;
  for (int d = 0; d < 4; ++d){
    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
      const Entity *e = gmsh_m->entities[d][ie];
      m->n_elements += e->n_elements;
      n_nodes_by_element = e->n_nodes_by_element;
    }
  }
  printf("n_elements : %d\n", m->n_elements);
  m->elements = malloc(sizeof(int)*m->n_elements*n_nodes_by_element);
  int i_element = 0;
  for (int d = 0; d < 4; d ++){
    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
      const Entity *e = gmsh_m->entities[d][ie];
      for(size_t i = 0; i < e->n_elements; ++i){
        for(size_t j = 0; j < e->n_nodes_by_element; ++j){
          m->elements[i_element] = e->elements[i*e->n_nodes_by_element+j]; // FIXME: the gmsh_free does not work with this line
          printf("i_element = %d\n", i_element);
          i_element++;
        }
      }
    }
  }
  printf("Elements Stored\n");
  gmsh_mesh_free(gmsh_m);
  printf("Free GMSH Mesh Structure\n");
  return m;
}
/*
280
281
282
Mesh *mesh_load_msh(const char *filename)
{
  Mesh *m = malloc(sizeof(Mesh));
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
283
  m->interfaces = NULL;
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
  FILE *f = fopen(filename, "r");
  if (!f){
    printf("Cannot open file \"%s\".\n", filename);
    return NULL;
  }
  check_word_in_file(f,"$MeshFormat");
  float v0;
  int v1, v2;
  if(fscanf(f, "%g %d %d", &v0, &v1, &v2) != 3){
    printf("Cannot parse mesh file version\n");
    return NULL;
  }
  if ((int)10*v0 != 22){
    printf("Mesh file format is %g and not msh 2.2\n", v0);
    return NULL;
  }
  check_word_in_file(f,"$EndMeshFormat");
  check_word_in_file(f,"$PhysicalNames");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
302
303
  int nphysfile;
  if (fscanf(f, "%d", &nphysfile) != 1){
304
305
306
    printf("Cannot read physical names\n");
    return NULL;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
307
308
309
310
311
312
  m->boundary_names = malloc(sizeof(char*)*nphysfile);
  int *physid = malloc(sizeof(int)*nphysfile);
  int nphys = 0;
  for (int i = 0; i < nphysfile; ++i) {
    int physdim, pid;
    if (fscanf(f, "%i %i", &physdim, &pid) != 2){
313
314
315
316
317
      printf("Cannot read physical names\n");
      return NULL;
    }
    char pname[256];
    if(read_quoted_string(f, pname, 256) < 0) return NULL;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
318
319
320
321
322
    if (physdim == DIMENSION-1) {
      physid[nphys] = pid;
      m->boundary_names[nphys] = strdup(pname);
      nphys++;
    }
323
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
324
  m->n_boundaries = nphys;
325
326
327
328
329
330
331
  check_word_in_file(f,"$EndPhysicalNames");
  check_word_in_file(f,"$Nodes");
  if (fscanf(f, "%d", &m->n_nodes) != 1){
      printf("Cannot read physical nodes.\n");
      return NULL;
  }
  m->x = malloc(sizeof(double)*m->n_nodes*3);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
332

333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
  for (int i = 0; i < m->n_nodes; ++i) {
    int nid;
    if (fscanf(f, "%d %le %le %le", &nid, &m->x[i*3], &m->x[i*3+1], &m->x[i*3+2]) != 4){
      printf("Cannot read nodes.\n");
      return NULL;
    }
    if (nid != i+1){
      printf("Nodes are not sequentialy numbered.\n");
      return NULL;
    }
  }
  check_word_in_file(f,"$EndNodes");
  check_word_in_file(f,"$Elements");
  int n_el;
  if (fscanf(f, "%d", &n_el) != 1){
    printf("Cannot read elements.\n");
    return NULL;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
351
352
  int *bnd_elements = malloc(sizeof(int)*n_el*DIMENSION);
  int *bnd_tags = malloc(sizeof(int)*n_el);
353
  size_t n_bnd_elements = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
354
355
  m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_el);
  m->n_elements = 0;
356
357
  for (int i = 0; i < n_el; ++i) {
    int eid, etype, nflags, flag;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
358
    int physfile=0;
359
360
361
362
363
364
365
366
367
368
    if (fscanf(f, "%d %d %d", &eid, &etype, &nflags) != 3){
      printf("Cannot read elements.\n");
      return NULL;
    }
    for (int j = 0; j < nflags; ++j){
      if (fscanf(f, "%d", &flag) != 1){
        printf("Cannot read elements.\n");
        return NULL;
      }
      if (j == 0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
369
        physfile = flag;
370
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
371
    int n[DIMENSION+1];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
372
    int phys = -1;
373
374
    switch (etype) {
      case(1):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
375
        if (fscanf(f, "%d %d", &n[0], &n[1]) != 2){
376
377
378
          printf("Cannot read elements.\n");
          return NULL;
        }
379
#if DIMENSION==2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
380
381
382
383
384
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
385
386
387
          bnd_tags[n_bnd_elements] = phys;
          for (int iD = 0; iD < 2; ++iD) bnd_elements[n_bnd_elements*2+iD] = n[iD]-1;
          sortDi(&bnd_elements[n_bnd_elements*2]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
389
          n_bnd_elements++;
        }
390
#endif
391
392
        break;
      case(15):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
393
        if (fscanf(f, "%d", &n[0]) != 1){
394
395
396
397
398
          printf("Cannot read elements.\n");
          return NULL;
        }
        break;
      case(2):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
399
        if (fscanf(f, "%d %d %d", &n[0], &n[1], &n[2]) != 3){
400
401
402
          printf("Cannot read elements.\n");
          return NULL;
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
403
404
405
406
407
#if DIMENSION == 2
        m->elements[m->n_elements*3+0] = n[0]-1;
        m->elements[m->n_elements*3+1] = n[1]-1;
        m->elements[m->n_elements*3+2] = n[2]-1;
        m->n_elements++;
408
        break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
409
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
410
411
412
413
414
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
415
416
417
          bnd_tags[n_bnd_elements] = phys;
          for (int iD = 0; iD < 3; ++iD) bnd_elements[n_bnd_elements*3+iD] = n[iD]-1;
          sortDi(&bnd_elements[n_bnd_elements*3]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
418
419
          n_bnd_elements++;
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
420
421
422
423
424
425
426
427
428
429
430
431
432
        break;
      case(4):
        if (fscanf(f, "%d %d %d %d", &n[0], &n[1], &n[2], &n[3]) != 4){
          printf("Cannot read elements.\n");
          return NULL;
        }
        m->elements[m->n_elements*4+0] = n[0]-1;
        m->elements[m->n_elements*4+1] = n[1]-1;
        m->elements[m->n_elements*4+2] = n[2]-1;
        m->elements[m->n_elements*4+3] = n[3]-1;
        m->n_elements++;
        break;
#endif
433
434
435
436
      default:
        break;
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
437
  m->elements = realloc(m->elements, (DIMENSION+1)*sizeof(int)*m->n_elements);
438
  check_word_in_file(f,"$EndElements");
Michel Henry's avatar
Michel Henry committed
439

440
  fclose(f);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
441
  free(physid);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
442
  mesh_gen_edges(m,n_bnd_elements,bnd_elements,bnd_tags);
443
  free(bnd_elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
444
  free(bnd_tags);
445
446
  return m;
}
Michel Henry's avatar
Michel Henry committed
447
*/
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
448
449
450
451
452
453
454
455
456
457
int mesh_write_msh(Mesh *mesh, FILE *f)
{
  fprintf(f,"$MeshFormat\n2.2 0 0\n$EndMeshFormat\n");
  fprintf(f,"$Nodes\n");
  fprintf(f, "%d\n", mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%i %.16g %.16g %.16g\n", i+1, mesh->x[i*3], mesh->x[i*3+1], mesh->x[i*3+2]);
  }
  fprintf(f, "$EndNodes\n");
  fprintf(f, "$Elements\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
458
459
460
461
  fprintf(f, "%i\n", mesh->n_elements);
  for (int i= 0; i < mesh->n_elements; ++i) {
#if DIMENSION==2
    int *tri = mesh->elements+i*3;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
462
    fprintf(f, "%i 2 2 1 1 %i %i %i\n", i+1, tri[0]+1, tri[1]+1, tri[2]+1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
463
464
465
466
#else
    int *tet = mesh->elements+i*4;
    fprintf(f, "%i 4 2 1 1 %i %i %i %i\n", i+1, tet[0]+1, tet[1]+1, tet[2]+1, tet[3]+1);
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
468
469
470
471
  }
  fprintf(f, "$EndElements\n");
  return 0;
}

472
473
void mesh_free(Mesh *m)
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
474
  free(m->elements);
475
  free(m->x);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
476
477
  for (int i = 0; i < m->n_boundaries; ++i) {
    free(m->boundary_names[i]);
478
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
479
  free(m->boundary_names);
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
480
  free(m->interfaces);
481
482
483
  free(m);
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
484
int mesh_write_msh_scalar(const Mesh *mesh, FILE *f, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id)
485
486
487
488
489
490
491
492
493
494
495
496
{
  fprintf(f,"$NodeData\n");
  fprintf(f, "1\n\"%s\"\n",field_name);
  fprintf(f, "1\n%.16g\n",t);
  fprintf(f, "3\n%i\n%i\n%i\n", iter, 1, mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%i %.16g\n", i+1, solution[i*n_fields+field_id]);
  }
  fprintf(f, "$EndNodeData\n");
  return 0;
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
497
int mesh_write_msh_vector(const Mesh *mesh, FILE *f, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id[DIMENSION])
498
499
500
501
502
503
{
  fprintf(f,"$NodeData\n");
  fprintf(f, "1\n\"%s\"\n",field_name);
  fprintf(f, "1\n%.16g\n",t);
  fprintf(f, "3\n%i\n%i\n%i\n", iter, 3, mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
504
505
506
507
508
#if DIMENSION==2
    fprintf(f, "%i %.16g %.16g 0.\n", i+1, solution[i*n_fields+field_id[0]], solution[i*n_fields+field_id[1]]);
#else
    fprintf(f, "%i %.16g %.16g %.16g\n", i+1, solution[i*n_fields+field_id[0]], solution[i*n_fields+field_id[1]], solution[i*n_fields+field_id[2]]);
#endif
509
510
  }
  fprintf(f, "$EndNodeData\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
511
512
513
514
515
516
517
518
519
520
521
  return 0;
}

int mesh_write_pos_scalar(const Mesh *mesh, const char *dir_name, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id)
{
  char file_name[1024];
  sprintf(file_name,"%s/%s_%05i.pos",dir_name, field_name, iter);
  FILE *f = fopen(file_name, "w");
  if (!f)
    return printf("Cannot open file \"%s\" for writing.\n", file_name);
  fprintf(f,"View\"%s\"{\nTIME {%g};\n", field_name, t);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
522
523
524
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION==2
    int *tri = mesh->elements+iel*3;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
525
526
527
528
529
530
531
532
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    double V[3] = {solution[tri[0]*n_fields+field_id], solution[tri[1]*n_fields+field_id], solution[tri[2]*n_fields+field_id]};
    fprintf(f, "ST(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
        V[0],V[1],V[2]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
533
534
535
536
537
538
539
540
541
542
543
544
545
#else
    int *tet = mesh->elements+iel*4;
    const double x[4][3] = {
      {mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
      {mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
      {mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
      {mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]},
    };
    double V[4] = {solution[tet[0]*n_fields+field_id], solution[tet[1]*n_fields+field_id], solution[tet[2]*n_fields+field_id], solution[tet[3]*n_fields+field_id]};
    fprintf(f, "ST(%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g){%.8g, %.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],x[2][1],x[2][2],x[3][0],x[3][1],x[3][2],
        V[0],V[1],V[2],V[3]);
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
546
547
548
549
550
551
  }
  fprintf(f, "};\n");
  fclose(f);
  return 0;
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
552
int mesh_write_pos_vector(const Mesh *mesh, const char *dir_name, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id[DIMENSION])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
553
554
555
556
557
558
559
{
  char file_name[1024];
  sprintf(file_name,"%s/%s_%05i.pos",dir_name, field_name, iter);
  FILE *f = fopen(file_name, "w");
  if (!f)
    return printf("Cannot open file \"%s\" for writing.\n", file_name);
  fprintf(f,"View\"%s\"{\nTIME {%g};\n", field_name, t);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
560
561
562
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION ==2
    int *tri = mesh->elements+iel*3;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
563
564
565
566
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
567
568
    double V0[3] = {solution[tri[0]*n_fields+field_id[0]], solution[tri[1]*n_fields+field_id[0]], solution[tri[2]*n_fields+field_id[0]]};
    double V1[3] = {solution[tri[0]*n_fields+field_id[1]], solution[tri[1]*n_fields+field_id[1]], solution[tri[2]*n_fields+field_id[1]]};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
569
570
571
    fprintf(f, "VT(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0};\n",
        x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
        V0[0],V1[0], V0[1], V1[1], V0[2], V1[2]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
#else
    int *tet = mesh->elements+iel*4;
    const double x[4][3] = {
      {mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
      {mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
      {mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
      {mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]},
    };
    double V0[4] = {solution[tet[0]*n_fields+field_id[0]], solution[tet[1]*n_fields+field_id[0]], solution[tet[2]*n_fields+field_id[0]], solution[tet[3]*n_fields+field_id[0]]};
    double V1[4] = {solution[tet[0]*n_fields+field_id[1]], solution[tet[1]*n_fields+field_id[1]], solution[tet[2]*n_fields+field_id[1]], solution[tet[3]*n_fields+field_id[1]]};
    double V2[4] = {solution[tet[0]*n_fields+field_id[2]], solution[tet[1]*n_fields+field_id[2]], solution[tet[2]*n_fields+field_id[2]], solution[tet[3]*n_fields+field_id[2]]};
    fprintf(f, "VT(%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g){%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],x[2][1],x[2][2],x[3][0],x[3][1],x[3][2],
        V0[0],V1[0],V2[0], V0[1], V1[1],V2[1], V0[2], V1[2],V2[2], V0[3], V1[3],V2[3]);
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
587
588
  }
  fprintf(f, "};\n");
589
590
591
  fclose(f);
  return 0;
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650

static int intcmp(const void *p0, const void *p1) {
  return *(int*)p0 - *(int*)p1;
}

MeshBoundary *mesh_boundaries_create(const Mesh *m) {
  MeshBoundary *b = malloc(sizeof(MeshBoundary)*m->n_boundaries);
  for (int i = 0; i < m->n_boundaries; ++i) {
    b[i].n_nodes = 0;
    b[i].n_interfaces = 0;
  }
  for (int i = 0; i < m->n_interfaces; ++i) {
    int *interface = m->interfaces + i*4;
    if (interface[2] == -1) {
      b[interface[3]].n_interfaces++;
    }
  }
  for (int i = 0; i < m->n_boundaries; ++i) {
    b[i].interfaces = malloc(b[i].n_interfaces*sizeof(int));
    b[i].nodes = malloc(b[i].n_interfaces*sizeof(int)*DIMENSION);
    b[i].n_interfaces = 0;
  }
  for (int i = 0; i < m->n_interfaces; ++i) {
    int *interface = m->interfaces + i*4;
    if (interface[2] == -1) {
      int tag = interface[3];
      int *el = m->elements + interface[0]*(DIMENSION+1);
      int *cl = elbnd[interface[1]];
      b[tag].interfaces[b[tag].n_interfaces++] = i;
      for (int i = 0; i < DIMENSION; ++i) {
        b[tag].nodes[b[tag].n_nodes++] = el[cl[i]];
      }
    }
  }
  for (int i = 0; i < m->n_boundaries; ++i) {
    qsort(b[i].nodes,b[i].n_nodes,sizeof(int),intcmp);
    int end = b[i].n_nodes;
    b[i].n_nodes = 0;
    int p = -1;
    for (int j = 0; j < end; ++j) {
      if (b[i].nodes[j] == p){
        continue;
      }
      b[i].nodes[b[i].n_nodes++] = b[i].nodes[j];
      p = b[i].nodes[j];
    }
    b[i].nodes = realloc(b[i].nodes, b[i].n_nodes*sizeof(int));
  }
  return b;
}

void mesh_boundaries_free(MeshBoundary *bs, int n)
{
  for (int i =0; i < n; ++i) {
    free(bs[i].nodes);
    free(bs[i].interfaces);
  }
  free(bs);
}
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668


void mesh_set_elements(Mesh *m, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals) {
  m->n_elements = n_elements;
  m->elements = malloc(n_elements*sizeof(int)*(DIMENSION+1));
  memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);
  m->n_nodes = n_nodes;
  m->n_elements = n_elements;
  m->x = malloc(n_nodes*3*sizeof(double));
  memcpy(m->x,x,sizeof(double)*3*n_nodes);
  m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_elements);
  memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);
  m->n_boundaries = n_physicals;
  m->boundary_names = malloc(sizeof(char*)*n_physicals);
  for (int i = 0; i < n_physicals; ++i) {
    m->boundary_names[i] = strdup(physicals[i]);
  }
}