mesh.c

/*
 * MigFlow - Copyright (C) <2010-2018>
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
 *
 * List of the contributors to the development of MigFlow: see AUTHORS file.
 * Description and complete License: see LICENSE file.
 *
 * This program (MigFlow) is free software:
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
 *
 * 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
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program (see COPYING and COPYING.LESSER files).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include "mesh.h"
#include "gmsh_mesh.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "tools.h"

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;
  if (e0[0] != e1[0]) return e0[0]-e1[0];
#if DIMENSION==2
  return e0[1]-e1[1];
#else
  if (e0[1] != e1[1]) return e0[1]-e1[1];
  return e0[2]-e1[2];
#endif
}

static void swapint(int *a, int *b) {
  int c = *a;
  *a = *b;
  *b = c;
}

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
}

void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
  int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements;
  HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges);
  HalfEdge *he = halfedges; // current
  int id = 0;
  for (int i=0; i<m->n_elements; ++i) {
    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]];
      }
      sort_edge_nodes(he->n);
      he->element = i;
      he->pos = j;
      he++;
    }
  }
  for (int i=0; i<n_bnd_elements; ++i) {
    for (int k = 0; k < DIMENSION; ++k){
      he->n[k] = bnd_elements[i*DIMENSION+k];
    }
    sort_edge_nodes(he->n);
    he->element = -1;
    he->pos = bnd_tags[i];
    he ++;
  }
  qsort(halfedges,n_half_edges,sizeof(HalfEdge),hedgecmp);
  int n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)) n_edges++;
  }
  int *edges = malloc(n_edges*4*sizeof(int));
  n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)){
      n_edges += 1;
      edges[(n_edges-1)*4+0] = halfedges[i].element;
      edges[(n_edges-1)*4+1] = halfedges[i].pos;
      edges[(n_edges-1)*4+2] = -1;
      edges[(n_edges-1)*4+3] = -1;
    }
    else {
      edges[(n_edges-1)*4+2] = halfedges[i].element;
      edges[(n_edges-1)*4+3] = halfedges[i].pos;
    }
  }
  m->n_interfaces = n_edges;
  m->interfaces = edges;
  free(halfedges);
  for (int i = 0; i < n_edges; ++i) {
    int *e = edges + i*4;
    // ensure boundary tag is on the second part of the closure
    if (e[0] < 0) {
      swapint(&e[0],&e[2]);
      swapint(&e[1],&e[3]);
    }
    // compute the closure id of the second element
    if (e[2] >= 0) {
#if DIMENSION==2
      e[3] += 3;
#else
      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;
#endif
    }
  }
}

static int check_word_in_file(FILE *f, const char *w) {
  char word[256];
  if(fscanf(f,"%255s", word) != 1){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    return -1;
  }
  if(strcmp(word, w) != 0){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    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;
}

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]) {
      i[1]=n[0];
      i[0]=n[1];
   } else {
      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];
         i[0]=n[2];
      } else i[1]=n[2];
   } else i[2]=n[2];
#endif
}
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;
}
/*
Mesh *mesh_load_msh(const char *filename)
{
  Mesh *m = malloc(sizeof(Mesh));
  m->interfaces = NULL;
  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");
  int nphysfile;
  if (fscanf(f, "%d", &nphysfile) != 1){
    printf("Cannot read physical names\n");
    return NULL;
  }
  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){
      printf("Cannot read physical names\n");
      return NULL;
    }
    char pname[256];
    if(read_quoted_string(f, pname, 256) < 0) return NULL;
    if (physdim == DIMENSION-1) {
      physid[nphys] = pid;
      m->boundary_names[nphys] = strdup(pname);
      nphys++;
    }
  }
  m->n_boundaries = nphys;
  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);

  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;
  }
  int *bnd_elements = malloc(sizeof(int)*n_el*DIMENSION);
  int *bnd_tags = malloc(sizeof(int)*n_el);
  size_t n_bnd_elements = 0;
  m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_el);
  m->n_elements = 0;
  for (int i = 0; i < n_el; ++i) {
    int eid, etype, nflags, flag;
    int physfile=0;
    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)
        physfile = flag;
    }
    int n[DIMENSION+1];
    int phys = -1;
    switch (etype) {
      case(1):
        if (fscanf(f, "%d %d", &n[0], &n[1]) != 2){
          printf("Cannot read elements.\n");
          return NULL;
        }
#if DIMENSION==2
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
          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]);
          n_bnd_elements++;
        }
#endif
        break;
      case(15):
        if (fscanf(f, "%d", &n[0]) != 1){
          printf("Cannot read elements.\n");
          return NULL;
        }
        break;
      case(2):
        if (fscanf(f, "%d %d %d", &n[0], &n[1], &n[2]) != 3){
          printf("Cannot read elements.\n");
          return NULL;
        }
#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++;
        break;
#else
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
          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]);
          n_bnd_elements++;
        }
        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
      default:
        break;
    }
  }
  m->elements = realloc(m->elements, (DIMENSION+1)*sizeof(int)*m->n_elements);
  check_word_in_file(f,"$EndElements");

  fclose(f);
  free(physid);
  mesh_gen_edges(m,n_bnd_elements,bnd_elements,bnd_tags);
  free(bnd_elements);
  free(bnd_tags);
  return m;
}
*/
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");
  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;
    fprintf(f, "%i 2 2 1 1 %i %i %i\n", i+1, tri[0]+1, tri[1]+1, tri[2]+1);
#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
  }
  fprintf(f, "$EndElements\n");
  return 0;
}

void mesh_free(Mesh *m)
{
  free(m->elements);
  free(m->x);
  for (int i = 0; i < m->n_boundaries; ++i) {
    free(m->boundary_names[i]);
  }
  free(m->boundary_names);
  free(m->interfaces);
  free(m);
}

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)
{
  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;
}

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])
{
  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) {
#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
  }
  fprintf(f, "$EndNodeData\n");
  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);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION==2
    int *tri = mesh->elements+iel*3;
    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]);
#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
  }
  fprintf(f, "};\n");
  fclose(f);
  return 0;
}

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])
{
  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);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION ==2
    int *tri = mesh->elements+iel*3;
    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 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]]};
    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]);
#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
  }
  fprintf(f, "};\n");
  fclose(f);
  return 0;
}

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);
}


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]);
  }
}