star-cad

scad_geometry.c (94751B)

---

 /* Copyright (C) 2022-2024 |Méso|Star> (contact@meso-star.com)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU 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 General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "scad.h"
 #include "scad_c.h"
 #include "scad_device.h"
 #include "scad_geometry.h"
 
 #include <rsys/rsys.h>
 #include <rsys/mem_allocator.h>
 #include <rsys/str.h>
 #include <rsys/math.h>
 #include <rsys/double2.h>
 #include <rsys/double3.h>
 #include <rsys/logger.h>
 #include <rsys/hash_table.h>
 
 #include <stdlib.h>
 #include <string.h>
 #include <gmsh/gmshc.h>
 
 /* A type used to deduplicate tags */
 #define HTABLE_NAME tags
 #define HTABLE_KEY int
 #define HTABLE_DATA char
 #include <rsys/hash_table.h>
 
 /* A type used for mappings */
 #define HTABLE_NAME mappings
 #define HTABLE_KEY int
 #define HTABLE_DATA size_t
 #include <rsys/hash_table.h>
 
 #include <limits.h>
 
 /*******************************************************************************
  * Utility functions
  ******************************************************************************/
 struct coord_pair {
   const double* x;
   const double* y;
 };
 
 static res_T
 geom_set_name
   (struct scad_geometry* geom,
    const struct str* name) /* Can be NULL but not "" */
 {
   struct scad_device* dev = get_device();
   int same_name, name_isnt_null, previous_name_is_null;
   res_T res = RES_OK;
 
   ASSERT(geom);
   ASSERT(!name || !str_is_empty(name));
 
   previous_name_is_null = str_is_empty(&geom->name);
   if(name) {
     name_isnt_null = 1;
     same_name = (!previous_name_is_null) && (0 == str_cmp(name, &geom->name));
   } else {
     name_isnt_null = 0;
     same_name = previous_name_is_null;
   }
 
   if(same_name) /* No change needed */
     goto exit;
 
   if(name_isnt_null) {
     struct scad_geometry** g
       = htable_names_find(&dev->geometry_names, name);
     if(g) { /* if defined, names must be unique */
       res = RES_BAD_ARG;
       log_error(get_device(), "Geometry name '%s' is allready in use.\n",
         str_cget(name));
       goto error;
     }
   }
 
   if(!previous_name_is_null) {
     size_t n = htable_names_erase(&dev->geometry_names, &geom->name);
     ASSERT((n == 1)); (void)n;
   }
   if(name_isnt_null) {
     ERR(str_copy(&geom->name, name));
     ERR(htable_names_set(&dev->geometry_names, name, &geom));
   } else {
     str_clear(&geom->name);
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static int
 mixed_dim_err_msg
   (const char* op,
    struct scad_geometry** geometries,
    const size_t geometries_count,
    int* dim)
 {
   size_t i, j;
   struct scad_device* dev = get_device();
   ASSERT(op && (geometries || geometries_count == 0) && dim);
   for(i = 0; i < geometries_count; i++) {
     for(j = 0; j < geometries[i]->gmsh_dimTags_n; j += 2) {
       int d = geometries[i]->gmsh_dimTags[j];
       if(*dim == INT_MAX) *dim = d;
       else if (*dim != d) {
         log_error(dev,
             "Dimension mismatch in %s operation creating geometry.\n", op);
         return 1;
       }
     }
   }
   return 0;
 }
 
 static res_T
 geometry_create
   (struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = dev->allocator;
   char one = 1;
 
   ASSERT(out_geometry);
 
   geom = (struct scad_geometry*)MEM_CALLOC(allocator, 1, sizeof(*geom));
   if(!geom) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   ref_init(&geom->ref);
   str_init(allocator, &geom->name);
   ERR(htable_geometries_set(&dev->allgeom, &geom, &one));
 
   if(dev->options.Misc.LogRefCounting & SCAD_LOG_GEOMETRY) {
     if(str_is_empty(&geom->name)) {
       logger_print(dev->logger, dev->log_type,
           "Creating unnamed geometry %p (count set to 1).\n", (void*)geom);
     } else {
       logger_print(dev->logger, dev->log_type,
           "Creating geometry '%s' (%p, count set to 1).\n",
           str_cget(&geom->name), (void*)geom);
     }
   }
 
 end:
   *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto end;
 }
 
 enum tag_operation {
   COMMON_TAGS,
   UNIQUE_TAGS
 };
 
 static res_T
 process_tag_list
   (const int* dt1,
    const size_t c1,
    const int* dt2,
    const size_t c2,
    const enum tag_operation op,
    int** out_dimTags,
    size_t* out_dimTags_n)
 {
   res_T res = RES_OK;
   size_t i, d, c = 0, sz = 0;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = dev->allocator;
   struct htable_tags t2, t3;
   struct htable_tags_iterator it, end;
   int* dimTags = NULL;
   const char set1 = BIT(0), set2 = BIT(1), both = set1 | set2;
 
   ASSERT((dt1 || c1 == 0) && (dt2 || c2 == 0) && out_dimTags && out_dimTags_n);
 
   htable_tags_init(allocator, &t2);
   htable_tags_init(allocator, &t3);
 
   /* list tags and flags them whith set(s) they are from */
   for(i = 0; i < c1; i += 2) {
     int dim = dt1[i];
     int tag = dt1[i+1];
     struct htable_tags* tn = (dim == 2) ? &t2 : &t3;
     ASSERT(dim == 2 || dim == 3);
     ERR(htable_tags_set(tn, &tag, &set1));
   }
   for(i = 0; i < c2; i += 2) {
     char *fl, tmp;
     int dim = dt2[i];
     int tag = dt2[i+1];
     struct htable_tags* tn = (dim == 2) ? &t2 : &t3;
     ASSERT(dim == 2 || dim == 3);
     fl = htable_tags_find(tn, &tag);
     if(fl == NULL) {
       tmp = set2;
       fl = &tmp;
     } else {
       *fl |= set2;
     }
     ERR(htable_tags_set(tn, &tag, fl));
   }
 
   switch(op) {
     case COMMON_TAGS:
       for(d = 2; d <= 3; d++) {
         struct htable_tags* tn = (d == 2) ? &t2 : &t3;
         htable_tags_begin(tn, &it);
         htable_tags_end(tn, &end);
         while(!htable_tags_iterator_eq(&it, &end)) {
           if(*htable_tags_iterator_data_get(&it) == both) sz++;
           htable_tags_iterator_next(&it);
         }
       }
       break;
     case UNIQUE_TAGS:
       sz = htable_tags_size_get(&t2) + htable_tags_size_get(&t3);
       break;
     default: FATAL("Invalid enum value");
   }
 
   /* Build result */
   if(sz > 0) {
     dimTags = MEM_ALLOC(allocator, sz * 2 * sizeof(*dimTags));
     if(!dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
 
     for(d = 2; d <= 3; d++) {
       struct htable_tags* tn = (d == 2) ? &t2 : &t3;
       htable_tags_begin(tn, &it);
       htable_tags_end(tn, &end);
       while(!htable_tags_iterator_eq(&it, &end)) {
         int add;
         switch(op) {
           case COMMON_TAGS:
             add = (*htable_tags_iterator_data_get(&it) == both);
             break;
           case UNIQUE_TAGS:
             add = 1;
             break;
           default: FATAL("Invalid enum value");
         }
         if(add) {
           dimTags[c++] = 2;
           dimTags[c++] = *htable_tags_iterator_key_get(&it);
         }
         htable_tags_iterator_next(&it);
       }
     }
     ASSERT(sz*2 == c);
   }
 
 exit:
   *out_dimTags_n = c;
   *out_dimTags = dimTags;
   htable_tags_release(&t2);
   htable_tags_release(&t3);
   return res;
 error:
   MEM_RM(allocator, dimTags);
   dimTags = NULL;
   c = 0;
   goto exit;
 }
 
 static res_T
 gather_tags
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    int** out_dimTags,
    size_t* out_dimTags_n)
 {
   res_T res = RES_OK;
   int* dimTags = NULL;
   size_t i, j, c, sz;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = dev->allocator;
   struct htable_tags t2, t3;
   struct htable_tags_iterator it, end;
 
   ASSERT((geometries || !geometries_count) && out_dimTags && out_dimTags_n);
 
   htable_tags_init(allocator, &t2);
   htable_tags_init(allocator, &t3);
 
   /* list tags and remove duplicates */
   for(i = 0; i < geometries_count; i++) {
     if(!geometries[i]) {
       res = RES_BAD_ARG;
       goto error;
     }
     for(j = 0; j < geometries[i]->gmsh_dimTags_n; j += 2) {
       char one = 1;
       int dim = geometries[i]->gmsh_dimTags[j];
       int tag = geometries[i]->gmsh_dimTags[j+1];
       struct htable_tags* tn = (dim == 2) ? &t2 : &t3;
       ASSERT(dim == 2 || dim == 3);
       ERR(htable_tags_set(tn, &tag, &one));
     }
   }
 
   /* Build result */
   sz = htable_tags_size_get(&t2) + htable_tags_size_get(&t3);
   dimTags = MEM_ALLOC(allocator, sz * 2 * sizeof(*dimTags));
   if(!dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   c = 0;
   htable_tags_begin(&t2, &it);
   htable_tags_end(&t2, &end);
   while(!htable_tags_iterator_eq(&it, &end)) {
     dimTags[c++] = 2;
     dimTags[c++] = *htable_tags_iterator_key_get(&it);
     htable_tags_iterator_next(&it);
   }
   htable_tags_begin(&t3, &it);
   htable_tags_end(&t3, &end);
   while(!htable_tags_iterator_eq(&it, &end)) {
     dimTags[c++] = 3;
     dimTags[c++] = *htable_tags_iterator_key_get(&it);
     htable_tags_iterator_next(&it);
   }
   ASSERT(sz*2 == c);
 
   *out_dimTags_n = c;
   *out_dimTags = dimTags;
 
 exit:
   htable_tags_release(&t2);
   htable_tags_release(&t3);
   return res;
 error:
   MEM_RM(allocator, dimTags);
   goto exit;
 }
 
 static res_T
 store_tags
   (int* dimTags,
    size_t count,
    const int down_to_dim,
    struct htable_tags t[4])
 {
   res_T res = RES_OK;
   char one = 1;
   size_t i;
   ASSERT((dimTags || count == 0) && t);
   for(i = 0; i < count; i += 2) {
     int dim = dimTags[i];
     int tag = dimTags[i+1];
     ASSERT(dim >= 0 && dim <= 3);
     if(down_to_dim > dim) continue;
     ERR(htable_tags_set(t+dim, &tag, &one));
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* Recursivelly get the tags of geometries, possibly down to dim 0.
  * Trigger a call to sync_device() */
 static res_T
 gather_tags_recursive
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    const int down_to_dim,
    int* out_dimTags[4],
    size_t out_dimTags_n[4])
 {
   res_T res = RES_OK;
   int* dimTags[4] = { NULL, NULL, NULL, NULL }, *sub = NULL;
   size_t i, sz[4];
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = dev->allocator;
   struct htable_tags t[4];
   struct htable_tags_iterator it, end;
   int dim;
 
   ASSERT((geometries || geometries_count == 0) && out_dimTags && out_dimTags_n);
   ASSERT(0 <= down_to_dim && down_to_dim <= 3);
 
   for(i = (size_t)down_to_dim; i < 4; i++) {
     htable_tags_init(allocator, t+i);
   }
 
   /* list tags by dimension and remove duplicates */
   for(i = 0; i < geometries_count; i++) {
     if(!geometries[i]) {
       res = RES_BAD_ARG;
       goto error;
     }
     ERR(store_tags(geometries[i]->gmsh_dimTags, geometries[i]->gmsh_dimTags_n,
          down_to_dim, t));
   }
 
   /* Recursively build result by dimension and list constituents,
    * begining with dim==3 */
   ERR(sync_device());
   for(dim = 3; dim >= down_to_dim; dim--) {
     size_t c = 0;
     sz[dim] = 2 * htable_tags_size_get(t+dim);
     if(sz[dim] == 0)
       continue;
     dimTags[dim] = MEM_ALLOC(allocator, sz[dim] * sizeof(*dimTags));
     if(!dimTags[dim]) {
       res = RES_MEM_ERR;
       goto error;
     }
     htable_tags_begin(t+dim, &it);
     htable_tags_end(t+dim, &end);
     while(!htable_tags_iterator_eq(&it, &end)) {
       dimTags[dim][c++] = dim;
       dimTags[dim][c++] = *htable_tags_iterator_key_get(&it);
       htable_tags_iterator_next(&it);
     }
     if(dim > down_to_dim) {
       int ierr;
       size_t subn;
       gmshModelGetBoundary(dimTags[dim], sz[dim], &sub, &subn, 0, 0, 0, &ierr);
       ERR(gmsh_err_to_res_T(ierr));
       ERR(store_tags(sub, subn, down_to_dim, t));
       gmshFree(sub); sub = NULL;
     }
     ASSERT(sz[dim] == c);
   }
 
   for(i = (size_t)down_to_dim; i < 4; i++) {
     out_dimTags_n[i] = sz[i];
     out_dimTags[i] = dimTags[i];
   }
 
 exit:
   if(sub) gmshFree(sub);
   for(i = (size_t)down_to_dim; i < 4; i++) {
     htable_tags_release(t+i);
   }
   return res;
 error:
   for(i = (size_t)down_to_dim; i < 4; i++) {
     MEM_RM(allocator, dimTags[i]);
   }
   goto exit;
 }
 
 static double size_callback
   (int dim, int tag, double x, double y, double z, double lc, void* data_)
 {
   struct scad_device* dev = get_device();
   double *modifier =
     htable_size_modifiers_find(dev->size_modifiers_by_dim + dim, &tag);
   (void)x;(void)y;(void)z;(void)data_;
   if(!modifier)
     return lc; /* No modifier defined */
   if(*modifier < 0)
     return -(*modifier); /* Absolute size modifier */
   return (*modifier) * lc; /* Size factor modifier */
 }
 
 /* gmsh documentation states that memory allocated by gmsh should be freed using
  * gmshFree.
  * According to valgrind map results as allocated by gmsh are not fully freed if
  * simply freed using gmshFree.
  * A code review shows that maps content is recursivelly allocated.
  * We end writting this helper function that free maps content and makes
  * valgrind happy. */
 static void
 free_gmsh_map
   (int** map,
    size_t mapnn)
 {
   size_t i;
   if(!mapnn) return;
   for(i = 0; i < mapnn; i++) {
     gmshFree(map[i]);
   }
   gmshFree(map);
 }
 
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
 static void
 geometry_release(ref_T* ref)
 {
   struct scad_geometry* geom;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = dev->allocator;
   size_t n;
   ASSERT(ref);
 
   geom = CONTAINER_OF(ref, struct scad_geometry, ref);
   CHK(RES_OK == device_unregister_tags(dev->log, dev->log_type, geom));
   MEM_RM(allocator, geom->gmsh_dimTags);
   if(!str_is_empty(&geom->name)) {
     n = htable_names_erase(&dev->geometry_names, &geom->name);
     ASSERT(n == 1);
   }
   str_release(&geom->name);
   n = htable_geometries_erase(&dev->allgeom, &geom);
   ASSERT(n == 1); (void)n;
   if(geom->release && geom->custom) {
     (*geom->release)(geom->custom);
   }
   MEM_RM(allocator, geom);
 }
 
 /******************************************************************************
  * Exported functions
  *****************************************************************************/
 res_T
 scad_geometry_ref_get
   (struct scad_geometry* geom)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
 
   if(!geom) return RES_BAD_ARG;
   ERR(check_device(FUNC_NAME));
 
   ref_get(&geom->ref);
   if(dev->options.Misc.LogRefCounting & SCAD_LOG_GEOMETRY) {
     if(str_is_empty(&geom->name)) {
       logger_print(dev->logger, dev->log_type,
           "Getting a reference on unnamed geometry %p (count set to %lu).\n",
           (void*)geom, (long unsigned)geom->ref);
     } else {
       logger_print(dev->logger, dev->log_type,
           "Getting a reference on geometry '%s' (count set to %lu).\n",
           str_cget(&geom->name), (long unsigned)geom->ref);
     }
   }
 
 end:
   return res;
 error:
   goto end;
 }
 
 res_T
 scad_geometry_ref_put
   (struct scad_geometry* geom)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
 
   if(!geom) return RES_BAD_ARG;
   ERR(check_device(FUNC_NAME));
 
   if(dev->options.Misc.LogRefCounting & SCAD_LOG_GEOMETRY) {
     if(str_is_empty(&geom->name)) {
       logger_print(dev->logger, dev->log_type,
           "Putting a reference on unnamed geometry %p (count set to %lu).\n",
           (void*)geom, (long unsigned)geom->ref - 1);
     } else {
       logger_print(dev->logger, dev->log_type,
           "Putting a reference on geometry '%s' (count set to %lu).\n",
           str_cget(&geom->name), (long unsigned)geom->ref - 1);
     }
   }
   ref_put(&geom->ref, geometry_release);
 
 end:
   return res;
 error:
   goto end;
 }
 
 res_T
 scad_scene_count
   (size_t* count)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
 
   if(!count) {
     res =RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   *count = htable_geometries_size_get(&dev->allgeom);
 
 end:
   return res;
 error:
   goto end;
 }
 
 res_T
 scad_scene_clear
   (void)
 {
   res_T res = RES_OK;
   int ierr;
   struct htable_geometries tmp;
   struct htable_geometries_iterator it, end;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
   int log, empty;
   enum scad_log_refcounting option;
   enum log_type log_type;
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
   option = dev->options.Misc.LogRefCounting;
 
   htable_geometries_init(allocator, &tmp);
   ERR(htable_geometries_copy(&tmp, &dev->allgeom));
   htable_geometries_begin(&tmp, &it);
   htable_geometries_end(&tmp, &end);
   empty = htable_geometries_is_empty(&dev->allgeom);
   log_type = empty ? LOG_OUTPUT : LOG_WARNING;
   log = (option & SCAD_LOG_DIMTAGS_ALL)
     || (!empty && (option & SCAD_LOG_DIMTAGS_ONLY_UNDELETED));
   SWAP(int, dev->log, log);
   SWAP(enum log_type, dev->log_type, log_type);
   if(dev->log) {
     logger_print(dev->logger, dev->log_type, "Clearing scene.\n");
     if(empty) {
       logger_print(dev->logger, dev->log_type, "scene is empty.\n");
     }
   }
   while(!htable_geometries_iterator_eq(&it, &end)) {
     struct scad_geometry* geom = *htable_geometries_iterator_key_get(&it);
     long cpt = geom->ref;
     while(cpt-- > 0) {
       ERR(scad_geometry_ref_put(geom));
     }
     htable_geometries_iterator_next(&it);
   }
   if(dev->log) {
     logger_print(dev->logger, dev->log_type, "End clearing scene.\n");
   }
   SWAP(int, log, dev->log);
   SWAP(enum log_type, log_type, dev->log_type);
 
   /* Ensure clear is complete (not scad-registered tags) */
   gmshClear(&ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 end:
   if(allocator) htable_geometries_release(&tmp);
   return res;
 error:
   goto end;
 }
 
 res_T
 scad_geometry_get_count
   (const struct scad_geometry* geom,
    size_t* count)
 {
   res_T res = RES_OK;
 
   if(!geom || !count) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ASSERT(geom->gmsh_dimTags_n % 2 == 0);
   *count = geom->gmsh_dimTags_n / 2;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_is_empty
   (const struct scad_geometry* geom,
    int* empty)
 {
   res_T res = RES_OK;
 
   if(!geom || !empty) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   *empty = (geom->gmsh_dimTags_n == 0);
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_set_custom_data
   (struct scad_geometry* geom,
    void (*release) (void* data),
    void* data)
 {
   res_T res = RES_OK;
 
   if(!geom) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   geom->custom = data;
   geom->release = release;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_get_custom_data
   (struct scad_geometry* geom,
    void** data)
 {
   res_T res = RES_OK;
 
   if(!geom || !data) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   if(geom->release && geom->custom) {
     (*geom->release)(geom->custom);
   }
   *data = geom->custom;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_set_visibility
   (const struct scad_geometry* geom,
    int visible,
    int recursive)
 {
   res_T res = RES_OK;
   int ierr;
   int* data = NULL;
   size_t sz;
 
   if(!geom) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   data = geom->gmsh_dimTags;
   sz = geom->gmsh_dimTags_n;
 
   ERR(sync_device());
   gmshModelSetVisibility(data, sz, visible, recursive, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometries_clear_mesh
   (struct scad_geometry** geometries,
    const size_t geometries_count)
 {
   res_T res = RES_OK;
   int ierr, dim;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
   size_t sz[4];
   int *data[4] = { NULL, NULL, NULL, NULL};
 
   if(!geometries || geometries_count == 0) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags_recursive(geometries, geometries_count, 0, data, sz));
   for(dim = 3; dim >= 0; dim--) {
     /* Cannot clear the mesh of lower dim entities if linked to higher dim
      * entities with uncleared mesh: start from the higher dim down to 0 */
     gmshModelMeshClear(data[dim], sz[dim], &ierr);
     ERR(gmsh_err_to_res_T(ierr));
   }
 
 exit:
   if(allocator) {
     for(dim = 3; dim >= 0; dim--) {
       MEM_RM(allocator, data[dim]);
     }
   }
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_get_name
   (const struct scad_geometry* geom,
    const char** name)
 {
   res_T res = RES_OK;
 
   if(!geom || !name) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   if(str_is_empty(&geom->name)) {
     *name = NULL;
   } else {
     *name = str_cget(&geom->name);
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_get_mass
   (struct scad_geometry* geom,
    double* mass)
 {
   res_T res = RES_OK;
   int ref_dim = 0;
   size_t i;
   int* data = NULL;
   size_t sz = 0;
   double m = 0;
 
   if(!geom || !mass) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   data = geom->gmsh_dimTags;
   sz = geom->gmsh_dimTags_n;
 
   ref_dim = data[0];
   for(i = 0; i < sz; i += 2) {
     double geom_mass;
     int dim = data[i], tag = data[i+1];
     int ierr = 0;
     if(ref_dim != dim) goto error;
     gmshModelOccGetMass(dim, tag, &geom_mass, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
     m += geom_mass;
   }
   *mass = m;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_get_closest_point
   (struct scad_geometry* geom,
    const double from[3],
    double closest[3],
    double* closest_distance,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   size_t i = 0;
   double* coord = NULL;
   double* pcoord = NULL;
   double tmp[3], min[3], min_d = DBL_MAX;
   struct scad_device* dev = get_device();
   struct darray_int tags;
   const int* data = NULL;
   size_t sz = 0;
   int initialized = 0, min_tag = -1;
   struct scad_geometry* out = NULL;
 
   if(!geom || !from || !closest || !closest_distance) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   darray_int_init(dev->allocator, &tags);
   initialized = 1;
 
   ERR(get_2d_tags(geom, &tags));
   data = darray_int_cdata_get(&tags);
   sz = darray_int_size_get(&tags);
 
   for(i = 0; i < sz; i++) {
     double d;
     int ierr = 0;
     int dim = 2;
     int tag = data[i];
     size_t pcoord_n;
     size_t coord_n;
 
     gmshModelGetClosestPoint(dim, tag, from, 3, &coord, &coord_n, &pcoord,
         &pcoord_n, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
     ASSERT(pcoord_n == (size_t)dim);
     ASSERT(coord_n == 3);
     d = d3_len(d3_sub(tmp, from, coord));
     if(d < min_d) {
       min_tag = tag;
       min_d = d;
       d3_set(min, coord);
     }
     gmshFree(coord);
     gmshFree(pcoord);
     coord = pcoord = NULL;
   }
   ASSERT(min_tag != -1);
   d3_set(closest, min);
   *closest_distance = min_d;
   if(out_geometry) {
     ERR(geometry_create(&out));
     out->gmsh_dimTags =
       MEM_ALLOC(dev->allocator, 2 * sizeof(*out->gmsh_dimTags));
     if(!out->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     out->gmsh_dimTags_n = 2;
     out->gmsh_dimTags[0] = 2;
     out->gmsh_dimTags[1] = min_tag;
     ERR(device_register_tags(out));
   }
 
 exit:
   if(initialized) darray_int_release(&tags);
   if(out_geometry) *out_geometry = out;
   gmshFree(coord);
   gmshFree(pcoord);
   return res;
 error:
   if(out) {
     SCAD(geometry_ref_put(out));
     out = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_geometry_get_centerofmass
   (struct scad_geometry* geom,
    double center[3])
 {
   res_T res = RES_OK;
   double *centers = NULL, *masses = NULL;
   size_t i = 0, count;
   struct mem_allocator* allocator;
   int ierr;
 
   if(!geom || !center) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ASSERT(geom->gmsh_dimTags_n % 2 == 0);
   count = geom->gmsh_dimTags_n / 2;
   if(count == 1) {
     int dim = geom->gmsh_dimTags[i];
     int tag = geom->gmsh_dimTags[i + 1];
     gmshModelOccGetCenterOfMass(dim, tag, center, center+1, center+2, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
   } else {
     double c[3] = { 0, 0, 0 }, tmp[3], m = 0;
     allocator = get_device()->allocator;
     centers = MEM_ALLOC(allocator, count * 3 * sizeof(*centers));
     masses = MEM_ALLOC(allocator, count * sizeof(*masses));
     if(!centers || !masses) {
       res = RES_MEM_ERR;
       goto error;
     }
     for(i = 0; i < count; i++) {
       int dim = geom->gmsh_dimTags[2*i];
       int tag = geom->gmsh_dimTags[2*i + 1];
       gmshModelOccGetCenterOfMass(dim, tag,
           centers+3*i, centers+3*i+1, centers+3*i+2, &ierr);
       ERR(gmsh_err_to_res_T(ierr));
       gmshModelOccGetMass(dim, tag, masses+i, &ierr);
       ERR(gmsh_err_to_res_T(ierr));
       d3_add(c, c, d3_muld(tmp, centers+3*i, masses[i]));
       m += masses[i];
     }
     d3_muld(center, c, 1/m);
   }
 
 exit:
   if(centers) MEM_RM(allocator, centers);
   if(masses) MEM_RM(allocator, masses);
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_get_bounding_box
   (struct scad_geometry* geom,
    double min[3],
    double max[3])
 {
   res_T res = RES_OK;
   size_t n = 0;
   double mi[3], ma[3];
 
   if(!geom || !min || !max) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   d3_splat(mi, DBL_MAX);
   d3_splat(ma, -DBL_MAX);
   for(n = 0; n < geom->gmsh_dimTags_n; n += 2) {
     double i[3], a[3];
     int ierr = 0;
     int dim = geom->gmsh_dimTags[n];
     int tag = geom->gmsh_dimTags[n + 1];
     gmshModelOccGetBoundingBox(dim, tag, i, i+1, i+2, a, a+1, a+2, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
     d3_min(mi, mi, i);
     d3_max(ma, ma, a);
   }
   d3_set(min, mi);
   d3_set(max, ma);
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_add_rectangle
   (const double xyz[3],
    const double dxdy[2],
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!xyz || !dxdy || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID = gmshModelOccAddRectangle(SPLIT3(xyz), SPLIT2(dxdy), -1, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags
     = MEM_ALLOC(allocator, geom->gmsh_dimTags_n * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 2;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_add_disk
   (const double xyz[3],
    const double radius,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!xyz || radius <= 0 || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID = gmshModelOccAddDisk(SPLIT3(xyz), radius, radius, -1, NULL, 0,
       NULL, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags
     = MEM_ALLOC(allocator, geom->gmsh_dimTags_n * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 2;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_add_polygon
   (void (*get_position)(const size_t ivert, double pos[2], void* data),
    void* data,
    const double z,
    const size_t count,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   size_t i;
   int* points = NULL;
   int* lines = NULL;
   int loop;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!get_position || count < 3 || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   points = MEM_ALLOC(allocator, count * sizeof(*points));
   lines = MEM_ALLOC(allocator, count * sizeof(*lines));
   if(!points || !lines) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   for(i=0; i<count; ++i) {
     double pos[2];
     get_position(i, pos, data);
     points[i] = gmshModelOccAddPoint(pos[0], pos[1], z, -1, -1, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
   }
 
   for(i=0; i<count; i++) {
     size_t end = (i == count-1) ? 0 : i+1;
     lines[i] = gmshModelOccAddLine(points[i], points[end], -1, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
   }
 
   loop = gmshModelOccAddCurveLoop(lines, count, -1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   gmsh_ID = gmshModelOccAddPlaneSurface(&loop, 1, -1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 2;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   if(allocator) {
     MEM_RM(allocator, points);
     MEM_RM(allocator, lines);
   }
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_add_box
   (const double xyz[3],
    const double dxdydz[3],
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!xyz || !dxdydz || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID = gmshModelOccAddBox(SPLIT3(xyz), SPLIT3(dxdydz), -1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 3;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_add_cylinder
   (const double xyz[3],
    const double axis[3],
    const double radius,
    const double angle,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!xyz || !axis || radius <= 0 || angle < 0 || angle > 2*PI || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID = gmshModelOccAddCylinder(SPLIT3(xyz), SPLIT3(axis), radius, -1,
       angle, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 3;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 int
 scad_add_sphere
   (const double xyz[3],
    const double radius,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!xyz || radius <= 0 || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID =
     gmshModelOccAddSphere(SPLIT3(xyz), radius, -1, -PI/2, PI/2, 2*PI, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 3;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 SCAD_API res_T
 scad_add_cone
   (const double center[3],
    const double axis[3],
    const double radius1,
    const double radius2,
    const double angle,
    struct scad_geometry** cone)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!center || !axis || radius1 < 0 || radius2 < 0 || (radius1 == radius2)
       || angle < 0 || angle > 2*PI || !cone)
   {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID =
     gmshModelOccAddCone(SPLIT3(center), SPLIT3(axis), radius1, radius2,
         -1, 2*PI, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 3;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(cone) *cone = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 SCAD_API res_T
 scad_add_torus
   (const double center[3],
    const double radius1,
    const double radius2,
    const double angle,
    const double z_axis[3], /* Can be NULL */
    struct scad_geometry** torus)
 {
   res_T res = RES_OK;
   int ierr, gmsh_ID;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!center || radius1 <= 0 || radius2 <= 0 || angle < 0 || angle > 2*PI
       || !torus)
   {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmsh_ID =
     gmshModelOccAddTorus(SPLIT3(center), radius1, radius2,
         -1, 2*PI, z_axis, (z_axis == NULL ? 0 : 3), &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = 2;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*geom->gmsh_dimTags));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   geom->gmsh_dimTags[0] = 3;
   geom->gmsh_dimTags[1] = gmsh_ID;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(torus) *torus = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 SCAD_API res_T
 scad_geometry_equal
   (const struct scad_geometry* geom1,
    const struct scad_geometry* geom2,
    int* equal)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
   struct htable_tags t2, t3;
   int  eq = 1, initialized = 0;
 
   if(!geom1 || !geom2 || !equal) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   /* Trivial cases */
   if(geom1 == geom2) {
     eq = 1;
   }
   else if(geom1->gmsh_dimTags_n != geom2->gmsh_dimTags_n) {
     eq = 0;
   } else {
     size_t i;
     htable_tags_init(allocator, &t2);
     htable_tags_init(allocator, &t3);
     initialized = 1;
     /* Create tables from geom1 tags */
     for(i = 0; i < geom1->gmsh_dimTags_n; i += 2) {
       char one = 1;
       int d = geom1->gmsh_dimTags[i];
       int t = geom1->gmsh_dimTags[i+1];
       struct htable_tags* tn = (d == 2) ? &t2 : &t3;
       ERR(htable_tags_set(tn, &t, &one));
     }
     ASSERT((htable_tags_size_get(&t2) + htable_tags_size_get(&t3)) * 2
         == geom1->gmsh_dimTags_n);
     /* Check if tags from geom2 are included */
     for(i = 0; i < geom2->gmsh_dimTags_n; i += 2) {
       char* found;
       int d = geom2->gmsh_dimTags[i];
       int t = geom2->gmsh_dimTags[i+1];
       struct htable_tags* tn = (d == 2) ? &t2 : &t3;
       found = htable_tags_find(tn, &t);
       if(!found) {
         eq = 0;
         break;
       }
     }
   }
 
   *equal = eq;
 
 exit:
   if(initialized) {
     htable_tags_release(&t2);
     htable_tags_release(&t3);
   }
   return res;
 error:
   goto exit;
 }
 
 SCAD_API res_T
 scad_geometry_is_included
   (const struct scad_geometry* geometry,
    struct scad_geometry** geometries,
    const size_t geometries_count,
    int* included)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
   struct htable_tags t2, t3;
   int initialized = 0;
   size_t i, n;
 
   if(!geometry || !geometries || geometries_count == 0 || !included) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   /* Trivial case */
   for(i = 0; i < geometries_count; i++) {
     if(geometry == geometries[i]) {
       *included = 1;
       goto exit;
     }
   }
 
   /* Create tables from geometries tags */
   htable_tags_init(allocator, &t2);
   htable_tags_init(allocator, &t3);
   initialized = 1;
   for(n = 0; n < geometries_count; n++) {
     const struct scad_geometry* geom = geometries[n];
     for(i = 0; i < geometries[n]->gmsh_dimTags_n; i += 2) {
       char one = 1;
       int d = geom->gmsh_dimTags[i];
       int t = geom->gmsh_dimTags[i+1];
       struct htable_tags* tn = (d == 2) ? &t2 : &t3;
       ERR(htable_tags_set(tn, &t, &one));
     }
   }
 
   /* Check if tags from geometry are included */
   for(i = 0; i < geometry->gmsh_dimTags_n; i += 2) {
     char* found;
     int d = geometry->gmsh_dimTags[i];
     int t = geometry->gmsh_dimTags[i+1];
     struct htable_tags* tn = (d == 2) ? &t2 : &t3;
     found = htable_tags_find(tn, &t);
     if(!found) {
       *included = 0;
       goto exit;
     }
   }
 
   /* If here, no not-included tag was found */
   *included = 1;
 
 exit:
   if(initialized) {
     htable_tags_release(&t2);
     htable_tags_release(&t3);
   }
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometries_fuse
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** tools,
    const size_t tools_count,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   int** map = NULL;
   size_t* mapn = NULL;
   size_t tagoutn, mapnn = 0, sz1, sz2;
   int* data1 = NULL;
   int* data2 = NULL;
   int ierr = 0;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   int *tmp1 = NULL, *tmp2 = NULL;
   size_t c1, c2;
 #endif
 
   if(!geometries || !geometries_count || !tools || !tools_count || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data1, &sz1));
   ERR(gather_tags(tools, tools_count, &data2, &sz2));
 
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshModelOccCopy(data1, sz1, &tmp1, &c1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   gmshModelOccCopy(data2, sz2, &tmp2, &c2, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   /* We do remove gmsh objects here as they are temporary copies created to
    * allow the remove flag to be used (gmsh is broken here if remove=0) */
   gmshModelOccFuse(tmp1, c1, tmp2, c2, &tagout, &tagoutn, &map, &mapn, &mapnn,
       -1, 1, 1, &ierr);
 #else
   /* We don't remove gmsh objects here; they are only removed when their tags are
    * no longer used by any star-cad geometry */
   gmshModelOccFuse(data1, sz1, data2, sz2, &tagout, &tagoutn, &map, &mapn,
       &mapnn, -1, 0, 0, &ierr);
 #endif
 
   ASSERT(tagoutn % 2 == 0);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = tagoutn;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   memcpy(geom->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
 
   ERR(device_register_tags(geom));
 
 exit:
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshFree(tmp1);
   gmshFree(tmp2);
 #endif
   if(out_geometry) *out_geometry = geom;
   if(allocator) {
     MEM_RM(allocator, data1);
     MEM_RM(allocator, data2);
   }
   gmshFree(mapn);
   gmshFree(tagout);
   free_gmsh_map(map, mapnn);
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometries_collect
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   size_t sz;
   int* data = NULL;
   struct scad_geometry* geom = NULL;
 
   if(!geometries || !geometries_count || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(gather_tags(geometries, geometries_count, &data, &sz));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = sz;
   geom->gmsh_dimTags = data;
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_geometries_cut
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** tools,
    const size_t tools_count,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   int** map = NULL;
   size_t* mapn = NULL;
   size_t tagoutn, mapnn = 0, sz1, sz2;
   int* data1 = NULL;
   int* data2 = NULL;
   int ierr = 0;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   int *tmp1 = NULL, *tmp2 = NULL;
   size_t c1, c2;
 #endif
 
   if(!geometries || !geometries_count || !tools || !tools_count || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data1, &sz1));
   ERR(gather_tags(tools, tools_count, &data2, &sz2));
 
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshModelOccCopy(data1, sz1, &tmp1, &c1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   gmshModelOccCopy(data2, sz2, &tmp2, &c2, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   /* We do remove gmsh objects here as they are temporary copies created to
    * allow the remove flag to be used (gmsh is broken here if remove=0) */
   gmshModelOccCut(tmp1, c1, tmp2, c2, &tagout, &tagoutn, &map, &mapn,
     &mapnn, -1, 1, 1, &ierr);
 #else
   /* We don't remove gmsh objects here; they are only removed when their tags are
    * no longer used by any star-cad geometry */
   gmshModelOccCut(data1, sz1, data2, sz2, &tagout, &tagoutn, &map, &mapn,
       &mapnn, -1, 0, 0, &ierr);
 #endif
 
   ASSERT(tagoutn % 2 == 0);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = tagoutn;
   geom->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
   if(!geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   memcpy(geom->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   if(allocator) {
     MEM_RM(allocator, data1);
     MEM_RM(allocator, data2);
   }
   gmshFree(mapn);
   gmshFree(tagout);
   free_gmsh_map(map, mapnn);
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshFree(tmp1);
   gmshFree(tmp2);
 #endif
   return res;
 error:
   if(ierr) {
     int dim = INT_MAX;
     if(!mixed_dim_err_msg("cut", geometries, geometries_count, &dim))
       mixed_dim_err_msg("cut", tools, tools_count, &dim);
     if(geom) {
       SCAD(geometry_ref_put(geom));
       geom = NULL;
     }
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometries_intersect
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** tools,
    const size_t tools_count,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   int** map = NULL;
   size_t* mapn = NULL;
   size_t tagoutn, mapnn = 0, sz1, sz2;
   int* data1 = NULL;
   int* data2 = NULL;
   int ierr = 0;
   struct scad_geometry* geom = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   int *tmp1 = NULL, *tmp2 = NULL;
   size_t c1, c2;
 #endif
 
   if(!geometries || !geometries_count || !tools || !tools_count || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data1, &sz1));
   ERR(gather_tags(tools, tools_count, &data2, &sz2));
 
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshModelOccCopy(data1, sz1, &tmp1, &c1, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   gmshModelOccCopy(data2, sz2, &tmp2, &c2, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   /* We do remove gmsh objects here as they are temporary copies created to
    * allow the remove flag to be used (gmsh is broken here if remove=0) */
   gmshModelOccIntersect(tmp1, c1, tmp2, c2, &tagout, &tagoutn, &map,
     &mapn, &mapnn, -1, 1, 1, &ierr);
 #else
   /* We don't remove gmsh objects here; they are only removed when their tags are
    * no longer used by any star-cad geometry */
   gmshModelOccIntersect(data1, sz1, data2, sz2, &tagout, &tagoutn, &map, &mapn,
       &mapnn, -1, 0, 0, &ierr);
 #endif
 
   ASSERT(tagoutn % 2 == 0);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&geom));
   geom->gmsh_dimTags_n = tagoutn;
   if(tagoutn != 0){
     geom->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
     if(!geom->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     memcpy(geom->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
   }
 
   ERR(device_register_tags(geom));
 
 exit:
   if(out_geometry) *out_geometry = geom;
   if(allocator) {
     MEM_RM(allocator, data1);
     MEM_RM(allocator, data2);
   }
   gmshFree(mapn);
   gmshFree(tagout);
   free_gmsh_map(map, mapnn);
 #ifdef FIX_GMSH_BOOLEAN_OUTPUTS
   gmshFree(tmp1);
   gmshFree(tmp2);
 #endif
   return res;
 error:
   if(geom) {
     SCAD(geometry_ref_put(geom));
     geom = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometries_common_boundaries
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** tools,
    const size_t tools_count,
    struct scad_geometry*** out_boundaries,
    size_t *out_count)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   size_t tagoutn, sz1, sz2, u_sz;
   int* data1 = NULL;
   int* data2 = NULL;
   int* unique = NULL;
   int ierr = 0;
   int* bound1 = NULL;
   int* bound2 = NULL;
   size_t n1, n2;
   struct scad_geometry** out_geom = NULL;
   struct mem_allocator* allocator = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   size_t i, c = 0, n;
   struct str msg;
   int msg_initialized = 0;
 
   if(!geometries || !geometries_count || !tools || !tools_count
       || !out_boundaries || !out_count)
   {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data1, &sz1));
   ERR(gather_tags(tools, tools_count, &data2, &sz2));
 
   /* data1 and data2 can have tags in common: deduplicate them!
    * (even if the refcounting stuff can manage duplicates) */
   ERR(process_tag_list(data1, sz1, data2, sz2, UNIQUE_TAGS, &unique, &u_sz));
 
   ERR(sync_device());
   gmshModelGetBoundary(data1, sz1, &bound1, &n1, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   gmshModelGetBoundary(data2, sz2, &bound2, &n2, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(process_tag_list(bound1, n1, bound2, n2, COMMON_TAGS, &tagout, &tagoutn));
   ASSERT(tagoutn % 2 == 0);
   c = tagoutn / 2;
 
   if(tagoutn == 0) {
     log_message(dev, "Common boundaries list is empty.\n");
     goto exit;
   }
 
   log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
   log_type = dev->log_type;
   if(log) {
     str_init(allocator, &msg);
     msg_initialized = 1;
     logger_print(dev->logger, log_type,
         "Common boundaries specific tag management:\n");
     ERR(str_printf(&msg, "  tags ["));
     for(i = 0; i < tagoutn; i += 2) {
       const int dim = tagout[i];
       const int tag = tagout[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     ERR(str_append_printf(&msg, "] getting a ref to tags ["));
     for(i = 0; i < u_sz; i += 2) {
       const int dim = unique[i];
       const int tag = unique[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     logger_print(dev->logger, log_type, "%s].\n", str_cget(&msg));
   }
 
   out_geom = MEM_CALLOC(allocator, c, sizeof(*out_geom));
   if(!out_geom) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   for(i = 0, n = 0; i < tagoutn; i += 2, n++) {
     int dim = tagout[i];
     int tag = tagout[i+1];
     struct tag_desc* desc;
     ASSERT(dim == 2);
     ERR(geometry_create(out_geom+n));
     out_geom[n]->gmsh_dimTags_n = 2;
     out_geom[n]->gmsh_dimTags = MEM_ALLOC(allocator,
         2 * sizeof(*out_geom[n]->gmsh_dimTags));
     if(!out_geom[n]->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     out_geom[n]->gmsh_dimTags[0] = dim;
     out_geom[n]->gmsh_dimTags[1] = tag;
     ERR(device_register_tags(out_geom[n]));
     /* Need to protect out_geometry's tags by getting a ref on input tags or
      * deleting input geometry will possibly delete them */
     ERR(device_register_ref_to_tags(dim, tag, unique, u_sz));
     /* As the 2D tags will be deleted when the 3D tag they are part of are
      * deleted, they shouldn't be deleted when the geometry they belongs to are
      * released. */
     desc = device_get_description(dim, tag);
     ASSERT(desc);
     desc->delete_policy = Scad_do_not_delete;
   }
 
 exit:
   if(msg_initialized) str_release(&msg);
   if(out_boundaries) *out_boundaries = out_geom;
   if(out_count) *out_count = c;
   if(allocator) {
     MEM_RM(allocator, data1);
     MEM_RM(allocator, data2);
     MEM_RM(allocator, unique);
     MEM_RM(allocator, tagout);
   }
   gmshFree(bound1);
   gmshFree(bound2);
   return res;
 error:
   if(ierr) {
     int dim = INT_MAX;
     if(!mixed_dim_err_msg("common boundaries", geometries, geometries_count, &dim))
       mixed_dim_err_msg("common boundaries", tools, tools_count, &dim);
   }
   if(out_geom) {
     for(i = 0; i < c; i++) {
       if(out_geom[i]) SCAD(geometry_ref_put(out_geom[i]));
     }
     MEM_RM(allocator, out_geom);
     out_geom = NULL;
     c = 0;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometries_common_boundary
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** tools,
    const size_t tools_count,
    struct scad_geometry** out_boundary)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   size_t tagoutn, sz1, sz2, u_sz;
   int* data1 = NULL;
   int* data2 = NULL;
   int* unique = NULL;
   int ierr = 0;
   int* bound1 = NULL;
   int* bound2 = NULL;
   size_t n1, n2;
   struct scad_geometry* out_geom = NULL;
   struct mem_allocator* allocator = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   size_t i, n;
   struct str msg;
   int msg_initialized = 0;
 
   if(!geometries || !geometries_count || !tools || !tools_count || !out_boundary) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data1, &sz1));
   ERR(gather_tags(tools, tools_count, &data2, &sz2));
 
   /* data1 and data2 can have tags in common: deduplicate them!
    * (even if the refcounting stuff can manage duplicates) */
   ERR(process_tag_list(data1, sz1, data2, sz2, UNIQUE_TAGS, &unique, &u_sz));
 
   ERR(sync_device());
   gmshModelGetBoundary(data1, sz1, &bound1, &n1, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   gmshModelGetBoundary(data2, sz2, &bound2, &n2, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(process_tag_list(bound1, n1, bound2, n2, COMMON_TAGS, &tagout, &tagoutn));
   ASSERT(tagoutn % 2 == 0);
 
   if(tagoutn == 0) {
     log_message(dev, "Common boundary list is empty.\n");
   } else {
     log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
     log_type = dev->log_type;
     if(log) {
       str_init(allocator, &msg);
       msg_initialized = 1;
       logger_print(dev->logger, log_type,
           "Common boundary specific tag management:\n");
       ERR(str_printf(&msg, "  tags ["));
       for(i = 0; i < tagoutn; i += 2) {
         const int dim = tagout[i];
         const int tag = tagout[i+1];
         ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
       }
       ERR(str_append_printf(&msg, "] getting a ref to tags ["));
       for(i = 0; i < u_sz; i += 2) {
         const int dim = unique[i];
         const int tag = unique[i+1];
         ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
       }
       logger_print(dev->logger, log_type, "%s].\n", str_cget(&msg));
     }
   }
 
   ERR(geometry_create(&out_geom));
   out_geom->gmsh_dimTags_n = tagoutn;
   out_geom->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
   if(!out_geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   memcpy(out_geom->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
   ERR(device_register_tags(out_geom));
 
   for(i = 0, n = 0; i < tagoutn; i += 2, n++) {
     int dim = tagout[i];
     int tag = tagout[i+1];
     struct tag_desc* desc = device_get_description(dim, tag);
     ASSERT(desc);
     ASSERT(dim == 2);
     /* Need to protect out_geometry's tags by getting a ref on input tags or
      * deleting input geometry will possibly delete them */
     ERR(device_register_ref_to_tags(dim, tag, unique, u_sz));
     /* As the 2D tags will be deleted when the 3D tag they are part of are
      * deleted, they shouldn't be deleted when the geometry they belongs to are
      * released. */
     desc->delete_policy = Scad_do_not_delete;
   }
 
 exit:
   if(msg_initialized) str_release(&msg);
   if(out_boundary) *out_boundary = out_geom;
   if(allocator) {
     MEM_RM(allocator, data1);
     MEM_RM(allocator, data2);
     MEM_RM(allocator, unique);
     MEM_RM(allocator, tagout);
   }
   gmshFree(bound1);
   gmshFree(bound2);
   return res;
 error:
   if(ierr) {
     int dim = INT_MAX;
     if(!mixed_dim_err_msg("common boundary", geometries, geometries_count, &dim))
       mixed_dim_err_msg("common boundary", tools, tools_count, &dim);
   }
   if(out_geom) {
     SCAD(geometry_ref_put(out_geom));
     out_geom = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometry_rotate
   (const struct scad_geometry* geom,
    const double pt[3],
    const double axis[3],
    const double angle,
    struct scad_geometry** out_geometry)
 {
   int ierr = 0;
   struct scad_geometry* out = NULL;
   res_T res = RES_OK;
 
   if(!geom || !pt || !axis || ! out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(scad_geometry_copy(geom, &out));
 
   gmshModelOccRotate(out->gmsh_dimTags, out->gmsh_dimTags_n, SPLIT3(pt),
       SPLIT3(axis), angle, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   if(out_geometry) *out_geometry = out;
   return res;
 error:
   if(out) SCAD(geometry_ref_put(out));
   out = NULL;
   goto exit;
 }
 
 res_T
 scad_geometry_extrude
   (const struct scad_geometry* geom,
    const double dxdydz[3],
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int *tagout = NULL;
   size_t tagoutn;
   size_t i;
   int *ed, *extrude_data = NULL;
   size_t extrude_sz = 0;
   int ierr = 0;
   struct scad_geometry* extrude_geom = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   struct mem_allocator* allocator = NULL;
   struct str msg;
   int init = 0;
 
   if(!geom || !dxdydz || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   /* OCC cannot extrude 3D geometries */
   for(i = 0; i < geom->gmsh_dimTags_n; i += 2) {
     const int dim = geom->gmsh_dimTags[i];
     if(dim != 2) {
       if(str_is_empty(&geom->name)) {
         log_error(get_device(),
             "Cannot extrude unnamed 3D geometry '%p'.\n",
             (void*)geom);
       } else {
         log_error(get_device(),
             "Cannot extrude 3D geometry '%s'.\n",
             str_cget(&geom->name));
       }
       res = RES_BAD_ARG;
       goto error;
     }
   }
 
   gmshModelOccExtrude(geom->gmsh_dimTags, geom->gmsh_dimTags_n, SPLIT3(dxdydz),
       &tagout, &tagoutn, NULL, 0, NULL, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   /* Output includes both the 3D result and its 2D constituents.
    * Keep only 3D entities. */
   for(i = 0; i < tagoutn; i += 2) {
     int dim = tagout[i];
     if(dim == 3) extrude_sz += 2;
   }
   extrude_data = MEM_ALLOC(allocator, extrude_sz * sizeof(*extrude_data));
   if(!extrude_data) {
     res = RES_MEM_ERR;
     goto error;
   }
   ed = extrude_data;
   for(i = 0; i < tagoutn; i += 2) {
     int dim = tagout[i];
     int tag = tagout[i+1];
     if(dim == 3) {
       *ed++ = dim;
       *ed++ = tag;
     }
   }
 
   ERR(geometry_create(&extrude_geom));
   extrude_geom->gmsh_dimTags_n = extrude_sz;
   extrude_geom->gmsh_dimTags = extrude_data;
 
   ERR(device_register_tags(extrude_geom));
 
   log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
   log_type = dev->log_type;
   if(log) {
     str_init(allocator, &msg);
     init = 1;
     logger_print(dev->logger, log_type, "Extrude specific tag management:\n");
     ERR(str_printf(&msg, "  tags ["));
     for(i = 0; i < geom->gmsh_dimTags_n; i += 2) {
       const int dim = geom->gmsh_dimTags[i];
       const int tag = geom->gmsh_dimTags[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     ERR(str_append_printf(&msg, "] getting a ref to tags ["));
     for(i = 0; i < extrude_sz; i += 2) {
       const int dim = extrude_data[i];
       const int tag = extrude_data[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     logger_print(dev->logger, log_type, "%s].\n", str_cget(&msg));
   }
   for(i = 0; i < geom->gmsh_dimTags_n; i += 2) {
     const int dim = geom->gmsh_dimTags[i];
     const int tag = geom->gmsh_dimTags[i+1];
     struct tag_desc* desc = device_get_description(dim, tag);
     ASSERT(dim == 2);
     ASSERT(desc != NULL);
     /* Need to protect input geometry's tags by getting a ref on output tags or
      * deleting out_geometry will possibly delete them */
     ERR(device_register_ref_to_tags(dim, tag, extrude_data, extrude_sz));
     /* As the 2D tags will be deleted when the 3D tag they are part of are
      * deleted, they shouldn't be deleted when the geometry they belongs to is
      * released. */
     desc->delete_policy = Scad_do_not_delete;
   }
 
 exit:
   if(init) str_release(&msg);
   if(out_geometry) *out_geometry = extrude_geom;
   gmshFree(tagout);
   return res;
 error:
   if(extrude_geom) {
     SCAD(geometry_ref_put(extrude_geom));
     extrude_geom = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometry_explode
   (const struct scad_geometry* geom,
    struct scad_geometry*** out_geometry,
    size_t* out_n)
 {
   res_T res = RES_OK;
   int* data = NULL;
   size_t i, n, sz = 0;
   struct scad_geometry** geom_array = NULL;
   struct str name;
   int name_initialized = 0;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!geom || !out_geometry || !out_n) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   allocator = dev->allocator;
   data = geom->gmsh_dimTags;
   sz = geom->gmsh_dimTags_n;
 
   ASSERT(sz % 2 == 0);
   geom_array = MEM_CALLOC(allocator, sz/2, sizeof(*geom_array));
   if(!geom_array) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   for(i = 0, n = 0; i < sz; i += 2, n++) {
     const int dim = data[i];
     const int tag = data[i+1];
     ERR(geometry_create(geom_array+n));
     geom_array[n]->gmsh_dimTags_n = 2;
     geom_array[n]->gmsh_dimTags
       = MEM_ALLOC(allocator, 2 * sizeof(*geom_array[0]->gmsh_dimTags));
     if(!geom_array[n]->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     geom_array[n]->gmsh_dimTags[0] = dim;
     geom_array[n]->gmsh_dimTags[1] = tag;
 
     ERR(device_register_tags(geom_array[n]));
   }
 
 exit:
   if(out_n) *out_n = sz/2 ;
   if(out_geometry) *out_geometry = geom_array;
   if(name_initialized) str_release(&name);
   return res;
 error:
   if(geom_array) {
     for(i = 0, n = 0; i < sz; i += 2, n++) {
       if(geom_array[n]) SCAD(geometry_ref_put(geom_array[n]));
     }
     MEM_RM(allocator, geom_array);
     geom_array = NULL;
   }
   goto exit;
 }
 
 
 res_T
 scad_geometry_copy
   (const struct scad_geometry* geom,
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int* data1;
   int* tagout = NULL;
   size_t sz1, tagoutn;
   int ierr = 0;
   struct scad_geometry* copy = NULL;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
 
   if(!geom || !out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   sz1 = geom->gmsh_dimTags_n;
   data1 = geom->gmsh_dimTags;
   gmshModelOccCopy(data1, sz1, &tagout, &tagoutn, &ierr);
 
   ASSERT(tagoutn % 2 == 0);
   ERR(gmsh_err_to_res_T(ierr));
 
   ERR(geometry_create(&copy));
   copy->gmsh_dimTags_n = tagoutn;
   copy->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
   if(!copy->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   memcpy(copy->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
 
   ERR(device_register_tags(copy));
 
 exit:
   if(out_geometry) *out_geometry = copy;
   gmshFree(tagout);
   return res;
 error:
   if(copy) {
     SCAD(geometry_ref_put(copy));
     copy = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_geometries_set_name
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    const char* prefix_name, /* Can be NULL */
    const size_t from_rank) /* Can be NULL */
 {
   res_T res = RES_OK;
   struct str tmp;
   int initialized = 0;
   struct mem_allocator* allocator = NULL;
   size_t i, cpt;
 
   if(!geometries || geometries_count == 0) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   if(!prefix_name) {
     for(i = 0; i < geometries_count; i++) {
       ERR(geom_set_name(geometries[i], NULL));
     }
   }
   else if(strlen(prefix_name) == 0) {
     log_error(get_device(), "Geometry name \"\" is invalid.\n");
     res = RES_BAD_ARG;
     goto error;
   } else {
     str_init(allocator, &tmp);
     initialized = 1;
     cpt = from_rank;
     for(i = 0; i < geometries_count; i++) {
       ERR(str_printf(&tmp, "%s_%ld", prefix_name, cpt++));
       ERR(geom_set_name(geometries[i], &tmp));
     }
   }
 
 exit:
   if(initialized) str_release(&tmp);
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_set_name
   (struct scad_geometry* geom,
    const char* name) /* Can be NULL */
 {
   res_T res = RES_OK;
   struct str tmp;
   int initialized = 0;
   struct mem_allocator* allocator = NULL;
 
   if(!geom) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   if(!name) {
     ERR(geom_set_name(geom, NULL));
   }
   else if(strlen(name) == 0) {
     log_error(get_device(), "Geometry name \"\" is invalid.\n");
     res = RES_BAD_ARG;
     goto error;
   } else {
     str_init(allocator, &tmp);
     initialized = 1;
     ERR(str_set(&tmp, name));
     ERR(geom_set_name(geom, &tmp));
   }
 
 exit:
   if(initialized) str_release(&tmp);
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometry_translate
   (const struct scad_geometry* geom,
    const double dxdydz[3],
    struct scad_geometry** out_geometry)
 {
   int ierr = 0;
   struct scad_geometry* out = NULL;
   res_T res = RES_OK;
 
   if(!geom || !dxdydz || ! out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(scad_geometry_copy(geom, &out));
 
   gmshModelOccTranslate(out->gmsh_dimTags, out->gmsh_dimTags_n, SPLIT3(dxdydz),
       &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   if(out_geometry) *out_geometry = out;
   return res;
 error:
   if(out) SCAD(geometry_ref_put(out));
   out = NULL;
   goto exit;
 }
 
 static char
 geom_uses_dimTag
   (const struct scad_geometry* g,
    const int dim,
    const int tag)
 {
   size_t i;
   ASSERT(g);
   for(i = 0; i < g->gmsh_dimTags_n; i+=2) {
     int d = g->gmsh_dimTags[i];
     int t = g->gmsh_dimTags[i+1];
     if(dim == d && tag == t) return 1;
   }
   return 0;
 }
 
 res_T
 scad_geometries_partition
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    const int flags,
    struct scad_geometry** out_geometries)
 {
   res_T res = RES_OK;
   size_t i;
   int* tagout = NULL;
   int** map = NULL;
   size_t* mapn = NULL;
   size_t tagoutn = 0, mapnn = 0, sz;
   int* data = NULL;
   int ierr = 0;
   struct scad_geometry** geoms = NULL;
   struct htable_mappings m2, m3;
   struct scad_device* dev = get_device();
   struct htable_tags t2, t3;
   struct htable_tags_iterator it, end;
   int ht_initialized = 0, tmp_initialized = 0;
   struct mem_allocator* allocator = NULL;
   int dont_call_fragment = 0;
   char* overlap = NULL;
   const int invalid_flags =
     ~(SCAD_ALLOW_OVERLAPPING | SCAD_DUMP_ON_OVERLAPPING_ERROR);
   const int dump_overlapping_err = flags & SCAD_DUMP_ON_OVERLAPPING_ERROR;
   static size_t err_cpt = 0;
   struct str tmp;
 
   if(!geometries || !geometries_count || !out_geometries || flags & invalid_flags) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
   str_init(allocator, &tmp);
   tmp_initialized = 1;
 
   ERR(gather_tags(geometries, geometries_count, &data, &sz));
 
   /* Create output geometries */
   geoms = MEM_CALLOC(allocator, geometries_count, sizeof(*geoms));
   if(!geoms) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   dont_call_fragment = (sz == 2);
   if(dont_call_fragment) {
     /* gmshModelOccFragment doesn't allow to process a single entity: need to
      * create the result (same as input) here.
      * Not so simple as one could have provided more than 1 input geometry with
      * identical or no dim.tag(s). */
     mapnn = 1;
     map = gmshMalloc(sizeof(*map));
     map[0] = gmshMalloc(2 * sizeof(**map));
     mapn = gmshMalloc(sizeof(*mapn));
     if(!map || !map[0] || !mapn) {
       res = RES_MEM_ERR;
       goto error;
     }
     mapn[0] = 2;
     map[0][0] = data[0];
     map[0][1] = data[1];
   } else {
     /* As a general principle, we don't remove gmsh objects directly; they are
      * only removed from scad_geometry_ref_put when their tags are no longer
      * used by any star-cad geometry. */
     gmshModelOccFragment(data, sz, NULL, 0, &tagout, &tagoutn, &map, &mapn,
         &mapnn, -1, 0, 0, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
     ASSERT(sz == 2*mapnn); /* Because input tags where deduplicated */
 
     /* Check first if there was an overlapping problem */
     if(!(flags & SCAD_ALLOW_OVERLAPPING)) {
       /* No overlapping means that each tag in geometries is translated into a
        * single tag in map */
       size_t ov = 0;
       overlap = MEM_CALLOC(allocator, geometries_count, sizeof(*overlap));
       if(!overlap) {
         res = RES_MEM_ERR;
         goto error;
       }
       for(i = 0; i < mapnn; i++) {
         if(mapn[i] != 2) {
           int dim = data[2*i];
           int tag = data[2*i+1];
           size_t k;
           ov++;
           res = RES_BAD_ARG;
           /* dim.tag #i overlaps: search geometries using it in input */
           for(k = 0; k < geometries_count; k++) {
             overlap[k] |= geom_uses_dimTag(geometries[k], dim, tag);
           }
         }
       }
       if(ov) {
         size_t k, item_cpt = 0;
         res_T tmp_err;
         if(dump_overlapping_err) {
           ERR(scad_scene_mesh());
         }
         for(k = 0; k < geometries_count; k++) {
           struct scad_geometry* g = geometries[k];
           if(!overlap[k])
             continue;
           tmp_err = RES_BAD_OP;
           if(dump_overlapping_err) {
             if(str_is_empty(&g->name)) {
               str_printf(&tmp, "unamed_partition_error_%lu_%lu",
                   err_cpt, (long unsigned)item_cpt++);
               tmp_err =
                 scad_stl_export(g, str_cget(&tmp), SCAD_KEEP_NORMALS_UNCHANGED, 0);
               if(tmp_err == RES_OK) {
                 log_error(get_device(),
                     "Unnamed geometry '%p' overlapping (dumped in '%s).\n",
                     (void*)g, str_cget(&tmp));
               } else {
                 log_error(get_device(), "Could not dump geometry.\n");
               }
             } else {
               str_printf(&tmp, "%s_partition_error_%lu_%lu",
                   str_cget(&g->name), err_cpt, (long unsigned)item_cpt++);
               tmp_err =
                 scad_stl_export(g, str_cget(&tmp), SCAD_KEEP_NORMALS_UNCHANGED, 0);
               if(tmp_err == RES_OK) {
                 log_error(get_device(),
                     "Geometry '%s' overlapping (dumped in '%s).\n",
                     str_cget(&g->name), str_cget(&tmp));
               } else {
                 log_error(get_device(), "Could not dump geometry.\n");
               }
             }
           }
           if(tmp_err != RES_OK) { /* not dumped */
             if(str_is_empty(&g->name)) {
               log_error(get_device(), "Unnamed geometry '%p' overlapping.\n",
                   (void*)g);
             } else {
               log_error(get_device(), "Geometry '%s' overlapping.\n",
                   str_cget(&g->name));
             }
           }
         }
         err_cpt++;
       }
       if(ov) {
         res = RES_BAD_ARG;
         goto error;
       }
     }
   }
 
   /* Create htables (mappings are used to ease access) */
   htable_mappings_init(allocator, &m2);
   htable_mappings_init(allocator, &m3);
   htable_tags_init(allocator, &t2);
   htable_tags_init(allocator, &t3);
   ht_initialized = 1;
   for(i = 0; i < sz; i += 2) {
     int dim = data[i];
     int tag = data[i+1];
     size_t mapping = i/2;
     struct htable_mappings* mn = (dim == 2) ? &m2 : &m3;
     ASSERT(dim == 2 || dim == 3);
     ERR(htable_mappings_set(mn, &tag, &mapping));
   }
 
   for(i = 0; i < geometries_count; i++) {
     struct scad_geometry* geom = geometries[i];
     size_t c, n, j;
     int* dt = NULL;
     /* For each tag in geometries[i] out_geometries[i] includes the mapped tags.
      * Because of the overlapping case, the resulting tags need to be
      * deduplicated though. */
     htable_tags_clear(&t2);
     htable_tags_clear(&t3);
     for(j = 0; j < geom->gmsh_dimTags_n; j += 2) {
       int dim = geom->gmsh_dimTags[j];
       int tag = geom->gmsh_dimTags[j+1];
       struct htable_mappings* mn = (dim == 2) ? &m2 : &m3;
       size_t k;
       size_t* mapping = htable_mappings_find(mn, &tag);
       ASSERT(dim == 2 || dim == 3);
       ASSERT(mapping && *mapping < mapnn);
       for(k = 0; k < mapn[*mapping]; k += 2) {
         char one = 1;
         int d = map[*mapping][k];
         int t = map[*mapping][k+1];
         struct htable_tags* tn = (d == 2) ? &t2 : &t3;
         ERR(htable_tags_set(tn, &t, &one));
       }
     }
     /* Allocate result */
     n = htable_tags_size_get(&t2) + htable_tags_size_get(&t3);
     dt = MEM_ALLOC(allocator, sizeof(*dt) * 2 * n);
     if(!dt) {
       res = RES_MEM_ERR;
       goto error;
     }
     /* Copy tags */
     c = 0;
     htable_tags_begin(&t2, &it);
     htable_tags_end(&t2, &end);
     while(!htable_tags_iterator_eq(&it, &end)) {
       dt[c++] = 2;
       dt[c++] = *htable_tags_iterator_key_get(&it);
       htable_tags_iterator_next(&it);
     }
     htable_tags_begin(&t3, &it);
     htable_tags_end(&t3, &end);
     while(!htable_tags_iterator_eq(&it, &end)) {
       dt[c++] = 3;
       dt[c++] = *htable_tags_iterator_key_get(&it);
       htable_tags_iterator_next(&it);
     }
     ASSERT(c == 2*n);
 
     /* Create geometry */
     ERR(geometry_create(geoms+i));
     geoms[i]->gmsh_dimTags_n = c;
     geoms[i]->gmsh_dimTags = dt;
     ERR(device_register_tags(geoms[i]));
   }
   memcpy(out_geometries, geoms, geometries_count * sizeof(*geoms));
 
 exit:
   gmshFree(mapn);
   free_gmsh_map(map, mapnn);
   if(ht_initialized) {
     htable_mappings_release(&m2);
     htable_mappings_release(&m3);
     htable_tags_release(&t2);
     htable_tags_release(&t3);
   }
   if(tmp_initialized) str_release(&tmp);
   if(allocator) {
     MEM_RM(allocator, data);
     MEM_RM(allocator, geoms);
     MEM_RM(allocator, overlap);
   }
   gmshFree(tagout);
   return res;
 error:
   if(geoms) {
     for(i = 0; i < geometries_count; i++) {
       if(geoms[i]) {
         SCAD(geometry_ref_put(geoms[i]));
       }
     }
   }
   if(tagout) {
      gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
      /* Do not check ierr, first because we are in error handler already, but
       * also because an error is expected if these OCC entities have allready
       * been freed through the call to ref_put above */
   }
   goto exit;
 }
 
 res_T
 scad_geometries_swap
   (struct scad_geometry** pool1,
    struct scad_geometry** pool2,
    const size_t count,
    const int flags)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   size_t i;
   struct mem_allocator* allocator = NULL;
   struct str tmp, msg;
 
   if(!pool1 || !pool2 || !flags || count == 0) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
   if(flags & SCAD_SWAP_NAME) str_init(allocator, &tmp);
   if(flags & SCAD_SWAP_GEOMETRY && dev->log) str_init(allocator, &msg);
   for(i = 0; i < count; i++) {
     struct scad_geometry *g1 = pool1[i], *g2 = pool2[i];
     size_t c1 = g1->gmsh_dimTags_n, c2 = g2->gmsh_dimTags_n;
     int *dt1 = g1->gmsh_dimTags, *dt2 = g2->gmsh_dimTags;
     if(pool1[i] == pool2[i])
       continue;
     /* Swap content according to flags. Don't swap  refcount! */
     if(flags & SCAD_SWAP_NAME) {
       if(dev->log) {
         if(str_is_empty(&g1->name) && str_is_empty(&g2->name)) {
           /* Do nothing */
         }
         else if(str_is_empty(&g1->name)) {
           logger_print(dev->logger, dev->log_type,
               "Swapping names for geometry %p and geometry '%s'.\n",
               (void*)g1, str_cget(&g2->name));
           logger_print(dev->logger, dev->log_type,
               "Geometry '%s' is now unnamed geometry %p.\n",
               str_cget(&g2->name), (void*)g2);
         }
         else if(str_is_empty(&g2->name)) {
           logger_print(dev->logger, dev->log_type,
               "Swapping names for geometry %p and geometry '%s'.\n",
               (void*)g2, str_cget(&g1->name));
           logger_print(dev->logger, dev->log_type,
               "Geometry '%s' is now unnamed geometry %p.\n",
               str_cget(&g1->name), (void*)g1);
         } else { /* Both named */
           logger_print(dev->logger, dev->log_type,
               "Swapping names for geometries '%s' and '%s'.\n",
               str_cget(&g1->name), str_cget(&g1->name));
         }
       }
       if(!str_is_empty(&g1->name)) {
         ERR(htable_names_set(&dev->geometry_names, &g1->name, &g2));
       }
       if(!str_is_empty(&g2->name)) {
         ERR(htable_names_set(&dev->geometry_names, &g2->name, &g1));
       }
       if(!str_is_empty(&g1->name) || !str_is_empty(&g2->name)) {
         ERR(str_copy(&tmp, &g1->name));
         ERR(str_copy(&g1->name, &g2->name));
         ERR(str_copy(&g2->name, &tmp));
       }
     }
     if(flags & SCAD_SWAP_GEOMETRY) {
       /* Swap in tag2geom tables */
       size_t n;
       if(dev->log) {
         if(str_is_empty(&g1->name) && str_is_empty(&g2->name)) {
           logger_print(dev->logger, dev->log_type,
               "Swapping tags for unnamed geometries %p and %p.\n",
               (void*)g1, (void*)g2);
         }
         else if(str_is_empty(&g1->name)) {
           logger_print(dev->logger, dev->log_type,
               "Swapping tags for unnamed geometry %p and geometry '%s'.\n",
               (void*)g1, str_cget(&g2->name));
         }
         else if(str_is_empty(&g2->name)) {
           logger_print(dev->logger, dev->log_type,
               "Swapping tags for unnamed geometry %p and geometry '%s'.\n",
               (void*)g2, str_cget(&g1->name));
         }
         else {
           logger_print(dev->logger, dev->log_type,
               "Swapping tags for geometries '%s' and '%s'.\n",
               str_cget(&g1->name), str_cget(&g1->name));
         }
         if(str_is_empty(&g1->name)) {
           ERR(str_printf(&msg,
                 "Tags now registered against unnamed geometry '%p': ",
                 (void*)g1));
         } else {
           ERR(str_printf(&msg, "Tags now registered against geometry '%s': ",
                 str_cget(&g1->name)));
         }
         for(n = 0; n < c2; n += 2) {
           int dim = dt2[n];
           int tag = dt2[n+1];
           ERR(str_append_printf(&msg, (n ? ", %d.%d" : "%d.%d"), dim, tag));
         }
         logger_print(dev->logger, dev->log_type, "%s.\n", str_cget(&msg));
         if(str_is_empty(&g2->name)) {
           ERR(str_printf(&msg,
                 "Tags now registered against unnamed geometry '%p': ",
                 (void*)g2));
         } else {
           ERR(str_printf(&msg, "Tags now registered against geometry '%s': ",
                 str_cget(&g2->name)));
         }
         for(n = 0; n < c1; n += 2) {
           int dim = dt1[n];
           int tag = dt1[n+1];
           ERR(str_append_printf(&msg, (n ? ", %d.%d" : "%d.%d"), dim, tag));
         }
         logger_print(dev->logger, dev->log_type, "%s.\n", str_cget(&msg));
       }
       /* Swap tags */
       SWAP(int*, g1->gmsh_dimTags, g2->gmsh_dimTags);
       SWAP(size_t, g1->gmsh_dimTags_n, g2->gmsh_dimTags_n);
     }
   }
 
 exit:
   if(flags & SCAD_SWAP_NAME) str_release(&tmp);
   if(flags & SCAD_SWAP_GEOMETRY && dev->log) str_release(&msg);
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometries_boundaries
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry*** out_boundaries,
    size_t *out_count)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   size_t tagoutn, sz;
   int* data = NULL;
   int ierr = 0;
   struct scad_geometry** out_geom = NULL;
   struct mem_allocator* allocator = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   size_t i, c = 0, n;
   struct str msg;
   int msg_initialized = 0;
 
   if(!geometries || geometries_count == 0 || !out_boundaries || !out_count) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data, &sz));
   ERR(sync_device());
   gmshModelGetBoundary(data, sz, &tagout, &tagoutn, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ASSERT(tagoutn % 2 == 0);
   c = tagoutn / 2;
   out_geom = MEM_CALLOC(allocator, c, sizeof(*out_geom));
   if(!out_geom) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
   log_type = dev->log_type;
   if(log) {
     str_init(allocator, &msg);
     msg_initialized = 1;
     logger_print(dev->logger, log_type, "Boundary specific tag management:\n");
     ERR(str_printf(&msg, "  tags ["));
     for(i = 0; i < tagoutn; i += 2) {
       const int dim = tagout[i];
       const int tag = tagout[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     ERR(str_append_printf(&msg, "] getting a ref to tags ["));
     for(i = 0; i < sz; i += 2) {
       const int dim = data[i];
       const int tag = data[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     logger_print(dev->logger, log_type, "%s].\n", str_cget(&msg));
   }
 
   for(i = 0, n = 0; i < tagoutn; i += 2, n++) {
     int dim = tagout[i];
     int tag = tagout[i+1];
     struct tag_desc* desc = NULL;
     ASSERT(dim == 2);
     ERR(geometry_create(out_geom+n));
     out_geom[n]->gmsh_dimTags_n = 2;
     out_geom[n]->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*tagout));
     if(!out_geom[n]->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     out_geom[n]->gmsh_dimTags[0] = dim;
     out_geom[n]->gmsh_dimTags[1] = tag;
 
     ERR(device_register_tags(out_geom[n]));
 
     /* Need to protect out_geometry's tags by getting a ref on input tags or
      * deleting input geometry will possibly delete them */
     ERR(device_register_ref_to_tags(dim, tag, data, sz));
     /* As the 2D tags will be deleted when the 3D tag they are part of are
      * deleted, they shouldn't be deleted when the geometry they belongs to are
      * released. */
     desc = device_get_description(dim, tag);
     ASSERT(desc);
     desc->delete_policy = Scad_do_not_delete;
   }
 
 exit:
   if(msg_initialized) str_release(&msg);
   if(allocator) MEM_RM(allocator, data);
   if(out_boundaries) *out_boundaries = out_geom;
   if(out_count) *out_count = c;
   gmshFree(tagout);
   return res;
 error:
   if(out_geom) {
     for(i = 0; i < c; i++) {
       if(out_geom[i]) SCAD(geometry_ref_put(out_geom[i]));
     }
     MEM_RM(allocator, out_geom);
     out_geom = NULL;
     c = 0;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 
 res_T
 scad_geometries_boundary
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    struct scad_geometry** out_boundary)
 {
   res_T res = RES_OK;
   int* tagout = NULL;
   size_t tagoutn, sz;
   int* data = NULL;
   int ierr = 0;
   struct scad_geometry* out_geom = NULL;
   struct mem_allocator* allocator = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   size_t i;
   struct str msg;
   int msg_initialized = 0;
 
   if(!geometries || geometries_count == 0 || !out_boundary) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   ERR(gather_tags(geometries, geometries_count, &data, &sz));
   ERR(sync_device());
   gmshModelGetBoundary(data, sz, &tagout, &tagoutn, 1, 0, 0, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   ASSERT(tagoutn % 2 == 0);
 
   log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
   log_type = dev->log_type;
   if(log) {
     str_init(allocator, &msg);
     msg_initialized = 1;
     logger_print(dev->logger, log_type, "Boundary specific tag management:\n");
     ERR(str_printf(&msg, "  tags ["));
     for(i = 0; i < tagoutn; i += 2) {
       const int dim = tagout[i];
       const int tag = tagout[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     ERR(str_append_printf(&msg, "] getting a ref to tags ["));
     for(i = 0; i < sz; i += 2) {
       const int dim = data[i];
       const int tag = data[i+1];
       ERR(str_append_printf(&msg, (i ? ", %d.%d" : "%d.%d"), dim, tag));
     }
     logger_print(dev->logger, log_type, "%s].\n", str_cget(&msg));
   }
 
   ERR(geometry_create(&out_geom));
   out_geom->gmsh_dimTags_n = tagoutn;
   out_geom->gmsh_dimTags = MEM_ALLOC(allocator, tagoutn * sizeof(*tagout));
   if(!out_geom->gmsh_dimTags) {
     res = RES_MEM_ERR;
     goto error;
   }
   memcpy(out_geom->gmsh_dimTags, tagout, tagoutn * sizeof(*tagout));
   ERR(device_register_tags(out_geom));
   for(i = 0; i < out_geom->gmsh_dimTags_n; i += 2) {
     const int dim = out_geom->gmsh_dimTags[i];
     const int tag = out_geom->gmsh_dimTags[i+1];
     struct tag_desc* desc = device_get_description(dim, tag);
     ASSERT(dim == 2);
     ASSERT(desc != NULL);
     /* Need to protect input geometry's tags by getting a ref on output tags or
      * deleting out_geometry will possibly delete them */
     ERR(device_register_ref_to_tags(dim, tag, data, sz));
     /* As the 2D tags will be deleted when the 3D tag they are part of are
      * deleted, they shouldn't be deleted when the geometry they belongs to is
      * released. */
     desc->delete_policy = Scad_do_not_delete;
   }
 
 exit:
   if(msg_initialized) str_release(&msg);
   if(allocator) MEM_RM(allocator, data);
   if(out_boundary) *out_boundary = out_geom;
   gmshFree(tagout);
   return res;
 error:
   if(out_geom) {
     SCAD(geometry_ref_put(out_geom));
     out_geom = NULL;
   }
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   goto exit;
 }
 
 res_T
 scad_step_import
   (const char* filename,
    struct scad_geometry*** out_geometry,
    size_t* out_n)
 {
   int ierr;
   int* tagout = NULL;
   size_t tagoutn, i, ga_sz;
   struct str strname;
   int name_initialized = 0;
   struct scad_geometry** geom_array = NULL;
   struct mem_allocator* allocator = NULL;
   struct scad_device* dev = get_device();
   res_T res = RES_OK;
 
   if(!filename || !out_geometry || !out_n) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   gmshModelOccImportShapes(filename, &tagout, &tagoutn, 1, "step", &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
   ASSERT(tagoutn % 2 == 0);
   ga_sz = tagoutn / 2;
   allocator = get_device()->allocator;
   geom_array = MEM_CALLOC(allocator, ga_sz, sizeof(*geom_array));
   if(!geom_array) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   str_init(allocator, &strname);
   name_initialized = 1;
   for(i=0; i<ga_sz; ++i) {
     ERR(geometry_create(geom_array+i));
     geom_array[i]->gmsh_dimTags_n = 2;
     geom_array[i]->gmsh_dimTags
       = MEM_ALLOC(allocator, 2 * sizeof(*geom_array[i]->gmsh_dimTags));
     if(!geom_array[i]->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     geom_array[i]->gmsh_dimTags[0] = tagout[2*i];
     geom_array[i]->gmsh_dimTags[1] = tagout[2*i+1];
 
     ERR(device_register_tags(geom_array[i]));
   }
 
 exit:
   gmshFree(tagout);
   if(out_n) *out_n = ga_sz ;
   if(out_geometry) *out_geometry = geom_array;
   if(name_initialized) str_release(&strname);
   return res;
 error:
   if(tagout) gmshModelOccRemove(tagout, tagoutn, 1, &ierr);
   ga_sz = 0;
   if(geom_array) {
     for(i=0; i<ga_sz; ++i) {
       if(geom_array[i]) SCAD(geometry_ref_put(geom_array[i]));
     }
     MEM_RM(allocator, geom_array);
     geom_array = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_geometry_get_normal
   (struct scad_geometry* geom,
    const double p[3],
    double N[3],
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr = 0;
   size_t i;
   const int* data = NULL;
   size_t sz = 0;
   struct scad_geometry* out = NULL;
   struct scad_device* dev = get_device();
   int log;
   enum log_type log_type;
   struct mem_allocator* allocator = NULL;
   double* coord = NULL;
   double* pcoord = NULL;
   double* normals = NULL;
   struct darray_int tags;
   int initialized = 0;
   double d, min_d = DBL_MAX, pcoord_min[2];
   int min_tag = -1;
   size_t normals_n;
 
   if(!geom || !p || !N) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
   allocator = dev->allocator;
 
   darray_int_init(dev->allocator, &tags);
   initialized = 1;
 
   ERR(get_2d_tags(geom, &tags));
   data = darray_int_cdata_get(&tags);
   sz = darray_int_size_get(&tags);
 
   log = (dev->options.Misc.LogRefCounting & SCAD_LOG_DIMTAGS_ALL);
   log_type = dev->log_type;
 
   /* Find the closest point on tags of geom */
   for(i = 0; i < sz; ++i) {
     size_t pcoord_n;
     size_t coord_n;
     const int dim = 2;
     int tag = data[i];
     double tmp[3];
 
     gmshModelGetClosestPoint(dim, tag, p, 3, &coord, &coord_n, &pcoord,
         &pcoord_n, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
     ASSERT(pcoord_n == (size_t)dim);
     ASSERT(coord_n == 3);
     d = d3_len(d3_sub(tmp, p, coord));
     if(d < min_d) {
       min_d = d;
       min_tag = tag;
       d2_set(pcoord_min, pcoord);
     }
     gmshFree(coord);
     gmshFree(pcoord);
     coord = pcoord = NULL;
   }
   if(min_d == DBL_MAX) { /* At least if sz==0 */
     goto exit;
   }
 
   /* Get the normal at the selected point */
   gmshModelGetNormal(min_tag, pcoord_min, 2, &normals, &normals_n, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
   ASSERT(normals_n == 3);
   d3_set(N, normals);
 
   /* Create a geometry if required */
   if(out_geometry) {
     ERR(geometry_create(&out));
     out->gmsh_dimTags = MEM_ALLOC(allocator, 2 * sizeof(*out->gmsh_dimTags));
     if(!out->gmsh_dimTags) {
       res = RES_MEM_ERR;
       goto error;
     }
     out->gmsh_dimTags_n = 2;
     out->gmsh_dimTags[0] = 2;
     out->gmsh_dimTags[1] = min_tag;
     ERR(device_register_tags(out));
 
     /* Need to protect geometries' tags or deleting out geometry will possibly
      * delete them */
     if(log) {
       logger_print(dev->logger, log_type,
           "Tag %d.%d getting a reference to other tags.\n", 2, min_tag);
     }
     ERR(device_register_ref_to_tags(2, min_tag, geom->gmsh_dimTags,
           geom->gmsh_dimTags_n));
     if(log) {
       logger_print(dev->logger, log_type,
           "Tag %d.%d getting a reference to other tags done.\n", 2, min_tag);
     }
   }
 
 exit:
   gmshFree(coord);
   gmshFree(pcoord);
   gmshFree(normals);
   if(initialized) darray_int_release(&tags);
   if(out_geometry) *out_geometry = out;
   return res;
 error:
   if(out) {
     SCAD(geometry_ref_put(out));
     out = NULL;
   }
   goto exit;
 }
 
 res_T
 scad_geometry_dilate
   (const struct scad_geometry* geom,
    const double center[3],
    const double scale[3],
    struct scad_geometry** out_geometry)
 {
   res_T res = RES_OK;
   int ierr = 0;
   struct scad_geometry* out = NULL;
 
   if(!geom || !scale|| !center || ! out_geometry) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(scad_geometry_copy(geom, &out));
 
   gmshModelOccDilate(out->gmsh_dimTags, out->gmsh_dimTags_n, SPLIT3(center),
       SPLIT3(scale), &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   if(out_geometry) *out_geometry = out;
   return res;
 error:
   if(out) SCAD(geometry_ref_put(out));
   out = NULL;
   goto exit;
 }
 
 res_T
 scad_geometries_set_mesh_size_modifier
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    enum scad_size_modifier_type type,
    double modifier)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   int ierr, dim, some = 0;
   int* tagout[4] = { NULL, NULL, NULL, NULL };
   size_t tagoutn[4] = { 0, 0, 0, 0}, i;
 
   if(!geometries || geometries_count == 0 || modifier <= 0
       || (type != SCAD_ABSOLUTE_SIZE && type != SCAD_SIZE_FACTOR))
   {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   if(type == SCAD_ABSOLUTE_SIZE) modifier = -modifier;
   ERR(gather_tags_recursive(geometries, geometries_count, 0, tagout, tagoutn));
   for(dim = 0; dim < 4; dim++) {
     for(i = 0; i < tagoutn[dim]; i += 2) {
       int d = tagout[dim][i];
       int tag = tagout[dim][i+1];
       struct htable_size_modifiers* modifiers = dev->size_modifiers_by_dim + dim;
       double* f_ptr = htable_size_modifiers_find(modifiers, &tag);
       ASSERT(d == dim);(void)d;
       if(f_ptr && modifier == 1) { /* Size factor of 1 => no modifier */
         size_t c = htable_size_modifiers_erase(modifiers, &tag);
         ASSERT(c == 1);(void)c;
       } else if(f_ptr) {
         *f_ptr = modifier;
       } else {
         ERR(htable_size_modifiers_set(modifiers, &tag, &modifier));
       }
     }
   }
   for(dim = 0; dim < 4; dim++) {
     if(htable_size_modifiers_size_get(dev->size_modifiers_by_dim + dim) > 0) {
       some = 1;
       break;
     }
   }
 
   gmshModelMeshSetSizeCallback((some ? size_callback : NULL), NULL, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   if(dev) {
     for(dim = 0; dim < 4; dim++) {
       MEM_RM(dev->allocator, tagout[dim]);
     }
   }
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometries_set_mesh_algorithm
   (struct scad_geometry** geometries,
    const size_t geometries_count,
    enum scad_mesh_algorithm algorithm)
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   int ierr;
   int* tagout[4] = { NULL, NULL, NULL, NULL };
   size_t tagoutn[4], i;
 
   if(!geometries || geometries_count == 0 ||
       (algorithm != SCAD_MESHADAPT &&
        algorithm != SCAD_AUTOMATIC &&
        algorithm != SCAD_INITIAL_MESH_ONLY &&
        algorithm != SCAD_DELAUNAY &&
        algorithm != SCAD_FRONTAL_DELAUNAY &&
        algorithm != SCAD_QUASI_STRUCTURED))
   {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(gather_tags_recursive(geometries, geometries_count, 2, tagout, tagoutn));
   for(i = 0; i < tagoutn[2]; i += 2) {
     int dim = tagout[2][i];
     int tag = tagout[2][i+1];
     gmshModelMeshSetAlgorithm(dim, tag, (int)algorithm, &ierr);
     ERR(gmsh_err_to_res_T(ierr));
   }
 
 exit:
   if(dev) {
     for(i = 2; i < 4; i++) {
       MEM_RM(dev->allocator, tagout[i]);
     }
   }
   return res;
 error:
   goto exit;
 }
 
 res_T
 scad_geometries_set_periodic
   (struct scad_geometry** source,
    const size_t source_count,
    struct scad_geometry** target,
    const size_t target_count,
    double affine[16])
 {
   res_T res = RES_OK;
   struct scad_device* dev = get_device();
   struct mem_allocator* allocator = NULL;
   int ierr;
   int* src_dimTagout[4] = { NULL, NULL, NULL, NULL };
   int* tgt_dimTagout[4] = { NULL, NULL, NULL, NULL };
   int* src_tags = NULL, *tgt_tags = NULL;
   size_t src_dimTagoutn[4] = { 0,0,0,0 }, tgt_dimTagoutn[4] = { 0,0,0,0 };
   size_t src_tagsn = 0, tgt_tagsn = 0,  i;
 
   if(!source || source_count == 0 || !target || target_count == 0 || !affine) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   ERR(check_device(FUNC_NAME));
 
   ERR(gather_tags_recursive(source, source_count, 2, src_dimTagout, src_dimTagoutn));
   ERR(gather_tags_recursive(target, target_count, 2, tgt_dimTagout, tgt_dimTagoutn));
   ASSERT(src_dimTagoutn[2] % 2 == 0 && tgt_dimTagoutn[2] % 2 == 0);
   src_tagsn = src_dimTagoutn[2] / 2;
   tgt_tagsn = tgt_dimTagoutn[2] / 2;
   if(src_tagsn == 0 || tgt_tagsn == 0) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   allocator = dev->allocator;
   src_tags = MEM_ALLOC(allocator, src_tagsn * sizeof(*src_tags));
   tgt_tags = MEM_ALLOC(allocator, tgt_tagsn * sizeof(*tgt_tags));
   if(!src_tags || !tgt_tags) {
     res = RES_MEM_ERR;
     goto error;
   }
   for(i = 0; i < src_tagsn; i += 2) {
     src_tags[i] = src_dimTagout[2][i+1];
   }
   for(i = 0; i < tgt_tagsn; i += 2) {
     tgt_tags[i] = tgt_dimTagout[2][i+1];
   }
   gmshModelMeshSetPeriodic(2, tgt_tags, tgt_tagsn, src_tags, src_tagsn,
       affine, 16, &ierr);
   ERR(gmsh_err_to_res_T(ierr));
 
 exit:
   if(dev) {
     for(i = 2; i < 4; i++) {
       MEM_RM(dev->allocator, src_dimTagout[i]);
       MEM_RM(dev->allocator, tgt_dimTagout[i]);
     }
     MEM_RM(dev->allocator, src_tags);
     MEM_RM(dev->allocator, tgt_tags);
   }
   return res;
 error:
   goto exit;
 }